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