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