ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GeneratorPrecompoundInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GeneratorPrecompoundInterface.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -----------------------------------------------------------------------------
28 // GEANT 4 class file
29 //
30 // History: first implementation
31 // HPW, 10DEC 98, the decay part originally written by Gunter Folger
32 // in his FTF-test-program.
33 //
34 // M.Kelsey, 28 Jul 2011 -- Replace loop to decay input secondaries
35 // with new utility class, simplify cleanup loops
36 // -----------------------------------------------------------------------------
37 
38 #include <algorithm>
39 #include <vector>
40 
41 #include "G4GeneratorPrecompoundInterface.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4DynamicParticleVector.hh"
45 #include "G4KineticTrackVector.hh"
46 #include "G4Proton.hh"
47 #include "G4Neutron.hh"
48 
49 #include "G4Deuteron.hh"
50 #include "G4Triton.hh"
51 #include "G4He3.hh"
52 #include "G4Alpha.hh"
53 
54 #include "G4V3DNucleus.hh"
55 #include "G4Nucleon.hh"
56 
57 #include "G4AntiProton.hh"
58 #include "G4AntiNeutron.hh"
59 #include "G4AntiDeuteron.hh"
60 #include "G4AntiTriton.hh"
61 #include "G4AntiHe3.hh"
62 #include "G4AntiAlpha.hh"
63 
64 #include "G4FragmentVector.hh"
65 #include "G4ReactionProduct.hh"
66 #include "G4ReactionProductVector.hh"
67 #include "G4PreCompoundModel.hh"
68 #include "G4ExcitationHandler.hh"
69 #include "G4DecayKineticTracks.hh"
70 #include "G4HadronicInteractionRegistry.hh"
71 
72 //---------------------------------------------------------------------
73 #include "Randomize.hh"
74 #include "G4Log.hh"
75 
76 //#define debugPrecoInt
77 
78 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* preModel)
79 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), DeltaR(0.0)
80 {
81  proton = G4Proton::Proton();
82  neutron = G4Neutron::Neutron();
83 
84  deuteron=G4Deuteron::Deuteron();
85  triton =G4Triton::Triton();
86  He3 =G4He3::He3();
87  He4 =G4Alpha::Alpha();
88 
89  ANTIproton=G4AntiProton::AntiProton();
90  ANTIneutron=G4AntiNeutron::AntiNeutron();
91 
92  ANTIdeuteron=G4AntiDeuteron::AntiDeuteron();
93  ANTItriton =G4AntiTriton::AntiTriton();
94  ANTIHe3 =G4AntiHe3::AntiHe3();
95  ANTIHe4 =G4AntiAlpha::AntiAlpha();
96 
97  if(preModel) { SetDeExcitation(preModel); }
98  else {
99  G4HadronicInteraction* hadi =
100  G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
101  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
102  if(!pre) { pre = new G4PreCompoundModel(); }
103  SetDeExcitation(pre);
104  }
105 }
106 
107 G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
108 {
109 }
110 
111 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
112 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
113 {
114  #ifdef debugPrecoInt
115  G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::Propagate"<<G4endl;
116  G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
117  G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
118  #endif
119 
120  G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
121 
122  // decay the strong resonances
123  G4DecayKineticTracks decay(theSecondaries);
124  #ifdef debugPrecoInt
125  G4cout<<"Final stable particles number "<<theSecondaries->size()<<G4endl;
126  #endif
127 
128  // prepare the fragment
129  G4int anA=theNucleus->GetMassNumber();
130  G4int aZ=theNucleus->GetCharge();
131 // G4double TargetNucleusMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
132 
133  G4int numberOfEx = 0;
134  G4int numberOfCh = 0;
135  G4int numberOfHoles = 0;
136 
137  G4double R = theNucleus->GetNuclearRadius();
138 
139  G4LorentzVector captured4Momentum(0.,0.,0.,0.);
140  G4LorentzVector Residual4Momentum(0.,0.,0.,0.); // TargetNucleusMass is not need at the moment
141  G4LorentzVector Secondary4Momentum(0.,0.,0.,0.);
142 
143  // loop over secondaries
144  G4KineticTrackVector::iterator iter;
145  for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
146  {
147  const G4ParticleDefinition* part = (*iter)->GetDefinition();
148  G4double e = (*iter)->Get4Momentum().e();
149  G4double mass = (*iter)->Get4Momentum().mag();
150  G4ThreeVector mom = (*iter)->Get4Momentum().vect();
151  if((part != proton && part != neutron) ||
152  ((*iter)->GetPosition().mag() > R)) {
153  G4ReactionProduct * theNew = new G4ReactionProduct(part);
154  theNew->SetMomentum(mom);
155  theNew->SetTotalEnergy(e);
156  theTotalResult->push_back(theNew);
157  Secondary4Momentum += (*iter)->Get4Momentum();
158  #ifdef debugPrecoInt
159  G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
160  <<(*iter)->Get4Momentum().mag()<<G4endl;
161  #endif
162  } else {
163  if( e-mass > -CaptureThreshold*G4Log( G4UniformRand()) ) {
164  G4ReactionProduct * theNew = new G4ReactionProduct(part);
165  theNew->SetMomentum(mom);
166  theNew->SetTotalEnergy(e);
167  theTotalResult->push_back(theNew);
168  Secondary4Momentum += (*iter)->Get4Momentum();
169  #ifdef debugPrecoInt
170  G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
171  <<(*iter)->Get4Momentum().mag()<<G4endl;
172  #endif
173  } else {
174  // within the nucleus, neutron or proton
175  // now calculate A, Z of the fragment, momentum, number of exciton states
176  ++anA;
177  ++numberOfEx;
178  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
179  aZ += Z;
180  numberOfCh += Z;
181  captured4Momentum += (*iter)->Get4Momentum();
182  #ifdef debugPrecoInt
183  G4cout<<"Captured 4Mom "<<part->GetParticleName()<<(*iter)->Get4Momentum()<<G4endl;
184  #endif
185  }
186  }
187  delete (*iter);
188  }
189  delete theSecondaries;
190 
191  // loop over wounded nucleus
192  G4Nucleon * theCurrentNucleon =
193  theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
194  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
195  {
196  if(theCurrentNucleon->AreYouHit()) {
197  ++numberOfHoles;
198  ++numberOfEx;
199  --anA;
200  aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
201 
202  Residual4Momentum -= theCurrentNucleon->Get4Momentum();
203  }
204  theCurrentNucleon = theNucleus->GetNextNucleon();
205  }
206 
207  #ifdef debugPrecoInt
208  G4cout<<G4endl;
209  G4cout<<"Secondary 4Mom "<<Secondary4Momentum<<G4endl;
210  G4cout<<"Captured 4Mom "<<captured4Momentum<<G4endl;
211  G4cout<<"Sec + Captured "<<Secondary4Momentum+captured4Momentum<<G4endl;
212  G4cout<<"Residual4Mom "<<Residual4Momentum<<G4endl;
213  G4cout<<"Sum 4 momenta "
214  <<Secondary4Momentum + captured4Momentum + Residual4Momentum <<G4endl;
215  #endif
216 
217  // Check that we use QGS model; loop over wounded nucleus
218  G4bool QGSM(false);
219  theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
220  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
221  {
222  if(theCurrentNucleon->AreYouHit())
223  {
224  if(theCurrentNucleon->Get4Momentum().mag() <
225  theCurrentNucleon->GetDefinition()->GetPDGMass()) QGSM=true;
226  }
227  theCurrentNucleon = theNucleus->GetNextNucleon();
228  }
229 
230  #ifdef debugPrecoInt
231  if(!QGSM){
232  G4cout<<G4endl;
233  G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
234  G4cout<<"Residual 4Mom "<<Residual4Momentum<<G4endl;
235  if(numberOfEx == 0)
236  {G4cout<<"Residual 4Mom = 0 means that there were not wounded and captured nucleons"<<G4endl;}
237  }
238  #endif
239 
240  if(anA == 0) return theTotalResult;
241 
242  G4LorentzVector exciton4Momentum(0.,0.,0.,0.);
243  if(anA >= aZ)
244  {
245  if(!QGSM)
246  { // FTF model was used
247  G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
248 
249 // G4LorentzVector exciton4Momentum = Residual4Momentum + captured4Momentum;
250  exciton4Momentum = Residual4Momentum + captured4Momentum;
251 //exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
252  G4double ActualMass = exciton4Momentum.mag();
253  if(ActualMass <= fMass ) {
254  exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
255  }
256 
257  #ifdef debugPrecoInt
258  G4double exEnergy = 0.0;
259  if(ActualMass <= fMass ) {exEnergy = 0.;}
260  else {exEnergy = ActualMass - fMass;}
261  G4cout<<"Ground state residual Mass "<<fMass<<" E* "<<exEnergy<<G4endl;
262  #endif
263  }
264  else
265  { // QGS model was used
266  G4double InitialTargetMass =
267  G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(), theNucleus->GetCharge());
268 
269  exciton4Momentum =
270  GetPrimaryProjectile()->Get4Momentum() + G4LorentzVector(0.,0.,0.,InitialTargetMass)
271  -Secondary4Momentum;
272 
273  G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
274  G4double ActualMass = exciton4Momentum.mag();
275 
276  #ifdef debugPrecoInt
277  G4cout<<G4endl;
278  G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
279  G4cout<<"Residual4Momentum "<<exciton4Momentum<<G4endl;
280  G4cout<<"ResidualMass, GroundStateMass and E* "<<ActualMass<<" "<<fMass<<" "
281  <<ActualMass - fMass<<G4endl;
282  #endif
283 
284  if(ActualMass - fMass < 0.)
285  {
286  G4double ResE = std::sqrt(exciton4Momentum.vect().mag2() + sqr(fMass+10*MeV));
287  exciton4Momentum.setE(ResE);
288  #ifdef debugPrecoInt
289  G4cout<<"ActualMass - fMass < 0. "<<ActualMass<<" "<<fMass<<" "<<ActualMass - fMass<<G4endl;
290  G4int Uzhi; G4cin>>Uzhi;
291  #endif
292  }
293  }
294  // Need to de-excite the remnant nucleus only if excitation energy > 0.
295  G4Fragment anInitialState(anA, aZ, exciton4Momentum);
296  anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
297  anInitialState.SetNumberOfCharged(numberOfCh);
298  anInitialState.SetNumberOfHoles(numberOfHoles);
299 
300  G4ReactionProductVector * aPrecoResult =
301  theDeExcitation->DeExcite(anInitialState);
302  // fill pre-compound part into the result, and return
303  #ifdef debugPrecoInt
304  G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
305  #endif
306  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
307  {
308  theTotalResult->push_back(aPrecoResult->operator[](ll));
309  #ifdef debugPrecoInt
310  G4cout<<"Fragment "<<ll<<" "
311  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
312  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
313  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
314  <<aPrecoResult->operator[](ll)->GetDefinition()->GetPDGMass()<<G4endl;
315  #endif
316  }
317  delete aPrecoResult;
318  }
319 
320  return theTotalResult;
321 }
322 
323 G4HadFinalState* G4GeneratorPrecompoundInterface::
324 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
325 {
326  G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
327  << G4endl;
328  G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
329  G4cout << "Please remove from your physics list."<<G4endl;
330  throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
331  return new G4HadFinalState;
332 }
333 void G4GeneratorPrecompoundInterface::PropagateModelDescription(std::ostream& outFile) const
334 {
335  outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
336  << "energy model through the wounded nucleus to precompound de-excitation.\n"
337  << "Low energy protons and neutron present among secondaries produced by \n"
338  << "the high energy generator and within the nucleus are captured. The wounded\n"
339  << "nucleus and the captured particles form an excited nuclear fragment. This\n"
340  << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
341  << "Nuclear de-excitation:\n";
342  // preco
343 
344 }
345 
346 
347 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
348 PropagateNuclNucl(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus,
349  G4V3DNucleus* theProjectileNucleus)
350 {
351 #ifdef debugPrecoInt
352  G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
353  G4cout<<"Projectile A and Z "<<theProjectileNucleus->GetMassNumber()<<" "
354  <<theProjectileNucleus->GetCharge()<<G4endl;
355  G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "
356  <<theNucleus->GetCharge()<<G4endl;
357  G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
358  G4cout<<"Projectile 4Mom and mass "<<GetPrimaryProjectile()->Get4Momentum()<<" "
359  <<GetPrimaryProjectile()->Get4Momentum().mag()<<G4endl<<G4endl;
360 #endif
361 
362  // prepare the target residual
363  G4int anA=theNucleus->GetMassNumber();
364  G4int aZ=theNucleus->GetCharge();
365  G4int numberOfEx = 0;
366  G4int numberOfCh = 0;
367  G4int numberOfHoles = 0;
368  G4double exEnergy = 0.0;
369  G4double R = theNucleus->GetNuclearRadius();
370  G4LorentzVector Target4Momentum(0.,0.,0.,0.);
371 
372  // loop over wounded target nucleus
373  G4Nucleon * theCurrentNucleon =
374  theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
375  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
376  {
377  if(theCurrentNucleon->AreYouHit()) {
378  ++numberOfHoles;
379  ++numberOfEx;
380  --anA;
381  aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
382  eplus + 0.1);
383  exEnergy += theCurrentNucleon->GetBindingEnergy();
384  Target4Momentum -=theCurrentNucleon->Get4Momentum();
385  }
386  theCurrentNucleon = theNucleus->GetNextNucleon();
387  }
388 
389 #ifdef debugPrecoInt
390  G4cout<<"Residual Target A Z E* 4mom "<<anA<<" "<<aZ<<" "<<exEnergy<<" "
391  <<Target4Momentum<<G4endl;
392 #endif
393 
394  // prepare the projectile residual
395 
396  G4bool ProjectileIsAntiNucleus=
397  GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber() < -1;
398 
399  G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();
400 
401  G4int anAb=theProjectileNucleus->GetMassNumber();
402  G4int aZb=theProjectileNucleus->GetCharge();
403  G4int numberOfExB = 0;
404  G4int numberOfChB = 0;
405  G4int numberOfHolesB = 0;
406  G4double exEnergyB = 0.0;
407  G4double Rb = theProjectileNucleus->GetNuclearRadius();
408  G4LorentzVector Projectile4Momentum(0.,0.,0.,0.);
409 
410  // loop over wounded projectile nucleus
411  theCurrentNucleon =
412  theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
413  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
414  {
415  if(theCurrentNucleon->AreYouHit()) {
416  ++numberOfHolesB;
417  ++numberOfExB;
418  --anAb;
419  if(!ProjectileIsAntiNucleus) {
420  aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
421  eplus + 0.1);
422  } else {
423  aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
424  eplus - 0.1);
425  }
426  exEnergyB += theCurrentNucleon->GetBindingEnergy();
427  Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
428  }
429  theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
430  }
431 
432  G4bool ExistTargetRemnant = G4double (numberOfHoles) <
433  0.3* G4double (numberOfHoles + anA);
434  G4bool ExistProjectileRemnant= G4double (numberOfHolesB) <
435  0.3*G4double (numberOfHolesB + anAb);
436 
437 #ifdef debugPrecoInt
438  G4cout<<"Projectile residual A Z E* 4mom "<<anAb<<" "<<aZb<<" "<<exEnergyB<<" "
439  <<Projectile4Momentum<<G4endl;
440  G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
441  <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
442 #endif
443  //-----------------------------------------------------------------------------
444  // decay the strong resonances
445  G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
446  G4DecayKineticTracks decay(theSecondaries);
447 
448  MakeCoalescence(theSecondaries);
449 
450  #ifdef debugPrecoInt
451  G4cout<<"Secondary stable particles number "<<theSecondaries->size()<<G4endl;
452  #endif
453 
454 #ifdef debugPrecoInt
455  G4LorentzVector secondary4Momemtum(0,0,0,0);
456  G4int SecondrNum(0);
457 #endif
458 
459  // loop over secondaries
460  G4KineticTrackVector::iterator iter;
461  for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
462  {
463  const G4ParticleDefinition* part = (*iter)->GetDefinition();
464  G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();
465 
466  if( part != proton && part != neutron &&
467  (part != ANTIproton && ProjectileIsAntiNucleus) &&
468  (part != ANTIneutron && ProjectileIsAntiNucleus) )
469  {
470  G4ReactionProduct * theNew = new G4ReactionProduct(part);
471  theNew->SetMomentum(aTrack4Momentum.vect());
472  theNew->SetTotalEnergy(aTrack4Momentum.e());
473  theTotalResult->push_back(theNew);
474 #ifdef debugPrecoInt
475  SecondrNum++;
476  secondary4Momemtum += (*iter)->Get4Momentum();
477  G4cout<<"Secondary "<<SecondrNum<<" "
478  <<theNew->GetDefinition()->GetParticleName()<<" "
479  <<theNew->GetMomentum()<<" "<<theNew->GetTotalEnergy()<<G4endl;
480 
481 #endif
482  delete (*iter);
483  continue;
484  }
485 
486  G4bool CanBeCapturedByTarget = false;
487  if( part == proton || part == neutron)
488  {
489  CanBeCapturedByTarget = ExistTargetRemnant &&
490  (-CaptureThreshold*G4Log( G4UniformRand()) >
491  (aTrack4Momentum + Target4Momentum).mag() -
492  aTrack4Momentum.mag() - Target4Momentum.mag()) &&
493  ((*iter)->GetPosition().mag() < R);
494  }
495  // ---------------------------
496  G4LorentzVector Position((*iter)->GetPosition(), (*iter)->GetFormationTime());
497  Position.boost(bst);
498 
499  G4bool CanBeCapturedByProjectile = false;
500 
501  if( !ProjectileIsAntiNucleus &&
502  ( part == proton || part == neutron))
503  {
504  CanBeCapturedByProjectile = ExistProjectileRemnant &&
505  (-CaptureThreshold*G4Log( G4UniformRand()) >
506  (aTrack4Momentum + Projectile4Momentum).mag() -
507  aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
508  (Position.vect().mag() < Rb);
509  }
510 
511  if( ProjectileIsAntiNucleus &&
512  ( part == ANTIproton || part == ANTIneutron))
513  {
514  CanBeCapturedByProjectile = ExistProjectileRemnant &&
515  (-CaptureThreshold*G4Log( G4UniformRand()) >
516  (aTrack4Momentum + Projectile4Momentum).mag() -
517  aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
518  (Position.vect().mag() < Rb);
519  }
520 
521  if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
522  {
523  if(G4UniformRand() < 0.5)
524  { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;}
525  else
526  { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
527  }
528 
529  if(CanBeCapturedByTarget)
530  {
531  // within the target nucleus, neutron or proton
532  // now calculate A, Z of the fragment, momentum,
533  // number of exciton states
534 #ifdef debugPrecoInt
535  G4cout<<"Track is CapturedByTarget "<<" "<<part->GetParticleName()<<" "
536  <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
537 #endif
538  ++anA;
539  ++numberOfEx;
540  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
541  aZ += Z;
542  numberOfCh += Z;
543  Target4Momentum +=aTrack4Momentum;
544  delete (*iter);
545  } else if(CanBeCapturedByProjectile)
546  {
547  // within the projectile nucleus, neutron or proton
548  // now calculate A, Z of the fragment, momentum,
549  // number of exciton states
550 #ifdef debugPrecoInt
551  G4cout<<"Track is CapturedByProjectile"<<" "<<part->GetParticleName()<<" "
552  <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
553 #endif
554  ++anAb;
555  ++numberOfExB;
556  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
557  if( ProjectileIsAntiNucleus ) Z=-Z;
558  aZb += Z;
559  numberOfChB += Z;
560  Projectile4Momentum +=aTrack4Momentum;
561  delete (*iter);
562  } else
563  { // the track is not captured
564  G4ReactionProduct * theNew = new G4ReactionProduct(part);
565  theNew->SetMomentum(aTrack4Momentum.vect());
566  theNew->SetTotalEnergy(aTrack4Momentum.e());
567  theTotalResult->push_back(theNew);
568 
569 #ifdef debugPrecoInt
570  SecondrNum++;
571  secondary4Momemtum += (*iter)->Get4Momentum();
572 /*
573  G4cout<<"Secondary "<<SecondrNum<<" "
574  <<theNew->GetDefinition()->GetParticleName()<<" "
575  <<secondary4Momemtum<<G4endl;
576 */
577 #endif
578  delete (*iter);
579  continue;
580  }
581  }
582  delete theSecondaries;
583  //-----------------------------------------------------
584 
585  #ifdef debugPrecoInt
586  G4cout<<"Final target residual A Z E* 4mom "<<anA<<" "<<aZ<<" "
587  <<exEnergy<<" "<<Target4Momentum<<G4endl;
588  #endif
589 
590  if(0!=anA )
591  {
592  G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
593 
594  if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
595  {Target4Momentum.setE(fMass);}
596 
597  G4double RemnMass=Target4Momentum.mag();
598 
599  if(RemnMass < fMass)
600  {
601  RemnMass=fMass + exEnergy;
602  Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
603  RemnMass*RemnMass));
604  } else
605  { exEnergy=RemnMass-fMass;}
606 
607  if( exEnergy < 0.) exEnergy=0.;
608 
609  // Need to de-excite the remnant nucleus
610  G4Fragment anInitialState(anA, aZ, Target4Momentum);
611  anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
612  anInitialState.SetNumberOfCharged(numberOfCh);
613  anInitialState.SetNumberOfHoles(numberOfHoles);
614 
615  G4ReactionProductVector * aPrecoResult =
616  theDeExcitation->DeExcite(anInitialState);
617 
618  #ifdef debugPrecoInt
619  G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
620  #endif
621 
622  // fill pre-compound part into the result, and return
623  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
624  {
625  theTotalResult->push_back(aPrecoResult->operator[](ll));
626  #ifdef debugPrecoInt
627  G4cout<<"Target fragment "<<ll<<" "
628  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
629  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
630  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
631  <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
632  #endif
633  }
634  delete aPrecoResult;
635  }
636 
637  //-----------------------------------------------------
638  if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
639  {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
640 
641  #ifdef debugPrecoInt
642  G4cout<<"Final projectile residual A Z E* Pmom Pmag2 "<<anAb<<" "<<aZb<<" "
643  <<exEnergyB<<" "<<Projectile4Momentum<<" "
644  <<Projectile4Momentum.mag2()<<G4endl;
645  #endif
646 
647  if(0!=anAb)
648  {
649  G4double fMass = G4NucleiProperties::GetNuclearMass(anAb, aZb);
650  G4double RemnMass=Projectile4Momentum.mag();
651 
652  if(RemnMass < fMass)
653  {
654  RemnMass=fMass + exEnergyB;
655  Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() +
656  RemnMass*RemnMass));
657  } else
658  { exEnergyB=RemnMass-fMass;}
659 
660  if( exEnergyB < 0.) exEnergyB=0.;
661 
662  G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM();
663  Projectile4Momentum.boost(bstToCM);
664 
665  // Need to de-excite the remnant nucleus
666  G4Fragment anInitialState(anAb, aZb, Projectile4Momentum);
667  anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
668  anInitialState.SetNumberOfCharged(numberOfChB);
669  anInitialState.SetNumberOfHoles(numberOfHolesB);
670 
671  G4ReactionProductVector * aPrecoResult =
672  theDeExcitation->DeExcite(anInitialState);
673 
674  #ifdef debugPrecoInt
675  G4cout<<"Projectile fragment number "<<aPrecoResult->size()<<G4endl;
676  #endif
677 
678  // fill pre-compound part into the result, and return
679  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
680  {
681  G4LorentzVector tmp=G4LorentzVector(aPrecoResult->operator[](ll)->GetMomentum(),
682  aPrecoResult->operator[](ll)->GetTotalEnergy());
683  tmp.boost(-bstToCM); // Transformation to the system of original remnant
684  aPrecoResult->operator[](ll)->SetMomentum(tmp.vect());
685  aPrecoResult->operator[](ll)->SetTotalEnergy(tmp.e());
686 
687  if(ProjectileIsAntiNucleus)
688  {
689  const G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
690  const G4ParticleDefinition * LastFragment=aFragment;
691  if (aFragment == proton) {LastFragment=G4AntiProton::AntiProtonDefinition();}
692  else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
693  else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
694  else if(aFragment == triton) {LastFragment=G4AntiTriton::AntiTritonDefinition();}
695  else if(aFragment == He3) {LastFragment=G4AntiHe3::AntiHe3Definition();}
696  else if(aFragment == He4) {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
697  else {}
698 
699  aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
700  }
701 
702  #ifdef debugPrecoInt
703  G4cout<<"Projectile fragment "<<ll<<" "
704  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
705  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
706  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
707  <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
708  #endif
709 
710  theTotalResult->push_back(aPrecoResult->operator[](ll));
711  }
712 
713  delete aPrecoResult;
714  }
715 
716  return theTotalResult;
717 }
718 
719 
720 void G4GeneratorPrecompoundInterface::MakeCoalescence(G4KineticTrackVector *tracks) {
721 
722  if (!tracks) return;
723 
724  G4double MassCut = deuteron->GetPDGMass() + DeltaM; // In MeV
725 
726  for ( size_t i = 0; i < tracks->size(); ++i ) { // search for protons
727 
728  G4KineticTrack* trackP = (*tracks)[i];
729  if ( ! trackP ) continue;
730  if (trackP->GetDefinition() != proton) continue;
731 
732  G4LorentzVector Prot4Mom = trackP->Get4Momentum();
733  G4LorentzVector ProtSPposition = G4LorentzVector(trackP->GetPosition(), trackP->GetFormationTime());
734 
735  for ( size_t j = 0; j < tracks->size(); ++j ) { // search for neutron
736 
737  G4KineticTrack* trackN = (*tracks)[j];
738  if (! trackN ) continue;
739  if (trackN->GetDefinition() != neutron) continue;
740 
741  G4LorentzVector Neut4Mom = trackN->Get4Momentum();
742  G4LorentzVector NeutSPposition = G4LorentzVector( trackN->GetPosition(), trackN->GetFormationTime()*hbarc/fermi);
743  G4double EffMass = (Prot4Mom + Neut4Mom).mag();
744 
745  if ( EffMass <= MassCut ) { // && (EffDistance <= SpaceCut)) { // Create deuteron
746  G4KineticTrack* aDeuteron =
747  new G4KineticTrack( deuteron,
748  (trackP->GetFormationTime() + trackN->GetFormationTime())/2.0,
749  (trackP->GetPosition() + trackN->GetPosition() )/2.0,
750  ( Prot4Mom + Neut4Mom ));
751  tracks->push_back(aDeuteron);
752  delete trackP; delete trackN;
753  (*tracks)[i] = nullptr; (*tracks)[j] = nullptr;
754  break;
755  }
756  }
757  }
758 
759  // Find and remove null pointers created by decays above
760  for ( int jj = tracks->size()-1; jj >= 0; --jj ) {
761  if ( ! (*tracks)[jj] ) tracks->erase(tracks->begin()+jj);
762  }
763 
764 }
765