| Geant4 Cross Reference |
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 // A.Ribon, 27 Oct 2021 -- Extended the method PropagateNuclNucl
38 // to deal with projectile hypernuclei and anti-hypernuclei
39 //
40 // -----------------------------------------------------------------------------
41
42 #include <algorithm>
43 #include <vector>
44
45 #include "G4GeneratorPrecompoundInterface.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4DynamicParticleVector.hh"
49 #include "G4KineticTrackVector.hh"
50 #include "G4Proton.hh"
51 #include "G4Neutron.hh"
52 #include "G4Lambda.hh"
53
54 #include "G4Deuteron.hh"
55 #include "G4Triton.hh"
56 #include "G4He3.hh"
57 #include "G4Alpha.hh"
58
59 #include "G4V3DNucleus.hh"
60 #include "G4Nucleon.hh"
61
62 #include "G4AntiProton.hh"
63 #include "G4AntiNeutron.hh"
64 #include "G4AntiLambda.hh"
65 #include "G4AntiDeuteron.hh"
66 #include "G4AntiTriton.hh"
67 #include "G4AntiHe3.hh"
68 #include "G4AntiAlpha.hh"
69
70 #include "G4HyperTriton.hh"
71 #include "G4HyperH4.hh"
72 #include "G4HyperAlpha.hh"
73 #include "G4HyperHe5.hh"
74 #include "G4DoubleHyperH4.hh"
75 #include "G4DoubleHyperDoubleNeutron.hh"
76
77 #include "G4AntiHyperTriton.hh"
78 #include "G4AntiHyperH4.hh"
79 #include "G4AntiHyperAlpha.hh"
80 #include "G4AntiHyperHe5.hh"
81 #include "G4AntiDoubleHyperH4.hh"
82 #include "G4AntiDoubleHyperDoubleNeutron.hh"
83
84 #include "G4FragmentVector.hh"
85 #include "G4ReactionProduct.hh"
86 #include "G4ReactionProductVector.hh"
87 #include "G4PreCompoundModel.hh"
88 #include "G4ExcitationHandler.hh"
89 #include "G4DecayKineticTracks.hh"
90 #include "G4HadronicInteractionRegistry.hh"
91
92 #include "G4PhysicsModelCatalog.hh"
93 #include "G4HyperNucleiProperties.hh"
94 //---------------------------------------------------------------------
95 #include "Randomize.hh"
96 #include "G4Log.hh"
97
98 //#define debugPrecoInt
99
100 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* preModel)
101 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), DeltaR(0.0), secID(-1)
102 {
103 proton = G4Proton::Proton();
104 neutron = G4Neutron::Neutron();
105 lambda = G4Lambda::Lambda();
106
107 deuteron=G4Deuteron::Deuteron();
108 triton =G4Triton::Triton();
109 He3 =G4He3::He3();
110 He4 =G4Alpha::Alpha();
111
112 ANTIproton=G4AntiProton::AntiProton();
113 ANTIneutron=G4AntiNeutron::AntiNeutron();
114
115 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron();
116 ANTItriton =G4AntiTriton::AntiTriton();
117 ANTIHe3 =G4AntiHe3::AntiHe3();
118 ANTIHe4 =G4AntiAlpha::AntiAlpha();
119
120 if(preModel) { SetDeExcitation(preModel); }
121 else {
122 G4HadronicInteraction* hadi =
123 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
124 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
125 if(!pre) { pre = new G4PreCompoundModel(); }
126 SetDeExcitation(pre);
127 }
128
129 secID = G4PhysicsModelCatalog::GetModelID("model_PRECO");
130 }
131
132 G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
133 {
134 }
135
136 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
137 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
138 {
139 #ifdef debugPrecoInt
140 G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::Propagate"<<G4endl;
141 G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
142 G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
143 #endif
144
145 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
146
147 // decay the strong resonances
148 G4DecayKineticTracks decay(theSecondaries);
149 #ifdef debugPrecoInt
150 G4cout<<"Final stable particles number "<<theSecondaries->size()<<G4endl;
151 #endif
152
153 // prepare the fragment (it is assumed that target nuclei are never hypernuclei)
154 G4int anA=theNucleus->GetMassNumber();
155 G4int aZ=theNucleus->GetCharge();
156 // G4double TargetNucleusMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
157
158 G4int numberOfEx = 0;
159 G4int numberOfCh = 0;
160 G4int numberOfHoles = 0;
161
162 G4double R = theNucleus->GetNuclearRadius();
163
164 G4LorentzVector captured4Momentum(0.,0.,0.,0.);
165 G4LorentzVector Residual4Momentum(0.,0.,0.,0.); // TargetNucleusMass is not need at the moment
166 G4LorentzVector Secondary4Momentum(0.,0.,0.,0.);
167
168 // loop over secondaries
169 G4KineticTrackVector::iterator iter;
170 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
171 {
172 const G4ParticleDefinition* part = (*iter)->GetDefinition();
173 G4double e = (*iter)->Get4Momentum().e();
174 G4double mass = (*iter)->Get4Momentum().mag();
175 G4ThreeVector mom = (*iter)->Get4Momentum().vect();
176 if((part != proton && part != neutron) ||
177 ((*iter)->GetPosition().mag() > R)) {
178 G4ReactionProduct * theNew = new G4ReactionProduct(part);
179 theNew->SetMomentum(mom);
180 theNew->SetTotalEnergy(e);
181 theNew->SetCreatorModelID((*iter)->GetCreatorModelID());
182 theNew->SetParentResonanceDef((*iter)->GetParentResonanceDef());
183 theNew->SetParentResonanceID((*iter)->GetParentResonanceID());
184 theTotalResult->push_back(theNew);
185 Secondary4Momentum += (*iter)->Get4Momentum();
186 #ifdef debugPrecoInt
187 G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
188 <<(*iter)->Get4Momentum().mag()<<G4endl;
189 #endif
190 } else {
191 if( e-mass > -CaptureThreshold*G4Log( G4UniformRand()) ) {
192 G4ReactionProduct * theNew = new G4ReactionProduct(part);
193 theNew->SetMomentum(mom);
194 theNew->SetTotalEnergy(e);
195 theNew->SetCreatorModelID((*iter)->GetCreatorModelID());
196 theNew->SetParentResonanceDef((*iter)->GetParentResonanceDef());
197 theNew->SetParentResonanceID((*iter)->GetParentResonanceID());
198 theTotalResult->push_back(theNew);
199 Secondary4Momentum += (*iter)->Get4Momentum();
200 #ifdef debugPrecoInt
201 G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
202 <<(*iter)->Get4Momentum().mag()<<G4endl;
203 #endif
204 } else {
205 // within the nucleus, neutron or proton
206 // now calculate A, Z of the fragment, momentum, number of exciton states
207 ++anA;
208 ++numberOfEx;
209 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
210 aZ += Z;
211 numberOfCh += Z;
212 captured4Momentum += (*iter)->Get4Momentum();
213 #ifdef debugPrecoInt
214 G4cout<<"Captured 4Mom "<<part->GetParticleName()<<(*iter)->Get4Momentum()<<G4endl;
215 #endif
216 }
217 }
218 delete (*iter);
219 }
220 delete theSecondaries;
221
222 // loop over wounded nucleus
223 G4Nucleon * theCurrentNucleon =
224 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
225 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
226 {
227 if(theCurrentNucleon->AreYouHit()) {
228 ++numberOfHoles;
229 ++numberOfEx;
230 --anA;
231 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
232
233 Residual4Momentum -= theCurrentNucleon->Get4Momentum();
234 }
235 theCurrentNucleon = theNucleus->GetNextNucleon();
236 }
237
238 #ifdef debugPrecoInt
239 G4cout<<G4endl;
240 G4cout<<"Secondary 4Mom "<<Secondary4Momentum<<G4endl;
241 G4cout<<"Captured 4Mom "<<captured4Momentum<<G4endl;
242 G4cout<<"Sec + Captured "<<Secondary4Momentum+captured4Momentum<<G4endl;
243 G4cout<<"Residual4Mom "<<Residual4Momentum<<G4endl;
244 G4cout<<"Sum 4 momenta "
245 <<Secondary4Momentum + captured4Momentum + Residual4Momentum <<G4endl;
246 #endif
247
248 // Check that we use QGS model; loop over wounded nucleus
249 G4bool QGSM(false);
250 theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
251 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
252 {
253 if(theCurrentNucleon->AreYouHit())
254 {
255 if(theCurrentNucleon->Get4Momentum().mag() <
256 theCurrentNucleon->GetDefinition()->GetPDGMass()) QGSM=true;
257 }
258 theCurrentNucleon = theNucleus->GetNextNucleon();
259 }
260
261 #ifdef debugPrecoInt
262 if(!QGSM){
263 G4cout<<G4endl;
264 G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
265 G4cout<<"Residual 4Mom "<<Residual4Momentum<<G4endl;
266 if(numberOfEx == 0)
267 {G4cout<<"Residual 4Mom = 0 means that there were not wounded and captured nucleons"<<G4endl;}
268 }
269 #endif
270
271 if(anA == 0) return theTotalResult;
272
273 G4LorentzVector exciton4Momentum(0.,0.,0.,0.);
274 if(anA >= aZ)
275 {
276 if(!QGSM)
277 { // FTF model was used
278 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
279
280 // G4LorentzVector exciton4Momentum = Residual4Momentum + captured4Momentum;
281 exciton4Momentum = Residual4Momentum + captured4Momentum;
282 //exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
283 G4double ActualMass = exciton4Momentum.mag();
284 if(ActualMass <= fMass ) {
285 exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
286 }
287
288 #ifdef debugPrecoInt
289 G4double exEnergy = 0.0;
290 if(ActualMass <= fMass ) {exEnergy = 0.;}
291 else {exEnergy = ActualMass - fMass;}
292 G4cout<<"Ground state residual Mass "<<fMass<<" E* "<<exEnergy<<G4endl;
293 #endif
294 }
295 else
296 { // QGS model was used
297 G4double InitialTargetMass =
298 G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(), theNucleus->GetCharge());
299
300 exciton4Momentum =
301 GetPrimaryProjectile()->Get4Momentum() + G4LorentzVector(0.,0.,0.,InitialTargetMass)
302 -Secondary4Momentum;
303
304 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
305 G4double ActualMass = exciton4Momentum.mag();
306
307 #ifdef debugPrecoInt
308 G4cout<<G4endl;
309 G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
310 G4cout<<"Residual4Momentum "<<exciton4Momentum<<G4endl;
311 G4cout<<"ResidualMass, GroundStateMass and E* "<<ActualMass<<" "<<fMass<<" "
312 <<ActualMass - fMass<<G4endl;
313 #endif
314
315 if(ActualMass - fMass < 0.)
316 {
317 G4double ResE = std::sqrt(exciton4Momentum.vect().mag2() + sqr(fMass+10*MeV));
318 exciton4Momentum.setE(ResE);
319 #ifdef debugPrecoInt
320 G4cout<<"ActualMass - fMass < 0. "<<ActualMass<<" "<<fMass<<" "<<ActualMass - fMass<<G4endl;
321 #endif
322 }
323 }
324
325 // Need to de-excite the remnant nucleus only if excitation energy > 0.
326 G4Fragment anInitialState(anA, aZ, exciton4Momentum);
327 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
328 anInitialState.SetNumberOfCharged(numberOfCh);
329 anInitialState.SetNumberOfHoles(numberOfHoles);
330 anInitialState.SetCreatorModelID(secID);
331
332 G4ReactionProductVector * aPrecoResult =
333 theDeExcitation->DeExcite(anInitialState);
334 // fill pre-compound part into the result, and return
335 #ifdef debugPrecoInt
336 G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
337 #endif
338 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
339 {
340 theTotalResult->push_back(aPrecoResult->operator[](ll));
341 #ifdef debugPrecoInt
342 G4cout<<"Fragment "<<ll<<" "
343 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
344 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
345 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
346 <<aPrecoResult->operator[](ll)->GetDefinition()->GetPDGMass()<<G4endl;
347 #endif
348 }
349 delete aPrecoResult;
350 }
351
352 return theTotalResult;
353 }
354
355 G4HadFinalState* G4GeneratorPrecompoundInterface::
356 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
357 {
358 G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
359 << G4endl;
360 G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
361 G4cout << "Please remove from your physics list."<<G4endl;
362 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
363 return new G4HadFinalState;
364 }
365 void G4GeneratorPrecompoundInterface::PropagateModelDescription(std::ostream& outFile) const
366 {
367 outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
368 << "energy model through the wounded nucleus to precompound de-excitation.\n"
369 << "Low energy protons and neutron present among secondaries produced by \n"
370 << "the high energy generator and within the nucleus are captured. The wounded\n"
371 << "nucleus and the captured particles form an excited nuclear fragment. This\n"
372 << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
373 << "Nuclear de-excitation:\n";
374 // preco
375
376 }
377
378
379 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
380 PropagateNuclNucl(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus,
381 G4V3DNucleus* theProjectileNucleus)
382 {
383 #ifdef debugPrecoInt
384 G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
385 G4cout<<"Projectile A and Z (and numberOfLambdas) "<<theProjectileNucleus->GetMassNumber()<<" "
386 <<theProjectileNucleus->GetCharge()<<" ("
387 <<theProjectileNucleus->GetNumberOfLambdas()<<")"<<G4endl;
388 G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "
389 <<theNucleus->GetCharge()<<" ("
390 <<theNucleus->GetNumberOfLambdas()<<")"<<G4endl;
391 G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
392 G4cout<<"Projectile 4Mom and mass "<<GetPrimaryProjectile()->Get4Momentum()<<" "
393 <<GetPrimaryProjectile()->Get4Momentum().mag()<<G4endl<<G4endl;
394 #endif
395
396 // prepare the target residual (assumed to be never a hypernucleus)
397 G4int anA=theNucleus->GetMassNumber();
398 G4int aZ=theNucleus->GetCharge();
399 //G4int aL=theNucleus->GetNumberOfLambdas(); // Should be 0
400 G4int numberOfEx = 0;
401 G4int numberOfCh = 0;
402 G4int numberOfHoles = 0;
403 G4double exEnergy = 0.0;
404 G4double R = theNucleus->GetNuclearRadius();
405 G4LorentzVector Target4Momentum(0.,0.,0.,0.);
406
407 // loop over the wounded target nucleus
408 G4Nucleon * theCurrentNucleon =
409 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
410 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
411 {
412 if(theCurrentNucleon->AreYouHit()) {
413 ++numberOfHoles;
414 ++numberOfEx;
415 --anA;
416 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
417 eplus + 0.1);
418 exEnergy += theCurrentNucleon->GetBindingEnergy();
419 Target4Momentum -=theCurrentNucleon->Get4Momentum();
420 }
421 theCurrentNucleon = theNucleus->GetNextNucleon();
422 }
423
424 #ifdef debugPrecoInt
425 G4cout<<"Residual Target A Z (numberOfLambdas) E* 4mom "<<anA<<" "<<aZ<<" (0"//<<aL
426 <<") "<<exEnergy<<" "<<Target4Momentum<<G4endl;
427 #endif
428
429 // prepare the projectile residual - which can be a hypernucleus or anti-hypernucleus
430
431 G4bool ProjectileIsAntiNucleus=
432 GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber() < -1;
433
434 G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();
435
436 G4int anAb=theProjectileNucleus->GetMassNumber();
437 G4int aZb=theProjectileNucleus->GetCharge();
438 G4int aLb=theProjectileNucleus->GetNumberOfLambdas(); // Non negative number of (anti-)lambdas in (anti-)nucleus
439 G4int numberOfExB = 0;
440 G4int numberOfChB = 0;
441 G4int numberOfHolesB = 0;
442 G4double exEnergyB = 0.0;
443 G4double Rb = theProjectileNucleus->GetNuclearRadius();
444 G4LorentzVector Projectile4Momentum(0.,0.,0.,0.);
445
446 // loop over the wounded projectile nucleus or anti-nucleus
447 theCurrentNucleon =
448 theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
449 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
450 {
451 if(theCurrentNucleon->AreYouHit()) {
452 ++numberOfHolesB;
453 ++numberOfExB;
454 --anAb;
455 if(!ProjectileIsAntiNucleus) {
456 aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
457 eplus + 0.1);
458 if (theCurrentNucleon->GetParticleType()==G4Lambda::Definition()) --aLb;
459 } else {
460 aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
461 eplus - 0.1);
462 if (theCurrentNucleon->GetParticleType()==G4AntiLambda::Definition()) --aLb;
463 }
464 exEnergyB += theCurrentNucleon->GetBindingEnergy();
465 Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
466 }
467 theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
468 }
469
470 G4bool ExistTargetRemnant = G4double (numberOfHoles) <
471 0.3* G4double (numberOfHoles + anA);
472 G4bool ExistProjectileRemnant= G4double (numberOfHolesB) <
473 0.3*G4double (numberOfHolesB + anAb);
474
475 #ifdef debugPrecoInt
476 G4cout<<"Projectile residual A Z (numberOfLambdas) E* 4mom "<<anAb<<" "<<aZb<<" ("<<aLb
477 <<") "<<exEnergyB<<" "<<Projectile4Momentum<<G4endl;
478 G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
479 <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
480 #endif
481 //-----------------------------------------------------------------------------
482 // decay the strong resonances
483 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
484 G4DecayKineticTracks decay(theSecondaries);
485
486 MakeCoalescence(theSecondaries);
487
488 #ifdef debugPrecoInt
489 G4cout<<"Secondary stable particles number "<<theSecondaries->size()<<G4endl;
490 #endif
491
492 #ifdef debugPrecoInt
493 G4LorentzVector secondary4Momemtum(0,0,0,0);
494 G4int SecondrNum(0);
495 #endif
496
497 // Loop over secondaries.
498 // We are assuming that only protons and neutrons - for nuclei -
499 // and only antiprotons and antineutrons - for antinuclei - can be absorbed,
500 // not instead lambdas (or hyperons more generally) - for nuclei - or anti-lambdas
501 // (or anti-hyperons more generally) - for antinuclei. This is a simplification,
502 // to be eventually reviewed later on, in particular when generic hypernuclei and
503 // anti-hypernuclei are introduced, instead of the few light hypernuclei and
504 // anti-hypernuclei which currently exist in Geant4.
505 G4KineticTrackVector::iterator iter;
506 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
507 {
508 const G4ParticleDefinition* part = (*iter)->GetDefinition();
509 G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();
510
511 if( part != proton && part != neutron &&
512 (part != ANTIproton && ProjectileIsAntiNucleus) &&
513 (part != ANTIneutron && ProjectileIsAntiNucleus) )
514 {
515 G4ReactionProduct * theNew = new G4ReactionProduct(part);
516 theNew->SetMomentum(aTrack4Momentum.vect());
517 theNew->SetTotalEnergy(aTrack4Momentum.e());
518 theNew->SetCreatorModelID((*iter)->GetCreatorModelID());
519 theNew->SetParentResonanceDef((*iter)->GetParentResonanceDef());
520 theNew->SetParentResonanceID((*iter)->GetParentResonanceID());
521 theTotalResult->push_back(theNew);
522 #ifdef debugPrecoInt
523 SecondrNum++;
524 secondary4Momemtum += (*iter)->Get4Momentum();
525 G4cout<<"Secondary "<<SecondrNum<<" "
526 <<theNew->GetDefinition()->GetParticleName()<<" "
527 <<theNew->GetMomentum()<<" "<<theNew->GetTotalEnergy()<<G4endl;
528
529 #endif
530 delete (*iter);
531 continue;
532 }
533
534 G4bool CanBeCapturedByTarget = false;
535 if( part == proton || part == neutron)
536 {
537 CanBeCapturedByTarget = ExistTargetRemnant &&
538 (-CaptureThreshold*G4Log( G4UniformRand()) >
539 (aTrack4Momentum + Target4Momentum).mag() -
540 aTrack4Momentum.mag() - Target4Momentum.mag()) &&
541 ((*iter)->GetPosition().mag() < R);
542 }
543 // ---------------------------
544 G4LorentzVector Position((*iter)->GetPosition(), (*iter)->GetFormationTime());
545 Position.boost(bst);
546
547 G4bool CanBeCapturedByProjectile = false;
548
549 if( !ProjectileIsAntiNucleus &&
550 ( part == proton || part == neutron))
551 {
552 CanBeCapturedByProjectile = ExistProjectileRemnant &&
553 (-CaptureThreshold*G4Log( G4UniformRand()) >
554 (aTrack4Momentum + Projectile4Momentum).mag() -
555 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
556 (Position.vect().mag() < Rb);
557 }
558
559 if( ProjectileIsAntiNucleus &&
560 ( part == ANTIproton || part == ANTIneutron))
561 {
562 CanBeCapturedByProjectile = ExistProjectileRemnant &&
563 (-CaptureThreshold*G4Log( G4UniformRand()) >
564 (aTrack4Momentum + Projectile4Momentum).mag() -
565 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
566 (Position.vect().mag() < Rb);
567 }
568
569 if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
570 {
571 if(G4UniformRand() < 0.5)
572 { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;}
573 else
574 { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
575 }
576
577 if(CanBeCapturedByTarget)
578 {
579 // within the target nucleus, neutron or proton
580 // now calculate A, Z of the fragment, momentum,
581 // number of exciton states
582 #ifdef debugPrecoInt
583 G4cout<<"Track is CapturedByTarget "<<" "<<part->GetParticleName()<<" "
584 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
585 #endif
586 ++anA;
587 ++numberOfEx;
588 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
589 aZ += Z;
590 numberOfCh += Z;
591 Target4Momentum +=aTrack4Momentum;
592 delete (*iter);
593 } else if(CanBeCapturedByProjectile)
594 {
595 // within the projectile nucleus, neutron or proton
596 // now calculate A, Z of the fragment, momentum,
597 // number of exciton states
598 #ifdef debugPrecoInt
599 G4cout<<"Track is CapturedByProjectile"<<" "<<part->GetParticleName()<<" "
600 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
601 #endif
602 ++anAb;
603 ++numberOfExB;
604 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
605 if( ProjectileIsAntiNucleus ) Z=-Z;
606 aZb += Z;
607 numberOfChB += Z;
608 Projectile4Momentum +=aTrack4Momentum;
609 delete (*iter);
610 } else
611 { // the track is not captured
612 G4ReactionProduct * theNew = new G4ReactionProduct(part);
613 theNew->SetMomentum(aTrack4Momentum.vect());
614 theNew->SetTotalEnergy(aTrack4Momentum.e());
615 theNew->SetCreatorModelID((*iter)->GetCreatorModelID());
616 theNew->SetParentResonanceDef((*iter)->GetParentResonanceDef());
617 theNew->SetParentResonanceID((*iter)->GetParentResonanceID());
618 theTotalResult->push_back(theNew);
619
620 #ifdef debugPrecoInt
621 SecondrNum++;
622 secondary4Momemtum += (*iter)->Get4Momentum();
623 /*
624 G4cout<<"Secondary "<<SecondrNum<<" "
625 <<theNew->GetDefinition()->GetParticleName()<<" "
626 <<secondary4Momemtum<<G4endl;
627 */
628 #endif
629 delete (*iter);
630 continue;
631 }
632 }
633 delete theSecondaries;
634 //-----------------------------------------------------
635
636 #ifdef debugPrecoInt
637 G4cout<<"Final target residual A Z (numberOfLambdas) E* 4mom "<<anA<<" "<<aZ<<" (0"//<<aL
638 <<") "<<exEnergy<<" "<<Target4Momentum<<G4endl;
639 #endif
640
641 if(0!=anA )
642 {
643 // We assume that the target residual is never a hypernucleus
644 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
645
646 if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
647 {Target4Momentum.setE(fMass);}
648
649 G4double RemnMass=Target4Momentum.mag();
650
651 if(RemnMass < fMass)
652 {
653 RemnMass=fMass + exEnergy;
654 Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
655 RemnMass*RemnMass));
656 } else
657 { exEnergy=RemnMass-fMass;}
658
659 if( exEnergy < 0.) exEnergy=0.;
660
661 // Need to de-excite the remnant nucleus
662 G4Fragment anInitialState(anA, aZ, Target4Momentum);
663 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
664 anInitialState.SetNumberOfCharged(numberOfCh);
665 anInitialState.SetNumberOfHoles(numberOfHoles);
666 anInitialState.SetCreatorModelID(secID);
667
668 G4ReactionProductVector * aPrecoResult =
669 theDeExcitation->DeExcite(anInitialState);
670
671 #ifdef debugPrecoInt
672 G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
673 #endif
674
675 // fill pre-compound part into the result, and return
676 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
677 {
678 theTotalResult->push_back(aPrecoResult->operator[](ll));
679 #ifdef debugPrecoInt
680 G4cout<<"Target fragment "<<ll<<" "
681 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
682 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
683 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
684 <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
685 #endif
686 }
687 delete aPrecoResult;
688 }
689
690 //-----------------------------------------------------
691 if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
692 {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
693
694 #ifdef debugPrecoInt
695 G4cout<<"Final projectile residual A Z (numberOfLambdas) E* Pmom Pmag2 "<<anAb<<" "<<aZb<<" ("
696 <<aLb<<") "<<exEnergyB<<" "<<Projectile4Momentum<<" "
697 <<Projectile4Momentum.mag2()<<G4endl;
698 #endif
699
700 if(0!=anAb)
701 {
702 // The projectile residual can be a hypernucleus or anti-hypernucleus
703 G4double fMass = 0.0;
704 if ( aLb > 0 ) {
705 fMass = G4HyperNucleiProperties::GetNuclearMass(anAb, aZb, aLb);
706 } else {
707 fMass = G4NucleiProperties::GetNuclearMass(anAb, aZb);
708 }
709 G4double RemnMass=Projectile4Momentum.mag();
710
711 if(RemnMass < fMass)
712 {
713 RemnMass=fMass + exEnergyB;
714 Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() +
715 RemnMass*RemnMass));
716 } else
717 { exEnergyB=RemnMass-fMass;}
718
719 if( exEnergyB < 0.) exEnergyB=0.;
720
721 G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM();
722 Projectile4Momentum.boost(bstToCM);
723
724 // Need to de-excite the remnant nucleus
725 G4Fragment anInitialState(anAb, aZb, aLb, Projectile4Momentum);
726 anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
727 anInitialState.SetNumberOfCharged(numberOfChB);
728 anInitialState.SetNumberOfHoles(numberOfHolesB);
729 anInitialState.SetCreatorModelID(secID);
730
731 G4ReactionProductVector * aPrecoResult =
732 theDeExcitation->DeExcite(anInitialState);
733
734 #ifdef debugPrecoInt
735 G4cout<<"Projectile fragment number "<<aPrecoResult->size()<<G4endl;
736 #endif
737
738 // fill pre-compound part into the result, and return
739 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
740 {
741 G4LorentzVector tmp=G4LorentzVector(aPrecoResult->operator[](ll)->GetMomentum(),
742 aPrecoResult->operator[](ll)->GetTotalEnergy());
743 tmp.boost(-bstToCM); // Transformation to the system of original remnant
744 aPrecoResult->operator[](ll)->SetMomentum(tmp.vect());
745 aPrecoResult->operator[](ll)->SetTotalEnergy(tmp.e());
746
747 if(ProjectileIsAntiNucleus)
748 {
749 const G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
750 const G4ParticleDefinition * LastFragment=aFragment;
751 if (aFragment == proton) {LastFragment=G4AntiProton::AntiProtonDefinition();}
752 else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
753 else if(aFragment == lambda) {LastFragment=G4AntiLambda::AntiLambdaDefinition();}
754 else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
755 else if(aFragment == triton) {LastFragment=G4AntiTriton::AntiTritonDefinition();}
756 else if(aFragment == He3) {LastFragment=G4AntiHe3::AntiHe3Definition();}
757 else if(aFragment == He4) {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
758 else {}
759
760 if (aLb > 0) { // Anti-hypernucleus
761 if (aFragment == G4HyperTriton::Definition()) {
762 LastFragment=G4AntiHyperTriton::Definition();
763 } else if (aFragment == G4HyperH4::Definition()) {
764 LastFragment=G4AntiHyperH4::Definition();
765 } else if (aFragment == G4HyperAlpha::Definition()) {
766 LastFragment=G4AntiHyperAlpha::Definition();
767 } else if (aFragment == G4HyperHe5::Definition()) {
768 LastFragment=G4AntiHyperHe5::Definition();
769 } else if (aFragment == G4DoubleHyperH4::Definition()) {
770 LastFragment=G4AntiDoubleHyperH4::Definition();
771 } else if (aFragment == G4DoubleHyperDoubleNeutron::Definition()) {
772 LastFragment=G4AntiDoubleHyperDoubleNeutron::Definition();
773 }
774 }
775
776 aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
777 }
778
779 #ifdef debugPrecoInt
780 G4cout<<"Projectile fragment "<<ll<<" "
781 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
782 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
783 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
784 <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
785 #endif
786
787 theTotalResult->push_back(aPrecoResult->operator[](ll));
788 }
789
790 delete aPrecoResult;
791 }
792
793 return theTotalResult;
794 }
795
796
797 void G4GeneratorPrecompoundInterface::MakeCoalescence(G4KineticTrackVector *tracks) {
798 // This method replaces pairs of proton-neutron - in the case of nuclei - or
799 // antiproton-antineutron - in the case of anti-nuclei - which are close in
800 // momentum, with, respectively, deuterons and anti-deuterons.
801 // Note that in the case of hypernuclei or anti-hypernuclei, lambdas or anti-lambdas
802 // are not considered for coalescence because hyper-deuteron or anti-hyper-deuteron
803 // are assumed not to exist.
804
805 if (!tracks) return;
806
807 G4double MassCut = deuteron->GetPDGMass() + DeltaM; // In MeV
808
809 for ( std::size_t i = 0; i < tracks->size(); ++i ) { // search for protons
810
811 G4KineticTrack* trackP = (*tracks)[i];
812 if ( ! trackP ) continue;
813 if (trackP->GetDefinition() != proton) continue;
814
815 G4LorentzVector Prot4Mom = trackP->Get4Momentum();
816 G4LorentzVector ProtSPposition = G4LorentzVector(trackP->GetPosition(), trackP->GetFormationTime());
817
818 for ( std::size_t j = 0; j < tracks->size(); ++j ) { // search for neutron
819
820 G4KineticTrack* trackN = (*tracks)[j];
821 if (! trackN ) continue;
822 if (trackN->GetDefinition() != neutron) continue;
823
824 G4LorentzVector Neut4Mom = trackN->Get4Momentum();
825 G4LorentzVector NeutSPposition = G4LorentzVector( trackN->GetPosition(), trackN->GetFormationTime()*hbarc/fermi);
826 G4double EffMass = (Prot4Mom + Neut4Mom).mag();
827
828 if ( EffMass <= MassCut ) { // && (EffDistance <= SpaceCut)) { // Create deuteron
829 G4KineticTrack* aDeuteron =
830 new G4KineticTrack( deuteron,
831 (trackP->GetFormationTime() + trackN->GetFormationTime())/2.0,
832 (trackP->GetPosition() + trackN->GetPosition() )/2.0,
833 ( Prot4Mom + Neut4Mom ));
834 aDeuteron->SetCreatorModelID(secID);
835 tracks->push_back(aDeuteron);
836 delete trackP; delete trackN;
837 (*tracks)[i] = nullptr; (*tracks)[j] = nullptr;
838 break;
839 }
840 }
841 }
842
843 // Find and remove null pointers created by decays above
844 for ( G4int jj = (G4int)tracks->size()-1; jj >= 0; --jj ) {
845 if ( ! (*tracks)[jj] ) tracks->erase(tracks->begin()+jj);
846 }
847 }
848