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