Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BinaryCascade.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 
27 #include "globals.hh"
28 #include "G4PhysicalConstants.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4Proton.hh"
31 #include "G4Neutron.hh"
32 #include "G4LorentzRotation.hh"
33 #include "G4BinaryCascade.hh"
34 #include "G4KineticTrackVector.hh"
35 #include "G4DecayKineticTracks.hh"
36 #include "G4ReactionProductVector.hh"
37 #include "G4Track.hh"
38 #include "G4V3DNucleus.hh"
39 #include "G4Fancy3DNucleus.hh"
40 #include "G4Scatterer.hh"
41 #include "G4MesonAbsorption.hh"
42 #include "G4ping.hh"
43 #include "G4Delete.hh"
44 
45 #include "G4CollisionManager.hh"
46 #include "G4Absorber.hh"
47 
48 #include "G4CollisionInitialState.hh"
49 #include "G4ListOfCollisions.hh"
50 #include "G4Fragment.hh"
51 #include "G4RKPropagation.hh"
52 
53 #include "G4NuclearShellModelDensity.hh"
54 #include "G4NuclearFermiDensity.hh"
55 #include "G4FermiMomentum.hh"
56 
57 #include "G4PreCompoundModel.hh"
58 #include "G4ExcitationHandler.hh"
59 #include "G4HadronicInteractionRegistry.hh"
60 
61 #include "G4FermiPhaseSpaceDecay.hh"
62 
63 #include "G4PreCompoundModel.hh"
64 
65 #include <algorithm>
66 #include "G4ShortLivedConstructor.hh"
67 #include <typeinfo>
68 
69 // turn on general debugging info, and consistency checks
70 
71 //#define debug_G4BinaryCascade 1
72 
73 // more detailed debugging -- deprecated
74 //#define debug_H1_BinaryCascade 1
75 
76 // specific debugging info per method or functionality
77 //#define debug_BIC_ApplyCollision 1
78 //#define debug_BIC_CheckPauli 1
79 //#define debug_BIC_CorrectFinalPandE 1
80 //#define debug_BIC_Propagate 1
81 //#define debug_BIC_Propagate_Excitation 1
82 //#define debug_BIC_Propagate_Collisions 1
83 //#define debug_BIC_Propagate_finals 1
84 //#define debug_BIC_DoTimeStep 1
85 //#define debug_BIC_CorrectBarionsOnBoundary 1
86 //#define debug_BIC_GetExcitationEnergy 1
87 //#define debug_BIC_DeexcitationProducts 1
88 //#define debug_BIC_FinalNucleusMomentum 1
89 //#define debug_BIC_Final4Momentum 1
90 //#define debug_BIC_FillVoidnucleus 1
91 //#define debug_BIC_FindFragments 1
92 //#define debug_BIC_BuildTargetList 1
93 //#define debug_BIC_FindCollision 1
94 //#define debug_BIC_return 1
95 
96 //-------
97 //#if defined(debug_G4BinaryCascade)
98 #if 0
99  #define _CheckChargeAndBaryonNumber_(val) CheckChargeAndBaryonNumber(val)
100  //#define debugCheckChargeAndBaryonNumberverbose 1
101 #else
102  #define _CheckChargeAndBaryonNumber_(val)
103 #endif
104 //#if defined(debug_G4BinaryCascade)
105 #if 0
106  #define _DebugEpConservation(val) DebugEpConservation(val)
107  //#define debugCheckChargeAndBaryonNumberverbose 1
108 #else
109  #define _DebugEpConservation(val)
110 #endif
111 //
112 // C O N S T R U C T O R S A N D D E S T R U C T O R S
113 //
114 
115 G4BinaryCascade::G4BinaryCascade(G4VPreCompoundModel* ptr) :
116 G4VIntraNuclearTransportModel("Binary Cascade", ptr)
117 {
118  // initialise the resonance sector
119  G4ShortLivedConstructor ShortLived;
120  ShortLived.ConstructParticle();
121 
122  theCollisionMgr = new G4CollisionManager;
123  theDecay=new G4BCDecay;
124  theImR.push_back(theDecay);
125  theLateParticle= new G4BCLateParticle;
126  G4MesonAbsorption * aAb=new G4MesonAbsorption;
127  theImR.push_back(aAb);
128  G4Scatterer * aSc=new G4Scatterer;
129  theH1Scatterer = new G4Scatterer;
130  theImR.push_back(aSc);
131 
132  thePropagator = new G4RKPropagation;
133  theCurrentTime = 0.;
134  theBCminP = 45*MeV;
135  theCutOnP = 90*MeV;
136  theCutOnPAbsorb= 0*MeV; // No Absorption of slow Mesons, other than above G4MesonAbsorption
137 
138  // reuse existing pre-compound model
139  if(!ptr) {
140  G4HadronicInteraction* p =
141  G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
142  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
143  if(!pre) { pre = new G4PreCompoundModel(); }
144  SetDeExcitation(pre);
145  }
146  theExcitationHandler = GetDeExcitation()->GetExcitationHandler();
147  SetMinEnergy(0.0*GeV);
148  SetMaxEnergy(10.1*GeV);
149  //PrintWelcomeMessage();
150  thePrimaryEscape = true;
151  thePrimaryType = 0;
152 
153  SetEnergyMomentumCheckLevels(1.0*perCent, 1.0*MeV);
154 
155  // init data members
156  currentA=currentZ=0;
157  lateA=lateZ=0;
158  initialA=initialZ=0;
159  projectileA=projectileZ=0;
160  currentInitialEnergy=initial_nuclear_mass=0.;
161  massInNucleus=0.;
162  theOuterRadius=0.;
163 }
164 
165 /*
166 G4BinaryCascade::G4BinaryCascade(const G4BinaryCascade& )
167 : G4VIntraNuclearTransportModel("Binary Cascade")
168 {
169 }
170  */
171 
172 G4BinaryCascade::~G4BinaryCascade()
173 {
174  ClearAndDestroy(&theTargetList);
175  ClearAndDestroy(&theSecondaryList);
176  ClearAndDestroy(&theCapturedList);
177  delete thePropagator;
178  delete theCollisionMgr;
179  std::for_each(theImR.begin(), theImR.end(), Delete<G4BCAction>());
180  delete theLateParticle;
181  //delete theExcitationHandler;
182  delete theH1Scatterer;
183 }
184 
185 void G4BinaryCascade::ModelDescription(std::ostream& outFile) const
186 {
187  outFile << "G4BinaryCascade is an intra-nuclear cascade model in which\n"
188  << "an incident hadron collides with a nucleon, forming two\n"
189  << "final-state particles, one or both of which may be resonances.\n"
190  << "The resonances then decay hadronically and the decay products\n"
191  << "are then propagated through the nuclear potential along curved\n"
192  << "trajectories until they re-interact or leave the nucleus.\n"
193  << "This model is valid for incident pions up to 1.5 GeV and\n"
194  << "nucleons up to 10 GeV.\n"
195  << "The remaining excited nucleus is handed on to ";
196  if (theDeExcitation) // pre-compound
197  {
198  outFile << theDeExcitation->GetModelName() << " : \n ";
199  theDeExcitation->DeExciteModelDescription(outFile);
200  }
201  else if (theExcitationHandler) // de-excitation
202  {
203  outFile << "G4ExcitationHandler"; //theExcitationHandler->GetModelName();
204  theExcitationHandler->ModelDescription(outFile);
205  }
206  else
207  {
208  outFile << "void.\n";
209  }
210  outFile<< " \n";
211 }
212 void G4BinaryCascade::PropagateModelDescription(std::ostream& outFile) const
213 {
214  outFile << "G4BinaryCascade propagtes secondaries produced by a high\n"
215  << "energy model through the wounded nucleus.\n"
216  << "Secondaries are followed after the formation time and if\n"
217  << "within the nucleus are propagated through the nuclear\n"
218  << "potential along curved trajectories until they interact\n"
219  << "with a nucleon, decay, or leave the nucleus.\n"
220  << "An interaction of a secondary with a nucleon produces two\n"
221  << "final-state particles, one or both of which may be resonances.\n"
222  << "Resonances decay hadronically and the decay products\n"
223  << "are in turn propagated through the nuclear potential along curved\n"
224  << "trajectories until they re-interact or leave the nucleus.\n"
225  << "This model is valid for pions up to 1.5 GeV and\n"
226  << "nucleons up to about 3.5 GeV.\n"
227  << "The remaining excited nucleus is handed on to ";
228  if (theDeExcitation) // pre-compound
229  {
230  outFile << theDeExcitation->GetModelName() << " : \n ";
231  theDeExcitation->DeExciteModelDescription(outFile);
232  }
233  else if (theExcitationHandler) // de-excitation
234  {
235  outFile << "G4ExcitationHandler"; //theExcitationHandler->GetModelName();
236  theExcitationHandler->ModelDescription(outFile);
237  }
238  else
239  {
240  outFile << "void.\n";
241  }
242 outFile<< " \n";
243 }
244 
245 //----------------------------------------------------------------------------
246 
247 //
248 // I M P L E M E N T A T I O N
249 //
250 
251 
252 //----------------------------------------------------------------------------
253 G4HadFinalState * G4BinaryCascade::ApplyYourself(const G4HadProjectile & aTrack,
254  G4Nucleus & aNucleus)
255 //----------------------------------------------------------------------------
256 {
257  if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction starts ######### "<< G4endl;
258 
259  G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
260  const G4ParticleDefinition * definition = aTrack.GetDefinition();
261 
262  if(initial4Momentum.e()-initial4Momentum.m()<theBCminP &&
263  ( definition==G4Neutron::NeutronDefinition() || definition==G4Proton::ProtonDefinition() ) )
264  {
265  return theDeExcitation->ApplyYourself(aTrack, aNucleus);
266  }
267 
268  theParticleChange.Clear();
269  // initialize the G4V3DNucleus from G4Nucleus
270  the3DNucleus = new G4Fancy3DNucleus;
271 
272  // Build a KineticTrackVector with the G4Track
273  G4KineticTrackVector * secondaries;// = new G4KineticTrackVector;
274  G4ThreeVector initialPosition(0., 0., 0.); // will be set later
275 
276  if(!getenv("I_Am_G4BinaryCascade_Developer") )
277  {
278  if(definition!=G4Neutron::NeutronDefinition() &&
279  definition!=G4Proton::ProtonDefinition() &&
280  definition!=G4PionPlus::PionPlusDefinition() &&
281  definition!=G4PionMinus::PionMinusDefinition() )
282  {
283  G4cerr << "You are trying to use G4BinaryCascade with " <<definition->GetParticleName()<<" as projectile."<<G4endl;
284  G4cerr << "G4BinaryCascade should not be used for projectiles other than nucleons or pions."<<G4endl;
285  G4cerr << "If you want to continue, please switch on the developer environment: "<<G4endl;
286  G4cerr << "setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<G4endl;
287  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade - used for unvalid particle type - Fatal");
288  }
289  }
290 
291  // keep primary
292  thePrimaryType = definition;
293  thePrimaryEscape = false;
294 
295  // try until an interaction will happen
296  G4ReactionProductVector * products=0;
297  G4int interactionCounter = 0,collisionLoopMaxCount;
298  do
299  {
300  // reset status that could be changed in previous loop event
301  theCollisionMgr->ClearAndDestroy();
302 
303  if(products != 0)
304  { // free memory from previous loop event
305  ClearAndDestroy(products);
306  delete products;
307  products=0;
308  }
309 
310  G4int massNumber=aNucleus.GetA_asInt();
311  the3DNucleus->Init(massNumber, aNucleus.GetZ_asInt());
312  thePropagator->Init(the3DNucleus);
313  G4KineticTrack * kt;
314  collisionLoopMaxCount = 200;
315  do // sample impact parameter until collisions are found
316  {
317  theCurrentTime=0;
318  G4double radius = the3DNucleus->GetOuterRadius()+3*fermi;
319  initialPosition=GetSpherePoint(1.1*radius, initial4Momentum); // get random position
320  kt = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
321  kt->SetState(G4KineticTrack::outside);
322  // secondaries has been cleared by Propagate() in the previous loop event
323  secondaries= new G4KineticTrackVector;
324  secondaries->push_back(kt);
325  if(massNumber > 1) // 1H1 is special case
326  {
327  products = Propagate(secondaries, the3DNucleus);
328  } else {
329  products = Propagate1H1(secondaries,the3DNucleus);
330  }
331  // until we FIND a collision ... or give up
332  } while(! products && --collisionLoopMaxCount>0); /* Loop checking, 31.08.2015, G.Folger */
333 
334  if(++interactionCounter>99) break;
335  // ...until we find an ALLOWED collision ... or give up
336  } while(products && products->size() == 0); /* Loop checking, 31.08.2015, G.Folger */
337 
338  if(products && products->size()>0)
339  {
340  // G4cout << "BIC Applyyourself: number of products " << products->size() << G4endl;
341 
342  // Fill the G4ParticleChange * with products
343  theParticleChange.SetStatusChange(stopAndKill);
344  G4ReactionProductVector::iterator iter;
345 
346  for(iter = products->begin(); iter != products->end(); ++iter)
347  {
348  G4DynamicParticle * aNew =
349  new G4DynamicParticle((*iter)->GetDefinition(),
350  (*iter)->GetTotalEnergy(),
351  (*iter)->GetMomentum());
352  theParticleChange.AddSecondary(aNew);
353  }
354 
355  //DebugFinalEpConservation(aTrack, products);
356 
357 
358  } else { // no interaction, return primary
359  if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction void, return intial state ######### "<< G4endl;
360  theParticleChange.SetStatusChange(isAlive);
361  theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
362  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
363  }
364 
365  if ( products )
366  {
367  ClearAndDestroy(products);
368  delete products;
369  }
370 
371  delete the3DNucleus;
372  the3DNucleus = NULL;
373 
374  if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction ends ######### "<< G4endl;
375 
376  return &theParticleChange;
377 }
378 //----------------------------------------------------------------------------
379 G4ReactionProductVector * G4BinaryCascade::Propagate(
380  G4KineticTrackVector * secondaries, G4V3DNucleus * aNucleus)
381 //----------------------------------------------------------------------------
382 {
383  G4ping debug("debug_G4BinaryCascade");
384 #ifdef debug_BIC_Propagate
385  G4cout << "G4BinaryCascade Propagate starting -------------------------------------------------------" <<G4endl;
386 #endif
387 
388  the3DNucleus=aNucleus;
389  G4ReactionProductVector * products = new G4ReactionProductVector;
390  theOuterRadius = the3DNucleus->GetOuterRadius();
391  theCurrentTime=0;
392  theProjectile4Momentum=G4LorentzVector(0,0,0,0);
393  theMomentumTransfer=G4ThreeVector(0,0,0);
394  // build theSecondaryList, theProjectileList and theCapturedList
395  ClearAndDestroy(&theCapturedList);
396  ClearAndDestroy(&theSecondaryList);
397  theSecondaryList.clear();
398  ClearAndDestroy(&theFinalState);
399  std::vector<G4KineticTrack *>::iterator iter;
400  theCollisionMgr->ClearAndDestroy();
401 
402  theCutOnP=90*MeV;
403  if(the3DNucleus->GetMass()>30) theCutOnP = 70*MeV;
404  if(the3DNucleus->GetMass()>60) theCutOnP = 50*MeV;
405  if(the3DNucleus->GetMass()>120) theCutOnP = 45*MeV;
406 
407 
408  BuildTargetList();
409 
410 #ifdef debug_BIC_GetExcitationEnergy
411  G4cout << "ExcitationEnergy0 " << GetExcitationEnergy() << G4endl;
412 #endif
413 
414  thePropagator->Init(the3DNucleus);
415 
416  G4bool success = BuildLateParticleCollisions(secondaries);
417  if (! success ) // fails if no excitation energy left....
418  {
419  products=HighEnergyModelFSProducts(products, secondaries);
420  ClearAndDestroy(secondaries);
421  delete secondaries;
422 
423 #ifdef debug_G4BinaryCascade
424  G4cout << "G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" << G4endl;
425 #endif
426 
427  return products;
428  }
429  // check baryon and charge ...
430 
431  _CheckChargeAndBaryonNumber_("lateparticles");
432  _DebugEpConservation(" be4 findcollisions");
433 
434  // if called stand alone find first collisions
435  FindCollisions(&theSecondaryList);
436 
437 
438  if(theCollisionMgr->Entries() == 0 ) //late particles ALWAYS create Entries
439  {
440  //G4cout << " no collsions -> return 0" << G4endl;
441  delete products;
442 #ifdef debug_BIC_return
443  G4cout << "return @ begin2, no collisions "<< G4endl;
444 #endif
445  return 0;
446  }
447 
448  // end of initialization: do the job now
449  // loop until there are no more collisions
450 
451 
452  G4bool haveProducts = false;
453  G4int collisionCount=0;
454  G4int collisionLoopMaxCount=1000000;
455  while(theCollisionMgr->Entries() > 0 && currentZ && --collisionLoopMaxCount>0) /* Loop checking, 31.08.2015, G.Folger */
456  {
457  if(Absorb()) { // absorb secondaries, pions only
458  haveProducts = true;
459  }
460  if(Capture()) { // capture secondaries, nucleons only
461  haveProducts = true;
462  }
463 
464  // propagate to the next collision if any (collisions could have been deleted
465  // by previous absorption or capture)
466  if(theCollisionMgr->Entries() > 0)
467  {
468  G4CollisionInitialState *
469  nextCollision = theCollisionMgr->GetNextCollision();
470 #ifdef debug_BIC_Propagate_Collisions
471  G4cout << " NextCollision * , Time, curtime = " << nextCollision << " "
472  <<nextCollision->GetCollisionTime()<< " " <<
473  theCurrentTime<< G4endl;
474 #endif
475  if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) )
476  {
477  // Check if nextCollision is still valid, ie. particle did not leave nucleus
478  if (theCollisionMgr->GetNextCollision() != nextCollision )
479  {
480  nextCollision = 0;
481  }
482  }
483  //_DebugEpConservation("Stepped");
484 
485  if( nextCollision )
486  {
487  if (ApplyCollision(nextCollision))
488  {
489  //G4cerr << "ApplyCollision success " << G4endl;
490  haveProducts = true;
491  collisionCount++;
492  //_CheckChargeAndBaryonNumber_("ApplyCollision");
493  //_DebugEpConservation("ApplyCollision");
494 
495  } else {
496  //G4cerr << "ApplyCollision failure " << G4endl;
497  theCollisionMgr->RemoveCollision(nextCollision);
498  }
499  }
500  }
501  }
502 
503  //--------- end of on Collisions
504  //G4cout << "currentZ @ end loop " << currentZ << G4endl;
505  G4int nProtons(0);
506  for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
507  {
508  if ( (*iter)->GetDefinition() == G4Proton::Proton() ) ++nProtons;
509  }
510  if ( ! theTargetList.size() || ! nProtons ){
511  // nucleus completely destroyed, fill in ReactionProductVector
512  products = FillVoidNucleusProducts(products);
513 #ifdef debug_BIC_return
514  G4cout << "return @ Z=0 after collision loop "<< G4endl;
515  PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
516  G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
517  PrintKTVector(&theTargetList,std::string(" theTargetList"));
518  PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
519 
520  G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
521  G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
522  PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
523  G4cout << "returned products: " << products->size() << G4endl;
524  _CheckChargeAndBaryonNumber_("destroyed Nucleus");
525  _DebugEpConservation("destroyed Nucleus");
526 #endif
527 
528  return products;
529  }
530 
531  // No more collisions: absorb, capture and propagate the secondaries out of the nucleus
532  if(Absorb()) {
533  haveProducts = true;
534  // G4cout << "Absorb sucess " << G4endl;
535  }
536 
537  if(Capture()) {
538  haveProducts = true;
539  // G4cout << "Capture sucess " << G4endl;
540  }
541 
542  if(!haveProducts) // no collisions happened. Return an empty vector.
543  {
544 #ifdef debug_BIC_return
545  G4cout << "return 3, no products "<< G4endl;
546 #endif
547  return products;
548  }
549 
550 
551 #ifdef debug_BIC_Propagate
552  G4cout << " Momentum transfer to Nucleus " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
553  G4cout << " Stepping particles out...... " << G4endl;
554 #endif
555 
556  StepParticlesOut();
557  _DebugEpConservation("stepped out");
558 
559 
560  if ( theSecondaryList.size() > 0 )
561  {
562 #ifdef debug_G4BinaryCascade
563  G4cerr << "G4BinaryCascade: Warning, have active particles at end" << G4endl;
564  PrintKTVector(&theSecondaryList, "active particles @ end added to theFinalState");
565 #endif
566  // add left secondaries to FinalSate
567  for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
568  {
569  theFinalState.push_back(*iter);
570  }
571  theSecondaryList.clear();
572 
573  }
574  while ( theCollisionMgr->Entries() > 0 ) /* Loop checking, 31.08.2015, G.Folger */
575  {
576 #ifdef debug_G4BinaryCascade
577  G4cerr << " Warning: remove left over collision(s) " << G4endl;
578 #endif
579  theCollisionMgr->RemoveCollision(theCollisionMgr->GetNextCollision());
580  }
581 
582 #ifdef debug_BIC_Propagate_Excitation
583 
584  PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
585  G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
586  // PrintKTVector(&theTargetList,std::string(" theTargetList"));
587  PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
588 
589  G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
590  G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
591  PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
592 #endif
593 
594  //
595 
596 
597  G4double ExcitationEnergy=GetExcitationEnergy();
598 
599 #ifdef debug_BIC_Propagate_finals
600  PrintKTVector(&theFinalState,std::string(" FinalState be4 corr"));
601  G4cout << " Excitation Energy prefinal, #collisions:, out, captured "
602  << ExcitationEnergy << " "
603  << collisionCount << " "
604  << theFinalState.size() << " "
605  << theCapturedList.size()<<G4endl;
606 #endif
607 
608  if (ExcitationEnergy < 0 )
609  {
610  G4int maxtry=5, ntry=0;
611  do {
612  CorrectFinalPandE();
613  ExcitationEnergy=GetExcitationEnergy();
614  } while ( ++ntry < maxtry && ExcitationEnergy < 0 ); /* Loop checking, 31.08.2015, G.Folger */
615  }
616  _DebugEpConservation("corrected");
617 
618 #ifdef debug_BIC_Propagate_finals
619  PrintKTVector(&theFinalState,std::string(" FinalState corrected"));
620  G4cout << " Excitation Energy final, #collisions:, out, captured "
621  << ExcitationEnergy << " "
622  << collisionCount << " "
623  << theFinalState.size() << " "
624  << theCapturedList.size()<<G4endl;
625 #endif
626 
627 
628  if ( ExcitationEnergy < 0. )
629  {
630  #ifdef debug_G4BinaryCascade
631  G4cerr << "G4BinaryCascade-Warning: negative excitation energy ";
632  G4cerr <<ExcitationEnergy<<G4endl;
633  PrintKTVector(&theFinalState,std::string("FinalState"));
634  PrintKTVector(&theCapturedList,std::string("captured"));
635  G4cout << "negative ExE:Final 4Momentum .mag: " << GetFinal4Momentum()
636  << " "<< GetFinal4Momentum().mag()<< G4endl
637  << "negative ExE:FinalNucleusMom .mag: " << GetFinalNucleusMomentum()
638  << " "<< GetFinalNucleusMomentum().mag()<< G4endl;
639  #endif
640  #ifdef debug_BIC_return
641  G4cout << " negative Excitation E return empty products " << products << G4endl;
642  G4cout << "return 4, excit < 0 "<< G4endl;
643  #endif
644 
645  ClearAndDestroy(products);
646  return products; // return empty products- FIXME
647  }
648 
649  G4ReactionProductVector * precompoundProducts=DeExcite();
650 
651 
652  G4DecayKineticTracks decay(&theFinalState);
653  _DebugEpConservation("decayed");
654 
655  products= ProductsAddFinalState(products, theFinalState);
656 
657  products= ProductsAddPrecompound(products, precompoundProducts);
658 
659 // products=ProductsAddFakeGamma(products);
660 
661 
662  thePrimaryEscape = true;
663 
664  #ifdef debug_BIC_return
665  G4cout << "BIC: return @end, all ok "<< G4endl;
666  //G4cout << " return products " << products << G4endl;
667  #endif
668 
669  return products;
670 }
671 
672 //----------------------------------------------------------------------------
673 G4double G4BinaryCascade::GetExcitationEnergy()
674 //----------------------------------------------------------------------------
675 {
676 
677  // get A and Z for the residual nucleus
678 #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy)
679  G4int finalA = theTargetList.size()+theCapturedList.size();
680  G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
681  if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 )
682  {
683  G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} "
684  << "("<< currentA << "," << finalA << ") ("<< currentZ << "," << finalZ << ")" << G4endl;
685  }
686 
687 #endif
688 
689  G4double excitationE(0);
690  G4double nucleusMass(0);
691  if(currentZ>.5)
692  {
693  nucleusMass = GetIonMass(currentZ,currentA);
694  }
695  else if (currentZ==0 ) // Uzhi && currentA==1 ) // Uzhi
696  { // Uzhi
697  if(currentA == 1) {nucleusMass = G4Neutron::Neutron()->GetPDGMass();}// Uzhi
698  else {nucleusMass = GetFinalNucleusMomentum().mag() // Uzhi
699  - 3.*MeV*currentA;} // Uzhi
700  } // Uzhi
701  else
702  {
703 #ifdef debug_G4BinaryCascade
704  G4cout << "G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=("
705  << currentA << "," << currentZ << ")" << G4endl;
706 #endif
707  return 0;
708  }
709 
710 #ifdef debug_BIC_GetExcitationEnergy
711  G4ping debug("debug_ExcitationEnergy");
712  debug.push_back("====> current A, Z");
713  debug.push_back(currentZ);
714  debug.push_back(currentA);
715  debug.push_back("====> final A, Z");
716  debug.push_back(finalZ);
717  debug.push_back(finalA);
718  debug.push_back(nucleusMass);
719  debug.push_back(GetFinalNucleusMomentum().mag());
720  debug.dump();
721  // PrintKTVector(&theTargetList, std::string(" current target list info"));
722  //PrintKTVector(&theCapturedList, std::string(" current captured list info"));
723 #endif
724 
725  excitationE = GetFinalNucleusMomentum().mag() - nucleusMass;
726 
727  //G4double exE2 = GetFinal4Momentum().mag() - nucleusMass;
728 
729  //G4cout << "old/new excitE " << excitationE << " / "<< exE2 << G4endl;
730 
731 #ifdef debug_BIC_GetExcitationEnergy
732  // ------ debug
733  if ( excitationE < 0 )
734  {
735  G4cout << "negative ExE final Ion mass " <<nucleusMass<< G4endl;
736  G4LorentzVector Nucl_mom=GetFinalNucleusMomentum();
737  if(finalZ>.5) G4cout << " Final nuclmom/mass " << Nucl_mom << " " << Nucl_mom.mag()
738  << " (A,Z)=("<< finalA <<","<<finalZ <<")"
739  << " mass " << nucleusMass << " "
740  << " excitE " << excitationE << G4endl;
741 
742 
743  G4int A = the3DNucleus->GetMassNumber();
744  G4int Z = the3DNucleus->GetCharge();
745  G4double initialExc(0);
746  if(Z>.5)
747  {
748  initialExc = theInitial4Mom.mag()- GetIonMass(Z, A);
749  G4cout << "GetExcitationEnergy: Initial nucleus A Z " << A << " " << Z << " " << initialExc << G4endl;
750  }
751  }
752 
753 #endif
754 
755  return excitationE;
756 }
757 
758 
759 //----------------------------------------------------------------------------
760 //
761 // P R I V A T E M E M B E R F U N C T I O N S
762 //
763 //----------------------------------------------------------------------------
764 
765 //----------------------------------------------------------------------------
766 void G4BinaryCascade::BuildTargetList()
767 //----------------------------------------------------------------------------
768 {
769 
770  if(!the3DNucleus->StartLoop())
771  {
772  // G4cerr << "G4BinaryCascade::BuildTargetList(): StartLoop() error!"
773  // << G4endl;
774  return;
775  }
776 
777  ClearAndDestroy(&theTargetList); // clear theTargetList before rebuilding
778 
779  G4Nucleon * nucleon;
780  const G4ParticleDefinition * definition;
781  G4ThreeVector pos;
782  G4LorentzVector mom;
783  // if there are nucleon hit by higher energy models, then SUM(momenta) != 0
784  initialZ=the3DNucleus->GetCharge();
785  initialA=the3DNucleus->GetMassNumber();
786  initial_nuclear_mass=GetIonMass(initialZ,initialA);
787  theInitial4Mom = G4LorentzVector(0,0,0,initial_nuclear_mass);
788  currentA=0;
789  currentZ=0;
790  while((nucleon = the3DNucleus->GetNextNucleon()) != NULL) /* Loop checking, 31.08.2015, G.Folger */
791  {
792  // check if nucleon is hit by higher energy model.
793  if ( ! nucleon->AreYouHit() )
794  {
795  definition = nucleon->GetDefinition();
796  pos = nucleon->GetPosition();
797  mom = nucleon->GetMomentum();
798  // G4cout << "Nucleus " << pos.mag()/fermi << " " << mom.e() << G4endl;
799  //theInitial4Mom += mom;
800  // the potential inside the nucleus is taken into account, and nucleons are on mass shell.
801  mom.setE( std::sqrt( mom.vect().mag2() + sqr(definition->GetPDGMass()) ) );
802  G4KineticTrack * kt = new G4KineticTrack(definition, 0., pos, mom);
803  kt->SetState(G4KineticTrack::inside);
804  kt->SetNucleon(nucleon);
805  theTargetList.push_back(kt);
806  ++currentA;
807  if (definition->GetPDGCharge() > .5 ) ++currentZ;
808  }
809 
810 #ifdef debug_BIC_BuildTargetList
811  else { G4cout << "nucleon is hit" << nucleon << G4endl;}
812 #endif
813  }
814  massInNucleus = 0;
815  if(currentZ>.5)
816  {
817  massInNucleus = GetIonMass(currentZ,currentA);
818  } else if (currentZ==0 && currentA>=1 )
819  {
820  massInNucleus = currentA * G4Neutron::Neutron()->GetPDGMass();
821  } else
822  {
823  G4cerr << "G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=("
824  << currentA << "," << currentZ << ")" << G4endl;
825  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::BuildTargetList()");
826  }
827  currentInitialEnergy= theInitial4Mom.e() + theProjectile4Momentum.e();
828 
829 #ifdef debug_BIC_BuildTargetList
830  G4cout << "G4BinaryCascade::BuildTargetList(): nucleus (A,Z)=("
831  << currentA << "," << currentZ << ") mass: " << massInNucleus <<
832  ", theInitial4Mom " << theInitial4Mom <<
833  ", currentInitialEnergy " << currentInitialEnergy << G4endl;
834 #endif
835 
836 }
837 
838 //----------------------------------------------------------------------------
839 G4bool G4BinaryCascade::BuildLateParticleCollisions(G4KineticTrackVector * secondaries)
840 //----------------------------------------------------------------------------
841 {
842  G4bool success(false);
843  std::vector<G4KineticTrack *>::iterator iter;
844 
845  lateA=lateZ=0;
846  projectileA=projectileZ=0;
847 
848  G4double StartingTime=DBL_MAX; // Search for minimal formation time
849  for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
850  {
851  if((*iter)->GetFormationTime() < StartingTime)
852  StartingTime = (*iter)->GetFormationTime();
853  }
854 
855  //PrintKTVector(secondaries, "initial late particles ");
856  G4LorentzVector lateParticles4Momentum(0,0,0,0);
857  for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
858  {
859  // G4cout << " Formation time : " << (*iter)->GetDefinition()->GetParticleName() << " "
860  // << (*iter)->GetFormationTime() << G4endl;
861  G4double FormTime = (*iter)->GetFormationTime() - StartingTime;
862  (*iter)->SetFormationTime(FormTime);
863  if( (*iter)->GetState() == G4KineticTrack::undefined ) // particles from high energy generator
864  {
865  FindLateParticleCollision(*iter);
866  lateParticles4Momentum += (*iter)->GetTrackingMomentum();
867  lateA += (*iter)->GetDefinition()->GetBaryonNumber();
868  lateZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
869  //PrintKTVector(*iter, "late particle ");
870  } else
871  {
872  theSecondaryList.push_back(*iter);
873  //PrintKTVector(*iter, "incoming particle ");
874  theProjectile4Momentum += (*iter)->GetTrackingMomentum();
875  projectileA += (*iter)->GetDefinition()->GetBaryonNumber();
876  projectileZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
877 #ifdef debug_BIC_Propagate
878  G4cout << " Adding initial secondary " << *iter
879  << " time" << (*iter)->GetFormationTime()
880  << ", state " << (*iter)->GetState() << G4endl;
881 #endif
882  }
883  }
884  //theCollisionMgr->Print();
885  const G4HadProjectile * primary = GetPrimaryProjectile(); // check for primary from TheoHE model
886 
887  if (primary){
888  G4LorentzVector mom=primary->Get4Momentum();
889  theProjectile4Momentum += mom;
890  projectileA = primary->GetDefinition()->GetBaryonNumber();
891  projectileZ = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
892  // now check if "excitation" energy left by TheoHE model
893  G4double excitation= theProjectile4Momentum.e() + initial_nuclear_mass - lateParticles4Momentum.e() - massInNucleus;
894 #ifdef debug_BIC_GetExcitationEnergy
895  G4cout << "BIC: Proj.e, nucl initial, nucl final, lateParticles"
896  << theProjectile4Momentum << ", "
897  << initial_nuclear_mass<< ", " << massInNucleus << ", "
898  << lateParticles4Momentum << G4endl;
899  G4cout << "BIC: Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl;
900 #endif
901  success = excitation > 0;
902 #ifdef debug_G4BinaryCascade
903  if ( ! success ) {
904  G4cout << "G4BinaryCascade::BuildLateParticleCollisions(): Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl;
905  //PrintKTVector(secondaries);
906  }
907 #endif
908  } else {
909  // no primary from HE model -> cascade
910  success=true;
911  }
912 
913  if (success) {
914  secondaries->clear(); // Don't leave "G4KineticTrack *"s in two vectors
915  delete secondaries;
916  }
917  return success;
918 }
919 
920 //----------------------------------------------------------------------------
921 G4ReactionProductVector * G4BinaryCascade::DeExcite()
922 //----------------------------------------------------------------------------
923 {
924  // find a fragment and call the precompound model.
925  G4Fragment * fragment = 0;
926  G4ReactionProductVector * precompoundProducts = 0;
927 
928  G4LorentzVector pFragment(0);
929  // G4cout << " final4mon " << GetFinal4Momentum() /MeV << G4endl;
930 
931  // if ( ExcitationEnergy >= 0 ) // closed by Uzhi
932  // { // closed by Uzhi
933  fragment = FindFragments();
934 
935 
936  if(fragment) // Uzhi
937  { // Uzhi
938  if(fragment->GetA_asInt() >1) // Uzhi
939  {
940  pFragment=fragment->GetMomentum();
941  // G4cout << " going to preco with fragment 4 mom " << pFragment << G4endl;
942  if (theDeExcitation) // pre-compound
943  {
944  precompoundProducts= theDeExcitation->DeExcite(*fragment);
945  }
946  else if (theExcitationHandler) // de-excitation
947  {
948  precompoundProducts=theExcitationHandler->BreakItUp(*fragment);
949  }
950 
951  } else
952  { // fragment->GetA_asInt() <= 1, so a single proton, as a fragment must have Z>0
953  if (theTargetList.size() + theCapturedList.size() > 1 ) {
954  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde:: Invalid Fragment");
955  }
956 
957  std::vector<G4KineticTrack *>::iterator i;
958  if ( theTargetList.size() == 1 ) {i=theTargetList.begin();}
959  if ( theCapturedList.size() == 1 ) {i=theCapturedList.begin();} // Uzhi
960  G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition());
961  aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass());
962  aNew->SetMomentum(G4ThreeVector(0));// see boost for preCompoundProducts below..
963  precompoundProducts = new G4ReactionProductVector();
964  precompoundProducts->push_back(aNew);
965  } // End of fragment->GetA() < 1.5
966  delete fragment;
967  fragment=0;
968 
969  } else // End of if(fragment)
970  { // No fragment, can be neutrons only // Uzhi
971 
972  precompoundProducts = DecayVoidNucleus();
973  }
974  #ifdef debug_BIC_DeexcitationProducts
975 
976  G4LorentzVector fragment_momentum=GetFinalNucleusMomentum();
977  G4LorentzVector Preco_momentum;
978  if ( precompoundProducts )
979  {
980  std::vector<G4ReactionProduct *>::iterator j;
981  for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
982  {
983  G4LorentzVector pProduct((*j)->GetMomentum(),(*j)->GetTotalEnergy());
984  Preco_momentum += pProduct;
985  }
986  }
987  G4cout << "finalNuclMom / sum preco products" << fragment_momentum << " / " << Preco_momentum
988  << " delta E "<< fragment_momentum.e() - Preco_momentum.e() << G4endl;
989 
990  #endif
991 
992  return precompoundProducts;
993 }
994 
995 //----------------------------------------------------------------------------
996 G4ReactionProductVector * G4BinaryCascade::DecayVoidNucleus()
997 //----------------------------------------------------------------------------
998 {
999  G4ReactionProductVector * result=0;
1000  if ( (theTargetList.size()+theCapturedList.size()) > 0 )
1001  {
1002  result = new G4ReactionProductVector;
1003  std::vector<G4KineticTrack *>::iterator aNuc;
1004  G4LorentzVector aVec;
1005  std::vector<G4double> masses;
1006  G4double sumMass(0);
1007 
1008  if ( theTargetList.size() != 0) // Uzhi
1009  {
1010  for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
1011  {
1012  G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
1013  masses.push_back(mass);
1014  sumMass += mass;
1015  }
1016  } // Uzhi
1017 
1018  if ( theCapturedList.size() != 0) // Uzhi
1019  { // Uzhi
1020  for(aNuc = theCapturedList.begin(); // Uzhi
1021  aNuc != theCapturedList.end(); aNuc++) // Uzhi
1022  { // Uzhi
1023  G4double mass=(*aNuc)->GetDefinition()->GetPDGMass(); // Uzhi
1024  masses.push_back(mass); // Uzhi
1025  sumMass += mass; // Uzhi
1026  }
1027  }
1028 
1029  G4LorentzVector finalP=GetFinal4Momentum();
1030  G4FermiPhaseSpaceDecay decay;
1031  // G4cout << " some neutrons? " << masses.size() <<" " ;
1032  // G4cout<< theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl;
1033 
1034  G4double eCMS=finalP.mag();
1035  if ( eCMS < sumMass ) // @@GF --- Cheat!!
1036  {
1037  eCMS=sumMass + 2*MeV*masses.size();
1038  finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS)));
1039  }
1040 
1041  precompoundLorentzboost.set(finalP.boostVector());
1042  std::vector<G4LorentzVector*> * momenta=decay.Decay(eCMS,masses);
1043  std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
1044 
1045 
1046  if ( theTargetList.size() != 0)
1047  {
1048  for ( aNuc=theTargetList.begin();
1049  (aNuc != theTargetList.end()) && (aMom!=momenta->end());
1050  aNuc++, aMom++ )
1051  {
1052  G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition());
1053  aNew->SetTotalEnergy((*aMom)->e());
1054  aNew->SetMomentum((*aMom)->vect());
1055  result->push_back(aNew);
1056 
1057  delete *aMom;
1058  }
1059  }
1060 
1061  if ( theCapturedList.size() != 0) // Uzhi
1062  { // Uzhi
1063  for ( aNuc=theCapturedList.begin(); // Uzhi
1064  (aNuc != theCapturedList.end()) && (aMom!=momenta->end()); // Uzhi
1065  aNuc++, aMom++ ) // Uzhi
1066  { // Uzhi
1067  G4ReactionProduct * aNew = new G4ReactionProduct( // Uzhi
1068  (*aNuc)->GetDefinition()); // Uzhi
1069  aNew->SetTotalEnergy((*aMom)->e()); // Uzhi
1070  aNew->SetMomentum((*aMom)->vect()); // Uzhi
1071  result->push_back(aNew); // Uzhi
1072  delete *aMom; // Uzhi
1073  } // Uzhi
1074  } // Uzhi
1075 
1076  delete momenta;
1077  }
1078  return result;
1079 } // End if(!fragment)
1080 
1081 //----------------------------------------------------------------------------
1082 G4ReactionProductVector * G4BinaryCascade::ProductsAddFinalState(G4ReactionProductVector * products, G4KineticTrackVector & fs)
1083 //----------------------------------------------------------------------------
1084 {
1085 // fill in products the outgoing particles
1086  size_t i(0);
1087 #ifdef debug_BIC_Propagate_finals
1088  G4LorentzVector mom_fs;
1089 #endif
1090  for(i = 0; i< fs.size(); i++)
1091  {
1092  G4KineticTrack * kt = fs[i];
1093  G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
1094  aNew->SetMomentum(kt->Get4Momentum().vect());
1095  aNew->SetTotalEnergy(kt->Get4Momentum().e());
1096  aNew->SetNewlyAdded(kt->IsParticipant());
1097  products->push_back(aNew);
1098 
1099 #ifdef debug_BIC_Propagate_finals
1100  mom_fs += kt->Get4Momentum();
1101  G4cout <<kt->GetDefinition()->GetParticleName();
1102  G4cout << " Particle Ekin " << aNew->GetKineticEnergy();
1103  G4cout << ", is " << (kt->GetDefinition()->GetPDGStable() ? "stable" :
1104  (kt->GetDefinition()->IsShortLived() ? "short lived " : "non stable")) ;
1105  G4cout << G4endl;
1106 #endif
1107 
1108  }
1109 #ifdef debug_BIC_Propagate_finals
1110  G4cout << " Final state momentum " << mom_fs << G4endl;
1111 #endif
1112 
1113  return products;
1114 }
1115 //----------------------------------------------------------------------------
1116 G4ReactionProductVector * G4BinaryCascade::ProductsAddPrecompound(G4ReactionProductVector * products, G4ReactionProductVector * precompoundProducts)
1117 //----------------------------------------------------------------------------
1118 {
1119  G4LorentzVector pSumPreco(0), pPreco(0);
1120 
1121  if ( precompoundProducts )
1122  {
1123  std::vector<G4ReactionProduct *>::iterator j;
1124  for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
1125  {
1126  // boost back to system of moving nucleus
1127  G4LorentzVector pProduct((*j)->GetMomentum(),(*j)->GetTotalEnergy());
1128  pPreco+= pProduct;
1129 #ifdef debug_BIC_Propagate_finals
1130  G4cout << "BIC: pProduct be4 boost " <<pProduct << G4endl;
1131 #endif
1132  pProduct *= precompoundLorentzboost;
1133 #ifdef debug_BIC_Propagate_finals
1134  G4cout << "BIC: pProduct aft boost " <<pProduct << G4endl;
1135 #endif
1136  pSumPreco += pProduct;
1137  (*j)->SetTotalEnergy(pProduct.e());
1138  (*j)->SetMomentum(pProduct.vect());
1139  (*j)->SetNewlyAdded(true);
1140  products->push_back(*j);
1141  }
1142  // G4cout << " unboosted preco result mom " << pPreco / MeV << " ..- fragmentMom " << (pPreco - pFragment)/MeV<< G4endl;
1143  // G4cout << " preco result mom " << pSumPreco / MeV << " ..-file4Mom " << (pSumPreco - GetFinal4Momentum())/MeV<< G4endl;
1144  precompoundProducts->clear();
1145  delete precompoundProducts;
1146  }
1147  return products;
1148 }
1149 //----------------------------------------------------------------------------
1150 void G4BinaryCascade::FindCollisions(G4KineticTrackVector * secondaries)
1151 //----------------------------------------------------------------------------
1152 {
1153  for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin();
1154  i != secondaries->end(); ++i)
1155  {
1156  for(std::vector<G4BCAction *>::iterator j = theImR.begin();
1157  j!=theImR.end(); j++)
1158  {
1159  // G4cout << "G4BinaryCascade::FindCollisions: at action " << *j << G4endl;
1160  const std::vector<G4CollisionInitialState *> & aCandList
1161  = (*j)->GetCollisions(*i, theTargetList, theCurrentTime);
1162  for(size_t count=0; count<aCandList.size(); count++)
1163  {
1164  theCollisionMgr->AddCollision(aCandList[count]);
1165  //4cout << "====================== New Collision ================="<<G4endl;
1166  //theCollisionMgr->Print();
1167  }
1168  }
1169  }
1170 }
1171 
1172 
1173 //----------------------------------------------------------------------------
1174 void G4BinaryCascade::FindDecayCollision(G4KineticTrack * secondary)
1175 //----------------------------------------------------------------------------
1176 {
1177  const std::vector<G4CollisionInitialState *> & aCandList
1178  = theDecay->GetCollisions(secondary, theTargetList, theCurrentTime);
1179  for(size_t count=0; count<aCandList.size(); count++)
1180  {
1181  theCollisionMgr->AddCollision(aCandList[count]);
1182  }
1183 }
1184 
1185 //----------------------------------------------------------------------------
1186 void G4BinaryCascade::FindLateParticleCollision(G4KineticTrack * secondary)
1187 //----------------------------------------------------------------------------
1188 {
1189 
1190  G4double tin=0., tout=0.;
1191  if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(secondary,tin,tout))
1192  {
1193  if ( tin > 0 )
1194  {
1195  secondary->SetState(G4KineticTrack::outside);
1196  } else if ( tout > 0 )
1197  {
1198  secondary->SetState(G4KineticTrack::inside);
1199  } else {
1200  //G4cout << "G4BC set miss , tin, tout " << tin << " , " << tout <<G4endl;
1201  secondary->SetState(G4KineticTrack::miss_nucleus);
1202  }
1203  } else {
1204  secondary->SetState(G4KineticTrack::miss_nucleus);
1205  //G4cout << "G4BC set miss ,no intersect tin, tout " << tin << " , " << tout <<G4endl;
1206  }
1207 
1208 
1209 #ifdef debug_BIC_FindCollision
1210  G4cout << "FindLateP Particle, 4-mom, times newState "
1211  << secondary->GetDefinition()->GetParticleName() << " "
1212  << secondary->Get4Momentum()
1213  << " times " << tin << " " << tout << " "
1214  << secondary->GetState() << G4endl;
1215 #endif
1216 
1217  const std::vector<G4CollisionInitialState *> & aCandList
1218  = theLateParticle->GetCollisions(secondary, theTargetList, theCurrentTime);
1219  for(size_t count=0; count<aCandList.size(); count++)
1220  {
1221 #ifdef debug_BIC_FindCollision
1222  G4cout << " Adding a late Col : " << aCandList[count] << G4endl;
1223 #endif
1224  theCollisionMgr->AddCollision(aCandList[count]);
1225  }
1226 }
1227 
1228 
1229 //----------------------------------------------------------------------------
1230 G4bool G4BinaryCascade::ApplyCollision(G4CollisionInitialState * collision)
1231 //----------------------------------------------------------------------------
1232 {
1233  G4KineticTrack * primary = collision->GetPrimary();
1234 
1235 #ifdef debug_BIC_ApplyCollision
1236  G4cerr << "G4BinaryCascade::ApplyCollision start"<<G4endl;
1237  theCollisionMgr->Print();
1238  G4cout << "ApplyCollisions : projte 4mom " << primary->GetTrackingMomentum()<< G4endl;
1239 #endif
1240 
1241  G4KineticTrackVector target_collection=collision->GetTargetCollection();
1242  G4bool haveTarget=target_collection.size()>0;
1243  if( haveTarget && (primary->GetState() != G4KineticTrack::inside) )
1244  {
1245 #ifdef debug_G4BinaryCascade
1246  G4cout << "G4BinaryCasacde::ApplyCollision(): StateError " << primary << G4endl;
1247  PrintKTVector(primary,std::string("primay- ..."));
1248  PrintKTVector(&target_collection,std::string("... targets"));
1249  collision->Print();
1250  G4cout << G4endl;
1251  theCollisionMgr->Print();
1252  //*GF* throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
1253 #endif
1254  return false;
1255 // } else {
1256 // G4cout << "G4BinaryCasacde::ApplyCollision(): decay " << G4endl;
1257 // PrintKTVector(primary,std::string("primay- ..."));
1258 // G4double tin=0., tout=0.;
1259 // if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(primary,tin,tout))
1260 // {
1261 // G4cout << "tin tout: " << tin << " " << tout << G4endl;
1262 // }
1263 
1264  }
1265 
1266  G4LorentzVector mom4Primary=primary->Get4Momentum();
1267 
1268  G4int initialBaryon(0);
1269  G4int initialCharge(0);
1270  if (primary->GetState() == G4KineticTrack::inside)
1271  {
1272  initialBaryon = primary->GetDefinition()->GetBaryonNumber();
1273  initialCharge = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
1274  }
1275 
1276  // for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
1277  G4double initial_Efermi=CorrectShortlivedPrimaryForFermi(primary,target_collection);
1278  //****************************************
1279 
1280 
1281  G4KineticTrackVector * products = collision->GetFinalState();
1282 
1283 #ifdef debug_BIC_ApplyCollision
1284  DebugApplyCollisionFail(collision, products);
1285 #endif
1286 
1287  // reset primary to initial state, in case there is a veto...
1288  primary->Set4Momentum(mom4Primary);
1289 
1290  G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
1291  G4bool decayCollision= (!haveTarget) && products && products->size() > 1;
1292  G4bool Success(true);
1293 
1294 
1295 #ifdef debug_G4BinaryCascade
1296  G4int lateBaryon(0), lateCharge(0);
1297 #endif
1298 
1299  if ( lateParticleCollision )
1300  { // for late particles, reset charges
1301  //G4cout << "lateP, initial B C state " << initialBaryon << " "
1302  // << initialCharge<< " " << primary->GetState() << " "<< primary->GetDefinition()->GetParticleName()<< G4endl;
1303 #ifdef debug_G4BinaryCascade
1304  lateBaryon = initialBaryon;
1305  lateCharge = initialCharge;
1306 #endif
1307  initialBaryon=initialCharge=0;
1308  lateA -= primary->GetDefinition()->GetBaryonNumber();
1309  lateZ -= G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
1310  }
1311 
1312  initialBaryon += collision->GetTargetBaryonNumber();
1313  initialCharge += G4lrint(collision->GetTargetCharge());
1314  if (!lateParticleCollision)
1315  {
1316  if( !products || products->size()==0 || !CheckPauliPrinciple(products) )
1317  {
1318 #ifdef debug_BIC_ApplyCollision
1319  if (products) G4cout << " ======Failed Pauli =====" << G4endl;
1320  G4cerr << "G4BinaryCascade::ApplyCollision blocked"<<G4endl;
1321 #endif
1322  Success=false;
1323  }
1324 
1325 
1326 
1327  if (Success && primary->GetState() == G4KineticTrack::inside ) { // if the primary was outside, nothing to correct
1328  if (! CorrectShortlivedFinalsForFermi(products, initial_Efermi)){
1329  Success=false;
1330  }
1331  }
1332  }
1333 
1334 #ifdef debug_BIC_ApplyCollision
1335  DebugApplyCollision(collision, products);
1336 #endif
1337 
1338  if ( ! Success ){
1339  if (products) ClearAndDestroy(products);
1340  if ( decayCollision ) FindDecayCollision(primary); // for decay, sample new decay
1341  delete products;
1342  products=0;
1343  return false;
1344  }
1345 
1346  G4int finalBaryon(0);
1347  G4int finalCharge(0);
1348  G4KineticTrackVector toFinalState;
1349  for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1350  {
1351  if ( ! lateParticleCollision )
1352  {
1353  (*i)->SetState(primary->GetState()); // decay may be anywhere!
1354  if ( (*i)->GetState() == G4KineticTrack::inside ){
1355  finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1356  finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
1357  } else {
1358  G4double tin=0., tout=0.;
1359  if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout) &&
1360  tin < 0 && tout > 0 )
1361  {
1362  PrintKTVector((*i),"particle inside marked not-inside");
1363  G4cout << "tin tout: " << tin << " " << tout << G4endl;
1364  }
1365  }
1366  } else {
1367  G4double tin=0., tout=0.;
1368  if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout))
1369  {
1370  //G4cout << "tin tout: " << tin << " " << tout << G4endl;
1371  if ( tin > 0 )
1372  {
1373  (*i)->SetState(G4KineticTrack::outside);
1374  }
1375  else if ( tout > 0 )
1376  {
1377  (*i)->SetState(G4KineticTrack::inside);
1378  finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
1379  finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
1380  }
1381  else
1382  {
1383  (*i)->SetState(G4KineticTrack::gone_out);
1384  toFinalState.push_back((*i));
1385  }
1386  } else
1387  {
1388  (*i)->SetState(G4KineticTrack::miss_nucleus);
1389  //G4cout << " G4BC - miss -late Part- no intersection found " << G4endl;
1390  toFinalState.push_back((*i));
1391  }
1392 
1393  }
1394  }
1395  if(!toFinalState.empty())
1396  {
1397  theFinalState.insert(theFinalState.end(),
1398  toFinalState.begin(),toFinalState.end());
1399  std::vector<G4KineticTrack *>::iterator iter1, iter2;
1400  for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
1401  ++iter1)
1402  {
1403  iter2 = std::find(products->begin(), products->end(),
1404  *iter1);
1405  if ( iter2 != products->end() ) products->erase(iter2);
1406  }
1407  theCollisionMgr->RemoveTracksCollisions(&toFinalState);
1408  }
1409 
1410  //G4cout << " currentA, Z be4: " << currentA << " " << currentZ << G4endl;
1411  currentA += finalBaryon-initialBaryon;
1412  currentZ += finalCharge-initialCharge;
1413  //G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl;
1414 
1415  G4KineticTrackVector oldSecondaries;
1416  oldSecondaries.push_back(primary);
1417  primary->Hit();
1418 
1419 #ifdef debug_G4BinaryCascade
1420  if ( (finalBaryon-initialBaryon-lateBaryon) != 0 || (finalCharge-initialCharge-lateCharge) != 0 )
1421  {
1422  G4cout << "G4BinaryCascade: Error in Balancing: " << G4endl;
1423  G4cout << "initial/final baryon number, initial/final Charge "
1424  << initialBaryon <<" "<< finalBaryon <<" "
1425  << initialCharge <<" "<< finalCharge <<" "
1426  << " in Collision type: "<< typeid(*collision->GetGenerator()).name()
1427  << ", with number of products: "<< products->size() <<G4endl;
1428  G4cout << G4endl<<"Initial condition are these:"<<G4endl;
1429  G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
1430  for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
1431  {
1432  G4cout << "targ: "
1433  <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
1434  }
1435  PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
1436  G4cout << G4endl<<G4endl;
1437  }
1438 #endif
1439 
1440  G4KineticTrackVector oldTarget = collision->GetTargetCollection();
1441  for(size_t ii=0; ii< oldTarget.size(); ii++)
1442  {
1443  oldTarget[ii]->Hit();
1444  }
1445 
1446  UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products);
1447  std::for_each(oldSecondaries.begin(), oldSecondaries.end(), Delete<G4KineticTrack>());
1448  std::for_each(oldTarget.begin(), oldTarget.end(), Delete<G4KineticTrack>());
1449 
1450  delete products;
1451  return true;
1452 }
1453 
1454 
1455 //----------------------------------------------------------------------------
1456 G4bool G4BinaryCascade::Absorb()
1457 //----------------------------------------------------------------------------
1458 {
1459  // Do it in two step: cannot change theSecondaryList inside the first loop.
1460  G4Absorber absorber(theCutOnPAbsorb);
1461 
1462  // Build the vector of G4KineticTracks that must be absorbed
1463  G4KineticTrackVector absorbList;
1464  std::vector<G4KineticTrack *>::iterator iter;
1465  // PrintKTVector(&theSecondaryList, " testing for Absorb" );
1466  for(iter = theSecondaryList.begin();
1467  iter != theSecondaryList.end(); ++iter)
1468  {
1469  G4KineticTrack * kt = *iter;
1470  if(kt->GetState() == G4KineticTrack::inside)// absorption happens only inside the nucleus
1471  {
1472  if(absorber.WillBeAbsorbed(*kt))
1473  {
1474  absorbList.push_back(kt);
1475  }
1476  }
1477  }
1478 
1479  if(absorbList.empty())
1480  return false;
1481 
1482  G4KineticTrackVector toDelete;
1483  for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
1484  {
1485  G4KineticTrack * kt = *iter;
1486  if(!absorber.FindAbsorbers(*kt, theTargetList))
1487  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1488 
1489  if(!absorber.FindProducts(*kt))
1490  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1491 
1492  G4KineticTrackVector * products = absorber.GetProducts();
1493  G4int maxLoopCount = 1000;
1494  while(!CheckPauliPrinciple(products) && --maxLoopCount>0) /* Loop checking, 31.08.2015, G.Folger */
1495  {
1496  ClearAndDestroy(products);
1497  if(!absorber.FindProducts(*kt))
1498  throw G4HadronicException(__FILE__, __LINE__,
1499  "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1500  }
1501  if ( --maxLoopCount < 0 ) throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
1502  // ------ debug
1503  // G4cerr << "Absorb CheckPauliPrinciple count= " << maxLoopCount << G4endl;
1504  // ------ end debug
1505  G4KineticTrackVector toRemove; // build a vector for UpdateTrack...
1506  toRemove.push_back(kt);
1507  toDelete.push_back(kt); // delete the track later
1508  G4KineticTrackVector * absorbers = absorber.GetAbsorbers();
1509  UpdateTracksAndCollisions(&toRemove, absorbers, products);
1510  ClearAndDestroy(absorbers);
1511  }
1512  ClearAndDestroy(&toDelete);
1513  return true;
1514 }
1515 
1516 
1517 
1518 // Capture all p and n with Energy < theCutOnP
1519 //----------------------------------------------------------------------------
1520 G4bool G4BinaryCascade::Capture(G4bool verbose)
1521 //----------------------------------------------------------------------------
1522 {
1523  G4KineticTrackVector captured;
1524  G4bool capture = false;
1525  std::vector<G4KineticTrack *>::iterator i;
1526 
1527  G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
1528 
1529  G4double capturedEnergy = 0;
1530  G4int particlesAboveCut=0;
1531  G4int particlesBelowCut=0;
1532  if ( verbose ) G4cout << " Capture: secondaries " << theSecondaryList.size() << G4endl;
1533  for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1534  {
1535  G4KineticTrack * kt = *i;
1536  if (verbose) G4cout << "Capture position, radius, state " <<kt->GetPosition().mag()<<" "<<theOuterRadius<<" "<<kt->GetState()<<G4endl;
1537  if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1538  {
1539  if((kt->GetDefinition() == G4Proton::Proton()) ||
1540  (kt->GetDefinition() == G4Neutron::Neutron()))
1541  {
1542  //GF cut on kinetic energy if(kt->Get4Momentum().vect().mag() >= theCutOnP)
1543  G4double field=RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
1544  -RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
1545  G4double energy= kt->Get4Momentum().e() - kt->GetActualMass() + field;
1546  if (verbose ) G4cout << "Capture: .e(), mass, field, energy" << kt->Get4Momentum().e() <<" "<<kt->GetActualMass()<<" "<<field<<" "<<energy<< G4endl;
1547  // if( energy < theCutOnP )
1548  // {
1549  capturedEnergy+=energy;
1550  ++particlesBelowCut;
1551  // } else
1552  // {
1553  // ++particlesAboveCut;
1554  // }
1555  }
1556  }
1557  }
1558  if (verbose) G4cout << "Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
1559  << particlesAboveCut << " " << particlesBelowCut << " " << capturedEnergy
1560  << " " << G4endl;
1561 // << " " << (particlesBelowCut>0) ? (capturedEnergy/particlesBelowCut) : (capturedEnergy) << " " << 0.2*theCutOnP << G4endl;
1562  // if(particlesAboveCut==0 && particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1563  if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
1564  {
1565  capture=true;
1566  for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1567  {
1568  G4KineticTrack * kt = *i;
1569  if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
1570  {
1571  if((kt->GetDefinition() == G4Proton::Proton()) ||
1572  (kt->GetDefinition() == G4Neutron::Neutron()))
1573  {
1574  captured.push_back(kt);
1575  kt->Hit(); //
1576  theCapturedList.push_back(kt);
1577  }
1578  }
1579  }
1580  UpdateTracksAndCollisions(&captured, NULL, NULL);
1581  }
1582 
1583  return capture;
1584 }
1585 
1586 
1587 //----------------------------------------------------------------------------
1588 G4bool G4BinaryCascade::CheckPauliPrinciple(G4KineticTrackVector * products)
1589 //----------------------------------------------------------------------------
1590 {
1591  G4int A = the3DNucleus->GetMassNumber();
1592  G4int Z = the3DNucleus->GetCharge();
1593 
1594  G4FermiMomentum fermiMom;
1595  fermiMom.Init(A, Z);
1596 
1597  const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
1598 
1599  G4KineticTrackVector::iterator i;
1600  const G4ParticleDefinition * definition;
1601 
1602  // ------ debug
1603  G4bool myflag = true;
1604  // ------ end debug
1605  // G4ThreeVector xpos(0);
1606  for(i = products->begin(); i != products->end(); ++i)
1607  {
1608  definition = (*i)->GetDefinition();
1609  if((definition == G4Proton::Proton()) ||
1610  (definition == G4Neutron::Neutron()))
1611  {
1612  G4ThreeVector pos = (*i)->GetPosition();
1613  G4double d = density->GetDensity(pos);
1614  // energy correspondiing to fermi momentum
1615  G4double eFermi = std::sqrt( sqr(fermiMom.GetFermiMomentum(d)) + (*i)->Get4Momentum().mag2() );
1616  if( definition == G4Proton::Proton() )
1617  {
1618  eFermi -= the3DNucleus->CoulombBarrier();
1619  }
1620  G4LorentzVector mom = (*i)->Get4Momentum();
1621  // ------ debug
1622  /*
1623  * G4cout << "p:[" << (1/MeV)*mom.x() << " " << (1/MeV)*mom.y() << " "
1624  * << (1/MeV)*mom.z() << "] |p3|:"
1625  * << (1/MeV)*mom.vect().mag() << " E:" << (1/MeV)*mom.t() << " ] m: "
1626  * << (1/MeV)*mom.mag() << " pos["
1627  * << (1/fermi)*pos.x() << " "<< (1/fermi)*pos.y() << " "
1628  * << (1/fermi)*pos.z() << "] |Dpos|: "
1629  * << (1/fermi)*(pos-xpos).mag() << " Pfermi:"
1630  * << (1/MeV)*p << G4endl;
1631  * xpos=pos;
1632  */
1633  // ------ end debug
1634  if(mom.e() < eFermi )
1635  {
1636  // ------ debug
1637  myflag = false;
1638  // ------ end debug
1639  // return false;
1640  }
1641  }
1642  }
1643 #ifdef debug_BIC_CheckPauli
1644  if ( myflag )
1645  {
1646  for(i = products->begin(); i != products->end(); ++i)
1647  {
1648  definition = (*i)->GetDefinition();
1649  if((definition == G4Proton::Proton()) ||
1650  (definition == G4Neutron::Neutron()))
1651  {
1652  G4ThreeVector pos = (*i)->GetPosition();
1653  G4double d = density->GetDensity(pos);
1654  G4double pFermi = fermiMom.GetFermiMomentum(d);
1655  G4LorentzVector mom = (*i)->Get4Momentum();
1656  G4double field =((G4RKPropagation*)thePropagator)->GetField(definition->GetPDGEncoding(),pos);
1657  if ( mom.e()-mom.mag()+field > 160*MeV )
1658  {
1659  G4cout << "momentum problem pFermi=" << pFermi
1660  << " mom, mom.m " << mom << " " << mom.mag()
1661  << " field " << field << G4endl;
1662  }
1663  }
1664  }
1665  }
1666 #endif
1667 
1668  return myflag;
1669 }
1670 
1671 //----------------------------------------------------------------------------
1672 void G4BinaryCascade::StepParticlesOut()
1673 //----------------------------------------------------------------------------
1674 {
1675  G4int counter=0;
1676  G4int countreset=0;
1677  //G4cout << " nucl. Radius " << radius << G4endl;
1678  // G4cerr <<"pre-while- theSecondaryList "<<G4endl;
1679  while( theSecondaryList.size() > 0 ) /* Loop checking, 31.08.2015, G.Folger */
1680  // if countreset reaches limit, there is a break from while, see below.
1681  {
1682  G4int nsec=0;
1683  G4double minTimeStep = 1.e-12*ns; // about 30*fermi/(0.1*c_light);1.e-12*ns
1684  // i.e. a big step
1685  std::vector<G4KineticTrack *>::iterator i;
1686  for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
1687  {
1688  G4KineticTrack * kt = *i;
1689  if( kt->GetState() == G4KineticTrack::inside )
1690  {
1691  nsec++;
1692  G4double tStep(0), tdummy(0);
1693  G4bool intersect =
1694  ((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep);
1695 #ifdef debug_BIC_StepParticlesOut
1696  G4cout << " minTimeStep, tStep Particle " <<minTimeStep << " " <<tStep
1697  << " " <<kt->GetDefinition()->GetParticleName()
1698  << " 4mom " << kt->GetTrackingMomentum()<<G4endl;
1699  if ( ! intersect );
1700  {
1701  PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
1702  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1703  }
1704 #endif
1705  if(intersect && tStep<minTimeStep && tStep> 0 )
1706  {
1707  minTimeStep = tStep;
1708  }
1709  } else if ( kt->GetState() != G4KineticTrack::outside ){
1710  PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
1711  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
1712  }
1713  }
1714  minTimeStep *= 1.2;
1715  // G4cerr << "CaptureCount = "<<counter<<" "<<nsec<<" "<<minTimeStep<<" "<<1*ns<<G4endl;
1716  G4double timeToCollision=DBL_MAX;
1717  G4CollisionInitialState * nextCollision=0;
1718  if(theCollisionMgr->Entries() > 0)
1719  {
1720  nextCollision = theCollisionMgr->GetNextCollision();
1721  timeToCollision = nextCollision->GetCollisionTime()-theCurrentTime;
1722  // G4cout << " NextCollision * , Time= " << nextCollision << " " <<timeToCollision<< G4endl;
1723  }
1724  if ( timeToCollision > minTimeStep )
1725  {
1726  DoTimeStep(minTimeStep);
1727  ++counter;
1728  } else
1729  {
1730  if (!DoTimeStep(timeToCollision) )
1731  {
1732  // Check if nextCollision is still valid, ie. partcile did not leave nucleus
1733  if (theCollisionMgr->GetNextCollision() != nextCollision )
1734  {
1735  nextCollision = 0;
1736  }
1737  }
1738  // G4cerr <<"post- DoTimeStep 3"<<G4endl;
1739 
1740  if(nextCollision)
1741  {
1742  if ( ApplyCollision(nextCollision))
1743  {
1744  // G4cout << "ApplyCollision sucess " << G4endl;
1745  } else
1746  {
1747  theCollisionMgr->RemoveCollision(nextCollision);
1748  }
1749  }
1750  }
1751 
1752  if(countreset>100)
1753  {
1754 #ifdef debug_G4BinaryCascade
1755  G4cerr << "G4BinaryCascade.cc: Warning - aborting looping particle(s)" << G4endl;
1756  PrintKTVector(&theSecondaryList," looping particles added to theFinalState");
1757 #endif
1758 
1759  // add left secondaries to FinalSate
1760  std::vector<G4KineticTrack *>::iterator iter;
1761  for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
1762  {
1763  theFinalState.push_back(*iter);
1764  }
1765  theSecondaryList.clear();
1766 
1767  break;
1768  }
1769 
1770  if(Absorb())
1771  {
1772  // haveProducts = true;
1773  // G4cout << "Absorb sucess " << G4endl;
1774  }
1775 
1776  if(Capture(false))
1777  {
1778  // haveProducts = true;
1779 #ifdef debug_BIC_StepParticlesOut
1780  G4cout << "Capture sucess " << G4endl;
1781 #endif
1782  }
1783  if ( counter > 100 && theCollisionMgr->Entries() == 0) // no collision, and stepping for some time....
1784  {
1785 #ifdef debug_BIC_StepParticlesOut
1786  PrintKTVector(&theSecondaryList,std::string("stepping 100 steps"));
1787 #endif
1788  FindCollisions(&theSecondaryList);
1789  counter=0;
1790  ++countreset;
1791  }
1792  //G4cout << "currentZ @ end loop " << currentZ << G4endl;
1793  if ( ! currentZ ){
1794  // nucleus completely destroyed, fill in ReactionProductVector
1795  // products = FillVoidNucleusProducts(products);
1796  #ifdef debug_BIC_return
1797  G4cout << "return @ Z=0 after collision loop "<< G4endl;
1798  PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
1799  G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
1800  PrintKTVector(&theTargetList,std::string(" theTargetList"));
1801  PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
1802 
1803  G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
1804  G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
1805  PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
1806  // G4cout << "returned products: " << products->size() << G4endl;
1807  #endif
1808  }
1809 
1810  }
1811  // G4cerr <<"Finished capture loop "<<G4endl;
1812 
1813  //G4cerr <<"pre- DoTimeStep 4"<<G4endl;
1814  DoTimeStep(DBL_MAX);
1815  //G4cerr <<"post- DoTimeStep 4"<<G4endl;
1816 
1817 
1818 }
1819 
1820 //----------------------------------------------------------------------------
1821 G4double G4BinaryCascade::CorrectShortlivedPrimaryForFermi(
1822  G4KineticTrack* primary,G4KineticTrackVector target_collection)
1823 //----------------------------------------------------------------------------
1824 {
1825  G4double Efermi(0);
1826  if (primary->GetState() == G4KineticTrack::inside ) {
1827  G4int PDGcode=primary->GetDefinition()->GetPDGEncoding();
1828  Efermi=((G4RKPropagation *)thePropagator)->GetField(PDGcode,primary->GetPosition());
1829 
1830  if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1831  {
1832  Efermi = ((G4RKPropagation *)thePropagator)->GetField(G4Neutron::Neutron()->GetPDGEncoding(),primary->GetPosition());
1833  G4LorentzVector mom4Primary=primary->Get4Momentum();
1834  primary->Update4Momentum(mom4Primary.e() - Efermi);
1835  }
1836 
1837  std::vector<G4KineticTrack *>::iterator titer;
1838  for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
1839  {
1840  const G4ParticleDefinition * aDef=(*titer)->GetDefinition();
1841  G4int aCode=aDef->GetPDGEncoding();
1842  G4ThreeVector aPos=(*titer)->GetPosition();
1843  Efermi+= ((G4RKPropagation *)thePropagator)->GetField(aCode, aPos);
1844  }
1845  }
1846  return Efermi;
1847 }
1848 
1849 //----------------------------------------------------------------------------
1850 G4bool G4BinaryCascade::CorrectShortlivedFinalsForFermi(G4KineticTrackVector * products,
1851  G4double initial_Efermi)
1852 //----------------------------------------------------------------------------
1853 {
1854  G4double final_Efermi(0);
1855  G4KineticTrackVector resonances;
1856  for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
1857  {
1858  G4int PDGcode=(*i)->GetDefinition()->GetPDGEncoding();
1859  // G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
1860  final_Efermi+=((G4RKPropagation *)thePropagator)->GetField(PDGcode,(*i)->GetPosition());
1861  if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
1862  {
1863  resonances.push_back(*i);
1864  }
1865  }
1866  if ( resonances.size() > 0 )
1867  {
1868  G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
1869  for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
1870  {
1871  G4LorentzVector mom=(*res)->Get4Momentum();
1872  G4double mass2=mom.mag2();
1873  G4double newEnergy=mom.e() + delta_Fermi;
1874  G4double newEnergy2= newEnergy*newEnergy;
1875  //G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
1876  if ( newEnergy2 < mass2 )
1877  {
1878  return false;
1879  }
1880  G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
1881  (*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
1882  //G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<<
1883  // " 3mom from/to " << mom.vect() << " / " << mom3 << G4endl;
1884  }
1885  }
1886  return true;
1887 }
1888 
1889 //----------------------------------------------------------------------------
1890 void G4BinaryCascade::CorrectFinalPandE()
1891 //----------------------------------------------------------------------------
1892 //
1893 // Modify momenta of outgoing particles.
1894 // Assume two body decay, nucleus(@nominal mass) + sum of final state particles(SFSP).
1895 // momentum of SFSP shall be less than momentum for two body decay.
1896 //
1897 {
1898 #ifdef debug_BIC_CorrectFinalPandE
1899  G4cerr << "BIC: -CorrectFinalPandE called" << G4endl;
1900 #endif
1901 
1902  if ( theFinalState.size() == 0 ) return;
1903 
1904  G4KineticTrackVector::iterator i;
1905  G4LorentzVector pNucleus=GetFinal4Momentum();
1906  if ( pNucleus.e() == 0 ) return; // check against explicit 0 from GetNucleus4Momentum()
1907 #ifdef debug_BIC_CorrectFinalPandE
1908  G4cerr << " -CorrectFinalPandE 3" << G4endl;
1909 #endif
1910  G4LorentzVector pFinals(0);
1911  G4int nFinals(0);
1912  for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1913  {
1914  pFinals += (*i)->Get4Momentum();
1915  ++nFinals;
1916 #ifdef debug_BIC_CorrectFinalPandE
1917  G4cout <<"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
1918  << " 4mom " << (*i)->Get4Momentum()<< G4endl;
1919 #endif
1920  }
1921 #ifdef debug_BIC_CorrectFinalPandE
1922  G4cout << "CorrectFinalPandE pN pF: " <<pNucleus << " " <<pFinals << G4endl;
1923 #endif
1924  G4LorentzVector pCM=pNucleus + pFinals;
1925 
1926  G4LorentzRotation toCMS(-pCM.boostVector());
1927  pFinals *=toCMS;
1928 #ifdef debug_BIC_CorrectFinalPandE
1929  G4cout << "CorrectFinalPandE pCM, CMS pCM " << pCM << " " <<toCMS*pCM<< G4endl;
1930  G4cout << "CorrectFinal CMS pN pF " <<toCMS*pNucleus << " "
1931  <<pFinals << G4endl
1932  << " nucleus initial mass : " <<GetFinal4Momentum().mag()
1933  <<" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus << " " <<pNucleus.mag()<< " "
1934  << pFinals.mag() << " " << pCM.mag() << G4endl;
1935 #endif
1936 
1937  G4LorentzRotation toLab = toCMS.inverse();
1938 
1939  G4double s0 = pCM.mag2();
1940  G4double m10 = GetIonMass(currentZ,currentA);
1941  G4double m20 = pFinals.mag();
1942  if( s0-(m10+m20)*(m10+m20) < 0 )
1943  {
1944 #ifdef debug_BIC_CorrectFinalPandE
1945  G4cout << "G4BinaryCascade::CorrectFinalPandE() : error! " << G4endl;
1946 
1947  G4cout << "not enough mass to correct: mass^2, A,Z, mass(nucl), mass(finals) "
1948  << (s0-(m10+m20)*(m10+m20)) << " "
1949  << currentA << " " << currentZ << " "
1950  << m10 << " " << m20
1951  << G4endl;
1952  G4cerr << " -CorrectFinalPandE 4" << G4endl;
1953 
1954  PrintKTVector(&theFinalState," mass problem");
1955 #endif
1956  return;
1957  }
1958 
1959  // Three momentum in cm system
1960  G4double pInCM = std::sqrt((s0-(m10+m20)*(m10+m20))*(s0-(m10-m20)*(m10-m20))/(4.*s0));
1961 #ifdef debug_BIC_CorrectFinalPandE
1962  G4cout <<" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM
1963  << " " << (pFinals).vect().mag()<< " " << pInCM/(pFinals).vect().mag() << G4endl;
1964 #endif
1965  if ( pFinals.vect().mag() > pInCM )
1966  {
1967  G4ThreeVector p3finals=pInCM*pFinals.vect().unit();
1968 
1969  // G4ThreeVector deltap=(p3finals - pFinals.vect() ) / nFinals;
1970  G4double factor=std::max(0.98,pInCM/pFinals.vect().mag()); // small correction
1971  G4LorentzVector qFinals(0);
1972  for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
1973  {
1974  // G4ThreeVector p3((toCMS*(*i)->Get4Momentum()).vect() + deltap);
1975  G4ThreeVector p3(factor*(toCMS*(*i)->Get4Momentum()).vect());
1976  G4LorentzVector p(p3,std::sqrt((*i)->Get4Momentum().mag2() + p3.mag2()));
1977  qFinals += p;
1978  p *= toLab;
1979 #ifdef debug_BIC_CorrectFinalPandE
1980  G4cout << " final p corrected: " << p << G4endl;
1981 #endif
1982  (*i)->Set4Momentum(p);
1983  }
1984 #ifdef debug_BIC_CorrectFinalPandE
1985  G4cout << "CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() << " "
1986  <<GetFinal4Momentum().mag() << G4endl
1987  << " CMS pFinals , mag, 3.mag : " << qFinals << " " << qFinals.mag() << " " << qFinals.vect().mag()<< G4endl;
1988  G4cerr << " -CorrectFinalPandE 5 " << factor << G4endl;
1989 #endif
1990  }
1991 #ifdef debug_BIC_CorrectFinalPandE
1992  else { G4cerr << " -CorrectFinalPandE 6 - no correction done" << G4endl; }
1993 #endif
1994 
1995 }
1996 
1997 //----------------------------------------------------------------------------
1998 void G4BinaryCascade::UpdateTracksAndCollisions(
1999  //----------------------------------------------------------------------------
2000  G4KineticTrackVector * oldSecondaries,
2001  G4KineticTrackVector * oldTarget,
2002  G4KineticTrackVector * newSecondaries)
2003 {
2004  std::vector<G4KineticTrack *>::iterator iter1, iter2;
2005 
2006  // remove old secondaries from the secondary list
2007  if(oldSecondaries)
2008  {
2009  if(!oldSecondaries->empty())
2010  {
2011  for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
2012  ++iter1)
2013  {
2014  iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(),
2015  *iter1);
2016  if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2);
2017  }
2018  theCollisionMgr->RemoveTracksCollisions(oldSecondaries);
2019  }
2020  }
2021 
2022  // remove old target from the target list
2023  if(oldTarget)
2024  {
2025  // G4cout << "################## Debugging 0 "<<G4endl;
2026  if(oldTarget->size()!=0)
2027  {
2028 
2029  // G4cout << "################## Debugging 1 "<<oldTarget->size()<<G4endl;
2030  for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
2031  {
2032  iter2 = std::find(theTargetList.begin(), theTargetList.end(),
2033  *iter1);
2034  theTargetList.erase(iter2);
2035  }
2036  theCollisionMgr->RemoveTracksCollisions(oldTarget);
2037  }
2038  }
2039 
2040  if(newSecondaries)
2041  {
2042  if(!newSecondaries->empty())
2043  {
2044  // insert new secondaries in the secondary list
2045  for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
2046  ++iter1)
2047  {
2048  theSecondaryList.push_back(*iter1);
2049  if ((*iter1)->GetState() == G4KineticTrack::undefined)
2050  {
2051  PrintKTVector(*iter1, "undefined in FindCollisions");
2052  }
2053 
2054 
2055  }
2056  // look for collisions of new secondaries
2057  FindCollisions(newSecondaries);
2058  }
2059  }
2060  // G4cout << "Exiting ... "<<oldTarget<<G4endl;
2061 }
2062 
2063 
2064 class SelectFromKTV
2065 {
2066 private:
2067  G4KineticTrackVector * ktv;
2068  G4KineticTrack::CascadeState wanted_state;
2069 public:
2070  SelectFromKTV(G4KineticTrackVector * out, G4KineticTrack::CascadeState astate)
2071  :
2072  ktv(out), wanted_state(astate)
2073  {};
2074  void operator() (G4KineticTrack *& kt) const
2075  {
2076  if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt);
2077  };
2078 };
2079 
2080 
2081 
2082 //----------------------------------------------------------------------------
2083 G4bool G4BinaryCascade::DoTimeStep(G4double theTimeStep)
2084 //----------------------------------------------------------------------------
2085 {
2086 
2087 #ifdef debug_BIC_DoTimeStep
2088  G4ping debug("debug_G4BinaryCascade");
2089  debug.push_back("======> DoTimeStep 1"); debug.dump();
2090  G4cerr <<"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
2091  << " , time="<<theCurrentTime << G4endl;
2092  PrintKTVector(&theSecondaryList, std::string("DoTimeStep - theSecondaryList"));
2093  //PrintKTVector(&theTargetList, std::string("DoTimeStep - theTargetList"));
2094 #endif
2095 
2096  G4bool success=true;
2097  std::vector<G4KineticTrack *>::iterator iter;
2098 
2099  G4KineticTrackVector * kt_outside = new G4KineticTrackVector;
2100  std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2101  SelectFromKTV(kt_outside,G4KineticTrack::outside));
2102  //PrintKTVector(kt_outside, std::string("DoTimeStep - found outside"));
2103 
2104  G4KineticTrackVector * kt_inside = new G4KineticTrackVector;
2105  std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2106  SelectFromKTV(kt_inside, G4KineticTrack::inside));
2107  // PrintKTVector(kt_inside, std::string("DoTimeStep - found inside"));
2108  //-----
2109  G4KineticTrackVector dummy; // needed for re-usability
2110 #ifdef debug_BIC_DoTimeStep
2111  G4cout << "NOW WE ARE ENTERING THE TRANSPORT"<<G4endl;
2112 #endif
2113 
2114  // =================== Here we move the particles ===================
2115 
2116  thePropagator->Transport(theSecondaryList, dummy, theTimeStep);
2117 
2118  // =================== Here we move the particles ===================
2119 
2120  //------
2121 
2122  theMomentumTransfer += thePropagator->GetMomentumTransfer();
2123 #ifdef debug_BIC_DoTimeStep
2124  G4cout << "DoTimeStep : theMomentumTransfer = " << theMomentumTransfer << G4endl;
2125  PrintKTVector(&theSecondaryList, std::string("DoTimeStep - secondaries aft trsprt"));
2126 #endif
2127 
2128  //_DebugEpConservation(" after stepping");
2129 
2130  // Partclies which went INTO nucleus
2131 
2132  G4KineticTrackVector * kt_gone_in = new G4KineticTrackVector;
2133  std::for_each( kt_outside->begin(),kt_outside->end(),
2134  SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
2135  // PrintKTVector(kt_gone_in, std::string("DoTimeStep - gone in"));
2136 
2137 
2138  // Partclies which went OUT OF nucleus
2139  G4KineticTrackVector * kt_gone_out = new G4KineticTrackVector;
2140  std::for_each( kt_inside->begin(),kt_inside->end(),
2141  SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
2142 
2143  // PrintKTVector(kt_gone_out, std::string("DoTimeStep - gone out"));
2144 
2145  G4KineticTrackVector *fail=CorrectBarionsOnBoundary(kt_gone_in,kt_gone_out);
2146 
2147  if ( fail )
2148  {
2149  // some particle(s) supposed to enter/leave were miss_nucleus/captured by the correction
2150  kt_gone_in->clear();
2151  std::for_each( kt_outside->begin(),kt_outside->end(),
2152  SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
2153 
2154  kt_gone_out->clear();
2155  std::for_each( kt_inside->begin(),kt_inside->end(),
2156  SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
2157 
2158 #ifdef debug_BIC_DoTimeStep
2159  PrintKTVector(fail,std::string(" Failed to go in/out -> miss_nucleus/captured"));
2160  PrintKTVector(kt_gone_in, std::string("recreated kt_gone_in"));
2161  PrintKTVector(kt_gone_out, std::string("recreated kt_gone_out"));
2162 #endif
2163  delete fail;
2164  }
2165 
2166  // Add tracks missing nucleus and tracks going straight though to addFinals
2167  std::for_each( kt_outside->begin(),kt_outside->end(),
2168  SelectFromKTV(kt_gone_out,G4KineticTrack::miss_nucleus));
2169  //PrintKTVector(kt_gone_out, std::string("miss to append to final state.."));
2170  std::for_each( kt_outside->begin(),kt_outside->end(),
2171  SelectFromKTV(kt_gone_out,G4KineticTrack::gone_out));
2172 
2173 #ifdef debug_BIC_DoTimeStep
2174  PrintKTVector(kt_gone_out, std::string("append gone_outs to final state.. theFinalState"));
2175 #endif
2176 
2177  theFinalState.insert(theFinalState.end(),
2178  kt_gone_out->begin(),kt_gone_out->end());
2179 
2180  // Partclies which could not leave nucleus, captured...
2181  G4KineticTrackVector * kt_captured = new G4KineticTrackVector;
2182  std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2183  SelectFromKTV(kt_captured, G4KineticTrack::captured));
2184 
2185  // Check no track is part in next collision, ie.
2186  // this step was to far, and collisions should not occur any more
2187 
2188  if ( theCollisionMgr->Entries()> 0 )
2189  {
2190  if (kt_gone_out->size() )
2191  {
2192  G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
2193  iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
2194  if ( iter != kt_gone_out->end() )
2195  {
2196  success=false;
2197 #ifdef debug_BIC_DoTimeStep
2198  G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
2199 #endif
2200  }
2201  }
2202  if ( kt_captured->size() )
2203  {
2204  G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
2205  iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
2206  if ( iter != kt_captured->end() )
2207  {
2208  success=false;
2209 #ifdef debug_BIC_DoTimeStep
2210  G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
2211 #endif
2212  }
2213  }
2214 
2215  }
2216  // PrintKTVector(kt_gone_out," kt_gone_out be4 updatetrack...");
2217  UpdateTracksAndCollisions(kt_gone_out,0 ,0);
2218 
2219 
2220  if ( kt_captured->size() )
2221  {
2222  theCapturedList.insert(theCapturedList.end(),
2223  kt_captured->begin(),kt_captured->end());
2224  //should be std::for_each(kt_captured->begin(),kt_captured->end(),
2225  // std::mem_fun(&G4KineticTrack::Hit));
2226  // but VC 6 requires:
2227  std::vector<G4KineticTrack *>::iterator i_captured;
2228  for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
2229  {
2230  (*i_captured)->Hit();
2231  }
2232  // PrintKTVector(kt_captured," kt_captured be4 updatetrack...");
2233  UpdateTracksAndCollisions(kt_captured, NULL, NULL);
2234  }
2235 
2236 #ifdef debug_G4BinaryCascade
2237  delete kt_inside;
2238  kt_inside = new G4KineticTrackVector;
2239  std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
2240  SelectFromKTV(kt_inside, G4KineticTrack::inside));
2241  if ( currentZ != (GetTotalCharge(theTargetList)
2242  + GetTotalCharge(theCapturedList)
2243  + GetTotalCharge(*kt_inside)) )
2244  {
2245  G4cout << " error-DoTimeStep aft, A, Z: " << currentA << " " << currentZ
2246  << " sum(tgt,capt,active) "
2247  << GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside)
2248  << " targets: " << GetTotalCharge(theTargetList)
2249  << " captured: " << GetTotalCharge(theCapturedList)
2250  << " active: " << GetTotalCharge(*kt_inside)
2251  << G4endl;
2252  }
2253 #endif
2254 
2255  delete kt_inside;
2256  delete kt_outside;
2257  delete kt_captured;
2258  delete kt_gone_in;
2259  delete kt_gone_out;
2260 
2261  // G4cerr <<"G4BinaryCascade::DoTimeStep: exit "<<G4endl;
2262  theCurrentTime += theTimeStep;
2263 
2264  //debug.push_back("======> DoTimeStep 2"); debug.dump();
2265  return success;
2266 
2267 }
2268 
2269 //----------------------------------------------------------------------------
2270 G4KineticTrackVector* G4BinaryCascade::CorrectBarionsOnBoundary(
2271  G4KineticTrackVector *in,
2272  G4KineticTrackVector *out)
2273 //----------------------------------------------------------------------------
2274 {
2275  G4KineticTrackVector * kt_fail(0);
2276  std::vector<G4KineticTrack *>::iterator iter;
2277  // G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2278  // << currentZ << " "<< currentA << G4endl;
2279  if (in->size())
2280  {
2281  G4int secondaries_in(0);
2282  G4int secondaryBarions_in(0);
2283  G4int secondaryCharge_in(0);
2284  G4double secondaryMass_in(0);
2285 
2286  for ( iter =in->begin(); iter != in->end(); ++iter)
2287  {
2288  ++secondaries_in;
2289  secondaryCharge_in += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
2290  if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
2291  {
2292  secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
2293  if((*iter)->GetDefinition() == G4Neutron::Neutron() ||
2294  (*iter)->GetDefinition() == G4Proton::Proton() )
2295  {
2296  secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
2297  } else {
2298  secondaryMass_in += G4Proton::Proton()->GetPDGMass();
2299  }
2300  }
2301  }
2302  G4double mass_initial= GetIonMass(currentZ,currentA);
2303 
2304  currentZ += secondaryCharge_in;
2305  currentA += secondaryBarions_in;
2306 
2307  // G4cout << "CorrectBarionsOnBoundary,secondaryCharge_in, secondaryBarions_in "
2308  // << secondaryCharge_in << " "<< secondaryBarions_in << G4endl;
2309 
2310  G4double mass_final= GetIonMass(currentZ,currentA);
2311 
2312  G4double correction= secondaryMass_in + mass_initial - mass_final;
2313  if (secondaries_in>1)
2314  {correction /= secondaries_in;}
2315 
2316 #ifdef debug_BIC_CorrectBarionsOnBoundary
2317  G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
2318  << "secondaryCharge_in,secondaryBarions_in,"
2319  << "energy correction,m_secondry,m_nucl_init,m_nucl_final "
2320  << currentZ << " "<< currentA <<" "
2321  << secondaryCharge_in<<" "<<secondaryBarions_in<<" "
2322  << correction << " "
2323  << secondaryMass_in << " "
2324  << mass_initial << " "
2325  << mass_final << " "
2326  << G4endl;
2327  PrintKTVector(in,std::string("in be4 correction"));
2328 #endif
2329  for ( iter = in->begin(); iter != in->end(); ++iter)
2330  {
2331  if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2332  {
2333  (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2334  } else {
2335  //particle cannot go in, put to miss_nucleus
2336  G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
2337  (*iter)->SetState(G4KineticTrack::miss_nucleus);
2338  // Undo correction for Colomb Barrier
2339  G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2340  (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + barrier);
2341  if ( ! kt_fail ) kt_fail=new G4KineticTrackVector;
2342  kt_fail->push_back(*iter);
2343  currentZ -= G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
2344  currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
2345 
2346  }
2347 
2348  }
2349 #ifdef debug_BIC_CorrectBarionsOnBoundary
2350  G4cout << " CorrectBarionsOnBoundary, aft, Z, A, sec-Z,A,m,m_in_nucleus "
2351  << currentZ << " " << currentA << " "
2352  << secondaryCharge_in << " " << secondaryBarions_in << " "
2353  << secondaryMass_in << " "
2354  << G4endl;
2355  PrintKTVector(in,std::string("in AFT correction"));
2356 #endif
2357 
2358  }
2359  //----------------------------------------------
2360  if (out->size())
2361  {
2362  G4int secondaries_out(0);
2363  G4int secondaryBarions_out(0);
2364  G4int secondaryCharge_out(0);
2365  G4double secondaryMass_out(0);
2366 
2367  for ( iter =out->begin(); iter != out->end(); ++iter)
2368  {
2369  ++secondaries_out;
2370  secondaryCharge_out += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
2371  if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
2372  {
2373  secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
2374  if((*iter)->GetDefinition() == G4Neutron::Neutron() ||
2375  (*iter)->GetDefinition() == G4Proton::Proton() )
2376  {
2377  secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
2378  } else {
2379  secondaryMass_out += G4Neutron::Neutron()->GetPDGMass();
2380  }
2381  }
2382  }
2383 
2384  G4double mass_initial= GetIonMass(currentZ,currentA);
2385  currentA -=secondaryBarions_out;
2386  currentZ -=secondaryCharge_out;
2387 
2388  // G4cout << "CorrectBarionsOnBoundary,secondaryCharge_out, secondaryBarions_out "
2389  // << secondaryCharge_out << " "<< secondaryBarions_out << G4endl;
2390 
2391  // a delta minus will do currentZ < 0 in light nuclei
2392  // if (currentA < 0 || currentZ < 0 )
2393  if (currentA < 0 )
2394  {
2395  G4cerr << "G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
2396  secondaryBarions_out << " " << secondaryCharge_out << G4endl;
2397  PrintKTVector(&theTargetList,"CorrectBarionsOnBoundary Target");
2398  PrintKTVector(&theCapturedList,"CorrectBarionsOnBoundary Captured");
2399  PrintKTVector(&theSecondaryList,"CorrectBarionsOnBoundary Secondaries");
2400  G4cerr << "G4BinaryCascade - currentA, currentZ " << currentA << " " << currentZ << G4endl;
2401  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
2402  }
2403  G4double mass_final=GetIonMass(currentZ,currentA);
2404  G4double correction= mass_initial - mass_final - secondaryMass_out;
2405  // G4cout << "G4BinaryCascade::CorrectBarionsOnBoundary() total out correction: " << correction << G4endl;
2406 
2407  if (secondaries_out>1) correction /= secondaries_out;
2408 #ifdef debug_BIC_CorrectBarionsOnBoundary
2409  G4cout << "DoTimeStep,(current Z,A),"
2410  << "(secondaries out,Charge,Barions),"
2411  <<"* energy correction,(m_secondry,m_nucl_init,m_nucl_final) "
2412  << "("<< currentZ << ","<< currentA <<") ("
2413  << secondaries_out << ","
2414  << secondaryCharge_out<<","<<secondaryBarions_out<<") * "
2415  << correction << " ("
2416  << secondaryMass_out << ", "
2417  << mass_initial << ", "
2418  << mass_final << ")"
2419  << G4endl;
2420  PrintKTVector(out,std::string("out be4 correction"));
2421 #endif
2422 
2423  for ( iter = out->begin(); iter != out->end(); ++iter)
2424  {
2425  if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
2426  {
2427  (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
2428  } else
2429  {
2430  // particle cannot go out due to change of nuclear potential!
2431  // capture protons and neutrons;
2432  if(((*iter)->GetDefinition() == G4Proton::Proton()) ||
2433  ((*iter)->GetDefinition() == G4Neutron::Neutron()))
2434  {
2435  (*iter)->SetState(G4KineticTrack::captured);
2436  // Undo correction for Colomb Barrier
2437  G4double barrier=((G4RKPropagation *)thePropagator)->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
2438  (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() - barrier);
2439  if ( kt_fail == 0 ) kt_fail=new G4KineticTrackVector;
2440  kt_fail->push_back(*iter);
2441  currentZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
2442  currentA += (*iter)->GetDefinition()->GetBaryonNumber();
2443  }
2444 #ifdef debug_BIC_CorrectBarionsOnBoundary
2445  else
2446  {
2447  G4cout << "Not correcting outgoing " << *iter << " "
2448  << (*iter)->GetDefinition()->GetPDGEncoding() << " "
2449  << (*iter)->GetDefinition()->GetParticleName() << G4endl;
2450  PrintKTVector(out,std::string("outgoing, one not corrected"));
2451  }
2452 #endif
2453  }
2454  }
2455 
2456 #ifdef debug_BIC_CorrectBarionsOnBoundary
2457  PrintKTVector(out,std::string("out AFTER correction"));
2458  G4cout << " DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
2459  << currentA << " "<< currentZ << " "
2460  << secondaryCharge_out << " "<< secondaryBarions_out << " "<<
2461  secondaryMass_out << " "
2462  << massInNucleus << " "
2463  << GetIonMass(currentZ,currentA)
2464  << " " << massInNucleus - GetIonMass(currentZ,currentA)
2465  << G4endl;
2466 #endif
2467  }
2468 
2469  return kt_fail;
2470 }
2471 
2472 
2473 //----------------------------------------------------------------------------
2474 
2475 G4Fragment * G4BinaryCascade::FindFragments()
2476 //----------------------------------------------------------------------------
2477 {
2478 
2479 #ifdef debug_BIC_FindFragments
2480  G4cout << "target, captured, secondary: "
2481  << theTargetList.size() << " "
2482  << theCapturedList.size()<< " "
2483  << theSecondaryList.size()
2484  << G4endl;
2485 #endif
2486 
2487  G4int a = theTargetList.size()+theCapturedList.size();
2488  G4int zTarget = 0;
2489  G4KineticTrackVector::iterator i;
2490  for(i = theTargetList.begin(); i != theTargetList.end(); ++i)
2491  {
2492  if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 )
2493  {
2494  zTarget++;
2495  }
2496  }
2497 
2498  G4int zCaptured = 0;
2499  G4LorentzVector CapturedMomentum(0.,0.,0.,0.);
2500  for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2501  {
2502  CapturedMomentum += (*i)->Get4Momentum();
2503  if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 )
2504  {
2505  zCaptured++;
2506  }
2507  }
2508 
2509  G4int z = zTarget+zCaptured;
2510 
2511 #ifdef debug_G4BinaryCascade
2512  if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) )
2513  {
2514  G4cout << " FindFragment Counting error z a " << z << " " <<a << " "
2515  << GetTotalCharge(theTargetList) << " " << GetTotalCharge(theCapturedList)<<
2516  G4endl;
2517  PrintKTVector(&theTargetList, std::string("theTargetList"));
2518  PrintKTVector(&theCapturedList, std::string("theCapturedList"));
2519  }
2520 #endif
2521  //debug
2522  /*
2523  * G4cout << " Fragment mass table / real "
2524  * << GetIonMass(z, a)
2525  * << " / " << GetFinal4Momentum().mag()
2526  * << " difference "
2527  * << GetFinal4Momentum().mag() - GetIonMass(z, a)
2528  * << G4endl;
2529  */
2530  //
2531  // if(getenv("BCDEBUG") ) G4cerr << "Fragment A, Z "<< a <<" "<< z<<G4endl;
2532  if ( z < 1 ) return 0;
2533 
2534  G4int holes = the3DNucleus->GetMassNumber() - theTargetList.size();
2535  G4int excitons = theCapturedList.size();
2536 #ifdef debug_BIC_FindFragments
2537  G4cout << "Fragment: a= " << a << " z= " << z << " particles= " << excitons
2538  << " Charged= " << zCaptured << " holes= " << holes
2539  << " excitE= " <<GetExcitationEnergy()
2540  << " Final4Momentum= " << GetFinalNucleusMomentum() << " capturMomentum= " << CapturedMomentum
2541  << G4endl;
2542 #endif
2543 
2544  G4Fragment * fragment = new G4Fragment(a,z,GetFinalNucleusMomentum());
2545  fragment->SetNumberOfHoles(holes);
2546 
2547  //GF fragment->SetNumberOfParticles(excitons-holes);
2548  fragment->SetNumberOfParticles(excitons);
2549  fragment->SetNumberOfCharged(zCaptured);
2550 
2551  return fragment;
2552 }
2553 
2554 //----------------------------------------------------------------------------
2555 
2556 G4LorentzVector G4BinaryCascade::GetFinal4Momentum()
2557 //----------------------------------------------------------------------------
2558 // Return momentum of reminder nulceus;
2559 // ie. difference of (initial state(primaries+nucleus) - final state) particles, ignoring remnant nucleus
2560 {
2561  G4LorentzVector final4Momentum = theInitial4Mom + theProjectile4Momentum;
2562  G4LorentzVector finals(0,0,0,0);
2563  for(G4KineticTrackVector::iterator i = theFinalState.begin(); i != theFinalState.end(); ++i)
2564  {
2565  final4Momentum -= (*i)->Get4Momentum();
2566  finals += (*i)->Get4Momentum();
2567  }
2568 
2569  if(final4Momentum.e()> 0 && (final4Momentum.vect()/final4Momentum.e()).mag()>1.0 && currentA > 0)
2570  {
2571 #ifdef debug_BIC_Final4Momentum
2572  G4cerr << G4endl;
2573  G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl;
2574  G4KineticTrackVector::iterator i;
2575  G4cerr <<"Total initial 4-momentum " << theProjectile4Momentum << G4endl;
2576  G4cerr <<" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<G4endl;
2577  for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
2578  {
2579  G4cerr <<" Final state: "<<(*i)->Get4Momentum()<<(*i)->GetDefinition()->GetParticleName()<<G4endl;
2580  }
2581  G4cerr << "Sum( 4-mom ) finals " << finals << G4endl;
2582  G4cerr<< " Final4Momentum = "<<final4Momentum <<" "<<final4Momentum.m()<<G4endl;
2583  G4cerr <<" current A, Z = "<< currentA<<", "<<currentZ<<G4endl;
2584  G4cerr << G4endl;
2585 #endif
2586 
2587  final4Momentum=G4LorentzVector(0,0,0,0);
2588  }
2589  return final4Momentum;
2590 }
2591 
2592 //----------------------------------------------------------------------------
2593 G4LorentzVector G4BinaryCascade::GetFinalNucleusMomentum()
2594 //----------------------------------------------------------------------------
2595 {
2596  // return momentum of nucleus for use with precompound model; also keep transformation to
2597  // apply to precompoud products.
2598 
2599  G4LorentzVector CapturedMomentum(0,0,0,0);
2600  G4KineticTrackVector::iterator i;
2601  // G4cout << "GetFinalNucleusMomentum Captured size: " <<theCapturedList.size() << G4endl;
2602  for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
2603  {
2604  CapturedMomentum += (*i)->Get4Momentum();
2605  }
2606  //G4cout << "GetFinalNucleusMomentum CapturedMomentum= " <<CapturedMomentum << G4endl;
2607  // G4cerr << "it 9"<<G4endl;
2608 
2609  G4LorentzVector NucleusMomentum = GetFinal4Momentum();
2610  if ( NucleusMomentum.e() > 0 )
2611  {
2612  // G4cout << "GetFinalNucleusMomentum GetFinal4Momentum= " <<NucleusMomentum <<" "<<NucleusMomentum.mag()<<G4endl;
2613  // boost nucleus to a frame such that the momentum of nucleus == momentum of Captured
2614  G4ThreeVector boost= (NucleusMomentum.vect() -CapturedMomentum.vect())/NucleusMomentum.e();
2615  if(boost.mag2()>1.0)
2616  {
2617 # ifdef debug_BIC_FinalNucleusMomentum
2618  G4cerr << "G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<G4endl;
2619  G4cerr << "it 0"<<boost <<G4endl;
2620  G4cerr << "it 01"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2621  G4cout <<" testing boost "<<boost<<" "<<boost.mag()<<G4endl;
2622 # endif
2623  boost=G4ThreeVector(0,0,0);
2624  NucleusMomentum=G4LorentzVector(0,0,0,0);
2625  }
2626  G4LorentzRotation nucleusBoost( -boost );
2627  precompoundLorentzboost.set( boost );
2628 #ifdef debug_debug_BIC_FinalNucleusMomentum
2629  G4cout << "GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
2630 #endif
2631  NucleusMomentum *= nucleusBoost;
2632 #ifdef debug_BIC_FinalNucleusMomentum
2633  G4cout << "GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<G4endl;
2634 #endif
2635  }
2636  return NucleusMomentum;
2637 }
2638 
2639 //----------------------------------------------------------------------------
2640 G4ReactionProductVector * G4BinaryCascade::Propagate1H1(
2641  //----------------------------------------------------------------------------
2642  G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus)
2643 {
2644  G4ReactionProductVector * products = new G4ReactionProductVector;
2645  const G4ParticleDefinition * aHTarg = G4Proton::ProtonDefinition();
2646  if (nucleus->GetCharge() == 0) aHTarg = G4Neutron::NeutronDefinition();
2647  G4double mass = aHTarg->GetPDGMass();
2648  G4KineticTrackVector * secs = 0;
2649  G4ThreeVector pos(0,0,0);
2650  G4LorentzVector mom(mass);
2651  G4KineticTrack aTarget(aHTarg, 0., pos, mom);
2652  G4bool done(false);
2653  std::vector<G4KineticTrack *>::iterator iter, jter;
2654  // data member static G4Scatterer theH1Scatterer;
2655  //G4cout << " start 1H1 for " << (*secondaries).front()->GetDefinition()->GetParticleName()
2656  // << " on " << aHTarg->GetParticleName() << G4endl;
2657  G4int tryCount(0);
2658  while(!done && tryCount++ <200) /* Loop checking, 31.08.2015, G.Folger */
2659  {
2660  if(secs)
2661  {
2662  std::for_each(secs->begin(), secs->end(), DeleteKineticTrack());
2663  delete secs;
2664  }
2665  secs = theH1Scatterer->Scatter(*(*secondaries).front(), aTarget);
2666 #ifdef debug_H1_BinaryCascade
2667  PrintKTVector(secs," From Scatter");
2668 #endif
2669  for(size_t ss=0; secs && ss<secs->size(); ss++)
2670  {
2671  // must have one resonance in final state, or it was elastic, not allowed here.
2672  if((*secs)[ss]->GetDefinition()->IsShortLived()) done = true;
2673  }
2674  }
2675 
2676  ClearAndDestroy(&theFinalState);
2677  ClearAndDestroy(secondaries);
2678  delete secondaries;
2679 
2680  for(size_t current=0; secs && current<secs->size(); current++)
2681  {
2682  if((*secs)[current]->GetDefinition()->IsShortLived())
2683  {
2684  done = true; // must have one resonance in final state, elastic not allowed here!
2685  G4KineticTrackVector * dec = (*secs)[current]->Decay();
2686  for(jter=dec->begin(); jter != dec->end(); jter++)
2687  {
2688  //G4cout << "Decay"<<G4endl;
2689  secs->push_back(*jter);
2690  //G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<<G4endl;
2691  }
2692  delete (*secs)[current];
2693  delete dec;
2694  }
2695  else
2696  {
2697  //G4cout << "FS "<<G4endl;
2698  //G4cout << "FS "<<(*secs)[current]->GetDefinition()->GetParticleName()<<G4endl;
2699  theFinalState.push_back((*secs)[current]);
2700  }
2701  }
2702 
2703  delete secs;
2704 #ifdef debug_H1_BinaryCascade
2705  PrintKTVector(&theFinalState," FinalState");
2706 #endif
2707  for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2708  {
2709  G4KineticTrack * kt = *iter;
2710  G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
2711  aNew->SetMomentum(kt->Get4Momentum().vect());
2712  aNew->SetTotalEnergy(kt->Get4Momentum().e());
2713  products->push_back(aNew);
2714 #ifdef debug_H1_BinaryCascade
2715  if (! kt->GetDefinition()->GetPDGStable() )
2716  {
2717  if (kt->GetDefinition()->IsShortLived())
2718  {
2719  G4cout << "final shortlived : ";
2720  } else
2721  {
2722  G4cout << "final un stable : ";
2723  }
2724  G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
2725  }
2726 #endif
2727  delete kt;
2728  }
2729  theFinalState.clear();
2730  return products;
2731 
2732 }
2733 
2734 //----------------------------------------------------------------------------
2735 G4ThreeVector G4BinaryCascade::GetSpherePoint(
2736  G4double r, const G4LorentzVector & mom4)
2737 //----------------------------------------------------------------------------
2738 {
2739  // Get a point outside radius.
2740  // point is random in plane (circle of radius r) orthogonal to mom,
2741  // plus -1*r*mom->vect()->unit();
2742  G4ThreeVector o1, o2;
2743  G4ThreeVector mom = mom4.vect();
2744 
2745  o1= mom.orthogonal(); // we simply need any vector non parallel
2746  o2= mom.cross(o1); // o2 is now orthogonal to mom and o1, ie. o1 and o2 define plane.
2747 
2748  G4double x2, x1;
2749 
2750  do
2751  {
2752  x1=(G4UniformRand()-.5)*2;
2753  x2=(G4UniformRand()-.5)*2;
2754  } while (sqr(x1) +sqr(x2) > 1.); /* Loop checking, 31.08.2015, G.Folger */ // or random is badly broken.....
2755 
2756  return G4ThreeVector(r*(x1*o1.unit() + x2*o2.unit() - 1.5* mom.unit()));
2757 
2758 
2759 
2760  /*
2761  * // Get a point uniformly distributed on the surface of a sphere,
2762  * // with z < 0.
2763  * G4double b = r*G4UniformRand(); // impact parameter
2764  * G4double phi = G4UniformRand()*2*pi;
2765  * G4double x = b*std::cos(phi);
2766  * G4double y = b*std::sin(phi);
2767  * G4double z = -std::sqrt(r*r-b*b);
2768  * z *= 1.001; // Get position a little bit out of the sphere...
2769  * point.setX(x);
2770  * point.setY(y);
2771  * point.setZ(z);
2772  */
2773 }
2774 
2775 //----------------------------------------------------------------------------
2776 
2777 void G4BinaryCascade::ClearAndDestroy(G4KineticTrackVector * ktv)
2778 //----------------------------------------------------------------------------
2779 {
2780  std::vector<G4KineticTrack *>::iterator i;
2781  for(i = ktv->begin(); i != ktv->end(); ++i)
2782  delete (*i);
2783  ktv->clear();
2784 }
2785 
2786 //----------------------------------------------------------------------------
2787 void G4BinaryCascade::ClearAndDestroy(G4ReactionProductVector * rpv)
2788 //----------------------------------------------------------------------------
2789 {
2790  std::vector<G4ReactionProduct *>::iterator i;
2791  for(i = rpv->begin(); i != rpv->end(); ++i)
2792  delete (*i);
2793  rpv->clear();
2794 }
2795 
2796 //----------------------------------------------------------------------------
2797 void G4BinaryCascade::PrintKTVector(G4KineticTrackVector * ktv, std::string comment)
2798 //----------------------------------------------------------------------------
2799 {
2800  if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() " << comment << G4endl;
2801  if (ktv) {
2802  G4cout << " vector: " << ktv << ", number of tracks: " << ktv->size()
2803  << G4endl;
2804  std::vector<G4KineticTrack *>::iterator i;
2805  G4int count;
2806 
2807  for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
2808  {
2809  G4KineticTrack * kt = *i;
2810  G4cout << " track n. " << count;
2811  PrintKTVector(kt);
2812  }
2813  } else {
2814  G4cout << "G4BinaryCascade::PrintKTVector():No KineticTrackVector given " << G4endl;
2815  }
2816 }
2817 //----------------------------------------------------------------------------
2818 void G4BinaryCascade::PrintKTVector(G4KineticTrack * kt, std::string comment)
2819 //----------------------------------------------------------------------------
2820 {
2821  if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() "<< comment << G4endl;
2822  if ( kt ){
2823  G4cout << ", id: " << kt << G4endl;
2824  G4ThreeVector pos = kt->GetPosition();
2825  G4LorentzVector mom = kt->Get4Momentum();
2826  G4LorentzVector tmom = kt->GetTrackingMomentum();
2827  const G4ParticleDefinition * definition = kt->GetDefinition();
2828  G4cout << " definition: " << definition->GetPDGEncoding() << " pos: "
2829  << 1/fermi*pos << " R: " << 1/fermi*pos.mag() << " 4mom: "
2830  << 1/MeV*mom <<"Tr_mom" << 1/MeV*tmom << " P: " << 1/MeV*mom.vect().mag()
2831  << " M: " << 1/MeV*mom.mag() << G4endl;
2832  G4cout <<" trackstatus: "<<kt->GetState() << " isParticipant " << (kt->IsParticipant()?"T":"F") <<G4endl;
2833  } else {
2834  G4cout << "G4BinaryCascade::PrintKTVector(): No Kinetictrack given" << G4endl;
2835  }
2836 }
2837 
2838 
2839 //----------------------------------------------------------------------------
2840 G4double G4BinaryCascade::GetIonMass(G4int Z, G4int A)
2841 //----------------------------------------------------------------------------
2842 {
2843  G4double mass(0);
2844  if ( Z > 0 && A >= Z )
2845  {
2846  mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
2847 
2848  } else if ( A > 0 && Z>0 )
2849  {
2850  // charge Z > A; will happen for light nuclei with pions involved.
2851  mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(A,A);
2852 
2853  } else if ( A >= 0 && Z<=0 )
2854  {
2855  // all neutral, or empty nucleus
2856  mass = A * G4Neutron::Neutron()->GetPDGMass();
2857 
2858  } else if ( A == 0 )
2859  {
2860  // empty nucleus, except maybe pions
2861  mass = 0;
2862 
2863  } else
2864  {
2865  G4cerr << "G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
2866  << A << "," << Z << ")" <<G4endl;
2867  throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::GetIonMass() - giving up");
2868 
2869  }
2870  // G4cout << "G4BinaryCascade::GetIonMass() Z, A, mass " << Z << " " << A << " " << mass << G4endl;
2871  return mass;
2872 }
2873 G4ReactionProductVector * G4BinaryCascade::FillVoidNucleusProducts(G4ReactionProductVector * products)
2874 {
2875  // return product when nucleus is destroyed, i.e. charge=0, or theTaregtList.size()=0
2876  G4double Esecondaries(0.);
2877  G4LorentzVector psecondaries;
2878  std::vector<G4KineticTrack *>::iterator iter;
2879  std::vector<G4ReactionProduct *>::iterator rpiter;
2880  decayKTV.Decay(&theFinalState);
2881 
2882  for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
2883  {
2884  G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
2885  aNew->SetMomentum((*iter)->Get4Momentum().vect());
2886  aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
2887  Esecondaries +=(*iter)->Get4Momentum().e();
2888  psecondaries +=(*iter)->Get4Momentum();
2889  aNew->SetNewlyAdded(true);
2890  //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
2891  products->push_back(aNew);
2892  }
2893 
2894  // pull out late particles from collisions
2895  //theCollisionMgr->Print();
2896  while(theCollisionMgr->Entries() > 0) /* Loop checking, 31.08.2015, G.Folger */
2897  {
2898  G4CollisionInitialState *
2899  collision = theCollisionMgr->GetNextCollision();
2900 
2901  if ( ! collision->GetTargetCollection().size() ){
2902  G4KineticTrackVector * lates = collision->GetFinalState();
2903  if ( lates->size() == 1 ) {
2904  G4KineticTrack * atrack=*(lates->begin());
2905  //PrintKTVector(atrack, " late particle @ void Nucl ");
2906 
2907  G4ReactionProduct * aNew = new G4ReactionProduct(atrack->GetDefinition());
2908  aNew->SetMomentum(atrack->Get4Momentum().vect());
2909  aNew->SetTotalEnergy(atrack->Get4Momentum().e());
2910  Esecondaries +=atrack->Get4Momentum().e();
2911  psecondaries +=atrack->Get4Momentum();
2912  aNew->SetNewlyAdded(true);
2913  products->push_back(aNew);
2914 
2915  }
2916  }
2917  theCollisionMgr->RemoveCollision(collision);
2918 
2919  }
2920 
2921  // decay must be after loop on Collisions, and Decay() will delete entries in theSecondaryList, refered
2922  // to by Collisions.
2923  decayKTV.Decay(&theSecondaryList);
2924 
2925  // Correct for momentum transfered to Nucleus
2926  G4ThreeVector transferCorrection(0);
2927  if ( (theSecondaryList.size() + theCapturedList.size()) > 0)
2928  {
2929  transferCorrection= theMomentumTransfer /(theSecondaryList.size() + theCapturedList.size());
2930  }
2931 
2932  for(iter = theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
2933  {
2934  G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
2935  (*iter)->Update4Momentum((*iter)->Get4Momentum().vect()+transferCorrection);
2936  aNew->SetMomentum((*iter)->Get4Momentum().vect());
2937  aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
2938  Esecondaries +=(*iter)->Get4Momentum().e();
2939  psecondaries +=(*iter)->Get4Momentum();
2940  if ( (*iter)->IsParticipant() ) aNew->SetNewlyAdded(true);
2941  products->push_back(aNew);
2942  }
2943 
2944 
2945  for(iter = theCapturedList.begin(); iter != theCapturedList.end(); ++iter)
2946  {
2947  G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
2948  (*iter)->Update4Momentum((*iter)->Get4Momentum().vect()+transferCorrection);
2949  aNew->SetMomentum((*iter)->Get4Momentum().vect());
2950  aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
2951  Esecondaries +=(*iter)->Get4Momentum().e();
2952  psecondaries +=(*iter)->Get4Momentum();
2953  aNew->SetNewlyAdded(true);
2954  products->push_back(aNew);
2955  }
2956 
2957  G4double SumMassNucleons(0.);
2958  G4LorentzVector pNucleons(0.);
2959  for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
2960  {
2961  SumMassNucleons += (*iter)->GetDefinition()->GetPDGMass();
2962  pNucleons += (*iter)->Get4Momentum();
2963  }
2964 
2965  G4double Ekinetic=theProjectile4Momentum.e() + initial_nuclear_mass - Esecondaries - SumMassNucleons;
2966  #ifdef debug_BIC_FillVoidnucleus
2967  G4LorentzVector deltaP=theProjectile4Momentum + G4LorentzVector(initial_nuclear_mass) -
2968  psecondaries - pNucleons;
2969  //G4cout << "BIC::FillVoidNucleus() nucleons : "<<theTargetList.size() << " , T: " << Ekinetic <<
2970  // ", deltaP " << deltaP << " deltaPNoNucl " << deltaP + pNucleons << G4endl;
2971  #endif
2972  if (Ekinetic > 0. && theTargetList.size()){
2973  Ekinetic /= theTargetList.size();
2974  } else {
2975  G4double Ekineticrdm(0);
2976  if (theTargetList.size()) Ekineticrdm = ( 0.1 + G4UniformRand()*5.) * MeV; // leave some Energy for Nucleons
2977  G4double TotalEkin(Ekineticrdm);
2978  for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
2979  TotalEkin+=(*rpiter)->GetKineticEnergy();
2980  }
2981  G4double correction(1.);
2982  if ( std::abs(Ekinetic) < 20*perCent * TotalEkin ){
2983  correction=1. + (Ekinetic-Ekineticrdm)/TotalEkin; // Ekinetic < 0 == IS < FS, need to reduce energies
2984  }
2985  #ifdef debug_G4BinaryCascade
2986  else {
2987  G4cout << "BLIC::FillVoidNucleus() fail correction, Ekinetic, TotalEkin " << Ekinetic << ""<< TotalEkin << G4endl;
2988  }
2989  #endif
2990 
2991  for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
2992  (*rpiter)->SetKineticEnergy((*rpiter)->GetKineticEnergy()*correction); // this sets kinetic & total energy
2993  (*rpiter)->SetMomentum((*rpiter)->GetTotalMomentum() * (*rpiter)->GetMomentum().unit());
2994 
2995  }
2996 
2997  Ekinetic=Ekineticrdm*correction;
2998  if (theTargetList.size())Ekinetic /= theTargetList.size();
2999 
3000  }
3001 
3002  for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter) {
3003  // set Nucleon it to be hit - as it is in fact
3004  (*iter)->Hit();
3005  G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
3006  aNew->SetKineticEnergy(Ekinetic);
3007  aNew->SetMomentum(aNew->GetTotalMomentum() * ((*iter)->Get4Momentum().vect().unit()));
3008  aNew->SetNewlyAdded(true);
3009  products->push_back(aNew);
3010  Esecondaries += aNew->GetTotalEnergy();
3011  psecondaries += G4LorentzVector(aNew->GetMomentum(),aNew->GetTotalEnergy() );
3012  }
3013  psecondaries=G4LorentzVector(0);
3014  for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3015  psecondaries += G4LorentzVector((*rpiter)->GetMomentum(),(*rpiter)->GetTotalEnergy() );
3016  }
3017 
3018 
3019 
3020  G4LorentzVector initial4Mom=theProjectile4Momentum + G4LorentzVector(initial_nuclear_mass);
3021 
3022  //G4cout << "::FillVoidNucleus()final e/p conservation initial" <<initial4Mom
3023  // << " final " << psecondaries << " delta " << initial4Mom-psecondaries << G4endl;
3024 
3025  G4ThreeVector SumMom=psecondaries.vect();
3026 
3027  SumMom=initial4Mom.vect()-SumMom;
3028  G4int loopcount(0);
3029 
3030  std::vector<G4ReactionProduct *>::reverse_iterator reverse; // start to correct last added first
3031  while ( SumMom.mag() > 0.1*MeV && loopcount++ < 10) /* Loop checking, 31.08.2015, G.Folger */
3032  {
3033  G4int index=products->size();
3034  for (reverse=products->rbegin(); reverse!=products->rend(); ++reverse, --index){
3035  SumMom=initial4Mom.vect();
3036  for (rpiter=products->begin(); rpiter!=products->end(); ++rpiter){
3037  SumMom-=(*rpiter)->GetMomentum();
3038  }
3039 
3040  G4double p=((*reverse)->GetMomentum()).mag();
3041  (*reverse)->SetMomentum( p*(((*reverse)->GetMomentum()+SumMom).unit()));
3042 
3043  }
3044  }
3045 
3046 
3047  return products;
3048 }
3049 G4ReactionProductVector * G4BinaryCascade::HighEnergyModelFSProducts(G4ReactionProductVector * products,
3050  G4KineticTrackVector * secondaries)
3051 {
3052  std::vector<G4KineticTrack *>::iterator iter;
3053  for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
3054  {
3055  G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
3056  aNew->SetMomentum((*iter)->Get4Momentum().vect());
3057  aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
3058  aNew->SetNewlyAdded(true);
3059  //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
3060  products->push_back(aNew);
3061  }
3062  const G4ParticleDefinition* fragment = 0;
3063  if (currentA == 1 && currentZ == 0) {
3064  fragment = G4Neutron::NeutronDefinition();
3065  } else if (currentA == 1 && currentZ == 1) {
3066  fragment = G4Proton::ProtonDefinition();
3067  } else if (currentA == 2 && currentZ == 1) {
3068  fragment = G4Deuteron::DeuteronDefinition();
3069  } else if (currentA == 3 && currentZ == 1) {
3070  fragment = G4Triton::TritonDefinition();
3071  } else if (currentA == 3 && currentZ == 2) {
3072  fragment = G4He3::He3Definition();
3073  } else if (currentA == 4 && currentZ == 2) {
3074  fragment = G4Alpha::AlphaDefinition();;
3075  } else {
3076  fragment =
3077  G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(currentZ,currentA,0.0);
3078  }
3079  if (fragment != 0) {
3080  G4ReactionProduct * theNew = new G4ReactionProduct(fragment);
3081  theNew->SetMomentum(G4ThreeVector(0,0,0));
3082  theNew->SetTotalEnergy(massInNucleus);
3083  //theNew->SetFormationTime(??0.??);
3084  //G4cout << " Nucleus (" << currentZ << ","<< currentA << "), mass "<< massInNucleus << G4endl;
3085  products->push_back(theNew);
3086  }
3087  return products;
3088 }
3089 
3090 void G4BinaryCascade::PrintWelcomeMessage()
3091 {
3092  G4cout <<"Thank you for using G4BinaryCascade. "<<G4endl;
3093 }
3094 
3095 //----------------------------------------------------------------------------
3096 void G4BinaryCascade::DebugApplyCollisionFail(G4CollisionInitialState * collision,
3097  G4KineticTrackVector * products)
3098 {
3099  G4bool havePion=false;
3100  if (products)
3101  {
3102  for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
3103  {
3104  G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
3105  if (std::abs(PDGcode)==211 || PDGcode==111 ) havePion=true;
3106  }
3107  }
3108  if ( !products || havePion)
3109  {
3110  const G4BCAction &action= *collision->GetGenerator();
3111  G4cout << " Collision " << collision << ", type: "<< typeid(action).name()
3112  << ", with NO products! " <<G4endl;
3113  G4cout << G4endl<<"Initial condition are these:"<<G4endl;
3114  G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
3115  PrintKTVector(collision->GetPrimary());
3116  for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
3117  {
3118  G4cout << "targ: "
3119  <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
3120  }
3121  PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
3122  }
3123  // if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl;
3124  // if ( lateParticleCollision && products ) PrintKTVector(products, " reaction products");
3125 }
3126 
3127 //----------------------------------------------------------------------------
3128 
3129 G4bool G4BinaryCascade::CheckChargeAndBaryonNumber(G4String where)
3130 {
3131  static G4int lastdA(0), lastdZ(0);
3132  G4int iStateA = the3DNucleus->GetMassNumber() + projectileA;
3133  G4int iStateZ = the3DNucleus->GetCharge() + projectileZ;
3134 
3135  G4int fStateA(0);
3136  G4int fStateZ(0);
3137 
3138  std::vector<G4KineticTrack *>::iterator i;
3139  G4int CapturedA(0), CapturedZ(0);
3140  G4int secsA(0), secsZ(0);
3141  for ( i=theCapturedList.begin(); i!=theCapturedList.end(); ++i) {
3142  CapturedA += (*i)->GetDefinition()->GetBaryonNumber();
3143  CapturedZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
3144  }
3145 
3146  for ( i=theSecondaryList.begin(); i!=theSecondaryList.end(); ++i) {
3147  if ( (*i)->GetState() != G4KineticTrack::inside ) {
3148  secsA += (*i)->GetDefinition()->GetBaryonNumber();
3149  secsZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
3150  }
3151  }
3152 
3153  for ( i=theFinalState.begin(); i!=theFinalState.end(); ++i) {
3154  fStateA += (*i)->GetDefinition()->GetBaryonNumber();
3155  fStateZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
3156  }
3157 
3158  G4int deltaA= iStateA - secsA - fStateA -currentA - lateA;
3159  G4int deltaZ= iStateZ - secsZ - fStateZ -currentZ - lateZ;
3160 
3161 #ifdef debugCheckChargeAndBaryonNumberverbose
3162  G4cout << where <<" A: iState= "<< iStateA<<", secs= "<< secsA<< ", fState= "<< fStateA<< ", current= "<<currentA<< ", late= " <<lateA << G4endl;
3163  G4cout << where <<" Z: iState= "<< iStateZ<<", secs= "<< secsZ<< ", fState= "<< fStateZ<< ", current= "<<currentZ<< ", late= " <<lateZ << G4endl;
3164 #endif
3165 
3166  if (deltaA != 0 || deltaZ!=0 ) {
3167  if (deltaA != lastdA || deltaZ != lastdZ ) {
3168  G4cout << "baryon/charge imbalance - " << where << G4endl
3169  << "deltaA " <<deltaA<<", iStateA "<<iStateA<< ", CapturedA "<<CapturedA <<", secsA "<<secsA
3170  << ", fStateA "<<fStateA << ", currentA "<<currentA << ", lateA "<<lateA << G4endl
3171  << "deltaZ "<<deltaZ<<", iStateZ "<<iStateZ<< ", CapturedZ "<<CapturedZ <<", secsZ "<<secsZ
3172  << ", fStateZ "<<fStateZ << ", currentZ "<<currentZ << ", lateZ "<<lateZ << G4endl<< G4endl;
3173  lastdA=deltaA;
3174  lastdZ=deltaZ;
3175  }
3176  } else { lastdA=lastdZ=0;}
3177 
3178  return true;
3179 }
3180 //----------------------------------------------------------------------------
3181 void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision,
3182  G4KineticTrackVector * products)
3183 {
3184 
3185  PrintKTVector(collision->GetPrimary(),std::string(" Primary particle"));
3186  PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
3187  PrintKTVector(products,std::string(" Scatterer products"));
3188 
3189 #ifdef dontUse
3190  G4double thisExcitation(0);
3191  // excitation energy from this collision
3192  // initial state:
3193  G4double initial(0);
3194  G4KineticTrack * kt=collision->GetPrimary();
3195  initial += kt->Get4Momentum().e();
3196 
3197  G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
3198 
3199  initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
3200  initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
3201  G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
3202  << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
3203  << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
3204  << " " << initial << G4endl;;
3205 
3206  G4KineticTrackVector ktv=collision->GetTargetCollection();
3207  for ( unsigned int it=0; it < ktv.size(); it++)
3208  {
3209  kt=ktv[it];
3210  initial += kt->Get4Momentum().e();
3211  thisExcitation += kt->GetDefinition()->GetPDGMass()
3212  - kt->Get4Momentum().e()
3213  - RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
3214  // initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
3215  // initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
3216  G4cout << "Targ. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
3217  << " " << kt->Get4Momentum().e()
3218  << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
3219  << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
3220  << " " << initial <<" Excit " << thisExcitation << G4endl;;
3221  }
3222 
3223  G4double final(0);
3224  G4double mass_out(0);
3225  G4int product_barions(0);
3226  if ( products )
3227  {
3228  for ( unsigned int it=0; it < products->size(); it++)
3229  {
3230  kt=(*products)[it];
3231  final += kt->Get4Momentum().e();
3232  final += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
3233  final += RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
3234  if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
3235  mass_out += kt->GetDefinition()->GetPDGMass();
3236  G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
3237  << " " << kt->Get4Momentum().e()
3238  << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
3239  << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
3240  << " " << final << G4endl;;
3241  }
3242  }
3243 
3244 
3245  G4int finalA = currentA;
3246  G4int finalZ = currentZ;
3247  if ( products )
3248  {
3249  finalA -= product_barions;
3250  finalZ -= GetTotalCharge(*products);
3251  }
3252  G4double delta = GetIonMass(currentZ,currentA) - (GetIonMass(finalZ,finalA) + mass_out);
3253  G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ
3254  << " delta-mass " << delta<<G4endl;
3255  final+=delta;
3256  mass_out = GetIonMass(finalZ,finalA);
3257  G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
3258  << " " << final << " "
3259  << mass_out<<" "
3260  << currentInitialEnergy - final - mass_out
3261  << G4endl;
3262  currentInitialEnergy-=final;
3263 #endif
3264 }
3265 
3266 //----------------------------------------------------------------------------
3267 G4bool G4BinaryCascade::DebugFinalEpConservation(const G4HadProjectile & aTrack,
3268  G4ReactionProductVector* products)
3269 //----------------------------------------------------------------------------
3270 {
3271  G4ReactionProductVector::iterator iter;
3272  G4double Efinal(0);
3273  G4ThreeVector pFinal(0);
3274  if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
3275  {
3276  G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
3277  }
3278 
3279  for(iter = products->begin(); iter != products->end(); ++iter)
3280  {
3281 
3282  G4cout << " Secondary E - Ekin / p " <<
3283  (*iter)->GetDefinition()->GetParticleName() << " " <<
3284  (*iter)->GetTotalEnergy() << " - " <<
3285  (*iter)->GetKineticEnergy()<< " / " <<
3286  (*iter)->GetMomentum().x() << " " <<
3287  (*iter)->GetMomentum().y() << " " <<
3288  (*iter)->GetMomentum().z() << G4endl;
3289  Efinal += (*iter)->GetTotalEnergy();
3290  pFinal += (*iter)->GetMomentum();
3291  }
3292 
3293  G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
3294  G4cout << "BIC E/p delta " <<
3295  (aTrack.Get4Momentum().e()+theInitial4Mom.e() - Efinal)/MeV <<
3296  " MeV / mom " << (aTrack.Get4Momentum() - pFinal ) /MeV << G4endl;
3297 
3298  return (aTrack.Get4Momentum().e() + theInitial4Mom.e() - Efinal)/aTrack.Get4Momentum().e() < perCent;
3299 
3300 }
3301 //----------------------------------------------------------------------------
3302 G4bool G4BinaryCascade::DebugEpConservation(const G4String where)
3303 //----------------------------------------------------------------------------
3304 {
3305  G4cout << where << G4endl;
3306  G4LorentzVector psecs, ptgts, pcpts, pfins;
3307  if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
3308  {
3309  G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
3310  }
3311 
3312  std::vector<G4KineticTrack *>::iterator ktiter;
3313  for(ktiter = theSecondaryList.begin(); ktiter != theSecondaryList.end(); ++ktiter)
3314  {
3315 
3316  G4cout << " Secondary E - Ekin / p " <<
3317  (*ktiter)->GetDefinition()->GetParticleName() << " " <<
3318  (*ktiter)->Get4Momentum().e() << " - " <<
3319  (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() << " / " <<
3320  (*ktiter)->Get4Momentum().vect() << G4endl;
3321  psecs += (*ktiter)->Get4Momentum();
3322  }
3323 
3324  for(ktiter = theTargetList.begin(); ktiter != theTargetList.end(); ++ktiter)
3325  {
3326 
3327  G4cout << " Target E - Ekin / p " <<
3328  (*ktiter)->GetDefinition()->GetParticleName() << " " <<
3329  (*ktiter)->Get4Momentum().e() << " - " <<
3330  (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() << " / " <<
3331  (*ktiter)->Get4Momentum().vect() << G4endl;
3332  ptgts += (*ktiter)->Get4Momentum();
3333  }
3334 
3335  for(ktiter = theCapturedList.begin(); ktiter != theCapturedList.end(); ++ktiter)
3336  {
3337 
3338  G4cout << " Captured E - Ekin / p " <<
3339  (*ktiter)->GetDefinition()->GetParticleName() << " " <<
3340  (*ktiter)->Get4Momentum().e() << " - " <<
3341  (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() << " / " <<
3342  (*ktiter)->Get4Momentum().vect() << G4endl;
3343  pcpts += (*ktiter)->Get4Momentum();
3344  }
3345 
3346  for(ktiter = theFinalState.begin(); ktiter != theFinalState.end(); ++ktiter)
3347  {
3348 
3349  G4cout << " Finals E - Ekin / p " <<
3350  (*ktiter)->GetDefinition()->GetParticleName() << " " <<
3351  (*ktiter)->Get4Momentum().e() << " - " <<
3352  (*ktiter)->Get4Momentum().e() - (*ktiter)->Get4Momentum().mag() << " / " <<
3353  (*ktiter)->Get4Momentum().vect() << G4endl;
3354  pfins += (*ktiter)->Get4Momentum();
3355  }
3356 
3357  G4cout << " Secondaries " << psecs << ", Targets " << ptgts << G4endl
3358  <<" Captured " << pcpts << ", Finals " << pfins << G4endl
3359  <<" Sum " << psecs + ptgts + pcpts + pfins << " PTransfer " << theMomentumTransfer
3360  <<" Sum+PTransfer " << psecs + ptgts + pcpts + pfins + theMomentumTransfer
3361  << G4endl<< G4endl;
3362 
3363 
3364  return true;
3365 
3366 }
3367 
3368 //----------------------------------------------------------------------------
3369 G4ReactionProductVector * G4BinaryCascade::ProductsAddFakeGamma(G4ReactionProductVector *products )
3370 //----------------------------------------------------------------------------
3371 {
3372  // else
3373 // {
3374 // G4ReactionProduct * aNew=0;
3375 // // return nucleus e and p
3376 // if (fragment != 0 ) {
3377 // aNew = new G4ReactionProduct(G4Gamma::GammaDefinition()); // we only want to pass e/p
3378 // aNew->SetMomentum(fragment->GetMomentum().vect());
3379 // aNew->SetTotalEnergy(fragment->GetMomentum().e());
3380 // delete fragment;
3381 // fragment=0;
3382 // } else if (products->size() == 0) {
3383 // // FixMe GF: for testing without precompound, return 1 gamma of 0.01 MeV in +x
3384 //#include "G4Gamma.hh"
3385 // aNew = new G4ReactionProduct(G4Gamma::GammaDefinition());
3386 // aNew->SetMomentum(G4ThreeVector(0.01*MeV,0,0));
3387 // aNew->SetTotalEnergy(0.01*MeV);
3388 // }
3389 // if ( aNew != 0 ) products->push_back(aNew);
3390 // }
3391  return products;
3392 }
3393 
3394 //----------------------------------------------------------------------------
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4HadFinalState::GetWeightChange
G4double GetWeightChange() const
Definition: G4HadFinalState.hh:78
G4ping
Definition: G4ping.hh:33
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
name
const XML_Char * name
Definition: expat.h:151
G4KineticTrack::IsParticipant
G4bool IsParticipant() const
Definition: G4KineticTrack.hh:395
G4Triton::TritonDefinition
static G4Triton * TritonDefinition()
Definition: G4Triton.cc:90
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
G4KineticTrack::Update4Momentum
void Update4Momentum(G4double aEnergy)
Definition: G4KineticTrack.hh:259
G4He3::He3Definition
static G4He3 * He3Definition()
Definition: G4He3.cc:89
G4ExcitationHandler::ModelDescription
void ModelDescription(std::ostream &outFile) const
Definition: G4ExcitationHandler.cc:505
stopAndKill
Definition: G4HadFinalState.hh:42
G4KineticTrack::GetState
CascadeState GetState() const
Definition: G4KineticTrack.hh:402
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4Nucleon::GetMomentum
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
G4V3DNucleus::GetCharge
virtual G4int GetCharge()=0
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Track.hh
G4VIntraNuclearTransportModel::GetPrimaryProjectile
const G4HadProjectile * GetPrimaryProjectile() const
Definition: G4VIntraNuclearTransportModel.hh:145
G4VPreCompoundModel::DeExcite
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:230
G4V3DNucleus::StartLoop
virtual G4bool StartLoop()=0
G4VIntraNuclearTransportModel::GetDeExcitation
G4VPreCompoundModel * GetDeExcitation() const
Definition: G4VIntraNuclearTransportModel.hh:132
G4ExcitationHandler.hh
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4VFieldPropagation::Transport
virtual void Transport(G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep)=0
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4CollisionManager::RemoveCollision
void RemoveCollision(G4CollisionInitialState *collision)
Definition: G4CollisionManager.cc:70
G4V3DNucleus::GetNuclearDensity
virtual const G4VNuclearDensity * GetNuclearDensity() const =0
G4CollisionInitialState::GetPrimary
G4KineticTrack * GetPrimary(void)
Definition: G4CollisionInitialState.hh:69
perCent
static constexpr double perCent
Definition: G4SIunits.hh:332
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:171
G4ParticleDefinition::GetPDGStable
G4bool GetPDGStable() const
p
const char * p
Definition: xmltok.h:285
G4CollisionManager::GetNextCollision
G4CollisionInitialState * GetNextCollision()
Definition: G4CollisionManager.cc:140
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4V3DNucleus::GetMassNumber
virtual G4int GetMassNumber()=0
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4FermiPhaseSpaceDecay.hh
G4VFieldPropagation::GetMomentumTransfer
virtual G4ThreeVector GetMomentumTransfer() const =0
_DebugEpConservation
#define _DebugEpConservation(val)
Definition: G4BinaryCascade.cc:109
G4BCLateParticle::GetCollisions
virtual const std::vector< G4CollisionInitialState * > & GetCollisions(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &, G4double theCurrentTime)
Definition: G4BCLateParticle.hh:43
G4Fragment::SetNumberOfHoles
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:375
G4ListOfCollisions.hh
G4KineticTrack::GetActualMass
G4double GetActualMass() const
Definition: G4KineticTrack.hh:319
G4HadronicInteractionRegistry.hh
G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:116
G4Nucleon::GetPosition
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
_CheckChargeAndBaryonNumber_
#define _CheckChargeAndBaryonNumber_(val)
Definition: G4BinaryCascade.cc:102
G4BinaryCascade::PropagateModelDescription
virtual void PropagateModelDescription(std::ostream &) const
Definition: G4BinaryCascade.cc:212
G4int
int G4int
Definition: G4Types.hh:78
G4ExcitationHandler::BreakItUp
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
Definition: G4ExcitationHandler.cc:202
G4V3DNucleus::CoulombBarrier
virtual G4double CoulombBarrier()=0
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4CollisionManager::RemoveTracksCollisions
void RemoveTracksCollisions(G4KineticTrackVector *ktv)
Definition: G4CollisionManager.cc:80
G4Nucleon::GetDefinition
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
G4BinaryCascade::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
Definition: G4BinaryCascade.cc:253
G4FermiMomentum.hh
G4DecayKineticTracks.hh
G4ReactionProduct::SetNewlyAdded
void SetNewlyAdded(const G4bool f)
Definition: G4ReactionProduct.hh:171
G4RKPropagation::GetBarrier
G4double GetBarrier(G4int encoding)
Definition: G4RKPropagation.hh:84
G4HadFinalState::SetStatusChange
void SetStatusChange(G4HadFinalStateStatus aS)
Definition: G4HadFinalState.hh:67
G4ReactionProductVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
Definition: G4ReactionProductVector.hh:43
G4V3DNucleus::GetOuterRadius
virtual G4double GetOuterRadius()=0
G4NuclearShellModelDensity.hh
G4CollisionInitialState::GetFinalState
G4KineticTrackVector * GetFinalState()
Definition: G4CollisionInitialState.cc:95
G4KineticTrackVector.hh
G4ReactionProductVector.hh
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
G4BinaryCascade::~G4BinaryCascade
virtual ~G4BinaryCascade()
Definition: G4BinaryCascade.cc:172
G4PreCompoundModel.hh
G4VPreCompoundModel::GetExcitationHandler
G4ExcitationHandler * GetExcitationHandler() const
Definition: G4VPreCompoundModel.hh:93
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:658
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4Fragment::GetA_asInt
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
G4KineticTrack::SetState
CascadeState SetState(const CascadeState new_state)
Definition: G4KineticTrack.hh:408
G4INCL::ClusterDecay::decay
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
Definition: G4INCLClusterDecay.cc:553
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4V3DNucleus::Init
virtual void Init(G4int theA, G4int theZ)=0
G4InuclParticleNames::nucleon
G4bool nucleon(G4int ityp)
Definition: G4InuclParticleNames.hh:95
G4Fragment.hh
G4Neutron.hh
CLHEP::HepLorentzVector::mag
double mag() const
test1.atrack
tuple atrack
Definition: test1.py:9
G4PionMinus::PionMinusDefinition
static G4PionMinus * PionMinusDefinition()
Definition: G4PionMinus.cc:93
G4bool
bool G4bool
Definition: G4Types.hh:79
G4IonTable::GetIonMass
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1335
G4Fragment::GetMomentum
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:307
G4KineticTrack::SetNucleon
void SetNucleon(G4Nucleon *aN)
Definition: G4KineticTrack.hh:106
G4Nucleon::AreYouHit
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
G4Proton.hh
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
G4FermiMomentum::GetFermiMomentum
G4double GetFermiMomentum(G4double density)
Definition: G4FermiMomentum.hh:45
G4PionPlus::PionPlusDefinition
static G4PionPlus * PionPlusDefinition()
Definition: G4PionPlus.cc:93
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
CLHEP::HepLorentzRotation::set
HepLorentzRotation & set(double bx, double by, double bz)
Definition: LorentzRotation.cc:27
G4VFieldPropagation::Init
virtual void Init(G4V3DNucleus *theNucleus)=0
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4KineticTrack::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4KineticTrack.hh:250
G4Fragment::SetNumberOfParticles
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:384
G4CollisionManager.hh
G4CollisionInitialState::GetTargetCollection
G4KineticTrackVector & GetTargetCollection(void)
Definition: G4CollisionInitialState.hh:81
G4BinaryCascade::Propagate
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
Definition: G4BinaryCascade.cc:379
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4Scatterer.hh
G4BCDecay::GetCollisions
virtual const std::vector< G4CollisionInitialState * > & GetCollisions(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &, G4double theCurrentTime)
Definition: G4BCDecay.hh:41
reverse
Definition: F04UserTrackInformation.hh:37
G4BinaryCascade::G4BinaryCascade
G4BinaryCascade(G4VPreCompoundModel *ptr=0)
Definition: G4BinaryCascade.cc:115
G4CollisionManager::AddCollision
void AddCollision(G4double time, G4KineticTrack *proj, G4KineticTrack *target=NULL)
Definition: G4CollisionManager.cc:48
G4ShortLivedConstructor.hh
G4PhysicalConstants.hh
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4HadronicInteractionRegistry::FindModel
G4HadronicInteraction * FindModel(const G4String &name)
Definition: G4HadronicInteractionRegistry.cc:117
SelectFromKTV::operator()
void operator()(G4KineticTrack *&kt) const
Definition: G4BinaryCascade.cc:2074
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:39
G4BinaryCascade::ModelDescription
virtual void ModelDescription(std::ostream &) const
Definition: G4BinaryCascade.cc:185
G4Absorber.hh
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4FermiPhaseSpaceDecay::Decay
std::vector< G4LorentzVector * > * Decay(G4double parent_mass, const std::vector< G4double > &fragment_masses) const
Definition: G4FermiPhaseSpaceDecay.hh:95
isAlive
Definition: G4HadFinalState.hh:42
G4FermiMomentum::Init
void Init(G4int anA, G4int aZ)
Definition: G4FermiMomentum.hh:43
G4DecayKineticTracks::Decay
void Decay(G4KineticTrackVector *tracks) const
Definition: G4DecayKineticTracks.cc:43
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4VNuclearDensity::GetDensity
G4double GetDensity(const G4ThreeVector &aPosition) const
Definition: G4VNuclearDensity.hh:43
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
debug
#define debug
Definition: CCalPrimaryGeneratorAction.cc:47
G4Scatterer::Scatter
virtual G4KineticTrackVector * Scatter(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4Scatterer.cc:284
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4KineticTrack::GetTrackingMomentum
const G4LorentzVector & GetTrackingMomentum() const
Definition: G4KineticTrack.hh:245
G4HadronicInteractionRegistry::Instance
static G4HadronicInteractionRegistry * Instance()
Definition: G4HadronicInteractionRegistry.cc:36
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::HepLorentzVector::mag2
double mag2() const
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
G4Delete.hh
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
CLHEP::Hep3Vector::mag2
double mag2() const
G4RKPropagation::GetField
G4double GetField(G4int encoding, G4ThreeVector pos)
Definition: G4RKPropagation.hh:92
test.z
tuple z
Definition: test.py:28
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4VIntraNuclearTransportModel::SetDeExcitation
void SetDeExcitation(G4VPreCompoundModel *ptr)
Definition: G4VIntraNuclearTransportModel.hh:138
G4V3DNucleus::GetMass
virtual G4double GetMass()=0
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4HadFinalState::AddSecondary
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
Definition: G4HadFinalState.hh:61
G4MesonAbsorption.hh
G4CollisionInitialState.hh
G4V3DNucleus::GetNextNucleon
virtual G4Nucleon * GetNextNucleon()=0
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
G4Fragment::SetNumberOfCharged
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:389
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4VPreCompoundModel::DeExciteModelDescription
virtual void DeExciteModelDescription(std::ostream &outFile) const =0
Definition: G4VPreCompoundModel.cc:49
G4SystemOfUnits.hh
G4Fancy3DNucleus.hh
fermi
static constexpr double fermi
Definition: G4SIunits.hh:103
G4BinaryCascade.hh
G4ping.hh
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:240
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4RKPropagation.hh
CLHEP::Hep3Vector::mag
double mag() const
G4HadronicInteraction::SetEnergyMomentumCheckLevels
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel)
Definition: G4HadronicInteraction.hh:150
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:210
G4Deuteron::DeuteronDefinition
static G4Deuteron * DeuteronDefinition()
Definition: G4Deuteron.cc:89
G4Alpha::AlphaDefinition
static G4Alpha * AlphaDefinition()
Definition: G4Alpha.cc:84
G4Neutron::NeutronDefinition
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.hh:56
ns
#define ns
Definition: xmlparse.cc:614
G4CollisionInitialState::GetGenerator
const G4BCAction * GetGenerator()
Definition: G4CollisionInitialState.hh:111
SelectFromKTV::SelectFromKTV
SelectFromKTV(G4KineticTrackVector *out, G4KineticTrack::CascadeState astate)
Definition: G4BinaryCascade.cc:2070
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66
G4HadronicInteraction::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4NuclearFermiDensity.hh
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4V3DNucleus.hh
G4cerr
G4GLOB_DLL std::ostream G4cerr
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh