| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4GeneratorPrecompoundInterface.cc 66812 2013-01-12 16:06:46Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -----------------------------------------------------------------------------
28 // GEANT 4 class file 29 // GEANT 4 class file
29 // 30 //
30 // History: first implementation 31 // History: first implementation
31 // HPW, 10DEC 98, the decay part original 32 // HPW, 10DEC 98, the decay part originally written by Gunter Folger
32 // in his FTF-test-program. 33 // in his FTF-test-program.
33 // 34 //
34 // M.Kelsey, 28 Jul 2011 -- Replace loop 35 // M.Kelsey, 28 Jul 2011 -- Replace loop to decay input secondaries
35 // with new utility class, simplify cleanup 36 // with new utility class, simplify cleanup loops
36 // <<
37 // A.Ribon, 27 Oct 2021 -- Extended the m <<
38 // to deal with projectile hyper <<
39 // <<
40 // ------------------------------------------- 37 // -----------------------------------------------------------------------------
41 38
42 #include <algorithm> 39 #include <algorithm>
43 #include <vector> 40 #include <vector>
44 41
45 #include "G4GeneratorPrecompoundInterface.hh" 42 #include "G4GeneratorPrecompoundInterface.hh"
46 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
48 #include "G4DynamicParticleVector.hh" 45 #include "G4DynamicParticleVector.hh"
49 #include "G4KineticTrackVector.hh" 46 #include "G4KineticTrackVector.hh"
50 #include "G4Proton.hh" 47 #include "G4Proton.hh"
51 #include "G4Neutron.hh" 48 #include "G4Neutron.hh"
52 #include "G4Lambda.hh" <<
53 49
54 #include "G4Deuteron.hh" 50 #include "G4Deuteron.hh"
55 #include "G4Triton.hh" 51 #include "G4Triton.hh"
56 #include "G4He3.hh" 52 #include "G4He3.hh"
57 #include "G4Alpha.hh" 53 #include "G4Alpha.hh"
58 54
59 #include "G4V3DNucleus.hh" 55 #include "G4V3DNucleus.hh"
60 #include "G4Nucleon.hh" 56 #include "G4Nucleon.hh"
61 57
62 #include "G4AntiProton.hh" 58 #include "G4AntiProton.hh"
63 #include "G4AntiNeutron.hh" 59 #include "G4AntiNeutron.hh"
64 #include "G4AntiLambda.hh" <<
65 #include "G4AntiDeuteron.hh" 60 #include "G4AntiDeuteron.hh"
66 #include "G4AntiTriton.hh" 61 #include "G4AntiTriton.hh"
67 #include "G4AntiHe3.hh" 62 #include "G4AntiHe3.hh"
68 #include "G4AntiAlpha.hh" 63 #include "G4AntiAlpha.hh"
69 64
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" 65 #include "G4FragmentVector.hh"
85 #include "G4ReactionProduct.hh" 66 #include "G4ReactionProduct.hh"
86 #include "G4ReactionProductVector.hh" 67 #include "G4ReactionProductVector.hh"
87 #include "G4PreCompoundModel.hh" 68 #include "G4PreCompoundModel.hh"
88 #include "G4ExcitationHandler.hh" 69 #include "G4ExcitationHandler.hh"
89 #include "G4DecayKineticTracks.hh" 70 #include "G4DecayKineticTracks.hh"
90 #include "G4HadronicInteractionRegistry.hh" 71 #include "G4HadronicInteractionRegistry.hh"
91 72
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 73 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* preModel)
101 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), << 74 : CaptureThreshold(10*MeV)
102 { 75 {
103 proton = G4Proton::Proton(); << 76 proton = G4Proton::Proton();
104 neutron = G4Neutron::Neutron(); << 77 neutron = G4Neutron::Neutron();
105 lambda = G4Lambda::Lambda(); << 78
106 << 79 deuteron=G4Deuteron::Deuteron();
107 deuteron=G4Deuteron::Deuteron(); << 80 triton =G4Triton::Triton();
108 triton =G4Triton::Triton(); << 81 He3 =G4He3::He3();
109 He3 =G4He3::He3(); << 82 He4 =G4Alpha::Alpha();
110 He4 =G4Alpha::Alpha(); << 83
111 << 84 ANTIproton=G4AntiProton::AntiProton();
112 ANTIproton=G4AntiProton::AntiProton(); << 85 ANTIneutron=G4AntiNeutron::AntiNeutron();
113 ANTIneutron=G4AntiNeutron::AntiNeutron(); << 86
114 << 87 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron();
115 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron() << 88 ANTItriton =G4AntiTriton::AntiTriton();
116 ANTItriton =G4AntiTriton::AntiTriton(); << 89 ANTIHe3 =G4AntiHe3::AntiHe3();
117 ANTIHe3 =G4AntiHe3::AntiHe3(); << 90 ANTIHe4 =G4AntiAlpha::AntiAlpha();
118 ANTIHe4 =G4AntiAlpha::AntiAlpha(); << 91
119 << 92 if(preModel) { SetDeExcitation(preModel); }
120 if(preModel) { SetDeExcitation(preModel); } << 93 else {
121 else { << 94 G4HadronicInteraction* hadi =
122 G4HadronicInteraction* hadi = << 95 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
123 G4HadronicInteractionRegistry::Ins <<
124 G4VPreCompoundModel* pre = static_cast<G 96 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
125 if(!pre) { pre = new G4PreCompoundModel( 97 if(!pre) { pre = new G4PreCompoundModel(); }
126 SetDeExcitation(pre); << 98 SetDeExcitation(pre);
127 } << 99 }
128 <<
129 secID = G4PhysicsModelCatalog::GetModelID(" <<
130 } 100 }
131 101
132 G4GeneratorPrecompoundInterface::~G4GeneratorP 102 G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
133 { 103 {
134 } 104 }
135 105
>> 106 //---------------------------------------------------------------------
>> 107 // choose to calculate excitation energy from energy balance
>> 108 #define exactExcitationEnergy
>> 109 //#define debugPrecoInt
>> 110 //#define G4GPI_debug_excitation
>> 111
136 G4ReactionProductVector* G4GeneratorPrecompoun 112 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
137 Propagate(G4KineticTrackVector* theSecondaries 113 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
138 { 114 {
139 #ifdef debugPrecoInt << 115 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
140 G4cout<<G4endl<<"G4GeneratorPrecompoundI << 116
141 G4cout<<"Target A and Z "<<theNucleus->G << 117 // decay the strong resonances
142 G4cout<<"Directly produced particles num << 118 G4DecayKineticTracks decay(theSecondaries);
143 #endif << 119
144 << 120 // prepare the fragment
145 G4ReactionProductVector * theTotalResult = << 121 G4int anA=theNucleus->GetMassNumber();
146 << 122 G4int aZ=theNucleus->GetCharge();
147 // decay the strong resonances << 123 G4int numberOfEx = 0;
148 G4DecayKineticTracks decay(theSecondaries); << 124 G4int numberOfCh = 0;
149 #ifdef debugPrecoInt << 125 G4int numberOfHoles = 0;
150 G4cout<<"Final stable particles number " << 126 G4double exEnergy = 0.0;
151 #endif << 127 G4double R = theNucleus->GetNuclearRadius();
152 << 128 G4ThreeVector exciton3Momentum(0.,0.,0.);
153 // prepare the fragment (it is assumed that << 129 G4ThreeVector captured3Momentum(0.,0.,0.);
154 G4int anA=theNucleus->GetMassNumber(); << 130 G4ThreeVector wounded3Momentum(0.,0.,0.);
155 G4int aZ=theNucleus->GetCharge(); << 131
156 // G4double TargetNucleusMass = G4NucleiPrope << 132 // loop over secondaries
157 << 133 #ifdef exactExcitationEnergy
158 G4int numberOfEx = 0; << 134 G4LorentzVector secondary4Momemtum(0,0,0,0);
159 G4int numberOfCh = 0; << 135 #endif
160 G4int numberOfHoles = 0; << 136 G4KineticTrackVector::iterator iter;
161 << 137 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
162 G4double R = theNucleus->GetNuclearRadius() << 138 {
163 << 139 G4ParticleDefinition* part = (*iter)->GetDefinition();
164 G4LorentzVector captured4Momentum(0.,0.,0., << 140 G4double e = (*iter)->Get4Momentum().e();
165 G4LorentzVector Residual4Momentum(0.,0.,0., << 141 G4double mass = (*iter)->Get4Momentum().mag();
166 G4LorentzVector Secondary4Momentum(0.,0.,0. << 142 G4ThreeVector mom = (*iter)->Get4Momentum().vect();
167 << 143 if((part != proton && part != neutron) ||
168 // loop over secondaries << 144 (e > mass + CaptureThreshold) ||
169 G4KineticTrackVector::iterator iter; << 145 ((*iter)->GetPosition().mag() > R)) {
170 for(iter=theSecondaries->begin(); iter !=th <<
171 { <<
172 const G4ParticleDefinition* part = (*ite <<
173 G4double e = (*iter)->Get4Momentum().e() <<
174 G4double mass = (*iter)->Get4Momentum(). <<
175 G4ThreeVector mom = (*iter)->Get4Momentu <<
176 if((part != proton && part != neutron) | <<
177 ((*iter)->GetPosition().mag() > R) <<
178 G4ReactionProduct * theNew = new G4Re <<
179 theNew->SetMomentum(mom); <<
180 theNew->SetTotalEnergy(e); <<
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 146 G4ReactionProduct * theNew = new G4ReactionProduct(part);
193 theNew->SetMomentum(mom); 147 theNew->SetMomentum(mom);
194 theNew->SetTotalEnergy(e); 148 theNew->SetTotalEnergy(e);
195 theNew->SetCreatorModelID((*iter)->GetCr <<
196 theNew->SetParentResonanceDef((*iter)->G <<
197 theNew->SetParentResonanceID((*iter)->Ge <<
198 theTotalResult->push_back(theNew); 149 theTotalResult->push_back(theNew);
199 Secondary4Momentum += (*iter)->Get << 150 #ifdef exactExcitationEnergy
200 #ifdef debugPrecoInt << 151 secondary4Momemtum += (*iter)->Get4Momentum();
201 G4cout<<"Secondary 4Mom "<<part << 152 #endif
202 <<(*iter)->Get4Momentum() << 153 } else {
203 #endif <<
204 } else { <<
205 // within the nucleus, neutron or 154 // within the nucleus, neutron or proton
206 // now calculate A, Z of the frag 155 // now calculate A, Z of the fragment, momentum, number of exciton states
207 ++anA; 156 ++anA;
208 ++numberOfEx; 157 ++numberOfEx;
209 G4int Z = G4int(part->GetPDGCharge 158 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
210 aZ += Z; 159 aZ += Z;
211 numberOfCh += Z; 160 numberOfCh += Z;
212 captured4Momentum += (*iter)->Get4 << 161 captured3Momentum += mom;
213 #ifdef debugPrecoInt << 162 exEnergy += (e - mass);
214 G4cout<<"Captured 4Mom "<<part << 163 }
215 #endif << 164 delete (*iter);
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 { <<
276 if(!QGSM) <<
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 } 165 }
295 else << 166 delete theSecondaries;
296 { // QGS model was used << 167
297 G4double InitialTargetMass = << 168 // loop over wounded nucleus
298 G4NucleiProperties::GetNuclearMa << 169 G4Nucleon * theCurrentNucleon =
299 << 170 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
300 exciton4Momentum = << 171 while(theCurrentNucleon) {
301 GetPrimaryProjectile()->Get4Mome << 172 if(theCurrentNucleon->AreYouHit()) {
302 -Secondary4Momentum; << 173 ++numberOfHoles;
303 << 174 ++numberOfEx;
304 G4double fMass = G4NucleiProperties::Get << 175 --anA;
305 G4double ActualMass = exciton4Momentum.ma << 176 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
306 << 177 wounded3Momentum += theCurrentNucleon->Get4Momentum().vect();
307 #ifdef debugPrecoInt << 178 //G4cout << "hit nucleon " << theCurrentNucleon->Get4Momentum() << G4endl;
308 G4cout<<G4endl; << 179 exEnergy += theCurrentNucleon->GetBindingEnergy();
309 G4cout<<"Residual A and Z "<<anA<<" "< << 180 }
310 G4cout<<"Residual4Momentum "<<exciton4 << 181 theCurrentNucleon = theNucleus->GetNextNucleon();
311 G4cout<<"ResidualMass, GroundStateMass <<
312 <<ActualMass - fMass<<G4endl; <<
313 #endif <<
314 <<
315 if(ActualMass - fMass < 0.) <<
316 { <<
317 G4double ResE = std::sqrt(exciton4Moment <<
318 exciton4Momentum.setE(ResE); <<
319 #ifdef debugPrecoInt <<
320 G4cout<<"ActualMass - fMass < 0. "<<A <<
321 #endif <<
322 } <<
323 } 182 }
>> 183 exciton3Momentum = captured3Momentum - wounded3Momentum;
>> 184
>> 185 if(anA == 0) return theTotalResult;
324 186
325 // Need to de-excite the remnant nucleus o << 187 if(anA >= aZ)
326 G4Fragment anInitialState(anA, aZ, exciton <<
327 anInitialState.SetNumberOfParticles(number <<
328 anInitialState.SetNumberOfCharged(numberOf <<
329 anInitialState.SetNumberOfHoles(numberOfHo <<
330 anInitialState.SetCreatorModelID(secID); <<
331 <<
332 G4ReactionProductVector * aPrecoResult = <<
333 theDeExcitation->DeExcite(anInitialState <<
334 // fill pre-compound part into the result, <<
335 #ifdef debugPrecoInt <<
336 G4cout<<"Target fragment number "<<aPrecoR <<
337 #endif <<
338 for(unsigned int ll=0; ll<aPrecoResult->si <<
339 { 188 {
340 theTotalResult->push_back(aPrecoResult-> << 189 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
341 #ifdef debugPrecoInt << 190
342 G4cout<<"Fragment "<<ll<<" " << 191 #ifdef exactExcitationEnergy
343 <<aPrecoResult->operator[](ll)->GetDefin << 192 // recalculate exEnergy from Energy balance....
344 <<aPrecoResult->operator[](ll)->GetMomen << 193 const G4HadProjectile * primary = GetPrimaryProjectile();
345 <<aPrecoResult->operator[](ll)->GetTotal << 194 G4double Einitial= primary->Get4Momentum().e()
346 <<aPrecoResult->operator[](ll)->GetDefin << 195 + G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(),
347 #endif << 196 theNucleus->GetCharge());
>> 197 // Uzhi G4double Efinal = fMass + secondary4Momemtum.e();
>> 198 G4double Efinal = std::sqrt(exciton3Momentum.mag2() + fMass*fMass)
>> 199 + secondary4Momemtum.e();
>> 200 if ( (Einitial - Efinal) > 0 ) {
>> 201 // G4cout << "G4GPI::Propagate() : positive exact excitation Energy "
>> 202 // << (Einitial - Efinal)/MeV << " MeV, exciton estimate "
>> 203 // << exEnergy/MeV << " MeV" << G4endl;
>> 204
>> 205 // exEnergy=Einitial - Efinal;
>> 206 G4LorentzVector PrimMom=primary->Get4Momentum(); PrimMom.setE(Einitial);
>> 207
>> 208 exEnergy=(PrimMom - secondary4Momemtum).mag() - fMass;
>> 209 }
>> 210 else {
>> 211 // G4cout << "G4GeneratorPrecompoundInterface::Propagate() : "
>> 212 // << "negative exact excitation Energy "
>> 213 // << (Einitial - Efinal)/MeV
>> 214 // << " MeV, setting excitation to 0 MeV" << G4endl;
>> 215 exEnergy=0.;
>> 216 }
>> 217 #endif
>> 218
>> 219 if(exEnergy < 0.) exEnergy=0.; // Uzhi 11 Dec. 2012
>> 220
>> 221 fMass += exEnergy;
>> 222
>> 223 G4ThreeVector balance=primary->Get4Momentum().vect() -
>> 224 secondary4Momemtum.vect() - exciton3Momentum;
>> 225
>> 226 #ifdef G4GPI_debug_excitation
>> 227 G4cout << "momentum balance" << balance
>> 228 << " value " << balance.mag() <<G4endl
>> 229 << "primary "<< primary->Get4Momentum() <<G4endl
>> 230 << "secondary "<< secondary4Momemtum <<G4endl
>> 231 << "captured "<< captured3Momentum <<G4endl
>> 232 << "wounded "<< wounded3Momentum <<G4endl
>> 233 << "exciton "<< exciton3Momentum <<G4endl
>> 234 << "second + exciton"
>> 235 << secondary4Momemtum.vect() + exciton3Momentum << G4endl;
>> 236 #endif
>> 237 //#ifdef exactExcitationEnergy
>> 238 // G4LorentzVector exciton4Momentum(exciton3Momentum, fMass);
>> 239 // G4LorentzVector exciton4Momentum(exciton3Momentum,
>> 240 // std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
>> 241 //#else
>> 242 G4LorentzVector exciton4Momentum(exciton3Momentum,
>> 243 std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
>> 244 //#endif
>> 245 //G4cout<<"exciton4Momentum "<<exciton4Momentum<<G4endl;
>> 246 // Need to de-excite the remnant nucleus only if excitation energy > 0.
>> 247 G4Fragment anInitialState(anA, aZ, exciton4Momentum);
>> 248 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
>> 249 anInitialState.SetNumberOfCharged(numberOfCh);
>> 250 anInitialState.SetNumberOfHoles(numberOfHoles);
>> 251
>> 252 G4ReactionProductVector * aPrecoResult =
>> 253 theDeExcitation->DeExcite(anInitialState);
>> 254 // fill pre-compound part into the result, and return
>> 255 theTotalResult->insert(theTotalResult->end(),aPrecoResult->begin(),
>> 256 aPrecoResult->end() );
>> 257 delete aPrecoResult;
348 } 258 }
349 delete aPrecoResult; <<
350 } <<
351 259
352 return theTotalResult; << 260 return theTotalResult;
353 } 261 }
354 262
355 G4HadFinalState* G4GeneratorPrecompoundInterfa 263 G4HadFinalState* G4GeneratorPrecompoundInterface::
356 ApplyYourself(const G4HadProjectile &, G4Nucle 264 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
357 { 265 {
358 G4cout << "G4GeneratorPrecompoundInterface: << 266 G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
359 << G4endl; << 267 << G4endl;
360 G4cout << "This class is only a mediator be << 268 G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
361 G4cout << "Please remove from your physics << 269 G4cout << "Please remove from your physics list."<<G4endl;
362 throw G4HadronicException(__FILE__, __LINE_ << 270 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
363 return new G4HadFinalState; << 271 return new G4HadFinalState;
364 } 272 }
365 void G4GeneratorPrecompoundInterface::Propagat 273 void G4GeneratorPrecompoundInterface::PropagateModelDescription(std::ostream& outFile) const
366 { 274 {
367 outFile << "G4GeneratorPrecompoundInterface << 275 outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
368 << "energy model through the wounded << 276 << "energy model through the wounded nucleus to precompound de-excition.\n"
369 << "Low energy protons and neutron pr << 277 << "Low energy protons and neutron present among secondaries produced by \n"
370 << "the high energy generator and wit << 278 << "the high energy generator and within the nucleus are captured. The wounded\n"
371 << "nucleus and the captured particle << 279 << "nucleus and the captured particles form an excited nuclear fragment. This\n"
372 << "fragment is passed to the Geant4 << 280 << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
373 << "Nuclear de-excitation:\n"; << 281 << "Nuclear de-excitation:\n";
374 // preco << 282 // preco
375 283
376 } 284 }
377 285
378 << 286 // Uzhi Nov. 2012 ------------------------------------------------
379 G4ReactionProductVector* G4GeneratorPrecompoun 287 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
380 PropagateNuclNucl(G4KineticTrackVector* theSec 288 PropagateNuclNucl(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus,
381 G4V3DNucleus* theProjectileNucleus) << 289 G4V3DNucleus* theProjectileNucleus)
382 { 290 {
383 #ifdef debugPrecoInt 291 #ifdef debugPrecoInt
384 G4cout<<G4endl<<"G4GeneratorPrecompoundInte << 292 G4cout<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
385 G4cout<<"Projectile A and Z (and numberOfLa << 293 #endif
386 <<theProjectileNucleus->GetCharge()<<" (" << 294
387 <<theProjectileNucleus->GetNumberOfLambdas( << 295 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
388 G4cout<<"Target A and Z "<<theNucleus-> << 296
389 <<theNucleus->GetCharge()<<" (" << 297 // prepare the target residual
390 <<theNucleus->GetNumberOfLambdas()<<")"<<G4 << 298 G4int anA=theNucleus->GetMassNumber();
391 G4cout<<"Directly produced particles number << 299 G4int aZ=theNucleus->GetCharge();
392 G4cout<<"Projectile 4Mom and mass "<<GetPri << 300 G4int numberOfEx = 0;
393 <<GetPri << 301 G4int numberOfCh = 0;
394 #endif << 302 G4int numberOfHoles = 0;
395 << 303 G4double exEnergy = 0.0;
396 // prepare the target residual (assumed to << 304 G4double R = theNucleus->GetNuclearRadius();
397 G4int anA=theNucleus->GetMassNumber(); << 305 G4LorentzVector Target4Momentum(0,0,0,0);
398 G4int aZ=theNucleus->GetCharge(); << 306
399 //G4int aL=theNucleus->GetNumberOfLambdas() << 307 #ifdef debugPrecoInt
400 G4int numberOfEx = 0; << 308 G4cout<<"Target A Z "<<anA<<" "<<aZ<<G4endl;
401 G4int numberOfCh = 0; << 309 #endif
402 G4int numberOfHoles = 0; << 310
403 G4double exEnergy = 0.0; << 311 // loop over wounded target nucleus
404 G4double R = theNucleus->GetNuclearRadius() << 312 G4Nucleon * theCurrentNucleon =
405 G4LorentzVector Target4Momentum(0.,0.,0.,0. << 313 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
406 << 314 while(theCurrentNucleon) {
407 // loop over the wounded target nucleus << 315 if(theCurrentNucleon->AreYouHit()) {
408 G4Nucleon * theCurrentNucleon = << 316 ++numberOfHoles;
409 theNucleus->StartLoop() ? theNucleus- << 317 ++numberOfEx;
410 while(theCurrentNucleon) /* Loop checking << 318 --anA;
411 { << 319 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
412 if(theCurrentNucleon->AreYouHit()) { << 320 eplus + 0.1);
413 ++numberOfHoles; << 321 exEnergy += theCurrentNucleon->GetBindingEnergy();
414 ++numberOfEx; << 322 Target4Momentum -=theCurrentNucleon->Get4Momentum();
415 --anA; << 323 }
416 aZ -= G4int(theCurrentNucleon->GetDef << 324 theCurrentNucleon = theNucleus->GetNextNucleon();
417 eplus + 0.1); << 325 }
418 exEnergy += theCurrentNucleon->GetBin << 326
419 Target4Momentum -=theCurrentNucleon-> << 327 #ifdef debugPrecoInt
420 } << 328 G4cout<<"Residual Target A Z E* 4mom "<<anA<<" "<<aZ<<" "<<exEnergy<<" "
421 theCurrentNucleon = theNucleus->GetNextN << 329 <<Target4Momentum<<G4endl;
422 } << 330 #endif
423 << 331
424 #ifdef debugPrecoInt << 332 // prepare the projectile residual
425 G4cout<<"Residual Target A Z (numberOfLambd << 333 #ifdef debugPrecoInt
426 <<") "<<exEnergy<<" "<<Target4Momentu << 334 G4cout<<"Primary BaryonNumber "
427 #endif << 335 <<GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber()<<G4endl;
428 << 336 #endif
429 // prepare the projectile residual - which << 337
430 << 338 G4bool ProjectileIsAntiNucleus=
431 G4bool ProjectileIsAntiNucleus= << 339 GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber() < -1;
432 GetPrimaryProjectile()->GetDefinition << 340
433 << 341 G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();
434 G4ThreeVector bst = GetPrimaryProjectile()- << 342
435 << 343 G4int anAb=theProjectileNucleus->GetMassNumber();
436 G4int anAb=theProjectileNucleus->GetMassNum << 344 G4int aZb=theProjectileNucleus->GetCharge();
437 G4int aZb=theProjectileNucleus->GetCharge() << 345 G4int numberOfExB = 0;
438 G4int aLb=theProjectileNucleus->GetNumberOf << 346 G4int numberOfChB = 0;
439 G4int numberOfExB = 0; << 347 G4int numberOfHolesB = 0;
440 G4int numberOfChB = 0; << 348 G4double exEnergyB = 0.0;
441 G4int numberOfHolesB = 0; << 349 G4double Rb = theProjectileNucleus->GetNuclearRadius();
442 G4double exEnergyB = 0.0; << 350 G4LorentzVector Projectile4Momentum(0,0,0,0);
443 G4double Rb = theProjectileNucleus->GetNucl << 351
444 G4LorentzVector Projectile4Momentum(0.,0.,0 << 352 #ifdef debugPrecoInt
445 << 353 G4cout<<"Projectile A Z "<<anAb<<" "<<aZb<<G4endl;
446 // loop over the wounded projectile nucleus << 354 #endif
447 theCurrentNucleon = << 355
448 theProjectileNucleus->StartLoop() ? t << 356 // loop over wounded projectile nucleus
449 while(theCurrentNucleon) /* Loop checkin << 357 theCurrentNucleon =
>> 358 theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
>> 359 while(theCurrentNucleon) {
>> 360 if(theCurrentNucleon->AreYouHit()) {
>> 361 ++numberOfHolesB;
>> 362 ++numberOfExB;
>> 363 --anAb;
>> 364 if(!ProjectileIsAntiNucleus)
>> 365 {
>> 366 aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
>> 367 eplus + 0.1);
>> 368 } else
>> 369 {
>> 370 aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
>> 371 eplus - 0.1);
>> 372 }
>> 373 exEnergyB += theCurrentNucleon->GetBindingEnergy();
>> 374 Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
>> 375 }
>> 376 theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
>> 377 }
>> 378
>> 379 G4bool ExistTargetRemnant = G4double (numberOfHoles) <
>> 380 0.3* G4double (numberOfHoles + anA);
>> 381 G4bool ExistProjectileRemnant= G4double (numberOfHolesB) <
>> 382 0.3*G4double (numberOfHolesB + anAb);
>> 383
>> 384 #ifdef debugPrecoInt
>> 385 G4cout<<"Projectile residual A Z E* 4mom "<<anAb<<" "<<aZb<<" "<<exEnergyB<<" "
>> 386 <<Projectile4Momentum<<G4endl;
>> 387 G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
>> 388 <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
>> 389 #endif
>> 390 //-----------------------------------------------------------------------------
>> 391 // decay the strong resonances
>> 392 G4DecayKineticTracks decay(theSecondaries);
>> 393
>> 394 #ifdef debugPrecoInt
>> 395 G4LorentzVector secondary4Momemtum(0,0,0,0);
>> 396 G4int SecondrNum(0);
>> 397 #endif
>> 398
>> 399 // loop over secondaries
>> 400 G4KineticTrackVector::iterator iter;
>> 401 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
450 { 402 {
451 if(theCurrentNucleon->AreYouHit()) { << 403 G4ParticleDefinition* part = (*iter)->GetDefinition();
452 ++numberOfHolesB; << 404 G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();
453 ++numberOfExB; << 405
454 --anAb; << 406 if( part != proton && part != neutron &&
455 if(!ProjectileIsAntiNucleus) { << 407 (part != ANTIproton && ProjectileIsAntiNucleus) &&
456 aZb -= G4int(theCurrentNucleon->Ge << 408 (part != ANTIneutron && ProjectileIsAntiNucleus) )
457 eplus + 0.1); << 409 {
458 if (theCurrentNucleon->GetParticleType() << 410 G4ReactionProduct * theNew = new G4ReactionProduct(part);
459 } else { << 411 theNew->SetMomentum(aTrack4Momentum.vect());
460 aZb += G4int(theCurrentNucleon->Ge << 412 theNew->SetTotalEnergy(aTrack4Momentum.e());
461 eplus - 0.1); << 413 theTotalResult->push_back(theNew);
462 if (theCurrentNucleon->GetParticleType() << 414 #ifdef debugPrecoInt
463 } << 415 SecondrNum++;
464 exEnergyB += theCurrentNucleon->GetBi << 416 secondary4Momemtum += (*iter)->Get4Momentum();
465 Projectile4Momentum -=theCurrentNucle << 417 G4cout<<"Secondary "<<SecondrNum<<" "
466 } << 418 <<theNew->GetDefinition()->GetParticleName()<<" "
467 theCurrentNucleon = theProjectileNucleus << 419 <<secondary4Momemtum<<G4endl;
468 } << 420 #endif
469 << 421 delete (*iter);
470 G4bool ExistTargetRemnant = G4double (nu << 422 continue;
471 0.3* G4double (numberOfHoles + anA); << 423 }
472 G4bool ExistProjectileRemnant= G4double (nu << 424
473 0.3*G4double (numberOfHolesB + anAb); << 425 G4bool CanBeCapturedByTarget = false;
474 << 426 if( part == proton || part == neutron)
475 #ifdef debugPrecoInt << 427 {
476 G4cout<<"Projectile residual A Z (numberOfL << 428 CanBeCapturedByTarget = ExistTargetRemnant &&
477 <<") "<<exEnergyB<<" "<<Projectile4Momentum << 429 (CaptureThreshold >
478 G4cout<<" ExistTargetRemnant ExistProjectil << 430 (aTrack4Momentum + Target4Momentum).mag() -
479 <<ExistTargetRemnant<<" "<< ExistProj << 431 aTrack4Momentum.mag() - Target4Momentum.mag()) &&
480 #endif << 432 ((*iter)->GetPosition().mag() < R);
481 //----------------------------------------- << 433 }
482 // decay the strong resonances << 434 // ---------------------------
483 G4ReactionProductVector * theTotalResult = << 435 G4LorentzVector Position((*iter)->GetPosition(),
484 G4DecayKineticTracks decay(theSecondaries); << 436 (*iter)->GetFormationTime());
485 << 437 Position.boost(bst);
486 MakeCoalescence(theSecondaries); << 438
487 << 439 G4bool CanBeCapturedByProjectile = false;
488 #ifdef debugPrecoInt << 440
489 G4cout<<"Secondary stable particles numb << 441 if( !ProjectileIsAntiNucleus &&
490 #endif << 442 ( part == proton || part == neutron))
491 << 443 {
492 #ifdef debugPrecoInt << 444 CanBeCapturedByProjectile = ExistProjectileRemnant &&
493 G4LorentzVector secondary4Momemtum(0,0,0,0) << 445 (CaptureThreshold >
494 G4int SecondrNum(0); << 446 (aTrack4Momentum + Projectile4Momentum).mag() -
495 #endif << 447 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
496 << 448 (Position.vect().mag() < Rb);
497 // Loop over secondaries. << 449 }
498 // We are assuming that only protons and ne << 450
499 // and only antiprotons and antineutrons - << 451 if( ProjectileIsAntiNucleus &&
500 // not instead lambdas (or hyperons more ge << 452 ( part == ANTIproton || part == ANTIneutron))
501 // (or anti-hyperons more generally) - for << 453 {
502 // to be eventually reviewed later on, in p << 454 CanBeCapturedByProjectile = ExistProjectileRemnant &&
503 // anti-hypernuclei are introduced, instead << 455 (CaptureThreshold >
504 // anti-hypernuclei which currently exist i << 456 (aTrack4Momentum + Projectile4Momentum).mag() -
505 G4KineticTrackVector::iterator iter; << 457 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
506 for(iter=theSecondaries->begin(); iter !=th << 458 (Position.vect().mag() < Rb);
507 { << 459 }
508 const G4ParticleDefinition* part = (*ite << 460
509 G4LorentzVector aTrack4Momentum=(*iter)- << 461 if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
510 << 462 {
511 if( part != proton && part != neutron && << 463 if(G4UniformRand() < 0.5)
512 (part != ANTIproton && Projectile << 464 { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;}
513 (part != ANTIneutron && Projectile << 465 else
514 { << 466 { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
515 G4ReactionProduct * theNew = new G4Re << 467 }
516 theNew->SetMomentum(aTrack4Momentum.v << 468
517 theNew->SetTotalEnergy(aTrack4Momentu << 469 if(CanBeCapturedByTarget)
518 theNew->SetCreatorModelID((*iter)->GetCreat << 470 {
519 theNew->SetParentResonanceDef((*iter)->GetP << 471 // within the target nucleus, neutron or proton
520 theNew->SetParentResonanceID((*iter)->GetPa << 472 // now calculate A, Z of the fragment, momentum,
521 theTotalResult->push_back(theNew); << 473 // number of exciton states
522 #ifdef debugPrecoInt << 474 #ifdef debugPrecoInt
523 SecondrNum++; << 475 G4cout<<"Track is CapturedByTarget "<<" "
524 secondary4Momemtum += (*iter)->Get4Mo << 476 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
525 G4cout<<"Secondary "<<SecondrNum<<" << 477 #endif
526 <<theNew->GetDefinition()->GetP << 478 ++anA;
527 <<theNew->GetMomentum()<<" "<<t << 479 ++numberOfEx;
528 << 480 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
529 #endif << 481 aZ += Z;
530 delete (*iter); << 482 numberOfCh += Z;
531 continue; << 483 Target4Momentum +=aTrack4Momentum;
532 } << 484 delete (*iter);
533 << 485 } else if(CanBeCapturedByProjectile)
534 G4bool CanBeCapturedByTarget = false; << 486 {
535 if( part == proton || part == neutron) << 487 // within the projectile nucleus, neutron or proton
536 { << 488 // now calculate A, Z of the fragment, momentum,
537 CanBeCapturedByTarget = ExistTargetRe << 489 // number of exciton states
538 (-CaptureThreshold*G4Log( G4Uni << 490 #ifdef debugPrecoInt
539 (aTrack4Momentum + Target4 << 491 G4cout<<"Track is CapturedByProjectile"<<" "
540 aTrack4Momentum.mag() - Target4 << 492 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
541 ((*iter)->GetPosition().mag << 493 #endif
542 } << 494 ++anAb;
543 // --------------------------- << 495 ++numberOfExB;
544 G4LorentzVector Position((*iter)->GetPos << 496 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
545 Position.boost(bst); << 497 if( ProjectileIsAntiNucleus ) Z=-Z;
546 << 498 aZb += Z;
547 G4bool CanBeCapturedByProjectile = false << 499 numberOfChB += Z;
548 << 500 Projectile4Momentum +=aTrack4Momentum;
549 if( !ProjectileIsAntiNucleus && << 501 delete (*iter);
550 ( part == proton || part == neutro << 502 } else
551 { << 503 { // the track is not captured
552 CanBeCapturedByProjectile = ExistProj << 504 G4ReactionProduct * theNew = new G4ReactionProduct(part);
553 (-CaptureThreshold*G4Log( G4Uni << 505 theNew->SetMomentum(aTrack4Momentum.vect());
554 (aTrack4Momentum + Project << 506 theNew->SetTotalEnergy(aTrack4Momentum.e());
555 aTrack4Momentum.mag() - Project << 507 theTotalResult->push_back(theNew);
556 (Position.vect().mag() < Rb << 508
557 } << 509 #ifdef debugPrecoInt
558 << 510 SecondrNum++;
559 if( ProjectileIsAntiNucleus && << 511 secondary4Momemtum += (*iter)->Get4Momentum();
560 ( part == ANTIproton || part == AN << 512 G4cout<<"Secondary "<<SecondrNum<<" "
561 { << 513 <<theNew->GetDefinition()->GetParticleName()<<" "
562 CanBeCapturedByProjectile = ExistProj << 514 <<secondary4Momemtum<<G4endl;
563 (-CaptureThreshold*G4Log( G4Uni << 515 #endif
564 (aTrack4Momentum + Project << 516 delete (*iter);
565 aTrack4Momentum.mag() - Project << 517 continue;
566 (Position.vect().mag() < Rb << 518 }
567 } << 519 }
568 << 520 delete theSecondaries;
569 if(CanBeCapturedByTarget && CanBeCapture << 521 //-----------------------------------------------------
570 { << 522
571 if(G4UniformRand() < 0.5) << 523 #ifdef debugPrecoInt
572 { CanBeCapturedByTarget = true; CanBe << 524 G4cout<<"Final target residual A Z E* 4mom "<<anA<<" "<<aZ<<" "
573 else << 525 <<exEnergy<<" "<<Target4Momentum<<G4endl;
574 { CanBeCapturedByTarget = false; CanB << 526 #endif
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 527
747 if(ProjectileIsAntiNucleus) << 528 if(0!=anA )
>> 529 {
>> 530 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
>> 531
>> 532 if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
>> 533 {Target4Momentum.setE(fMass);}
>> 534
>> 535 G4double RemnMass=Target4Momentum.mag();
>> 536 if(RemnMass < fMass)
748 { 537 {
749 const G4ParticleDefinition * aFrag << 538 RemnMass=fMass + exEnergy;
750 const G4ParticleDefinition * LastF << 539 Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
751 if (aFragment == proton) {Las << 540 RemnMass*RemnMass));
752 else if(aFragment == neutron) {Las << 541 } else
753 else if(aFragment == lambda) {Las << 542 { exEnergy=RemnMass-fMass;}
754 else if(aFragment == deuteron){Las << 543
755 else if(aFragment == triton) {Las << 544 if( exEnergy < 0.) exEnergy=0.;
756 else if(aFragment == He3) {Las << 545
757 else if(aFragment == He4) {Las << 546 // Need to de-excite the remnant nucleus
758 else {} << 547 G4Fragment anInitialState(anA, aZ, Target4Momentum);
759 << 548 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
760 if (aLb > 0) { // Anti-hypernucle << 549 anInitialState.SetNumberOfCharged(numberOfCh);
761 if (aFragment == G4HyperTriton: << 550 anInitialState.SetNumberOfHoles(numberOfHoles);
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 551
787 theTotalResult->push_back(aPrecoResul << 552 G4ReactionProductVector * aPrecoResult =
788 } << 553 theDeExcitation->DeExcite(anInitialState);
789 554
790 delete aPrecoResult; << 555 // fill pre-compound part into the result, and return
791 } << 556 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
>> 557 {
>> 558 theTotalResult->push_back(aPrecoResult->operator[](ll));
>> 559 #ifdef debugPrecoInt
>> 560 G4cout<<"Tr frag "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
>> 561 <<" "<<aPrecoResult->operator[](ll)->GetMomentum()<<G4endl;
>> 562 #endif
>> 563 }
>> 564 delete aPrecoResult;
>> 565 }
>> 566
>> 567 //-----------------------------------------------------
>> 568 if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
>> 569 {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
792 570
793 return theTotalResult; << 571 #ifdef debugPrecoInt
794 } << 572 G4cout<<"Final projectile residual A Z E* Pmom "<<anAb<<" "<<aZb<<" "
>> 573 <<exEnergyB<<" "<<Projectile4Momentum<<G4endl;
>> 574 #endif
795 575
>> 576 if(0!=anAb)
>> 577 {
>> 578 G4double fMass = G4NucleiProperties::GetNuclearMass(anAb, aZb);
>> 579 G4double RemnMass=Projectile4Momentum.mag();
796 580
797 void G4GeneratorPrecompoundInterface::MakeCoal << 581 if(RemnMass < fMass)
798 // This method replaces pairs of proton-neut << 582 {
799 // antiproton-antineutron - in the case of a << 583 RemnMass=fMass + exEnergyB;
800 // momentum, with, respectively, deuterons a << 584 Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() +
801 // Note that in the case of hypernuclei or a << 585 RemnMass*RemnMass));
802 // are not considered for coalescence becaus << 586 } else
803 // are assumed not to exist. << 587 { exEnergyB=RemnMass-fMass;}
804 << 588
805 if (!tracks) return; << 589 if( exEnergyB < 0.) exEnergyB=0.;
806 << 590
807 G4double MassCut = deuteron->GetPDGMass() + << 591 // Need to de-excite the remnant nucleus
808 << 592 G4Fragment anInitialState(anAb, aZb, Projectile4Momentum);
809 for ( std::size_t i = 0; i < tracks->size(); << 593 anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
810 << 594 anInitialState.SetNumberOfCharged(numberOfChB);
811 G4KineticTrack* trackP = (*tracks)[i]; << 595 anInitialState.SetNumberOfHoles(numberOfHolesB);
812 if ( ! trackP ) continue; << 596
813 if (trackP->GetDefinition() != proton) con << 597 G4ReactionProductVector * aPrecoResult =
814 << 598 theDeExcitation->DeExcite(anInitialState);
815 G4LorentzVector Prot4Mom = trackP->Get4Mom << 599
816 G4LorentzVector ProtSPposition = G4Lorentz << 600 // fill pre-compound part into the result, and return
817 << 601 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
818 for ( std::size_t j = 0; j < tracks->size( << 602 {
819 << 603 if(ProjectileIsAntiNucleus)
820 G4KineticTrack* trackN = (*tracks)[j]; << 604 {
821 if (! trackN ) continue; << 605
822 if (trackN->GetDefinition() != neutron) << 606 #ifdef debugPrecoInt
823 << 607 G4cout<<"aPrecoRes "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
824 G4LorentzVector Neut4Mom = trackN->Get4M << 608 <<" "<<aPrecoResult->operator[](ll)->GetMomentum()
825 G4LorentzVector NeutSPposition = G4Loren << 609 <<" "<<aPrecoResult->operator[](ll)->GetTotalEnergy()
826 G4double EffMass = (Prot4Mom + Neut4Mom) << 610 <<" "<<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
827 << 611 #endif
828 if ( EffMass <= MassCut ) { // && (EffD << 612
829 G4KineticTrack* aDeuteron = << 613 G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
830 new G4KineticTrack( deuteron, << 614 G4ParticleDefinition * LastFragment=aFragment;
831 (trackP->GetForm << 615 if (aFragment == proton) {LastFragment=G4AntiProton::AntiProtonDefinition();}
832 (trackP->GetPosi << 616 else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
833 ( Prot4Mom << 617 else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
834 aDeuteron->SetCreatorModelID(secID); << 618 else if(aFragment == triton) {LastFragment=G4AntiTriton::AntiTritonDefinition();}
835 tracks->push_back(aDeuteron); << 619 else if(aFragment == He3) {LastFragment=G4AntiHe3::AntiHe3Definition();}
836 delete trackP; delete trackN; << 620 else if(aFragment == He4) {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
837 (*tracks)[i] = nullptr; (*tracks)[j] = << 621 else {}
838 break; << 622
839 } << 623 aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
840 } << 624 }
841 } <<
842 625
843 // Find and remove null pointers created by << 626 #ifdef debugPrecoInt
844 for ( G4int jj = (G4int)tracks->size()-1; jj << 627 G4cout<<"aPrecoResA "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
845 if ( ! (*tracks)[jj] ) tracks->erase(track << 628 <<" "<<aPrecoResult->operator[](ll)->GetMomentum()
846 } << 629 <<" "<<aPrecoResult->operator[](ll)->GetTotalEnergy()
>> 630 <<" "<<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
>> 631 #endif
>> 632 theTotalResult->push_back(aPrecoResult->operator[](ll));
>> 633 }
>> 634
>> 635 delete aPrecoResult;
>> 636 }
>> 637
>> 638 return theTotalResult;
847 } 639 }
>> 640
>> 641 // Uzhi Nov. 2012 ------------------------------------------------
>> 642
848 643