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