| 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 92692 2015-09-14 07:06:19Z 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 //-------------------------------------------- 73 //---------------------------------------------------------------------
95 #include "Randomize.hh" 74 #include "Randomize.hh"
96 #include "G4Log.hh" 75 #include "G4Log.hh"
97 76
98 //#define debugPrecoInt 77 //#define debugPrecoInt
99 78
100 G4GeneratorPrecompoundInterface::G4GeneratorPr 79 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* preModel)
101 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), << 80 : CaptureThreshold(70*MeV) // Uzhi 1.05.2015 10 ->70
102 { 81 {
103 proton = G4Proton::Proton(); 82 proton = G4Proton::Proton();
104 neutron = G4Neutron::Neutron(); 83 neutron = G4Neutron::Neutron();
105 lambda = G4Lambda::Lambda(); <<
106 84
107 deuteron=G4Deuteron::Deuteron(); 85 deuteron=G4Deuteron::Deuteron();
108 triton =G4Triton::Triton(); 86 triton =G4Triton::Triton();
109 He3 =G4He3::He3(); 87 He3 =G4He3::He3();
110 He4 =G4Alpha::Alpha(); 88 He4 =G4Alpha::Alpha();
111 89
112 ANTIproton=G4AntiProton::AntiProton(); 90 ANTIproton=G4AntiProton::AntiProton();
113 ANTIneutron=G4AntiNeutron::AntiNeutron(); 91 ANTIneutron=G4AntiNeutron::AntiNeutron();
114 92
115 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron() 93 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron();
116 ANTItriton =G4AntiTriton::AntiTriton(); 94 ANTItriton =G4AntiTriton::AntiTriton();
117 ANTIHe3 =G4AntiHe3::AntiHe3(); 95 ANTIHe3 =G4AntiHe3::AntiHe3();
118 ANTIHe4 =G4AntiAlpha::AntiAlpha(); 96 ANTIHe4 =G4AntiAlpha::AntiAlpha();
119 97
120 if(preModel) { SetDeExcitation(preModel); } 98 if(preModel) { SetDeExcitation(preModel); }
121 else { 99 else {
122 G4HadronicInteraction* hadi = 100 G4HadronicInteraction* hadi =
123 G4HadronicInteractionRegistry::Ins 101 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
124 G4VPreCompoundModel* pre = static_cast<G 102 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
125 if(!pre) { pre = new G4PreCompoundModel( 103 if(!pre) { pre = new G4PreCompoundModel(); }
126 SetDeExcitation(pre); 104 SetDeExcitation(pre);
127 } 105 }
128 <<
129 secID = G4PhysicsModelCatalog::GetModelID(" <<
130 } 106 }
131 107
132 G4GeneratorPrecompoundInterface::~G4GeneratorP 108 G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
133 { 109 {
134 } 110 }
135 111
136 G4ReactionProductVector* G4GeneratorPrecompoun 112 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
137 Propagate(G4KineticTrackVector* theSecondaries 113 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
138 { 114 {
139 #ifdef debugPrecoInt 115 #ifdef debugPrecoInt
140 G4cout<<G4endl<<"G4GeneratorPrecompoundI 116 G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::Propagate"<<G4endl;
141 G4cout<<"Target A and Z "<<theNucleus->G 117 G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
142 G4cout<<"Directly produced particles num 118 G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
143 #endif 119 #endif
144 120
145 G4ReactionProductVector * theTotalResult = 121 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
146 122
147 // decay the strong resonances 123 // decay the strong resonances
148 G4DecayKineticTracks decay(theSecondaries); 124 G4DecayKineticTracks decay(theSecondaries);
149 #ifdef debugPrecoInt 125 #ifdef debugPrecoInt
150 G4cout<<"Final stable particles number " 126 G4cout<<"Final stable particles number "<<theSecondaries->size()<<G4endl;
151 #endif 127 #endif
152 128
153 // prepare the fragment (it is assumed that << 129 // prepare the fragment
154 G4int anA=theNucleus->GetMassNumber(); 130 G4int anA=theNucleus->GetMassNumber();
155 G4int aZ=theNucleus->GetCharge(); 131 G4int aZ=theNucleus->GetCharge();
156 // G4double TargetNucleusMass = G4NucleiPrope 132 // G4double TargetNucleusMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
157 133
158 G4int numberOfEx = 0; 134 G4int numberOfEx = 0;
159 G4int numberOfCh = 0; 135 G4int numberOfCh = 0;
160 G4int numberOfHoles = 0; 136 G4int numberOfHoles = 0;
161 137
162 G4double R = theNucleus->GetNuclearRadius() 138 G4double R = theNucleus->GetNuclearRadius();
163 139
164 G4LorentzVector captured4Momentum(0.,0.,0., 140 G4LorentzVector captured4Momentum(0.,0.,0.,0.);
165 G4LorentzVector Residual4Momentum(0.,0.,0., 141 G4LorentzVector Residual4Momentum(0.,0.,0.,0.); // TargetNucleusMass is not need at the moment
166 G4LorentzVector Secondary4Momentum(0.,0.,0. 142 G4LorentzVector Secondary4Momentum(0.,0.,0.,0.);
167 143
168 // loop over secondaries 144 // loop over secondaries
169 G4KineticTrackVector::iterator iter; 145 G4KineticTrackVector::iterator iter;
170 for(iter=theSecondaries->begin(); iter !=th 146 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
171 { 147 {
172 const G4ParticleDefinition* part = (*ite 148 const G4ParticleDefinition* part = (*iter)->GetDefinition();
173 G4double e = (*iter)->Get4Momentum().e() 149 G4double e = (*iter)->Get4Momentum().e();
174 G4double mass = (*iter)->Get4Momentum(). 150 G4double mass = (*iter)->Get4Momentum().mag();
175 G4ThreeVector mom = (*iter)->Get4Momentu 151 G4ThreeVector mom = (*iter)->Get4Momentum().vect();
176 if((part != proton && part != neutron) | 152 if((part != proton && part != neutron) ||
>> 153 // Uzhi 2.05.2015 (e > mass + CaptureThreshold) ||
177 ((*iter)->GetPosition().mag() > R) 154 ((*iter)->GetPosition().mag() > R)) {
178 G4ReactionProduct * theNew = new G4Re 155 G4ReactionProduct * theNew = new G4ReactionProduct(part);
179 theNew->SetMomentum(mom); 156 theNew->SetMomentum(mom);
180 theNew->SetTotalEnergy(e); 157 theNew->SetTotalEnergy(e);
181 theNew->SetCreatorModelID((*iter)->GetCreat <<
182 theNew->SetParentResonanceDef((*iter)->GetP <<
183 theNew->SetParentResonanceID((*iter)->GetPa <<
184 theTotalResult->push_back(theNew); 158 theTotalResult->push_back(theNew);
185 Secondary4Momentum += (*iter)->Get4Mo << 159 Secondary4Momentum += (*iter)->Get4Momentum(); // Uzhi 29 April
186 #ifdef debugPrecoInt 160 #ifdef debugPrecoInt
187 G4cout<<"Secondary 4Mom "<<part->G 161 G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
188 <<(*iter)->Get4Momentum().ma 162 <<(*iter)->Get4Momentum().mag()<<G4endl;
189 #endif 163 #endif
190 } else { 164 } else {
191 if( e-mass > -CaptureThreshold*G4Log( << 165 if( e-mass > -CaptureThreshold*G4Log( G4UniformRand()) ) { // Added by Uzhi 2.05.2015
192 G4ReactionProduct * theNew = new G 166 G4ReactionProduct * theNew = new G4ReactionProduct(part);
193 theNew->SetMomentum(mom); 167 theNew->SetMomentum(mom);
194 theNew->SetTotalEnergy(e); 168 theNew->SetTotalEnergy(e);
195 theNew->SetCreatorModelID((*iter)->GetCr <<
196 theNew->SetParentResonanceDef((*iter)->G <<
197 theNew->SetParentResonanceID((*iter)->Ge <<
198 theTotalResult->push_back(theNew); 169 theTotalResult->push_back(theNew);
199 Secondary4Momentum += (*iter)->Get << 170 Secondary4Momentum += (*iter)->Get4Momentum(); // Uzhi 29 April
200 #ifdef debugPrecoInt 171 #ifdef debugPrecoInt
201 G4cout<<"Secondary 4Mom "<<part 172 G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
202 <<(*iter)->Get4Momentum() 173 <<(*iter)->Get4Momentum().mag()<<G4endl;
203 #endif 174 #endif
204 } else { 175 } else {
205 // within the nucleus, neutron or 176 // within the nucleus, neutron or proton
206 // now calculate A, Z of the frag 177 // now calculate A, Z of the fragment, momentum, number of exciton states
207 ++anA; 178 ++anA;
208 ++numberOfEx; 179 ++numberOfEx;
209 G4int Z = G4int(part->GetPDGCharge 180 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
210 aZ += Z; 181 aZ += Z;
211 numberOfCh += Z; 182 numberOfCh += Z;
212 captured4Momentum += (*iter)->Get4 183 captured4Momentum += (*iter)->Get4Momentum();
213 #ifdef debugPrecoInt 184 #ifdef debugPrecoInt
214 G4cout<<"Captured 4Mom "<<part 185 G4cout<<"Captured 4Mom "<<part->GetParticleName()<<(*iter)->Get4Momentum()<<G4endl;
215 #endif 186 #endif
216 } 187 }
217 } 188 }
218 delete (*iter); 189 delete (*iter);
219 } 190 }
220 delete theSecondaries; 191 delete theSecondaries;
221 192
222 // loop over wounded nucleus 193 // loop over wounded nucleus
223 G4Nucleon * theCurrentNucleon = 194 G4Nucleon * theCurrentNucleon =
224 theNucleus->StartLoop() ? theNucleus- 195 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
225 while(theCurrentNucleon) /* Loop checkin 196 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
226 { 197 {
227 if(theCurrentNucleon->AreYouHit()) { 198 if(theCurrentNucleon->AreYouHit()) {
228 ++numberOfHoles; 199 ++numberOfHoles;
229 ++numberOfEx; 200 ++numberOfEx;
230 --anA; 201 --anA;
231 aZ -= G4int(theCurrentNucleon->GetDef 202 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
232 203
233 Residual4Momentum -= theCurrentNucleo 204 Residual4Momentum -= theCurrentNucleon->Get4Momentum();
234 } 205 }
235 theCurrentNucleon = theNucleus->GetNextN 206 theCurrentNucleon = theNucleus->GetNextNucleon();
236 } 207 }
237 208
238 #ifdef debugPrecoInt 209 #ifdef debugPrecoInt
239 G4cout<<G4endl; 210 G4cout<<G4endl;
240 G4cout<<"Secondary 4Mom "<<Secondary4Mom 211 G4cout<<"Secondary 4Mom "<<Secondary4Momentum<<G4endl;
241 G4cout<<"Captured 4Mom "<<captured4Mome 212 G4cout<<"Captured 4Mom "<<captured4Momentum<<G4endl;
242 G4cout<<"Sec + Captured "<<Secondary4Mom 213 G4cout<<"Sec + Captured "<<Secondary4Momentum+captured4Momentum<<G4endl;
243 G4cout<<"Residual4Mom "<<Residual4Mome 214 G4cout<<"Residual4Mom "<<Residual4Momentum<<G4endl;
244 G4cout<<"Sum 4 momenta " 215 G4cout<<"Sum 4 momenta "
245 <<Secondary4Momentum + captured4Mo 216 <<Secondary4Momentum + captured4Momentum + Residual4Momentum <<G4endl;
246 #endif 217 #endif
247 218
248 // Check that we use QGS model; loop over w 219 // Check that we use QGS model; loop over wounded nucleus
249 G4bool QGSM(false); 220 G4bool QGSM(false);
250 theCurrentNucleon = theNucleus->StartLoop() 221 theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
251 while(theCurrentNucleon) /* Loop checking 222 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
252 { 223 {
253 if(theCurrentNucleon->AreYouHit()) 224 if(theCurrentNucleon->AreYouHit())
254 { 225 {
255 if(theCurrentNucleon->Get4Momentum().ma 226 if(theCurrentNucleon->Get4Momentum().mag() <
256 theCurrentNucleon->GetDefinition()-> 227 theCurrentNucleon->GetDefinition()->GetPDGMass()) QGSM=true;
257 } 228 }
258 theCurrentNucleon = theNucleus->GetNextN 229 theCurrentNucleon = theNucleus->GetNextNucleon();
259 } 230 }
260 231
261 #ifdef debugPrecoInt 232 #ifdef debugPrecoInt
262 if(!QGSM){ 233 if(!QGSM){
263 G4cout<<G4endl; 234 G4cout<<G4endl;
264 G4cout<<"Residual A and Z "<<anA<<" "< 235 G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
265 G4cout<<"Residual 4Mom "<<Residual4Mo 236 G4cout<<"Residual 4Mom "<<Residual4Momentum<<G4endl;
266 if(numberOfEx == 0) 237 if(numberOfEx == 0)
267 {G4cout<<"Residual 4Mom = 0 means tha 238 {G4cout<<"Residual 4Mom = 0 means that there were not wounded and captured nucleons"<<G4endl;}
268 } 239 }
269 #endif 240 #endif
270 241
271 if(anA == 0) return theTotalResult; 242 if(anA == 0) return theTotalResult;
272 243
273 G4LorentzVector exciton4Momentum(0.,0.,0.,0 << 244 G4LorentzVector exciton4Momentum(0.,0.,0.,0.); // Uzhi 29 April
274 if(anA >= aZ) 245 if(anA >= aZ)
275 { 246 {
276 if(!QGSM) 247 if(!QGSM)
277 { // FTF model was used 248 { // FTF model was used
278 G4double fMass = G4NucleiProperties::Ge 249 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
279 250
280 // G4LorentzVector exciton4Momentum = Res 251 // G4LorentzVector exciton4Momentum = Residual4Momentum + captured4Momentum;
281 exciton4Momentum = Residual4Momentum + c 252 exciton4Momentum = Residual4Momentum + captured4Momentum;
282 //exciton4Momentum.setE(std::sqrt(exciton4Mome 253 //exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
283 G4double ActualMass = exciton4Momentum.m 254 G4double ActualMass = exciton4Momentum.mag();
284 if(ActualMass <= fMass ) { << 255 if(ActualMass <= fMass ) { //E*<=0, Uzhi 5.05.2015
285 exciton4Momentum.setE(std::sqrt(excito << 256 exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass))); // Uzhi 13.05.2015
286 } 257 }
287 258
288 #ifdef debugPrecoInt 259 #ifdef debugPrecoInt
289 G4double exEnergy = 0.0; 260 G4double exEnergy = 0.0;
290 if(ActualMass <= fMass ) {exEnergy = << 261 if(ActualMass <= fMass ) {exEnergy = 0.;} // Uzhi 5.05.2015
291 else {exEnergy = 262 else {exEnergy = ActualMass - fMass;}
292 G4cout<<"Ground state residual Mass " 263 G4cout<<"Ground state residual Mass "<<fMass<<" E* "<<exEnergy<<G4endl;
293 #endif 264 #endif
294 } 265 }
295 else 266 else
296 { // QGS model was used 267 { // QGS model was used
297 G4double InitialTargetMass = 268 G4double InitialTargetMass =
298 G4NucleiProperties::GetNuclearMa 269 G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(), theNucleus->GetCharge());
299 270
300 exciton4Momentum = 271 exciton4Momentum =
301 GetPrimaryProjectile()->Get4Mome 272 GetPrimaryProjectile()->Get4Momentum() + G4LorentzVector(0.,0.,0.,InitialTargetMass)
302 -Secondary4Momentum; 273 -Secondary4Momentum;
303 274
304 G4double fMass = G4NucleiProperties::Get 275 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
305 G4double ActualMass = exciton4Momentum.ma 276 G4double ActualMass = exciton4Momentum.mag();
306 277
307 #ifdef debugPrecoInt 278 #ifdef debugPrecoInt
308 G4cout<<G4endl; 279 G4cout<<G4endl;
309 G4cout<<"Residual A and Z "<<anA<<" "< 280 G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
310 G4cout<<"Residual4Momentum "<<exciton4 281 G4cout<<"Residual4Momentum "<<exciton4Momentum<<G4endl;
311 G4cout<<"ResidualMass, GroundStateMass 282 G4cout<<"ResidualMass, GroundStateMass and E* "<<ActualMass<<" "<<fMass<<" "
312 <<ActualMass - fMass<<G4endl; 283 <<ActualMass - fMass<<G4endl;
313 #endif 284 #endif
314 285
315 if(ActualMass - fMass < 0.) 286 if(ActualMass - fMass < 0.)
316 { 287 {
317 G4double ResE = std::sqrt(exciton4Moment 288 G4double ResE = std::sqrt(exciton4Momentum.vect().mag2() + sqr(fMass+10*MeV));
318 exciton4Momentum.setE(ResE); 289 exciton4Momentum.setE(ResE);
319 #ifdef debugPrecoInt 290 #ifdef debugPrecoInt
320 G4cout<<"ActualMass - fMass < 0. "<<A 291 G4cout<<"ActualMass - fMass < 0. "<<ActualMass<<" "<<fMass<<" "<<ActualMass - fMass<<G4endl;
>> 292 G4int Uzhi; G4cin>>Uzhi;
321 #endif 293 #endif
322 } 294 }
323 } 295 }
>> 296 // Need to de-excite the remnant nucleus only if excitation energy > 0.
>> 297 G4Fragment anInitialState(anA, aZ, exciton4Momentum);
>> 298 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
>> 299 anInitialState.SetNumberOfCharged(numberOfCh);
>> 300 anInitialState.SetNumberOfHoles(numberOfHoles);
324 301
325 // Need to de-excite the remnant nucleus o << 302 G4ReactionProductVector * aPrecoResult =
326 G4Fragment anInitialState(anA, aZ, exciton << 303 theDeExcitation->DeExcite(anInitialState);
327 anInitialState.SetNumberOfParticles(number << 304 // fill pre-compound part into the result, and return
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 { <<
340 theTotalResult->push_back(aPrecoResult-> <<
341 #ifdef debugPrecoInt 305 #ifdef debugPrecoInt
342 G4cout<<"Fragment "<<ll<<" " << 306 G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
343 <<aPrecoResult->operator[](ll)->GetDefin << 307 #endif
344 <<aPrecoResult->operator[](ll)->GetMomen << 308 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
345 <<aPrecoResult->operator[](ll)->GetTotal << 309 {
346 <<aPrecoResult->operator[](ll)->GetDefin << 310 theTotalResult->push_back(aPrecoResult->operator[](ll));
347 #endif << 311 #ifdef debugPrecoInt
348 } << 312 G4cout<<"Fragment "<<ll<<" "
349 delete aPrecoResult; << 313 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
>> 314 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
>> 315 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
>> 316 <<aPrecoResult->operator[](ll)->GetDefinition()->GetPDGMass()<<G4endl;
>> 317 #endif
>> 318 }
>> 319 delete aPrecoResult;
350 } 320 }
351 321
352 return theTotalResult; 322 return theTotalResult;
353 } 323 }
354 324
355 G4HadFinalState* G4GeneratorPrecompoundInterfa 325 G4HadFinalState* G4GeneratorPrecompoundInterface::
356 ApplyYourself(const G4HadProjectile &, G4Nucle 326 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
357 { 327 {
358 G4cout << "G4GeneratorPrecompoundInterface: 328 G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
359 << G4endl; 329 << G4endl;
360 G4cout << "This class is only a mediator be 330 G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
361 G4cout << "Please remove from your physics 331 G4cout << "Please remove from your physics list."<<G4endl;
362 throw G4HadronicException(__FILE__, __LINE_ 332 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
363 return new G4HadFinalState; 333 return new G4HadFinalState;
364 } 334 }
365 void G4GeneratorPrecompoundInterface::Propagat 335 void G4GeneratorPrecompoundInterface::PropagateModelDescription(std::ostream& outFile) const
366 { 336 {
367 outFile << "G4GeneratorPrecompoundInterface 337 outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
368 << "energy model through the wounded << 338 << "energy model through the wounded nucleus to precompound de-excition.\n"
369 << "Low energy protons and neutron pr 339 << "Low energy protons and neutron present among secondaries produced by \n"
370 << "the high energy generator and wit 340 << "the high energy generator and within the nucleus are captured. The wounded\n"
371 << "nucleus and the captured particle 341 << "nucleus and the captured particles form an excited nuclear fragment. This\n"
372 << "fragment is passed to the Geant4 342 << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
373 << "Nuclear de-excitation:\n"; 343 << "Nuclear de-excitation:\n";
374 // preco 344 // preco
375 345
376 } 346 }
377 347
378 << 348 // Uzhi Nov. 2012 ------------------------------------------------
379 G4ReactionProductVector* G4GeneratorPrecompoun 349 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
380 PropagateNuclNucl(G4KineticTrackVector* theSec 350 PropagateNuclNucl(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus,
381 G4V3DNucleus* theProjectileNucleus) 351 G4V3DNucleus* theProjectileNucleus)
382 { 352 {
383 #ifdef debugPrecoInt 353 #ifdef debugPrecoInt
384 G4cout<<G4endl<<"G4GeneratorPrecompoundInte 354 G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
385 G4cout<<"Projectile A and Z (and numberOfLa << 355 G4cout<<"Projectile A and Z "<<theProjectileNucleus->GetMassNumber()<<" "
386 <<theProjectileNucleus->GetCharge()<<" (" << 356 <<theProjectileNucleus->GetCharge()<<G4endl;
387 <<theProjectileNucleus->GetNumberOfLambdas( <<
388 G4cout<<"Target A and Z "<<theNucleus-> 357 G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "
389 <<theNucleus->GetCharge()<<" (" << 358 <<theNucleus->GetCharge()<<G4endl;
390 <<theNucleus->GetNumberOfLambdas()<<")"<<G4 <<
391 G4cout<<"Directly produced particles number 359 G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
392 G4cout<<"Projectile 4Mom and mass "<<GetPri 360 G4cout<<"Projectile 4Mom and mass "<<GetPrimaryProjectile()->Get4Momentum()<<" "
393 <<GetPri 361 <<GetPrimaryProjectile()->Get4Momentum().mag()<<G4endl<<G4endl;
394 #endif 362 #endif
395 363
396 // prepare the target residual (assumed to << 364 // prepare the target residual
397 G4int anA=theNucleus->GetMassNumber(); 365 G4int anA=theNucleus->GetMassNumber();
398 G4int aZ=theNucleus->GetCharge(); 366 G4int aZ=theNucleus->GetCharge();
399 //G4int aL=theNucleus->GetNumberOfLambdas() <<
400 G4int numberOfEx = 0; 367 G4int numberOfEx = 0;
401 G4int numberOfCh = 0; 368 G4int numberOfCh = 0;
402 G4int numberOfHoles = 0; 369 G4int numberOfHoles = 0;
403 G4double exEnergy = 0.0; 370 G4double exEnergy = 0.0;
404 G4double R = theNucleus->GetNuclearRadius() 371 G4double R = theNucleus->GetNuclearRadius();
405 G4LorentzVector Target4Momentum(0.,0.,0.,0. 372 G4LorentzVector Target4Momentum(0.,0.,0.,0.);
406 373
407 // loop over the wounded target nucleus << 374 // loop over wounded target nucleus
408 G4Nucleon * theCurrentNucleon = 375 G4Nucleon * theCurrentNucleon =
409 theNucleus->StartLoop() ? theNucleus- 376 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
410 while(theCurrentNucleon) /* Loop checking 377 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
411 { 378 {
412 if(theCurrentNucleon->AreYouHit()) { 379 if(theCurrentNucleon->AreYouHit()) {
413 ++numberOfHoles; 380 ++numberOfHoles;
414 ++numberOfEx; 381 ++numberOfEx;
415 --anA; 382 --anA;
416 aZ -= G4int(theCurrentNucleon->GetDef 383 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
417 eplus + 0.1); 384 eplus + 0.1);
418 exEnergy += theCurrentNucleon->GetBin 385 exEnergy += theCurrentNucleon->GetBindingEnergy();
419 Target4Momentum -=theCurrentNucleon-> 386 Target4Momentum -=theCurrentNucleon->Get4Momentum();
420 } 387 }
421 theCurrentNucleon = theNucleus->GetNextN 388 theCurrentNucleon = theNucleus->GetNextNucleon();
422 } 389 }
423 390
424 #ifdef debugPrecoInt 391 #ifdef debugPrecoInt
425 G4cout<<"Residual Target A Z (numberOfLambd << 392 G4cout<<"Residual Target A Z E* 4mom "<<anA<<" "<<aZ<<" "<<exEnergy<<" "
426 <<") "<<exEnergy<<" "<<Target4Momentu << 393 <<Target4Momentum<<G4endl;
427 #endif 394 #endif
428 395
429 // prepare the projectile residual - which << 396 // prepare the projectile residual
430 397
431 G4bool ProjectileIsAntiNucleus= 398 G4bool ProjectileIsAntiNucleus=
432 GetPrimaryProjectile()->GetDefinition 399 GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber() < -1;
433 400
434 G4ThreeVector bst = GetPrimaryProjectile()- 401 G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();
435 402
436 G4int anAb=theProjectileNucleus->GetMassNum 403 G4int anAb=theProjectileNucleus->GetMassNumber();
437 G4int aZb=theProjectileNucleus->GetCharge() 404 G4int aZb=theProjectileNucleus->GetCharge();
438 G4int aLb=theProjectileNucleus->GetNumberOf <<
439 G4int numberOfExB = 0; 405 G4int numberOfExB = 0;
440 G4int numberOfChB = 0; 406 G4int numberOfChB = 0;
441 G4int numberOfHolesB = 0; 407 G4int numberOfHolesB = 0;
442 G4double exEnergyB = 0.0; 408 G4double exEnergyB = 0.0;
443 G4double Rb = theProjectileNucleus->GetNucl 409 G4double Rb = theProjectileNucleus->GetNuclearRadius();
444 G4LorentzVector Projectile4Momentum(0.,0.,0 410 G4LorentzVector Projectile4Momentum(0.,0.,0.,0.);
445 411
446 // loop over the wounded projectile nucleus << 412 // loop over wounded projectile nucleus
447 theCurrentNucleon = 413 theCurrentNucleon =
448 theProjectileNucleus->StartLoop() ? t 414 theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
449 while(theCurrentNucleon) /* Loop checkin 415 while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
450 { 416 {
451 if(theCurrentNucleon->AreYouHit()) { 417 if(theCurrentNucleon->AreYouHit()) {
452 ++numberOfHolesB; 418 ++numberOfHolesB;
453 ++numberOfExB; 419 ++numberOfExB;
454 --anAb; 420 --anAb;
455 if(!ProjectileIsAntiNucleus) { 421 if(!ProjectileIsAntiNucleus) {
456 aZb -= G4int(theCurrentNucleon->Ge 422 aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
457 eplus + 0.1); 423 eplus + 0.1);
458 if (theCurrentNucleon->GetParticleType() <<
459 } else { 424 } else {
460 aZb += G4int(theCurrentNucleon->Ge 425 aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
461 eplus - 0.1); 426 eplus - 0.1);
462 if (theCurrentNucleon->GetParticleType() <<
463 } 427 }
464 exEnergyB += theCurrentNucleon->GetBi 428 exEnergyB += theCurrentNucleon->GetBindingEnergy();
465 Projectile4Momentum -=theCurrentNucle 429 Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
466 } 430 }
467 theCurrentNucleon = theProjectileNucleus 431 theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
468 } 432 }
469 433
470 G4bool ExistTargetRemnant = G4double (nu 434 G4bool ExistTargetRemnant = G4double (numberOfHoles) <
471 0.3* G4double (numberOfHoles + anA); 435 0.3* G4double (numberOfHoles + anA);
472 G4bool ExistProjectileRemnant= G4double (nu 436 G4bool ExistProjectileRemnant= G4double (numberOfHolesB) <
473 0.3*G4double (numberOfHolesB + anAb); 437 0.3*G4double (numberOfHolesB + anAb);
474 438
475 #ifdef debugPrecoInt 439 #ifdef debugPrecoInt
476 G4cout<<"Projectile residual A Z (numberOfL << 440 G4cout<<"Projectile residual A Z E* 4mom "<<anAb<<" "<<aZb<<" "<<exEnergyB<<" "
477 <<") "<<exEnergyB<<" "<<Projectile4Momentum << 441 <<Projectile4Momentum<<G4endl;
478 G4cout<<" ExistTargetRemnant ExistProjectil 442 G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
479 <<ExistTargetRemnant<<" "<< ExistProj 443 <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
480 #endif 444 #endif
481 //----------------------------------------- 445 //-----------------------------------------------------------------------------
482 // decay the strong resonances 446 // decay the strong resonances
483 G4ReactionProductVector * theTotalResult = 447 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
484 G4DecayKineticTracks decay(theSecondaries); 448 G4DecayKineticTracks decay(theSecondaries);
485 <<
486 MakeCoalescence(theSecondaries); <<
487 <<
488 #ifdef debugPrecoInt 449 #ifdef debugPrecoInt
489 G4cout<<"Secondary stable particles numb 450 G4cout<<"Secondary stable particles number "<<theSecondaries->size()<<G4endl;
490 #endif 451 #endif
491 452
492 #ifdef debugPrecoInt 453 #ifdef debugPrecoInt
493 G4LorentzVector secondary4Momemtum(0,0,0,0) 454 G4LorentzVector secondary4Momemtum(0,0,0,0);
494 G4int SecondrNum(0); 455 G4int SecondrNum(0);
495 #endif 456 #endif
496 457
497 // Loop over secondaries. << 458 // 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; 459 G4KineticTrackVector::iterator iter;
506 for(iter=theSecondaries->begin(); iter !=th 460 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
507 { 461 {
508 const G4ParticleDefinition* part = (*ite 462 const G4ParticleDefinition* part = (*iter)->GetDefinition();
509 G4LorentzVector aTrack4Momentum=(*iter)- 463 G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();
510 464
511 if( part != proton && part != neutron && 465 if( part != proton && part != neutron &&
512 (part != ANTIproton && Projectile 466 (part != ANTIproton && ProjectileIsAntiNucleus) &&
513 (part != ANTIneutron && Projectile 467 (part != ANTIneutron && ProjectileIsAntiNucleus) )
514 { 468 {
515 G4ReactionProduct * theNew = new G4Re 469 G4ReactionProduct * theNew = new G4ReactionProduct(part);
516 theNew->SetMomentum(aTrack4Momentum.v 470 theNew->SetMomentum(aTrack4Momentum.vect());
517 theNew->SetTotalEnergy(aTrack4Momentu 471 theNew->SetTotalEnergy(aTrack4Momentum.e());
518 theNew->SetCreatorModelID((*iter)->GetCreat <<
519 theNew->SetParentResonanceDef((*iter)->GetP <<
520 theNew->SetParentResonanceID((*iter)->GetPa <<
521 theTotalResult->push_back(theNew); 472 theTotalResult->push_back(theNew);
522 #ifdef debugPrecoInt 473 #ifdef debugPrecoInt
523 SecondrNum++; 474 SecondrNum++;
524 secondary4Momemtum += (*iter)->Get4Mo 475 secondary4Momemtum += (*iter)->Get4Momentum();
525 G4cout<<"Secondary "<<SecondrNum<<" 476 G4cout<<"Secondary "<<SecondrNum<<" "
526 <<theNew->GetDefinition()->GetP 477 <<theNew->GetDefinition()->GetParticleName()<<" "
527 <<theNew->GetMomentum()<<" "<<t 478 <<theNew->GetMomentum()<<" "<<theNew->GetTotalEnergy()<<G4endl;
528 479
529 #endif 480 #endif
530 delete (*iter); 481 delete (*iter);
531 continue; 482 continue;
532 } 483 }
533 484
534 G4bool CanBeCapturedByTarget = false; 485 G4bool CanBeCapturedByTarget = false;
535 if( part == proton || part == neutron) 486 if( part == proton || part == neutron)
536 { 487 {
537 CanBeCapturedByTarget = ExistTargetRe 488 CanBeCapturedByTarget = ExistTargetRemnant &&
538 (-CaptureThreshold*G4Log( G4Uni 489 (-CaptureThreshold*G4Log( G4UniformRand()) >
539 (aTrack4Momentum + Target4 490 (aTrack4Momentum + Target4Momentum).mag() -
540 aTrack4Momentum.mag() - Target4 491 aTrack4Momentum.mag() - Target4Momentum.mag()) &&
541 ((*iter)->GetPosition().mag 492 ((*iter)->GetPosition().mag() < R);
542 } 493 }
543 // --------------------------- 494 // ---------------------------
544 G4LorentzVector Position((*iter)->GetPos 495 G4LorentzVector Position((*iter)->GetPosition(), (*iter)->GetFormationTime());
545 Position.boost(bst); 496 Position.boost(bst);
546 497
547 G4bool CanBeCapturedByProjectile = false 498 G4bool CanBeCapturedByProjectile = false;
548 499
549 if( !ProjectileIsAntiNucleus && 500 if( !ProjectileIsAntiNucleus &&
550 ( part == proton || part == neutro 501 ( part == proton || part == neutron))
551 { 502 {
552 CanBeCapturedByProjectile = ExistProj 503 CanBeCapturedByProjectile = ExistProjectileRemnant &&
553 (-CaptureThreshold*G4Log( G4Uni 504 (-CaptureThreshold*G4Log( G4UniformRand()) >
554 (aTrack4Momentum + Project 505 (aTrack4Momentum + Projectile4Momentum).mag() -
555 aTrack4Momentum.mag() - Project 506 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
556 (Position.vect().mag() < Rb 507 (Position.vect().mag() < Rb);
557 } 508 }
558 509
559 if( ProjectileIsAntiNucleus && 510 if( ProjectileIsAntiNucleus &&
560 ( part == ANTIproton || part == AN 511 ( part == ANTIproton || part == ANTIneutron))
561 { 512 {
562 CanBeCapturedByProjectile = ExistProj 513 CanBeCapturedByProjectile = ExistProjectileRemnant &&
563 (-CaptureThreshold*G4Log( G4Uni 514 (-CaptureThreshold*G4Log( G4UniformRand()) >
564 (aTrack4Momentum + Project 515 (aTrack4Momentum + Projectile4Momentum).mag() -
565 aTrack4Momentum.mag() - Project 516 aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
566 (Position.vect().mag() < Rb 517 (Position.vect().mag() < Rb);
567 } 518 }
568 519
569 if(CanBeCapturedByTarget && CanBeCapture 520 if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
570 { 521 {
571 if(G4UniformRand() < 0.5) 522 if(G4UniformRand() < 0.5)
572 { CanBeCapturedByTarget = true; CanBe 523 { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;}
573 else 524 else
574 { CanBeCapturedByTarget = false; CanB 525 { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
575 } 526 }
576 527
577 if(CanBeCapturedByTarget) 528 if(CanBeCapturedByTarget)
578 { 529 {
579 // within the target nucleus, neutron 530 // within the target nucleus, neutron or proton
580 // now calculate A, Z of the fragmen 531 // now calculate A, Z of the fragment, momentum,
581 // number of exciton states 532 // number of exciton states
582 #ifdef debugPrecoInt 533 #ifdef debugPrecoInt
583 G4cout<<"Track is CapturedByTarget "< 534 G4cout<<"Track is CapturedByTarget "<<" "<<part->GetParticleName()<<" "
584 <<aTrack4Momentum<<" "<<aTrack4 535 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
585 #endif 536 #endif
586 ++anA; 537 ++anA;
587 ++numberOfEx; 538 ++numberOfEx;
588 G4int Z = G4int(part->GetPDGCharge()/ 539 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
589 aZ += Z; 540 aZ += Z;
590 numberOfCh += Z; 541 numberOfCh += Z;
591 Target4Momentum +=aTrack4Momentum; 542 Target4Momentum +=aTrack4Momentum;
592 delete (*iter); 543 delete (*iter);
593 } else if(CanBeCapturedByProjectile) 544 } else if(CanBeCapturedByProjectile)
594 { 545 {
595 // within the projectile nucleus, neu 546 // within the projectile nucleus, neutron or proton
596 // now calculate A, Z of the fragmen 547 // now calculate A, Z of the fragment, momentum,
597 // number of exciton states 548 // number of exciton states
598 #ifdef debugPrecoInt 549 #ifdef debugPrecoInt
599 G4cout<<"Track is CapturedByProjectil 550 G4cout<<"Track is CapturedByProjectile"<<" "<<part->GetParticleName()<<" "
600 <<aTrack4Momentum<<" "<<aTrack4 551 <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
601 #endif 552 #endif
602 ++anAb; 553 ++anAb;
603 ++numberOfExB; 554 ++numberOfExB;
604 G4int Z = G4int(part->GetPDGCharge()/ 555 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
605 if( ProjectileIsAntiNucleus ) Z=-Z; 556 if( ProjectileIsAntiNucleus ) Z=-Z;
606 aZb += Z; 557 aZb += Z;
607 numberOfChB += Z; 558 numberOfChB += Z;
608 Projectile4Momentum +=aTrack4Momentum 559 Projectile4Momentum +=aTrack4Momentum;
609 delete (*iter); 560 delete (*iter);
610 } else 561 } else
611 { // the track is not captured 562 { // the track is not captured
612 G4ReactionProduct * theNew = new G4Re 563 G4ReactionProduct * theNew = new G4ReactionProduct(part);
613 theNew->SetMomentum(aTrack4Momentum.v 564 theNew->SetMomentum(aTrack4Momentum.vect());
614 theNew->SetTotalEnergy(aTrack4Momentu 565 theNew->SetTotalEnergy(aTrack4Momentum.e());
615 theNew->SetCreatorModelID((*iter)->GetCreat <<
616 theNew->SetParentResonanceDef((*iter)->GetP <<
617 theNew->SetParentResonanceID((*iter)->GetPa <<
618 theTotalResult->push_back(theNew); 566 theTotalResult->push_back(theNew);
619 567
620 #ifdef debugPrecoInt 568 #ifdef debugPrecoInt
621 SecondrNum++; 569 SecondrNum++;
622 secondary4Momemtum += (*iter)->Get4Mo 570 secondary4Momemtum += (*iter)->Get4Momentum();
623 /* 571 /*
624 G4cout<<"Secondary "<<SecondrNum<<" 572 G4cout<<"Secondary "<<SecondrNum<<" "
625 <<theNew->GetDefinition()->GetP 573 <<theNew->GetDefinition()->GetParticleName()<<" "
626 <<secondary4Momemtum<<G4endl; 574 <<secondary4Momemtum<<G4endl;
627 */ 575 */
628 #endif 576 #endif
629 delete (*iter); 577 delete (*iter);
630 continue; 578 continue;
631 } 579 }
632 } 580 }
633 delete theSecondaries; 581 delete theSecondaries;
634 //----------------------------------------- 582 //-----------------------------------------------------
635 583
636 #ifdef debugPrecoInt 584 #ifdef debugPrecoInt
637 G4cout<<"Final target residual A Z (numberO << 585 G4cout<<"Final target residual A Z E* 4mom "<<anA<<" "<<aZ<<" "
638 <<") "<<exEnergy<<" "<<Target4Momentum<< << 586 <<exEnergy<<" "<<Target4Momentum<<G4endl;
639 #endif 587 #endif
640 588
641 if(0!=anA ) 589 if(0!=anA )
642 { 590 {
643 // We assume that the target residual is <<
644 G4double fMass = G4NucleiProperties::Ge 591 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
645 592
646 if((anA == theNucleus->GetMassNumber()) 593 if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
647 {Target4Momentum.setE(fMass);} 594 {Target4Momentum.setE(fMass);}
648 595
649 G4double RemnMass=Target4Momentum.mag(); 596 G4double RemnMass=Target4Momentum.mag();
650 597
651 if(RemnMass < fMass) 598 if(RemnMass < fMass)
652 { 599 {
653 RemnMass=fMass + exEnergy; 600 RemnMass=fMass + exEnergy;
654 Target4Momentum.setE(std::sqrt(Target 601 Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
655 RemnMass*RemnMass)); 602 RemnMass*RemnMass));
656 } else 603 } else
657 { exEnergy=RemnMass-fMass;} 604 { exEnergy=RemnMass-fMass;}
658 605
659 if( exEnergy < 0.) exEnergy=0.; 606 if( exEnergy < 0.) exEnergy=0.;
660 607
661 // Need to de-excite the remnant nucleus 608 // Need to de-excite the remnant nucleus
662 G4Fragment anInitialState(anA, aZ, Targe 609 G4Fragment anInitialState(anA, aZ, Target4Momentum);
663 anInitialState.SetNumberOfParticles(numb 610 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
664 anInitialState.SetNumberOfCharged(number 611 anInitialState.SetNumberOfCharged(numberOfCh);
665 anInitialState.SetNumberOfHoles(numberOf 612 anInitialState.SetNumberOfHoles(numberOfHoles);
666 anInitialState.SetCreatorModelID(secID); <<
667 613
668 G4ReactionProductVector * aPrecoResult = 614 G4ReactionProductVector * aPrecoResult =
669 theDeExcitation->DeExcite(anInitia 615 theDeExcitation->DeExcite(anInitialState);
670 616
671 #ifdef debugPrecoInt 617 #ifdef debugPrecoInt
672 G4cout<<"Target fragment number "<<aP 618 G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
673 #endif 619 #endif
674 620
675 // fill pre-compound part into the resul 621 // fill pre-compound part into the result, and return
676 for(unsigned int ll=0; ll<aPrecoResult-> 622 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
677 { 623 {
678 theTotalResult->push_back(aPrecoResul 624 theTotalResult->push_back(aPrecoResult->operator[](ll));
679 #ifdef debugPrecoInt 625 #ifdef debugPrecoInt
680 G4cout<<"Target fragment "<<ll<<" 626 G4cout<<"Target fragment "<<ll<<" "
681 <<aPrecoResult->operator[](l 627 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
682 <<aPrecoResult->operator[](l 628 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
683 <<aPrecoResult->operator[](l 629 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
684 <<aPrecoResult->operator[](l 630 <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
685 #endif 631 #endif
686 } 632 }
687 delete aPrecoResult; 633 delete aPrecoResult;
688 } 634 }
689 635
690 //----------------------------------------- 636 //-----------------------------------------------------
691 if((anAb == theProjectileNucleus->GetMassNu 637 if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
692 {Projectile4Momentum = GetPrimaryProjectile 638 {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
693 639
694 #ifdef debugPrecoInt 640 #ifdef debugPrecoInt
695 G4cout<<"Final projectile residual A Z (num << 641 G4cout<<"Final projectile residual A Z E* Pmom Pmag2 "<<anAb<<" "<<aZb<<" "
696 <<aLb<<") "<<exEnergyB<<" "<<Projectile4Mom << 642 <<exEnergyB<<" "<<Projectile4Momentum<<" "
697 <<Projectile4Momen 643 <<Projectile4Momentum.mag2()<<G4endl;
698 #endif 644 #endif
699 645
700 if(0!=anAb) 646 if(0!=anAb)
701 { 647 {
702 // The projectile residual can be a hype << 648 // G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM(); // Uzhi Apr. 2015
703 G4double fMass = 0.0; << 649 // Projectile4Momentum.boost(bstToCM); // Uzhi Apr. 2015
704 if ( aLb > 0 ) { << 650
705 fMass = G4HyperNucleiProperties::GetNuclearM << 651 G4double fMass = G4NucleiProperties::GetNuclearMass(anAb, aZb);
706 } else { <<
707 fMass = G4NucleiProperties::GetNuclearMass(a <<
708 } <<
709 G4double RemnMass=Projectile4Momentum.ma 652 G4double RemnMass=Projectile4Momentum.mag();
710 653
711 if(RemnMass < fMass) 654 if(RemnMass < fMass)
712 { 655 {
713 RemnMass=fMass + exEnergyB; 656 RemnMass=fMass + exEnergyB;
714 Projectile4Momentum.setE(std::sqrt(Pr << 657 Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() + // Uzhi 8.05.2015
715 RemnMass*RemnMass)); << 658 RemnMass*RemnMass)); // Uzhi 8.05.2015
716 } else 659 } else
717 { exEnergyB=RemnMass-fMass;} 660 { exEnergyB=RemnMass-fMass;}
718 661
719 if( exEnergyB < 0.) exEnergyB=0.; 662 if( exEnergyB < 0.) exEnergyB=0.;
720 663
721 G4ThreeVector bstToCM =Projectile4Moment << 664 G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM(); // Uzhi Apr. 2015
722 Projectile4Momentum.boost(bstToCM); << 665 Projectile4Momentum.boost(bstToCM); // Uzhi Apr. 2015
723 666
724 // Need to de-excite the remnant nucleus 667 // Need to de-excite the remnant nucleus
725 G4Fragment anInitialState(anAb, aZb, aLb << 668 G4Fragment anInitialState(anAb, aZb, Projectile4Momentum);
726 anInitialState.SetNumberOfParticles(numb 669 anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
727 anInitialState.SetNumberOfCharged(number 670 anInitialState.SetNumberOfCharged(numberOfChB);
728 anInitialState.SetNumberOfHoles(numberOf 671 anInitialState.SetNumberOfHoles(numberOfHolesB);
729 anInitialState.SetCreatorModelID(secID); <<
730 672
731 G4ReactionProductVector * aPrecoResult = 673 G4ReactionProductVector * aPrecoResult =
732 theDeExcitation->DeExcite(anInitia 674 theDeExcitation->DeExcite(anInitialState);
733 675
734 #ifdef debugPrecoInt 676 #ifdef debugPrecoInt
735 G4cout<<"Projectile fragment number " 677 G4cout<<"Projectile fragment number "<<aPrecoResult->size()<<G4endl;
736 #endif 678 #endif
737 679
738 // fill pre-compound part into the resul 680 // fill pre-compound part into the result, and return
739 for(unsigned int ll=0; ll<aPrecoResult-> 681 for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
740 { 682 {
741 G4LorentzVector tmp=G4LorentzVector(a << 683 G4LorentzVector tmp=G4LorentzVector(aPrecoResult->operator[](ll)->GetMomentum(), // Uzhi 2015
742 a << 684 aPrecoResult->operator[](ll)->GetTotalEnergy());// Uzhi 2015
743 tmp.boost(-bstToCM); // Transformatio << 685 tmp.boost(-bstToCM); // Transformation to the system of original remnant // Uzhi 2015
744 aPrecoResult->operator[](ll)->SetMome << 686 aPrecoResult->operator[](ll)->SetMomentum(tmp.vect()); // Uzhi 2015
745 aPrecoResult->operator[](ll)->SetTota << 687 aPrecoResult->operator[](ll)->SetTotalEnergy(tmp.e()); // Uzhi 2015
746 688
747 if(ProjectileIsAntiNucleus) 689 if(ProjectileIsAntiNucleus)
748 { 690 {
749 const G4ParticleDefinition * aFrag 691 const G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
750 const G4ParticleDefinition * LastF 692 const G4ParticleDefinition * LastFragment=aFragment;
751 if (aFragment == proton) {Las 693 if (aFragment == proton) {LastFragment=G4AntiProton::AntiProtonDefinition();}
752 else if(aFragment == neutron) {Las 694 else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
753 else if(aFragment == lambda) {Las <<
754 else if(aFragment == deuteron){Las 695 else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
755 else if(aFragment == triton) {Las 696 else if(aFragment == triton) {LastFragment=G4AntiTriton::AntiTritonDefinition();}
756 else if(aFragment == He3) {Las 697 else if(aFragment == He3) {LastFragment=G4AntiHe3::AntiHe3Definition();}
757 else if(aFragment == He4) {Las 698 else if(aFragment == He4) {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
758 else {} 699 else {}
759 700
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 701 aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
777 } 702 }
778 703
779 #ifdef debugPrecoInt 704 #ifdef debugPrecoInt
780 G4cout<<"Projectile fragment "<<ll 705 G4cout<<"Projectile fragment "<<ll<<" "
781 <<aPrecoResult->operator[](l 706 <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
782 <<aPrecoResult->operator[](l 707 <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
783 <<aPrecoResult->operator[](l 708 <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
784 <<aPrecoResult->operator[](l 709 <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
785 #endif 710 #endif
786 << 711 //Uzhi
787 theTotalResult->push_back(aPrecoResul 712 theTotalResult->push_back(aPrecoResult->operator[](ll));
788 } 713 }
789 714
790 delete aPrecoResult; 715 delete aPrecoResult;
791 } 716 }
792 717
793 return theTotalResult; 718 return theTotalResult;
794 } 719 }
795 720
796 721
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 <<
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 722