ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BinaryCascade.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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(std::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(!std::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(std::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(std::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;
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();
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 {
1024  G4ReactionProductVector * result=0;
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();
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 //----------------------------------------------------------------------------