| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // <<
27 // ------------------------------------------- <<
28 // GEANT 4 class file <<
29 // <<
30 // History: first implementation <<
31 // HPW, 10DEC 98, the decay part original <<
32 // in his FTF-test-program. <<
33 // <<
34 // M.Kelsey, 28 Jul 2011 -- Replace loop <<
35 // with new utility class, simplify cleanup <<
36 // <<
37 // A.Ribon, 27 Oct 2021 -- Extended the m <<
38 // to deal with projectile hyper <<
39 // <<
40 // ------------------------------------------- <<
41 <<
42 #include <algorithm> <<
43 #include <vector> <<
44 <<
45 #include "G4GeneratorPrecompoundInterface.hh" 23 #include "G4GeneratorPrecompoundInterface.hh"
46 #include "G4PhysicalConstants.hh" <<
47 #include "G4SystemOfUnits.hh" <<
48 #include "G4DynamicParticleVector.hh" 24 #include "G4DynamicParticleVector.hh"
49 #include "G4KineticTrackVector.hh" << 25 #include "G4IonTable.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::G4GeneratorPr <<
101 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), <<
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 26
112 ANTIproton=G4AntiProton::AntiProton(); << 27 //
113 ANTIneutron=G4AntiNeutron::AntiNeutron(); << 28 // HPW, 10DEC 98, the decay part originally written by Gunter Folger in his FTF-test-program.
114 << 29 //
115 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron() << 30
116 ANTItriton =G4AntiTriton::AntiTriton(); << 31 G4HadFinalState* G4GeneratorPrecompoundInterface::
117 ANTIHe3 =G4AntiHe3::AntiHe3(); << 32 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
118 ANTIHe4 =G4AntiAlpha::AntiAlpha(); <<
119 <<
120 if(preModel) { SetDeExcitation(preModel); } <<
121 else { <<
122 G4HadronicInteraction* hadi = <<
123 G4HadronicInteractionRegistry::Ins <<
124 G4VPreCompoundModel* pre = static_cast<G <<
125 if(!pre) { pre = new G4PreCompoundModel( <<
126 SetDeExcitation(pre); <<
127 } <<
128 <<
129 secID = G4PhysicsModelCatalog::GetModelID(" <<
130 } <<
131 <<
132 G4GeneratorPrecompoundInterface::~G4GeneratorP <<
133 { <<
134 } <<
135 <<
136 G4ReactionProductVector* G4GeneratorPrecompoun <<
137 Propagate(G4KineticTrackVector* theSecondaries <<
138 { <<
139 #ifdef debugPrecoInt <<
140 G4cout<<G4endl<<"G4GeneratorPrecompoundI <<
141 G4cout<<"Target A and Z "<<theNucleus->G <<
142 G4cout<<"Directly produced particles num <<
143 #endif <<
144 <<
145 G4ReactionProductVector * theTotalResult = <<
146 <<
147 // decay the strong resonances <<
148 G4DecayKineticTracks decay(theSecondaries); <<
149 #ifdef debugPrecoInt <<
150 G4cout<<"Final stable particles number " <<
151 #endif <<
152 <<
153 // prepare the fragment (it is assumed that <<
154 G4int anA=theNucleus->GetMassNumber(); <<
155 G4int aZ=theNucleus->GetCharge(); <<
156 // G4double TargetNucleusMass = G4NucleiPrope <<
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., <<
165 G4LorentzVector Residual4Momentum(0.,0.,0., <<
166 G4LorentzVector Secondary4Momentum(0.,0.,0. <<
167 <<
168 // loop over secondaries <<
169 G4KineticTrackVector::iterator iter; <<
170 for(iter=theSecondaries->begin(); iter !=th <<
171 { 33 {
172 const G4ParticleDefinition* part = (*ite << 34 std::cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."<< G4endl;
173 G4double e = (*iter)->Get4Momentum().e() << 35 std::cout << "This class is only a mediator between generator and precompound"<<G4endl;
174 G4double mass = (*iter)->Get4Momentum(). << 36 std::cout << "Please remove from your physics list."<<G4endl;
175 G4ThreeVector mom = (*iter)->Get4Momentu << 37 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
176 if((part != proton && part != neutron) | << 38 return new G4HadFinalState;
177 ((*iter)->GetPosition().mag() > R) << 39 }
178 G4ReactionProduct * theNew = new G4Re << 40
179 theNew->SetMomentum(mom); << 41 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
180 theNew->SetTotalEnergy(e); << 42 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
181 theNew->SetCreatorModelID((*iter)->GetCreat <<
182 theNew->SetParentResonanceDef((*iter)->GetP <<
183 theNew->SetParentResonanceID((*iter)->GetPa <<
184 theTotalResult->push_back(theNew); <<
185 Secondary4Momentum += (*iter)->Get4Mo <<
186 #ifdef debugPrecoInt <<
187 G4cout<<"Secondary 4Mom "<<part->G <<
188 <<(*iter)->Get4Momentum().ma <<
189 #endif <<
190 } else { <<
191 if( e-mass > -CaptureThreshold*G4Log( <<
192 G4ReactionProduct * theNew = new G <<
193 theNew->SetMomentum(mom); <<
194 theNew->SetTotalEnergy(e); <<
195 theNew->SetCreatorModelID((*iter)->GetCr <<
196 theNew->SetParentResonanceDef((*iter)->G <<
197 theNew->SetParentResonanceID((*iter)->Ge <<
198 theTotalResult->push_back(theNew); <<
199 Secondary4Momentum += (*iter)->Get <<
200 #ifdef debugPrecoInt <<
201 G4cout<<"Secondary 4Mom "<<part <<
202 <<(*iter)->Get4Momentum() <<
203 #endif <<
204 } else { <<
205 // within the nucleus, neutron or <<
206 // now calculate A, Z of the frag <<
207 ++anA; <<
208 ++numberOfEx; <<
209 G4int Z = G4int(part->GetPDGCharge <<
210 aZ += Z; <<
211 numberOfCh += Z; <<
212 captured4Momentum += (*iter)->Get4 <<
213 #ifdef debugPrecoInt <<
214 G4cout<<"Captured 4Mom "<<part <<
215 #endif <<
216 } <<
217 } <<
218 delete (*iter); <<
219 } <<
220 delete theSecondaries; <<
221 <<
222 // loop over wounded nucleus <<
223 G4Nucleon * theCurrentNucleon = <<
224 theNucleus->StartLoop() ? theNucleus- <<
225 while(theCurrentNucleon) /* Loop checkin <<
226 { <<
227 if(theCurrentNucleon->AreYouHit()) { <<
228 ++numberOfHoles; <<
229 ++numberOfEx; <<
230 --anA; <<
231 aZ -= G4int(theCurrentNucleon->GetDef <<
232 <<
233 Residual4Momentum -= theCurrentNucleo <<
234 } <<
235 theCurrentNucleon = theNucleus->GetNextN <<
236 } <<
237 <<
238 #ifdef debugPrecoInt <<
239 G4cout<<G4endl; <<
240 G4cout<<"Secondary 4Mom "<<Secondary4Mom <<
241 G4cout<<"Captured 4Mom "<<captured4Mome <<
242 G4cout<<"Sec + Captured "<<Secondary4Mom <<
243 G4cout<<"Residual4Mom "<<Residual4Mome <<
244 G4cout<<"Sum 4 momenta " <<
245 <<Secondary4Momentum + captured4Mo <<
246 #endif <<
247 <<
248 // Check that we use QGS model; loop over w <<
249 G4bool QGSM(false); <<
250 theCurrentNucleon = theNucleus->StartLoop() <<
251 while(theCurrentNucleon) /* Loop checking <<
252 { <<
253 if(theCurrentNucleon->AreYouHit()) <<
254 { <<
255 if(theCurrentNucleon->Get4Momentum().ma <<
256 theCurrentNucleon->GetDefinition()-> <<
257 } <<
258 theCurrentNucleon = theNucleus->GetNextN <<
259 } <<
260 <<
261 #ifdef debugPrecoInt <<
262 if(!QGSM){ <<
263 G4cout<<G4endl; <<
264 G4cout<<"Residual A and Z "<<anA<<" "< <<
265 G4cout<<"Residual 4Mom "<<Residual4Mo <<
266 if(numberOfEx == 0) <<
267 {G4cout<<"Residual 4Mom = 0 means tha <<
268 } <<
269 #endif <<
270 <<
271 if(anA == 0) return theTotalResult; <<
272 <<
273 G4LorentzVector exciton4Momentum(0.,0.,0.,0 <<
274 if(anA >= aZ) <<
275 { 43 {
276 if(!QGSM) << 44 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
277 { // FTF model was used <<
278 G4double fMass = G4NucleiProperties::Ge <<
279 <<
280 // G4LorentzVector exciton4Momentum = Res <<
281 exciton4Momentum = Residual4Momentum + c <<
282 //exciton4Momentum.setE(std::sqrt(exciton4Mome <<
283 G4double ActualMass = exciton4Momentum.m <<
284 if(ActualMass <= fMass ) { <<
285 exciton4Momentum.setE(std::sqrt(excito <<
286 } <<
287 <<
288 #ifdef debugPrecoInt <<
289 G4double exEnergy = 0.0; <<
290 if(ActualMass <= fMass ) {exEnergy = <<
291 else {exEnergy = <<
292 G4cout<<"Ground state residual Mass " <<
293 #endif <<
294 } <<
295 else <<
296 { // QGS model was used <<
297 G4double InitialTargetMass = <<
298 G4NucleiProperties::GetNuclearMa <<
299 <<
300 exciton4Momentum = <<
301 GetPrimaryProjectile()->Get4Mome <<
302 -Secondary4Momentum; <<
303 45
304 G4double fMass = G4NucleiProperties::Get << 46 // decay the strong resonances
305 G4double ActualMass = exciton4Momentum.ma << 47 G4KineticTrackVector *result1, *secondaries, *result;
306 << 48 result1=theSecondaries;
307 #ifdef debugPrecoInt << 49 result=new G4KineticTrackVector();
308 G4cout<<G4endl; << 50
309 G4cout<<"Residual A and Z "<<anA<<" "< << 51 for (unsigned int aResult=0; aResult < result1->size(); aResult++)
310 G4cout<<"Residual4Momentum "<<exciton4 <<
311 G4cout<<"ResidualMass, GroundStateMass <<
312 <<ActualMass - fMass<<G4endl; <<
313 #endif <<
314 <<
315 if(ActualMass - fMass < 0.) <<
316 { 52 {
317 G4double ResE = std::sqrt(exciton4Moment << 53 G4ParticleDefinition * pdef;
318 exciton4Momentum.setE(ResE); << 54 pdef=result1->operator[](aResult)->GetDefinition();
319 #ifdef debugPrecoInt << 55 secondaries=NULL;
320 G4cout<<"ActualMass - fMass < 0. "<<A << 56 if ( pdef->IsShortLived() )
321 #endif << 57 {
>> 58 secondaries = result1->operator[](aResult)->Decay();
>> 59 }
>> 60 if ( secondaries == NULL )
>> 61 {
>> 62 result->push_back(result1->operator[](aResult));
>> 63 result1->operator[](aResult)=NULL; //protect for clearAndDestroy
>> 64 }
>> 65 else
>> 66 {
>> 67 for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
>> 68 {
>> 69 result1->push_back(secondaries->operator[](aSecondary));
>> 70 }
>> 71 delete secondaries;
>> 72 }
322 } 73 }
323 } << 74 std::for_each(result1->begin(), result1->end(), DeleteKineticTrack());
324 << 75 delete result1;
325 // Need to de-excite the remnant nucleus o << 76
326 G4Fragment anInitialState(anA, aZ, exciton << 77
327 anInitialState.SetNumberOfParticles(number << 78
328 anInitialState.SetNumberOfCharged(numberOf << 79 // prepare the fragment
329 anInitialState.SetNumberOfHoles(numberOfHo << 80 G4Fragment anInitialState;
330 anInitialState.SetCreatorModelID(secID); << 81 G4int anA=theNucleus->GetMassNumber();
331 << 82 G4int aZ=theNucleus->GetCharge();
332 G4ReactionProductVector * aPrecoResult = << 83 G4int numberOfEx = 0;
333 theDeExcitation->DeExcite(anInitialState << 84 G4int numberOfCh = 0;
334 // fill pre-compound part into the result, << 85 G4int numberOfHoles = 0;
335 #ifdef debugPrecoInt << 86 G4double exEnergy = 0;
336 G4cout<<"Target fragment number "<<aPrecoR << 87 G4ThreeVector exciton3Momentum(0,0,0);
337 #endif << 88 // loop over secondaries
338 for(unsigned int ll=0; ll<aPrecoResult->si << 89 for(unsigned int list=0; list < result->size(); list++)
339 { << 90 {
340 theTotalResult->push_back(aPrecoResult-> << 91 G4KineticTrack *aTrack = result->operator[](list);
341 #ifdef debugPrecoInt << 92 if(aTrack->GetDefinition() != G4Proton::Proton() &&
342 G4cout<<"Fragment "<<ll<<" " << 93 aTrack->GetDefinition() != G4Neutron::Neutron())
343 <<aPrecoResult->operator[](ll)->GetDefin << 94 {
344 <<aPrecoResult->operator[](ll)->GetMomen << 95 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
345 <<aPrecoResult->operator[](ll)->GetTotal << 96 theNew->SetMomentum(aTrack->Get4Momentum().vect());
346 <<aPrecoResult->operator[](ll)->GetDefin << 97 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
347 #endif << 98 theTotalResult->push_back(theNew);
348 } << 99 }
349 delete aPrecoResult; << 100 else if(aTrack->Get4Momentum().t() - aTrack->Get4Momentum().mag()>80*MeV)
350 } << 101 {
351 << 102 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
352 return theTotalResult; << 103 theNew->SetMomentum(aTrack->Get4Momentum().vect());
353 } << 104 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
354 << 105 theTotalResult->push_back(theNew);
355 G4HadFinalState* G4GeneratorPrecompoundInterfa << 106 }
356 ApplyYourself(const G4HadProjectile &, G4Nucle << 107 else if(aTrack->GetPosition().mag() > theNucleus->GetNuclearRadius())
357 { << 108 {
358 G4cout << "G4GeneratorPrecompoundInterface: << 109 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
359 << G4endl; << 110 theNew->SetMomentum(aTrack->Get4Momentum().vect());
360 G4cout << "This class is only a mediator be << 111 theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
361 G4cout << "Please remove from your physics << 112 theTotalResult->push_back(theNew);
362 throw G4HadronicException(__FILE__, __LINE_ << 113 }
363 return new G4HadFinalState; << 114 else
364 } << 115 {
365 void G4GeneratorPrecompoundInterface::Propagat << 116 // within the nucleus, neutron or proton
366 { << 117 // now calculate A, Z of the fragment, momentum, number of exciton states
367 outFile << "G4GeneratorPrecompoundInterface << 118 anA++;;
368 << "energy model through the wounded << 119 numberOfEx++;
369 << "Low energy protons and neutron pr << 120 aZ += G4int(aTrack->GetDefinition()->GetPDGCharge());
370 << "the high energy generator and wit << 121 numberOfCh += G4int(aTrack->GetDefinition()->GetPDGCharge());
371 << "nucleus and the captured particle << 122 exciton3Momentum += aTrack->Get4Momentum().vect();
372 << "fragment is passed to the Geant4 << 123 exEnergy += (aTrack->Get4Momentum().t()-aTrack->Get4Momentum().m());
373 << "Nuclear de-excitation:\n"; << 124 }
374 // preco << 125 }
375 << 126
376 } << 127 // loop over wounded nucleus
377 << 128 G4Nucleon * theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : NULL;
378 << 129 while(theCurrentNucleon != NULL)
379 G4ReactionProductVector* G4GeneratorPrecompoun << 130 {
380 PropagateNuclNucl(G4KineticTrackVector* theSec << 131 if(theCurrentNucleon->AreYouHit())
381 G4V3DNucleus* theProjectileNucleus) << 132 {
382 { << 133 numberOfHoles++;
383 #ifdef debugPrecoInt << 134 numberOfEx++;
384 G4cout<<G4endl<<"G4GeneratorPrecompoundInte << 135 anA--;
385 G4cout<<"Projectile A and Z (and numberOfLa << 136 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge());
386 <<theProjectileNucleus->GetCharge()<<" (" << 137 exciton3Momentum -= theCurrentNucleon->Get4Momentum().vect();
387 <<theProjectileNucleus->GetNumberOfLambdas( << 138 exEnergy+=theCurrentNucleon->GetBindingEnergy();
388 G4cout<<"Target A and Z "<<theNucleus-> << 139 }
389 <<theNucleus->GetCharge()<<" (" << 140 theCurrentNucleon = theNucleus->GetNextNucleon();
390 <<theNucleus->GetNumberOfLambdas()<<")"<<G4 << 141 }
391 G4cout<<"Directly produced particles number << 142
392 G4cout<<"Projectile 4Mom and mass "<<GetPri << 143 if(!theDeExcitation)
393 <<GetPri << 144 {
394 #endif << 145 // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
395 << 146 }
396 // prepare the target residual (assumed to << 147 else if(0!=anA && 0!=aZ)
397 G4int anA=theNucleus->GetMassNumber(); << 148 {
398 G4int aZ=theNucleus->GetCharge(); << 149 G4double residualMass =
399 //G4int aL=theNucleus->GetNumberOfLambdas() << 150 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(aZ ,anA);
400 G4int numberOfEx = 0; << 151 residualMass += exEnergy;
401 G4int numberOfCh = 0; << 152 G4LorentzVector exciton4Momentum(exciton3Momentum,
402 G4int numberOfHoles = 0; << 153 std::sqrt(exciton3Momentum.mag2()+residualMass*residualMass));
403 G4double exEnergy = 0.0; << 154
404 G4double R = theNucleus->GetNuclearRadius() << 155 anInitialState.SetA(anA);
405 G4LorentzVector Target4Momentum(0.,0.,0.,0. << 156 anInitialState.SetZ(aZ);
406 << 157 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
407 // loop over the wounded target nucleus << 158 anInitialState.SetNumberOfCharged(numberOfCh);
408 G4Nucleon * theCurrentNucleon = << 159 anInitialState.SetNumberOfHoles(numberOfHoles);
409 theNucleus->StartLoop() ? theNucleus- << 160 anInitialState.SetMomentum(exciton4Momentum);
410 while(theCurrentNucleon) /* Loop checking << 161 // anInitialState.SetExcitationEnergy(exEnergy); // now a redundant call.
411 { << 162
412 if(theCurrentNucleon->AreYouHit()) { << 163 // call pre-compound
413 ++numberOfHoles; << 164 const G4Fragment aFragment(anInitialState);
414 ++numberOfEx; << 165 G4ReactionProductVector * aPreResult = theDeExcitation->DeExcite(aFragment);
415 --anA; << 166 // G4ReactionProductVector * aPreResult = new G4ReactionProductVector;
416 aZ -= G4int(theCurrentNucleon->GetDef << 167 // fill pre-compound part into the result, and return
417 eplus + 0.1); << 168 for(unsigned int ll=0; ll<aPreResult->size(); ll++)
418 exEnergy += theCurrentNucleon->GetBin << 169 {
419 Target4Momentum -=theCurrentNucleon-> << 170 theTotalResult->push_back(aPreResult->operator[](ll));
420 } << 171 }
421 theCurrentNucleon = theNucleus->GetNextN << 172 delete aPreResult;
422 } << 173 }
423 << 174 else
424 #ifdef debugPrecoInt << 175 {
425 G4cout<<"Residual Target A Z (numberOfLambd << 176 // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
426 <<") "<<exEnergy<<" "<<Target4Momentu << 177 }
427 #endif << 178 // now return
428 << 179
429 // prepare the projectile residual - which << 180 std::for_each(result->begin(), result->end(), DeleteKineticTrack());
430 << 181 delete result;
431 G4bool ProjectileIsAntiNucleus= << 182 return theTotalResult;
432 GetPrimaryProjectile()->GetDefinition <<
433 <<
434 G4ThreeVector bst = GetPrimaryProjectile()- <<
435 <<
436 G4int anAb=theProjectileNucleus->GetMassNum <<
437 G4int aZb=theProjectileNucleus->GetCharge() <<
438 G4int aLb=theProjectileNucleus->GetNumberOf <<
439 G4int numberOfExB = 0; <<
440 G4int numberOfChB = 0; <<
441 G4int numberOfHolesB = 0; <<
442 G4double exEnergyB = 0.0; <<
443 G4double Rb = theProjectileNucleus->GetNucl <<
444 G4LorentzVector Projectile4Momentum(0.,0.,0 <<
445 <<
446 // loop over the wounded projectile nucleus <<
447 theCurrentNucleon = <<
448 theProjectileNucleus->StartLoop() ? t <<
449 while(theCurrentNucleon) /* Loop checkin <<
450 { <<
451 if(theCurrentNucleon->AreYouHit()) { <<
452 ++numberOfHolesB; <<
453 ++numberOfExB; <<
454 --anAb; <<
455 if(!ProjectileIsAntiNucleus) { <<
456 aZb -= G4int(theCurrentNucleon->Ge <<
457 eplus + 0.1); <<
458 if (theCurrentNucleon->GetParticleType() <<
459 } else { <<
460 aZb += G4int(theCurrentNucleon->Ge <<
461 eplus - 0.1); <<
462 if (theCurrentNucleon->GetParticleType() <<
463 } <<
464 exEnergyB += theCurrentNucleon->GetBi <<
465 Projectile4Momentum -=theCurrentNucle <<
466 } <<
467 theCurrentNucleon = theProjectileNucleus <<
468 } <<
469 <<
470 G4bool ExistTargetRemnant = G4double (nu <<
471 0.3* G4double (numberOfHoles + anA); <<
472 G4bool ExistProjectileRemnant= G4double (nu <<
473 0.3*G4double (numberOfHolesB + anAb); <<
474 <<
475 #ifdef debugPrecoInt <<
476 G4cout<<"Projectile residual A Z (numberOfL <<
477 <<") "<<exEnergyB<<" "<<Projectile4Momentum <<
478 G4cout<<" ExistTargetRemnant ExistProjectil <<
479 <<ExistTargetRemnant<<" "<< ExistProj <<
480 #endif <<
481 //----------------------------------------- <<
482 // decay the strong resonances <<
483 G4ReactionProductVector * theTotalResult = <<
484 G4DecayKineticTracks decay(theSecondaries); <<
485 <<
486 MakeCoalescence(theSecondaries); <<
487 <<
488 #ifdef debugPrecoInt <<
489 G4cout<<"Secondary stable particles numb <<
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 ne <<
499 // and only antiprotons and antineutrons - <<
500 // not instead lambdas (or hyperons more ge <<
501 // (or anti-hyperons more generally) - for <<
502 // to be eventually reviewed later on, in p <<
503 // anti-hypernuclei are introduced, instead <<
504 // anti-hypernuclei which currently exist i <<
505 G4KineticTrackVector::iterator iter; <<
506 for(iter=theSecondaries->begin(); iter !=th <<
507 { <<
508 const G4ParticleDefinition* part = (*ite <<
509 G4LorentzVector aTrack4Momentum=(*iter)- <<
510 <<
511 if( part != proton && part != neutron && <<
512 (part != ANTIproton && Projectile <<
513 (part != ANTIneutron && Projectile <<
514 { <<
515 G4ReactionProduct * theNew = new G4Re <<
516 theNew->SetMomentum(aTrack4Momentum.v <<
517 theNew->SetTotalEnergy(aTrack4Momentu <<
518 theNew->SetCreatorModelID((*iter)->GetCreat <<
519 theNew->SetParentResonanceDef((*iter)->GetP <<
520 theNew->SetParentResonanceID((*iter)->GetPa <<
521 theTotalResult->push_back(theNew); <<
522 #ifdef debugPrecoInt <<
523 SecondrNum++; <<
524 secondary4Momemtum += (*iter)->Get4Mo <<
525 G4cout<<"Secondary "<<SecondrNum<<" <<
526 <<theNew->GetDefinition()->GetP <<
527 <<theNew->GetMomentum()<<" "<<t <<
528 <<
529 #endif <<
530 delete (*iter); <<
531 continue; <<
532 } <<
533 <<
534 G4bool CanBeCapturedByTarget = false; <<
535 if( part == proton || part == neutron) <<
536 { <<
537 CanBeCapturedByTarget = ExistTargetRe <<
538 (-CaptureThreshold*G4Log( G4Uni <<
539 (aTrack4Momentum + Target4 <<
540 aTrack4Momentum.mag() - Target4 <<
541 ((*iter)->GetPosition().mag <<
542 } <<
543 // --------------------------- <<
544 G4LorentzVector Position((*iter)->GetPos <<
545 Position.boost(bst); <<
546 <<
547 G4bool CanBeCapturedByProjectile = false <<
548 <<
549 if( !ProjectileIsAntiNucleus && <<
550 ( part == proton || part == neutro <<
551 { <<
552 CanBeCapturedByProjectile = ExistProj <<
553 (-CaptureThreshold*G4Log( G4Uni <<
554 (aTrack4Momentum + Project <<
555 aTrack4Momentum.mag() - Project <<
556 (Position.vect().mag() < Rb <<
557 } <<
558 <<
559 if( ProjectileIsAntiNucleus && <<
560 ( part == ANTIproton || part == AN <<
561 { <<
562 CanBeCapturedByProjectile = ExistProj <<
563 (-CaptureThreshold*G4Log( G4Uni <<
564 (aTrack4Momentum + Project <<
565 aTrack4Momentum.mag() - Project <<
566 (Position.vect().mag() < Rb <<
567 } <<
568 <<
569 if(CanBeCapturedByTarget && CanBeCapture <<
570 { <<
571 if(G4UniformRand() < 0.5) <<
572 { CanBeCapturedByTarget = true; CanBe <<
573 else <<
574 { CanBeCapturedByTarget = false; CanB <<
575 } <<
576 <<
577 if(CanBeCapturedByTarget) <<
578 { <<
579 // within the target nucleus, neutron <<
580 // now calculate A, Z of the fragmen <<
581 // number of exciton states <<
582 #ifdef debugPrecoInt <<
583 G4cout<<"Track is CapturedByTarget "< <<
584 <<aTrack4Momentum<<" "<<aTrack4 <<
585 #endif <<
586 ++anA; <<
587 ++numberOfEx; <<
588 G4int Z = G4int(part->GetPDGCharge()/ <<
589 aZ += Z; <<
590 numberOfCh += Z; <<
591 Target4Momentum +=aTrack4Momentum; <<
592 delete (*iter); <<
593 } else if(CanBeCapturedByProjectile) <<
594 { <<
595 // within the projectile nucleus, neu <<
596 // now calculate A, Z of the fragmen <<
597 // number of exciton states <<
598 #ifdef debugPrecoInt <<
599 G4cout<<"Track is CapturedByProjectil <<
600 <<aTrack4Momentum<<" "<<aTrack4 <<
601 #endif <<
602 ++anAb; <<
603 ++numberOfExB; <<
604 G4int Z = G4int(part->GetPDGCharge()/ <<
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 G4Re <<
613 theNew->SetMomentum(aTrack4Momentum.v <<
614 theNew->SetTotalEnergy(aTrack4Momentu <<
615 theNew->SetCreatorModelID((*iter)->GetCreat <<
616 theNew->SetParentResonanceDef((*iter)->GetP <<
617 theNew->SetParentResonanceID((*iter)->GetPa <<
618 theTotalResult->push_back(theNew); <<
619 <<
620 #ifdef debugPrecoInt <<
621 SecondrNum++; <<
622 secondary4Momemtum += (*iter)->Get4Mo <<
623 /* <<
624 G4cout<<"Secondary "<<SecondrNum<<" <<
625 <<theNew->GetDefinition()->GetP <<
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 (numberO <<
638 <<") "<<exEnergy<<" "<<Target4Momentum<< <<
639 #endif <<
640 <<
641 if(0!=anA ) <<
642 { <<
643 // We assume that the target residual is <<
644 G4double fMass = G4NucleiProperties::Ge <<
645 <<
646 if((anA == theNucleus->GetMassNumber()) <<
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(Target <<
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, Targe <<
663 anInitialState.SetNumberOfParticles(numb <<
664 anInitialState.SetNumberOfCharged(number <<
665 anInitialState.SetNumberOfHoles(numberOf <<
666 anInitialState.SetCreatorModelID(secID); <<
667 <<
668 G4ReactionProductVector * aPrecoResult = <<
669 theDeExcitation->DeExcite(anInitia <<
670 <<
671 #ifdef debugPrecoInt <<
672 G4cout<<"Target fragment number "<<aP <<
673 #endif <<
674 <<
675 // fill pre-compound part into the resul <<
676 for(unsigned int ll=0; ll<aPrecoResult-> <<
677 { <<
678 theTotalResult->push_back(aPrecoResul <<
679 #ifdef debugPrecoInt <<
680 G4cout<<"Target fragment "<<ll<<" <<
681 <<aPrecoResult->operator[](l <<
682 <<aPrecoResult->operator[](l <<
683 <<aPrecoResult->operator[](l <<
684 <<aPrecoResult->operator[](l <<
685 #endif <<
686 } <<
687 delete aPrecoResult; <<
688 } <<
689 <<
690 //----------------------------------------- <<
691 if((anAb == theProjectileNucleus->GetMassNu <<
692 {Projectile4Momentum = GetPrimaryProjectile <<
693 <<
694 #ifdef debugPrecoInt <<
695 G4cout<<"Final projectile residual A Z (num <<
696 <<aLb<<") "<<exEnergyB<<" "<<Projectile4Mom <<
697 <<Projectile4Momen <<
698 #endif <<
699 <<
700 if(0!=anAb) <<
701 { <<
702 // The projectile residual can be a hype <<
703 G4double fMass = 0.0; <<
704 if ( aLb > 0 ) { <<
705 fMass = G4HyperNucleiProperties::GetNuclearM <<
706 } else { <<
707 fMass = G4NucleiProperties::GetNuclearMass(a <<
708 } <<
709 G4double RemnMass=Projectile4Momentum.ma <<
710 <<
711 if(RemnMass < fMass) <<
712 { <<
713 RemnMass=fMass + exEnergyB; <<
714 Projectile4Momentum.setE(std::sqrt(Pr <<
715 RemnMass*RemnMass)); <<
716 } else <<
717 { exEnergyB=RemnMass-fMass;} <<
718 <<
719 if( exEnergyB < 0.) exEnergyB=0.; <<
720 <<
721 G4ThreeVector bstToCM =Projectile4Moment <<
722 Projectile4Momentum.boost(bstToCM); <<
723 <<
724 // Need to de-excite the remnant nucleus <<
725 G4Fragment anInitialState(anAb, aZb, aLb <<
726 anInitialState.SetNumberOfParticles(numb <<
727 anInitialState.SetNumberOfCharged(number <<
728 anInitialState.SetNumberOfHoles(numberOf <<
729 anInitialState.SetCreatorModelID(secID); <<
730 <<
731 G4ReactionProductVector * aPrecoResult = <<
732 theDeExcitation->DeExcite(anInitia <<
733 <<
734 #ifdef debugPrecoInt <<
735 G4cout<<"Projectile fragment number " <<
736 #endif <<
737 <<
738 // fill pre-compound part into the resul <<
739 for(unsigned int ll=0; ll<aPrecoResult-> <<
740 { <<
741 G4LorentzVector tmp=G4LorentzVector(a <<
742 a <<
743 tmp.boost(-bstToCM); // Transformatio <<
744 aPrecoResult->operator[](ll)->SetMome <<
745 aPrecoResult->operator[](ll)->SetTota <<
746 <<
747 if(ProjectileIsAntiNucleus) <<
748 { <<
749 const G4ParticleDefinition * aFrag <<
750 const G4ParticleDefinition * LastF <<
751 if (aFragment == proton) {Las <<
752 else if(aFragment == neutron) {Las <<
753 else if(aFragment == lambda) {Las <<
754 else if(aFragment == deuteron){Las <<
755 else if(aFragment == triton) {Las <<
756 else if(aFragment == He3) {Las <<
757 else if(aFragment == He4) {Las <<
758 else {} <<
759 <<
760 if (aLb > 0) { // Anti-hypernucle <<
761 if (aFragment == G4HyperTriton: <<
762 LastFragment=G4AntiHyperTriton::Definition <<
763 } else if (aFragment == G4HyperH4::Def <<
764 LastFragment=G4AntiHyperH4::Definition(); <<
765 } else if (aFragment == G4HyperAlpha:: <<
766 LastFragment=G4AntiHyperAlpha::Definition( <<
767 } else if (aFragment == G4HyperHe5::De <<
768 LastFragment=G4AntiHyperHe5::Definition(); <<
769 } else if (aFragment == G4DoubleHyperH <<
770 LastFragment=G4AntiDoubleHyperH4::Definiti <<
771 } else if (aFragment == G4DoubleHyperD <<
772 LastFragment=G4AntiDoubleHyperDoubleNeutro <<
773 } <<
774 } <<
775 <<
776 aPrecoResult->operator[](ll)->SetD <<
777 } <<
778 <<
779 #ifdef debugPrecoInt <<
780 G4cout<<"Projectile fragment "<<ll <<
781 <<aPrecoResult->operator[](l <<
782 <<aPrecoResult->operator[](l <<
783 <<aPrecoResult->operator[](l <<
784 <<aPrecoResult->operator[](l <<
785 #endif <<
786 <<
787 theTotalResult->push_back(aPrecoResul <<
788 } <<
789 <<
790 delete aPrecoResult; <<
791 } 183 }
792 <<
793 return theTotalResult; <<
794 } <<
795 <<
796 <<
797 void G4GeneratorPrecompoundInterface::MakeCoal <<
798 // This method replaces pairs of proton-neut <<
799 // antiproton-antineutron - in the case of a <<
800 // momentum, with, respectively, deuterons a <<
801 // Note that in the case of hypernuclei or a <<
802 // are not considered for coalescence becaus <<
803 // are assumed not to exist. <<
804 184
805 if (!tracks) return; <<
806 <<
807 G4double MassCut = deuteron->GetPDGMass() + <<
808 <<
809 for ( std::size_t i = 0; i < tracks->size(); <<
810 <<
811 G4KineticTrack* trackP = (*tracks)[i]; <<
812 if ( ! trackP ) continue; <<
813 if (trackP->GetDefinition() != proton) con <<
814 <<
815 G4LorentzVector Prot4Mom = trackP->Get4Mom <<
816 G4LorentzVector ProtSPposition = G4Lorentz <<
817 <<
818 for ( std::size_t j = 0; j < tracks->size( <<
819 <<
820 G4KineticTrack* trackN = (*tracks)[j]; <<
821 if (! trackN ) continue; <<
822 if (trackN->GetDefinition() != neutron) <<
823 <<
824 G4LorentzVector Neut4Mom = trackN->Get4M <<
825 G4LorentzVector NeutSPposition = G4Loren <<
826 G4double EffMass = (Prot4Mom + Neut4Mom) <<
827 <<
828 if ( EffMass <= MassCut ) { // && (EffD <<
829 G4KineticTrack* aDeuteron = <<
830 new G4KineticTrack( deuteron, <<
831 (trackP->GetForm <<
832 (trackP->GetPosi <<
833 ( Prot4Mom <<
834 aDeuteron->SetCreatorModelID(secID); <<
835 tracks->push_back(aDeuteron); <<
836 delete trackP; delete trackN; <<
837 (*tracks)[i] = nullptr; (*tracks)[j] = <<
838 break; <<
839 } <<
840 } <<
841 } <<
842 <<
843 // Find and remove null pointers created by <<
844 for ( G4int jj = (G4int)tracks->size()-1; jj <<
845 if ( ! (*tracks)[jj] ) tracks->erase(track <<
846 } <<
847 } <<
848 185