| 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 // 26 //
27 // ------------------------------------------- 27 // -----------------------------------------------------------------------------
28 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file
29 // 29 //
30 // History: first implementation, A. Feli 30 // History: first implementation, A. Feliciello, 20th May 1998
31 // ------------------------------------------- 31 // -----------------------------------------------------------------------------
32 32
33 #include "globals.hh" 33 #include "globals.hh"
34 #include "G4ios.hh" 34 #include "G4ios.hh"
35 //#include <cmath> 35 //#include <cmath>
36 36
37 #include "Randomize.hh" 37 #include "Randomize.hh"
38 #include "G4SimpleIntegration.hh" 38 #include "G4SimpleIntegration.hh"
39 #include "G4ThreeVector.hh" 39 #include "G4ThreeVector.hh"
40 #include "G4LorentzVector.hh" 40 #include "G4LorentzVector.hh"
41 #include "G4KineticTrack.hh" 41 #include "G4KineticTrack.hh"
42 #include "G4KineticTrackVector.hh" 42 #include "G4KineticTrackVector.hh"
43 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh"
44 #include "G4DecayTable.hh" 44 #include "G4DecayTable.hh"
45 #include "G4GeneralPhaseSpaceDecay.hh" 45 #include "G4GeneralPhaseSpaceDecay.hh"
46 #include "G4DecayProducts.hh" 46 #include "G4DecayProducts.hh"
47 #include "G4LorentzRotation.hh" 47 #include "G4LorentzRotation.hh"
48 #include "G4SampleResonance.hh" 48 #include "G4SampleResonance.hh"
49 #include "G4Integrator.hh" 49 #include "G4Integrator.hh"
50 #include "G4KaonZero.hh" 50 #include "G4KaonZero.hh"
51 #include "G4KaonZeroShort.hh" 51 #include "G4KaonZeroShort.hh"
52 #include "G4KaonZeroLong.hh" 52 #include "G4KaonZeroLong.hh"
53 #include "G4AntiKaonZero.hh" 53 #include "G4AntiKaonZero.hh"
54 54
55 #include "G4HadTmpUtil.hh" 55 #include "G4HadTmpUtil.hh"
56 56
57 // 57 //
58 // Some static clobal for integration 58 // Some static clobal for integration
59 // 59 //
60 60
61 static G4ThreadLocal G4double G4KineticTrack_ 61 static G4ThreadLocal G4double G4KineticTrack_Gmass, G4KineticTrack_xmass1;
62 62
63 // 63 //
64 // Default constructor 64 // Default constructor
65 // 65 //
66 66
67 G4KineticTrack::G4KineticTrack() : 67 G4KineticTrack::G4KineticTrack() :
68 theDefinition(0), 68 theDefinition(0),
69 theFormationTime(0), 69 theFormationTime(0),
70 thePosition(0), 70 thePosition(0),
71 the4Momentum(0), 71 the4Momentum(0),
72 theFermi3Momentum(0), 72 theFermi3Momentum(0),
73 theTotal4Momentum(0), 73 theTotal4Momentum(0),
74 theNucleon(0), 74 theNucleon(0),
75 nChannels(0), 75 nChannels(0),
76 theActualMass(0), 76 theActualMass(0),
77 theActualWidth(0), 77 theActualWidth(0),
78 theDaughterMass(0), 78 theDaughterMass(0),
79 theDaughterWidth(0), 79 theDaughterWidth(0),
80 theStateToNucleus(undefined), 80 theStateToNucleus(undefined),
81 theProjectilePotential(0), 81 theProjectilePotential(0),
82 theCreatorModel(-1), << 82 theCreatorModel(-1)
83 theParentResonanceDef(nullptr) <<
84 theParentResonanceID(0) <<
85 { 83 {
86 //////////////// 84 ////////////////
87 // DEBUG // 85 // DEBUG //
88 //////////////// 86 ////////////////
89 87
90 /* 88 /*
91 G4cerr << G4endl << G4endl << G4endl; 89 G4cerr << G4endl << G4endl << G4endl;
92 G4cerr << " G4KineticTrack default construc 90 G4cerr << " G4KineticTrack default constructor invoked! \n";
93 G4cerr << " =============================== 91 G4cerr << " =========================================== \n" << G4endl;
94 */ 92 */
95 } 93 }
96 94
97 95
98 96
99 // 97 //
100 // Copy constructor 98 // Copy constructor
101 // 99 //
102 100
103 G4KineticTrack::G4KineticTrack(const G4Kinetic 101 G4KineticTrack::G4KineticTrack(const G4KineticTrack &right) : G4VKineticNucleon()
104 { 102 {
105 theDefinition = right.GetDefinition(); 103 theDefinition = right.GetDefinition();
106 theFormationTime = right.GetFormationTime(); 104 theFormationTime = right.GetFormationTime();
107 thePosition = right.GetPosition(); 105 thePosition = right.GetPosition();
108 the4Momentum = right.GetTrackingMomentum(); 106 the4Momentum = right.GetTrackingMomentum();
109 theFermi3Momentum = right.theFermi3Momentum; 107 theFermi3Momentum = right.theFermi3Momentum;
110 theTotal4Momentum = right.theTotal4Momentum; 108 theTotal4Momentum = right.theTotal4Momentum;
111 theNucleon=right.theNucleon; 109 theNucleon=right.theNucleon;
112 nChannels = right.GetnChannels(); 110 nChannels = right.GetnChannels();
113 theActualMass = right.GetActualMass(); 111 theActualMass = right.GetActualMass();
114 theActualWidth = new G4double[nChannels]; 112 theActualWidth = new G4double[nChannels];
115 for (G4int i = 0; i < nChannels; i++) 113 for (G4int i = 0; i < nChannels; i++)
116 { 114 {
117 theActualWidth[i] = right.theActualWidth[i 115 theActualWidth[i] = right.theActualWidth[i];
118 } 116 }
119 theDaughterMass = 0; 117 theDaughterMass = 0;
120 theDaughterWidth = 0; 118 theDaughterWidth = 0;
121 theStateToNucleus = right.theStateToNucleus; << 119 theStateToNucleus=right.theStateToNucleus;
122 theProjectilePotential = right.theProjectileP << 120 theProjectilePotential=right.theProjectilePotential;
123 theCreatorModel = right.GetCreatorModelID(); 121 theCreatorModel = right.GetCreatorModelID();
124 theParentResonanceDef = right.GetParentResona <<
125 theParentResonanceID = right.GetParentResonan <<
126 <<
127 //////////////// 122 ////////////////
128 // DEBUG // 123 // DEBUG //
129 //////////////// 124 ////////////////
130 125
131 /* 126 /*
132 G4cerr << G4endl << G4endl << G4endl; 127 G4cerr << G4endl << G4endl << G4endl;
133 G4cerr << " G4KineticTrack copy constructor 128 G4cerr << " G4KineticTrack copy constructor invoked! \n";
134 G4cerr << " =============================== 129 G4cerr << " ======================================== \n" <<G4endl;
135 */ 130 */
136 } 131 }
137 132
138 133
139 // 134 //
140 // By argument constructor 135 // By argument constructor
141 // 136 //
142 137
143 G4KineticTrack::G4KineticTrack(const G4Particl 138 G4KineticTrack::G4KineticTrack(const G4ParticleDefinition* aDefinition,
144 G4double aForma 139 G4double aFormationTime,
145 const G4ThreeVe 140 const G4ThreeVector& aPosition,
146 const G4Lorentz 141 const G4LorentzVector& a4Momentum) :
147 theDefinition(aDefinition), 142 theDefinition(aDefinition),
148 theFormationTime(aFormationTime), 143 theFormationTime(aFormationTime),
149 thePosition(aPosition), 144 thePosition(aPosition),
150 the4Momentum(a4Momentum), 145 the4Momentum(a4Momentum),
151 theFermi3Momentum(0), 146 theFermi3Momentum(0),
152 theTotal4Momentum(a4Momentum), 147 theTotal4Momentum(a4Momentum),
153 theNucleon(0), 148 theNucleon(0),
154 theStateToNucleus(undefined), 149 theStateToNucleus(undefined),
155 theProjectilePotential(0), 150 theProjectilePotential(0),
156 theCreatorModel(-1), << 151 theCreatorModel(-1)
157 theParentResonanceDef(nullptr) <<
158 theParentResonanceID(0) <<
159 { 152 {
160 if(G4KaonZero::KaonZero() == theDefinition | 153 if(G4KaonZero::KaonZero() == theDefinition ||
161 G4AntiKaonZero::AntiKaonZero() == theDefin 154 G4AntiKaonZero::AntiKaonZero() == theDefinition)
162 { 155 {
163 if(G4UniformRand()<0.5) 156 if(G4UniformRand()<0.5)
164 { 157 {
165 theDefinition = G4KaonZeroShort::KaonZer 158 theDefinition = G4KaonZeroShort::KaonZeroShort();
166 } 159 }
167 else 160 else
168 { 161 {
169 theDefinition = G4KaonZeroLong::KaonZero 162 theDefinition = G4KaonZeroLong::KaonZeroLong();
170 } 163 }
171 } 164 }
172 165
173 // 166 //
174 // Get the number of decay channels 167 // Get the number of decay channels
175 // 168 //
176 169
177 G4DecayTable* theDecayTable = theDefinition-> 170 G4DecayTable* theDecayTable = theDefinition->GetDecayTable();
178 if (theDecayTable != nullptr) << 171 if (theDecayTable != 0)
179 { 172 {
180 nChannels = theDecayTable->entries(); 173 nChannels = theDecayTable->entries();
181 174
182 } 175 }
183 else 176 else
184 { 177 {
185 nChannels = 0; 178 nChannels = 0;
186 } 179 }
187 180
188 // 181 //
189 // Get the actual mass value 182 // Get the actual mass value
190 // 183 //
191 184
192 theActualMass = GetActualMass(); 185 theActualMass = GetActualMass();
193 186
194 // 187 //
195 // Create an array to Store the actual parti 188 // Create an array to Store the actual partial widths
196 // of the decay channels 189 // of the decay channels
197 // 190 //
198 191
199 theDaughterMass = 0; 192 theDaughterMass = 0;
200 theDaughterWidth = 0; 193 theDaughterWidth = 0;
201 theActualWidth = 0; 194 theActualWidth = 0;
202 G4bool * theDaughterIsShortLived = nullptr; << 195 G4bool * theDaughterIsShortLived = 0;
203 196
204 if(nChannels!=0) theActualWidth = new G4doub 197 if(nChannels!=0) theActualWidth = new G4double[nChannels];
205 198
206 // cout << " ****CONSTR*** ActualMass ***** 199 // cout << " ****CONSTR*** ActualMass ******* " << theActualMass << G4endl;
207 G4int index; 200 G4int index;
208 for (index = nChannels - 1; index >= 0; --in 201 for (index = nChannels - 1; index >= 0; --index)
209 { 202 {
210 G4VDecayChannel* theChannel = theDecayTa 203 G4VDecayChannel* theChannel = theDecayTable->GetDecayChannel(index);
211 G4int nDaughters = theChannel->GetNumber 204 G4int nDaughters = theChannel->GetNumberOfDaughters();
212 G4double theMotherWidth; 205 G4double theMotherWidth;
213 if (nDaughters == 2 || nDaughters == 3) 206 if (nDaughters == 2 || nDaughters == 3)
214 { 207 {
215 G4double thePoleMass = theDefinitio 208 G4double thePoleMass = theDefinition->GetPDGMass();
216 theMotherWidth = theDefinition->GetP 209 theMotherWidth = theDefinition->GetPDGWidth();
217 G4double thePoleWidth = theChannel-> 210 G4double thePoleWidth = theChannel->GetBR()*theMotherWidth;
218 const G4ParticleDefinition* aDaughte 211 const G4ParticleDefinition* aDaughter;
219 theDaughterMass = new G4double[nDaug 212 theDaughterMass = new G4double[nDaughters];
220 theDaughterWidth = new G4double[nDau 213 theDaughterWidth = new G4double[nDaughters];
221 theDaughterIsShortLived = new G4bool[nDaug 214 theDaughterIsShortLived = new G4bool[nDaughters];
222 for (G4int n = 0; n < nDaughters; ++ 215 for (G4int n = 0; n < nDaughters; ++n)
223 { 216 {
224 aDaughter = theChannel->GetDaugh 217 aDaughter = theChannel->GetDaughter(n);
225 theDaughterMass[n] = aDaughter-> 218 theDaughterMass[n] = aDaughter->GetPDGMass();
226 theDaughterWidth[n] = aDaughter- 219 theDaughterWidth[n] = aDaughter->GetPDGWidth();
227 theDaughterIsShortLived[n] = aDaughter 220 theDaughterIsShortLived[n] = aDaughter->IsShortLived();
228 } 221 }
229 222
230 // 223 //
231 // Check whether both the decay prod 224 // Check whether both the decay products are stable
232 // 225 //
233 226
234 G4double theActualMom = 0.0; 227 G4double theActualMom = 0.0;
235 G4double thePoleMom = 0.0; 228 G4double thePoleMom = 0.0;
236 G4SampleResonance aSampler; 229 G4SampleResonance aSampler;
237 if (nDaughters==2) 230 if (nDaughters==2)
238 { 231 {
239 if ( !theDaughterIsShortLived[0] && !t 232 if ( !theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
240 { 233 {
241 234
242 // G4cout << G4endl << "Bot 235 // G4cout << G4endl << "Both the " << nDaughters <<
243 // " decay 236 // " decay products are stable!";
244 // cout << " LB: Both de 237 // cout << " LB: Both decay products STABLE !" << G4endl;
245 // cout << " parent: 238 // cout << " parent: " << theChannel->GetParentName() << G4endl;
246 // cout << " particle1: 239 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
247 // cout << " particle2: 240 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
248 241
249 theActualMom = EvaluateCMMomentum(theAct 242 theActualMom = EvaluateCMMomentum(theActualMass,
250 theDaughterMass); 243 theDaughterMass);
251 thePoleMom = EvaluateCMMomentum(thePoleM 244 thePoleMom = EvaluateCMMomentum(thePoleMass,
252 theDaughterMass); 245 theDaughterMass);
253 // cout << G4endl; 246 // cout << G4endl;
254 // cout << " LB: ActualMass/Daugh 247 // cout << " LB: ActualMass/DaughterMass " << theActualMass << " " << theDaughterMass << G4endl;
255 // cout << " LB: ActualMom " << t 248 // cout << " LB: ActualMom " << theActualMom << G4endl;
256 // cout << " LB: PoleMom " << t 249 // cout << " LB: PoleMom " << thePoleMom << G4endl;
257 // cout << G4endl; 250 // cout << G4endl;
258 } 251 }
259 else if ( !theDaughterIsShortLived[0] 252 else if ( !theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
260 { 253 {
261 254
262 // G4cout << G4endl << "Onl 255 // G4cout << G4endl << "Only the first of the " << nDaughters <<" decay products is stable!";
263 // cout << " LB: only the first 256 // cout << " LB: only the first decay product is STABLE !" << G4endl;
264 // cout << " parent: " << th 257 // cout << " parent: " << theChannel->GetParentName() << G4endl;
265 // cout << " particle1: " << th 258 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
266 // cout << " particle2: " << th 259 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
267 260
268 // global variable definition 261 // global variable definition
269 G4double lowerLimit = aSampler.GetMinimu 262 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(1));
270 theActualMom = IntegrateCMMomentum(lower 263 theActualMom = IntegrateCMMomentum(lowerLimit);
271 thePoleMom = IntegrateCMMomentum(lowerLi 264 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
272 // cout << " LB Parent Mass = " < 265 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
273 // cout << " LB Actual Mass = " < 266 // cout << " LB Actual Mass = " << theActualMass << G4endl;
274 // cout << " LB Daughter1 Mass = 267 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
275 // cout << " LB Daughter2 Mass = 268 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
276 // cout << " The Actual Momentum 269 // cout << " The Actual Momentum = " << theActualMom << G4endl;
277 // cout << " The Pole Momentum 270 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
278 // cout << G4endl; 271 // cout << G4endl;
279 272
280 } 273 }
281 else if ( theDaughterIsShortLived[0] & 274 else if ( theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
282 { 275 {
283 276
284 // G4cout << G4endl << "Onl 277 // G4cout << G4endl << "Only the second of the " << nDaughters <<
285 // " decay 278 // " decay products is stable!";
286 // cout << " LB: only th 279 // cout << " LB: only the second decay product is STABLE !" << G4endl;
287 // cout << " parent: " << th 280 // cout << " parent: " << theChannel->GetParentName() << G4endl;
288 // cout << " particle1: " << th 281 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
289 // cout << " particle2: " << th 282 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
290 283
291 // 284 //
292 // Swap the content of the theDa 285 // Swap the content of the theDaughterMass and theDaughterWidth arrays!!!
293 // 286 //
294 287
295 G4SwapObj(theDaughterMass, theDaughterMa 288 G4SwapObj(theDaughterMass, theDaughterMass + 1);
296 G4SwapObj(theDaughterWidth, theDaughterW 289 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
297 290
298 // global variable definition 291 // global variable definition
299 G4double lowerLimit = aSampler.GetMinimu 292 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(0));
300 theActualMom = IntegrateCMMomentum(lower 293 theActualMom = IntegrateCMMomentum(lowerLimit);
301 thePoleMom = IntegrateCMMomentum(lowerLi 294 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
302 // cout << " LB Parent Mass = " < 295 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
303 // cout << " LB Actual Mass = " < 296 // cout << " LB Actual Mass = " << theActualMass << G4endl;
304 // cout << " LB Daughter1 Mass = 297 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
305 // cout << " LB Daughter2 Mass = 298 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
306 // cout << " The Actual Momentum 299 // cout << " The Actual Momentum = " << theActualMom << G4endl;
307 // cout << " The Pole Momentum 300 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
308 // cout << G4endl; 301 // cout << G4endl;
309 302
310 } 303 }
311 else if ( theDaughterIsShortLived[0] & 304 else if ( theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
312 { 305 {
313 306
314 // G4cout << G4endl << "Both the 307 // G4cout << G4endl << "Both the " << nDaughters <<
315 // " decay produc 308 // " decay products are resonances!";
316 // cout << " LB: both decay prod 309 // cout << " LB: both decay products are RESONANCES !" << G4endl;
317 // cout << " parent: " << th 310 // cout << " parent: " << theChannel->GetParentName() << G4endl;
318 // cout << " particle1: " << th 311 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
319 // cout << " particle2: " << th 312 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
320 313
321 // global variable definition 314 // global variable definition
322 G4KineticTrack_Gmass = theActualMass; 315 G4KineticTrack_Gmass = theActualMass;
323 theActualMom = IntegrateCMMomentum2(); 316 theActualMom = IntegrateCMMomentum2();
324 G4KineticTrack_Gmass = thePoleMass; 317 G4KineticTrack_Gmass = thePoleMass;
325 thePoleMom = IntegrateCMMomentum2(); 318 thePoleMom = IntegrateCMMomentum2();
326 // cout << " LB Parent Mass = " < 319 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
327 // cout << " LB Daughter1 Mass = 320 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
328 // cout << " LB Daughter2 Mass = 321 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
329 // cout << " The Actual Mom 322 // cout << " The Actual Momentum = " << theActualMom << G4endl;
330 // cout << " The Pole Momen 323 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
331 // cout << G4endl; 324 // cout << G4endl;
332 325
333 } 326 }
334 } 327 }
335 else // (nDaughter==3) 328 else // (nDaughter==3)
336 { 329 {
337 330
338 G4int nShortLived = 0; 331 G4int nShortLived = 0;
339 if ( theDaughterIsShortLived[0] ) 332 if ( theDaughterIsShortLived[0] )
340 { 333 {
341 ++nShortLived; 334 ++nShortLived;
342 } 335 }
343 if ( theDaughterIsShortLived[1] ) 336 if ( theDaughterIsShortLived[1] )
344 { 337 {
345 ++nShortLived; 338 ++nShortLived;
346 G4SwapObj(theDaughterMass, theDaughterMa 339 G4SwapObj(theDaughterMass, theDaughterMass + 1);
347 G4SwapObj(theDaughterWidth, theDaughterW 340 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
348 } 341 }
349 if ( theDaughterIsShortLived[2] ) 342 if ( theDaughterIsShortLived[2] )
350 { 343 {
351 ++nShortLived; 344 ++nShortLived;
352 G4SwapObj(theDaughterMass, theDaughterMa 345 G4SwapObj(theDaughterMass, theDaughterMass + 2);
353 G4SwapObj(theDaughterWidth, theDaughterW 346 G4SwapObj(theDaughterWidth, theDaughterWidth + 2);
354 } 347 }
355 if ( nShortLived == 0 ) 348 if ( nShortLived == 0 )
356 { 349 {
357 theDaughterMass[1]+=theDaughterMass[2]; 350 theDaughterMass[1]+=theDaughterMass[2];
358 theActualMom = EvaluateCMMomentum(theAct 351 theActualMom = EvaluateCMMomentum(theActualMass,
359 theDaughterMass); 352 theDaughterMass);
360 thePoleMom = EvaluateCMMomentum(thePoleM 353 thePoleMom = EvaluateCMMomentum(thePoleMass,
361 theDaughterMass); 354 theDaughterMass);
362 } 355 }
363 // else if ( nShortLived == 1 ) 356 // else if ( nShortLived == 1 )
364 else if ( nShortLived >= 1 ) 357 else if ( nShortLived >= 1 )
365 { 358 {
366 // need the shortlived particle in slot 359 // need the shortlived particle in slot 1! (very bad style...)
367 G4SwapObj(theDaughterMass, theDaughterMa 360 G4SwapObj(theDaughterMass, theDaughterMass + 1);
368 G4SwapObj(theDaughterWidth, theDaughterW 361 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
369 theDaughterMass[0] += theDaughterMass[2] 362 theDaughterMass[0] += theDaughterMass[2];
370 theActualMom = IntegrateCMMomentum(0.0); 363 theActualMom = IntegrateCMMomentum(0.0);
371 thePoleMom = IntegrateCMMomentum(0.0, th 364 thePoleMom = IntegrateCMMomentum(0.0, thePoleMass);
372 } 365 }
373 // else 366 // else
374 // { 367 // {
375 // throw G4HadronicException(__FILE__, __ 368 // throw G4HadronicException(__FILE__, __LINE__, ("can't handle more than one shortlived in 3 particle output channel");
376 // } 369 // }
377 370
378 } 371 }
379 372
380 //if(nDaughters<3) theChannel->GetAngularM 373 //if(nDaughters<3) theChannel->GetAngularMomentum();
381 G4double theMassRatio = thePoleMass / theA 374 G4double theMassRatio = thePoleMass / theActualMass;
382 G4double theMomRatio = theActualMom 375 G4double theMomRatio = theActualMom / thePoleMom;
383 // VI 11.06.2015: for l=0 one not need use 376 // VI 11.06.2015: for l=0 one not need use pow
384 //G4double l=0; 377 //G4double l=0;
385 //theActualWidth[index] = thePoleWid 378 //theActualWidth[index] = thePoleWidth * theMassRatio *
386 // std::pow(t 379 // std::pow(theMomRatio, (2 * l + 1)) *
387 // (1.2 / (1+ 380 // (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l))));
388 theActualWidth[index] = thePoleWidth 381 theActualWidth[index] = thePoleWidth * theMassRatio *
389 theMomRatio; 382 theMomRatio;
390 delete [] theDaughterMass; 383 delete [] theDaughterMass;
391 theDaughterMass = nullptr; << 384 theDaughterMass = 0;
392 delete [] theDaughterWidth; 385 delete [] theDaughterWidth;
393 theDaughterWidth = nullptr; << 386 theDaughterWidth = 0;
394 delete [] theDaughterIsShortLived; 387 delete [] theDaughterIsShortLived;
395 theDaughterIsShortLived = nullptr; << 388 theDaughterIsShortLived = 0;
396 } 389 }
397 390
398 else // nDaughter = 1 ( e.g. K0 decays 391 else // nDaughter = 1 ( e.g. K0 decays 50% to Kshort, 50% Klong
399 { 392 {
400 theMotherWidth = theDefinition->GetPDGWidt 393 theMotherWidth = theDefinition->GetPDGWidth();
401 theActualWidth[index] = theChannel->GetBR( 394 theActualWidth[index] = theChannel->GetBR()*theMotherWidth;
402 } 395 }
403 } 396 }
404 397
405 //////////////// 398 ////////////////
406 // DEBUG // 399 // DEBUG //
407 //////////////// 400 ////////////////
408 401
409 // for (G4int y = nChannels - 1; y >= 0; --y) 402 // for (G4int y = nChannels - 1; y >= 0; --y)
410 // { 403 // {
411 // G4cout << G4endl << theActualWidth[y]; 404 // G4cout << G4endl << theActualWidth[y];
412 // } 405 // }
413 // G4cout << G4endl << G4endl << G4endl; 406 // G4cout << G4endl << G4endl << G4endl;
414 407
415 /* 408 /*
416 G4cerr << G4endl << G4endl << G4endl; 409 G4cerr << G4endl << G4endl << G4endl;
417 G4cerr << " G4KineticTrack by argument cons 410 G4cerr << " G4KineticTrack by argument constructor invoked! \n";
418 G4cerr << " =============================== 411 G4cerr << " =============================================== \n" << G4endl;
419 */ 412 */
420 413
421 } 414 }
422 415
423 G4KineticTrack::G4KineticTrack(G4Nucleon * nuc 416 G4KineticTrack::G4KineticTrack(G4Nucleon * nucleon,
424 const G4ThreeVector& aPosition, 417 const G4ThreeVector& aPosition,
425 const G4Lorentz 418 const G4LorentzVector& a4Momentum)
426 : theDefinition(nucleon->GetDefinition() 419 : theDefinition(nucleon->GetDefinition()),
427 theFormationTime(0), 420 theFormationTime(0),
428 thePosition(aPosition), 421 thePosition(aPosition),
429 the4Momentum(a4Momentum), 422 the4Momentum(a4Momentum),
430 theFermi3Momentum(nucleon->GetMomentum()), 423 theFermi3Momentum(nucleon->GetMomentum()),
431 theNucleon(nucleon), 424 theNucleon(nucleon),
432 nChannels(0), 425 nChannels(0),
433 theActualMass(nucleon->GetDefinition()->GetP 426 theActualMass(nucleon->GetDefinition()->GetPDGMass()),
434 theActualWidth(0), 427 theActualWidth(0),
435 theDaughterMass(0), 428 theDaughterMass(0),
436 theDaughterWidth(0), 429 theDaughterWidth(0),
437 theStateToNucleus(undefined), 430 theStateToNucleus(undefined),
438 theProjectilePotential(0), 431 theProjectilePotential(0),
439 theCreatorModel(-1), << 432 theCreatorModel(-1)
440 theParentResonanceDef(nullptr), <<
441 theParentResonanceID(0) <<
442 { 433 {
443 theFermi3Momentum.setE(0); 434 theFermi3Momentum.setE(0);
444 Set4Momentum(a4Momentum); 435 Set4Momentum(a4Momentum);
445 } 436 }
446 437
447 438
448 G4KineticTrack::~G4KineticTrack() 439 G4KineticTrack::~G4KineticTrack()
449 { 440 {
450 if (theActualWidth != 0) delete [] theActualW 441 if (theActualWidth != 0) delete [] theActualWidth;
451 if (theDaughterMass != 0) delete [] theDaught 442 if (theDaughterMass != 0) delete [] theDaughterMass;
452 if (theDaughterWidth != 0) delete [] theDaugh 443 if (theDaughterWidth != 0) delete [] theDaughterWidth;
453 } 444 }
454 445
455 446
456 447
457 G4KineticTrack& G4KineticTrack::operator=(cons 448 G4KineticTrack& G4KineticTrack::operator=(const G4KineticTrack& right)
458 { 449 {
459 if (this != &right) 450 if (this != &right)
460 { 451 {
461 theDefinition = right.GetDefinition(); 452 theDefinition = right.GetDefinition();
462 theFormationTime = right.GetFormationTime 453 theFormationTime = right.GetFormationTime();
463 the4Momentum = right.the4Momentum; 454 the4Momentum = right.the4Momentum;
464 the4Momentum = right.GetTrackingMomentum( 455 the4Momentum = right.GetTrackingMomentum();
465 theFermi3Momentum = right.theFermi3Moment 456 theFermi3Momentum = right.theFermi3Momentum;
466 theTotal4Momentum = right.theTotal4Moment 457 theTotal4Momentum = right.theTotal4Momentum;
467 theNucleon = right.theNucleon; << 458 theNucleon=right.theNucleon;
468 theStateToNucleus = right.theStateToNucle << 459 theStateToNucleus=right.theStateToNucleus;
469 delete [] theActualWidth; << 460 if (theActualWidth != 0) delete [] theActualWidth;
470 nChannels = right.GetnChannels(); 461 nChannels = right.GetnChannels();
471 theActualWidth = new G4double[nChannels]; 462 theActualWidth = new G4double[nChannels];
472 for (G4int i = 0; i < nChannels; ++i) the 463 for (G4int i = 0; i < nChannels; ++i) theActualWidth[i] = right.theActualWidth[i];
473 theCreatorModel = right.GetCreatorModelID 464 theCreatorModel = right.GetCreatorModelID();
474 theParentResonanceDef = right.GetParentRe <<
475 theParentResonanceID = right.GetParentRes <<
476 } 465 }
477 return *this; 466 return *this;
478 } 467 }
479 468
480 469
481 470
482 G4bool G4KineticTrack::operator==(const G4Kine 471 G4bool G4KineticTrack::operator==(const G4KineticTrack& right) const
483 { 472 {
484 return (this == & right); 473 return (this == & right);
485 } 474 }
486 475
487 476
488 477
489 G4bool G4KineticTrack::operator!=(const G4Kine 478 G4bool G4KineticTrack::operator!=(const G4KineticTrack& right) const
490 { 479 {
491 return (this != & right); 480 return (this != & right);
492 } 481 }
493 482
494 483
495 484
496 G4KineticTrackVector* G4KineticTrack::Decay() 485 G4KineticTrackVector* G4KineticTrack::Decay()
497 { 486 {
498 // 487 //
499 // Select a possible decay channel 488 // Select a possible decay channel
500 // 489 //
501 /* 490 /*
502 G4int index1; 491 G4int index1;
503 for (index1 = nChannels - 1; index1 >= 0; 492 for (index1 = nChannels - 1; index1 >= 0; --index1)
504 G4cout << "DECAY Actual Width IND/Actual 493 G4cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl;
505 G4cout << "DECAY Actual Mass " << theAct 494 G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
506 */ 495 */
507 const G4ParticleDefinition* thisDefinition = 496 const G4ParticleDefinition* thisDefinition = this->GetDefinition();
508 if(!thisDefinition) 497 if(!thisDefinition)
509 { 498 {
510 G4cerr << "Error condition encountered in 499 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
511 G4cerr << " track has no particle definit 500 G4cerr << " track has no particle definition associated."<<G4endl;
512 return nullptr; << 501 return 0;
513 } 502 }
514 G4DecayTable* theDecayTable = thisDefinition 503 G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
515 if(!theDecayTable) 504 if(!theDecayTable)
516 { 505 {
517 G4cerr << "Error condition encountered in 506 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
518 G4cerr << " particle definition has no de 507 G4cerr << " particle definition has no decay table associated."<<G4endl;
519 G4cerr << " particle was "<<thisDefinitio 508 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
520 return nullptr; << 509 return 0;
521 } 510 }
522 511
523 G4int chargeBalance = G4lrint(theDefinition-> 512 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
524 G4int baryonBalance = G4lrint(theDefinition-> 513 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
525 G4LorentzVector energyMomentumBalance(Get4Mom 514 G4LorentzVector energyMomentumBalance(Get4Momentum());
526 G4double theTotalActualWidth = this->Evaluate 515 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
527 if (theTotalActualWidth !=0) 516 if (theTotalActualWidth !=0)
528 { 517 {
529 518
530 //AR-16Aug2016 : Repeat the sampling of t 519 //AR-16Aug2016 : Repeat the sampling of the decay channel until is
531 // kinematically above thre 520 // kinematically above threshold or a max number of attempts is reached
532 G4bool isChannelBelowThreshold = true; 521 G4bool isChannelBelowThreshold = true;
533 const G4int maxNumberOfLoops = 10000; 522 const G4int maxNumberOfLoops = 10000;
534 G4int loopCounter = 0; 523 G4int loopCounter = 0;
535 524
536 G4int chosench; 525 G4int chosench;
537 G4String theParentName; 526 G4String theParentName;
538 G4double theParentMass; 527 G4double theParentMass;
539 G4double theBR; 528 G4double theBR;
540 G4int theNumberOfDaughters; 529 G4int theNumberOfDaughters;
541 G4String theDaughtersName1; 530 G4String theDaughtersName1;
542 G4String theDaughtersName2; 531 G4String theDaughtersName2;
543 G4String theDaughtersName3; 532 G4String theDaughtersName3;
544 G4String theDaughtersName4; 533 G4String theDaughtersName4;
545 G4double masses[4]={0.,0.,0.,0.}; 534 G4double masses[4]={0.,0.,0.,0.};
546 535
547 do { 536 do {
548 537
549 G4double theSumActualWidth = 0.0; 538 G4double theSumActualWidth = 0.0;
550 G4double* theCumActualWidth = new G4dou 539 G4double* theCumActualWidth = new G4double[nChannels]{};
551 for (G4int index = nChannels - 1; index 540 for (G4int index = nChannels - 1; index >= 0; --index)
552 { 541 {
553 theSumActualWidth += theActualWidth 542 theSumActualWidth += theActualWidth[index];
554 theCumActualWidth[index] = theSumAc 543 theCumActualWidth[index] = theSumActualWidth;
555 // cout << "DECAY Cum. Width " << index 544 // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
556 } 545 }
557 // cout << "DECAY Total Width " << the 546 // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
558 // cout << "DECAY Total Width " << the 547 // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
559 G4double r = theTotalActualWidth * G4Un 548 G4double r = theTotalActualWidth * G4UniformRand();
560 G4VDecayChannel* theDecayChannel(nullpt << 549 G4VDecayChannel* theDecayChannel(0);
561 chosench=-1; 550 chosench=-1;
562 for (G4int index = nChannels - 1; index 551 for (G4int index = nChannels - 1; index >= 0; --index)
563 { 552 {
564 if (r < theCumActualWidth[index]) 553 if (r < theCumActualWidth[index])
565 { 554 {
566 theDecayChannel = theDecayTable 555 theDecayChannel = theDecayTable->GetDecayChannel(index);
567 // cout << "DECAY SELECTED CHANN 556 // cout << "DECAY SELECTED CHANNEL" << index << G4endl;
568 chosench=index; 557 chosench=index;
569 break; 558 break;
570 } 559 }
571 } 560 }
572 561
573 delete [] theCumActualWidth; 562 delete [] theCumActualWidth;
574 563
575 if(theDecayChannel == nullptr) << 564 if(!theDecayChannel)
576 { 565 {
577 G4cerr << "Error condition encountere 566 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
578 G4cerr << " decay channel has 0x0 ch 567 G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
579 G4cerr << " particle was "<<thisDefi 568 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
580 G4cerr << " channel index "<< chosen 569 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
581 return nullptr; << 570 return 0;
582 } 571 }
583 theParentName = theDecayChannel->GetPar 572 theParentName = theDecayChannel->GetParentName();
584 theParentMass = this->GetActualMass(); 573 theParentMass = this->GetActualMass();
585 theBR = theActualWidth[chosench]; 574 theBR = theActualWidth[chosench];
586 // cout << "**BR*** DECAYNEW " << 575 // cout << "**BR*** DECAYNEW " << theBR << G4endl;
587 theNumberOfDaughters = theDecayChannel- 576 theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
588 theDaughtersName1 = ""; 577 theDaughtersName1 = "";
589 theDaughtersName2 = ""; 578 theDaughtersName2 = "";
590 theDaughtersName3 = ""; 579 theDaughtersName3 = "";
591 theDaughtersName4 = ""; 580 theDaughtersName4 = "";
592 581
593 for (G4int i=0; i < 4; ++i) masses[i]=0 582 for (G4int i=0; i < 4; ++i) masses[i]=0.;
594 G4int shortlivedDaughters[4]; 583 G4int shortlivedDaughters[4];
595 G4int numberOfShortliveds(0); 584 G4int numberOfShortliveds(0);
596 G4double SumLongLivedMass(0); 585 G4double SumLongLivedMass(0);
597 for (G4int aD=0; aD < theNumberOfDaught 586 for (G4int aD=0; aD < theNumberOfDaughters ; ++aD)
598 { 587 {
599 const G4ParticleDefinition* aDaughte 588 const G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
600 masses[aD] = aDaughter->GetPDGMass() 589 masses[aD] = aDaughter->GetPDGMass();
601 if ( aDaughter->IsShortLived() ) 590 if ( aDaughter->IsShortLived() )
602 { 591 {
603 shortlivedDaughters[numberOfShortlived 592 shortlivedDaughters[numberOfShortliveds]=aD;
604 ++numberOfShortliveds; 593 ++numberOfShortliveds;
605 } else { 594 } else {
606 SumLongLivedMass += aDaughter->GetPDGM 595 SumLongLivedMass += aDaughter->GetPDGMass();
607 } 596 }
608 597
609 } 598 }
610 switch (theNumberOfDaughters) 599 switch (theNumberOfDaughters)
611 { 600 {
612 case 0: 601 case 0:
613 break; 602 break;
614 case 1: 603 case 1:
615 theDaughtersName1 = theDecayChan 604 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
616 theDaughtersName2 = ""; 605 theDaughtersName2 = "";
617 theDaughtersName3 = ""; 606 theDaughtersName3 = "";
618 theDaughtersName4 = ""; 607 theDaughtersName4 = "";
619 break; 608 break;
620 case 2: 609 case 2:
621 theDaughtersName1 = theDecayChan 610 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
622 theDaughtersName2 = theDecayChan 611 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
623 theDaughtersName3 = ""; 612 theDaughtersName3 = "";
624 theDaughtersName4 = ""; 613 theDaughtersName4 = "";
625 if ( numberOfShortliveds == 1) 614 if ( numberOfShortliveds == 1)
626 { G4SampleResonance aSampler; 615 { G4SampleResonance aSampler;
627 G4double massmax=theParentMa 616 G4double massmax=theParentMass - SumLongLivedMass;
628 const G4ParticleDefinition * aDaughter=t 617 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
629 masses[shortlivedDaughters[0]]= aS 618 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
630 } else if ( numberOfShortliveds == 2) 619 } else if ( numberOfShortliveds == 2) {
631 // choose masses one after the oth 620 // choose masses one after the other, start with randomly choosen
632 G4int zero= (G4UniformRand() > 0.5 621 G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
633 G4int one = 1-zero; 622 G4int one = 1-zero;
634 G4SampleResonance aSampler; 623 G4SampleResonance aSampler;
635 G4double massmax=theParentMass - aSample 624 G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
636 const G4ParticleDefinition * aDaughter=t 625 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
637 masses[shortlivedDaughters[zero]]=aSampl 626 masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
638 massmax=theParentMass - masses[shortlive 627 massmax=theParentMass - masses[shortlivedDaughters[zero]];
639 aDaughter=theDecayChannel->GetDaughter(s 628 aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
640 masses[shortlivedDaughters[one]]=aSample 629 masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
641 } 630 }
642 break; 631 break;
643 case 3: 632 case 3:
644 theDaughtersName1 = theDecayChan 633 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
645 theDaughtersName2 = theDecayChan 634 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
646 theDaughtersName3 = theDecayChan 635 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
647 theDaughtersName4 = ""; 636 theDaughtersName4 = "";
648 if ( numberOfShortliveds == 1) 637 if ( numberOfShortliveds == 1)
649 { G4SampleResonance aSampler; 638 { G4SampleResonance aSampler;
650 G4double massmax=theParentMa 639 G4double massmax=theParentMass - SumLongLivedMass;
651 const G4ParticleDefinition * aDaught 640 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
652 masses[shortlivedDaughters[0]]= aS 641 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
653 } 642 }
654 break; 643 break;
655 default: 644 default:
656 theDaughtersName1 = theDecayChan 645 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
657 theDaughtersName2 = theDecayChan 646 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
658 theDaughtersName3 = theDecayChan 647 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
659 theDaughtersName4 = theDecayChan 648 theDaughtersName4 = theDecayChannel->GetDaughterName(3);
660 if ( numberOfShortliveds == 1) 649 if ( numberOfShortliveds == 1)
661 { G4SampleResonance aSampler; 650 { G4SampleResonance aSampler;
662 G4double massmax=theParentMa 651 G4double massmax=theParentMass - SumLongLivedMass;
663 const G4ParticleDefinition * aDaught 652 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
664 masses[shortlivedDaughters[0]]= aS 653 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
665 } 654 }
666 if ( theNumberOfDaughters > 4 ) 655 if ( theNumberOfDaughters > 4 ) {
667 G4ExceptionDescription ed; 656 G4ExceptionDescription ed;
668 ed << "More than 4 decay daugh 657 ed << "More than 4 decay daughters: kept only the first 4" << G4endl;
669 G4Exception( "G4KineticTrack:: 658 G4Exception( "G4KineticTrack::Decay()", "KINTRK5", JustWarning, ed );
670 } 659 }
671 break; 660 break;
672 } 661 }
673 662
674 //AR-16Aug2016 : Check whether the s 663 //AR-16Aug2016 : Check whether the sum of the masses of the daughters is smaller than the parent mass.
675 // If this is still no 664 // If this is still not the case, but the max number of attempts has been reached,
676 // then the subsequent 665 // then the subsequent call thePhaseSpaceDecayChannel.DecayIt() will throw an exception.
677 G4double sumDaughterMasses = 0.0; 666 G4double sumDaughterMasses = 0.0;
678 for (G4int i=0; i < 4; ++i) sumDaugh 667 for (G4int i=0; i < 4; ++i) sumDaughterMasses += masses[i];
679 if ( theParentMass - sumDaughterMass 668 if ( theParentMass - sumDaughterMasses > 0.0 ) isChannelBelowThreshold = false;
680 669
681 } while ( isChannelBelowThreshold && ++lo 670 } while ( isChannelBelowThreshold && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 16.08.2016, A.Ribon */
682 671
683 // 672 //
684 // Get the decay products List 673 // Get the decay products List
685 // 674 //
686 675
687 G4GeneralPhaseSpaceDecay thePhaseSpaceDec 676 G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
688 677 theParentMass,
689 678 theBR,
690 679 theNumberOfDaughters,
691 680 theDaughtersName1,
692 th 681 theDaughtersName2,
693 th 682 theDaughtersName3,
694 th 683 theDaughtersName4,
695 masses); 684 masses);
696 G4DecayProducts* theDecayProducts = thePh 685 G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
697 if(theDecayProducts == nullptr) << 686 if(!theDecayProducts)
698 { 687 {
699 G4ExceptionDescription ed; 688 G4ExceptionDescription ed;
700 ed << "Error condition encountered: pha 689 ed << "Error condition encountered: phase-space decay failed." << G4endl
701 << "\t the decaying particle is: " < 690 << "\t the decaying particle is: " << thisDefinition->GetParticleName() << G4endl
702 << "\t the channel index is: "<< cho 691 << "\t the channel index is: "<< chosench << " of "<< nChannels << "channels" << G4endl
703 << "\t " << theNumberOfDaughters << 692 << "\t " << theNumberOfDaughters << " daughter particles: "
704 << theDaughtersName1 << " " << theDa 693 << theDaughtersName1 << " " << theDaughtersName2 << " " << theDaughtersName3 << " "
705 << theDaughtersName4 << G4endl; 694 << theDaughtersName4 << G4endl;
706 G4Exception( "G4KineticTrack::Decay ", 695 G4Exception( "G4KineticTrack::Decay ", "HAD_KINTRACK_001", JustWarning, ed );
707 return nullptr; << 696 return 0;
708 } 697 }
709 698
710 // 699 //
711 // Create the kinetic track List associat 700 // Create the kinetic track List associated to the decay products
712 // 701 //
713 // For the decay products of hadronic res 702 // For the decay products of hadronic resonances, we assign as creator model ID
714 // the same as their parent 703 // the same as their parent
715 G4LorentzRotation toMoving(Get4Momentum() 704 G4LorentzRotation toMoving(Get4Momentum().boostVector());
716 G4DynamicParticle* theDynamicParticle; 705 G4DynamicParticle* theDynamicParticle;
717 G4double formationTime = 0.0; 706 G4double formationTime = 0.0;
718 G4ThreeVector position = this->GetPositio 707 G4ThreeVector position = this->GetPosition();
719 G4LorentzVector momentum; 708 G4LorentzVector momentum;
720 G4LorentzVector momentumBalanceCMS(0); 709 G4LorentzVector momentumBalanceCMS(0);
721 G4KineticTrackVector* theDecayProductList 710 G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
722 G4int dEntries = theDecayProducts->entrie 711 G4int dEntries = theDecayProducts->entries();
723 const G4ParticleDefinition * aProduct = n << 712 const G4ParticleDefinition * aProduct = 0;
724 // Use the integer round mass in keV to g <<
725 G4int uniqueID = static_cast< G4int >( ro <<
726 for (G4int i=dEntries; i > 0; --i) 713 for (G4int i=dEntries; i > 0; --i)
727 { 714 {
728 theDynamicParticle = theDecayProducts->PopP 715 theDynamicParticle = theDecayProducts->PopProducts();
729 aProduct = theDynamicParticle->GetDef 716 aProduct = theDynamicParticle->GetDefinition();
730 chargeBalance -= G4lrint(aProduct->Ge 717 chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
731 baryonBalance -= G4lrint(aProduct->Ge 718 baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
732 momentumBalanceCMS += theDynamicParticle->G 719 momentumBalanceCMS += theDynamicParticle->Get4Momentum();
733 momentum = toMoving*theDynamicParticl 720 momentum = toMoving*theDynamicParticle->Get4Momentum();
734 energyMomentumBalance -= momentum; 721 energyMomentumBalance -= momentum;
735 G4KineticTrack* aDaughter = new G4Kin 722 G4KineticTrack* aDaughter = new G4KineticTrack (aProduct,
736 723 formationTime,
737 724 position,
738 725 momentum);
739 if (aDaughter != nullptr) << 726 if (aDaughter != nullptr) aDaughter->SetCreatorModelID(GetCreatorModelID());
740 { <<
741 aDaughter->SetCreatorModelID(GetC <<
742 aDaughter->SetParentResonanceDef( <<
743 aDaughter->SetParentResonanceID(u <<
744 } <<
745 theDecayProductList->push_back(aDaugh 727 theDecayProductList->push_back(aDaughter);
746 delete theDynamicParticle; 728 delete theDynamicParticle;
747 } 729 }
748 delete theDecayProducts; 730 delete theDecayProducts;
749 if(std::getenv("DecayEnergyBalanceCheck") 731 if(std::getenv("DecayEnergyBalanceCheck"))
750 std::cout << "DEBUGGING energy balance 732 std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
751 << momentumBalanceCMS << " " 733 << momentumBalanceCMS << " "
752 <<energyMomentumBalance << " " 734 <<energyMomentumBalance << " "
753 <<chargeBalance<<" " 735 <<chargeBalance<<" "
754 <<baryonBalance<<" " 736 <<baryonBalance<<" "
755 <<G4endl; 737 <<G4endl;
756 return theDecayProductList; 738 return theDecayProductList;
757 } 739 }
758 else 740 else
759 { 741 {
760 return nullptr; << 742 return 0;
761 } 743 }
762 } 744 }
763 745
764 G4double G4KineticTrack::IntegrandFunction1(G4 746 G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const
765 { 747 {
766 G4double mass = theActualMass; /* the actu 748 G4double mass = theActualMass; /* the actual mass value */
767 G4double mass1 = theDaughterMass[0]; 749 G4double mass1 = theDaughterMass[0];
768 G4double mass2 = theDaughterMass[1]; 750 G4double mass2 = theDaughterMass[1];
769 G4double gamma2 = theDaughterWidth[1]; 751 G4double gamma2 = theDaughterWidth[1];
770 752
771 G4double result = (1. / (2 * mass)) * 753 G4double result = (1. / (2 * mass)) *
772 std::sqrt(std::max((((mass * mass) - (mass 754 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
773 ((mass * mass) - (mass1 - xmass) * (mas 755 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
774 BrWig(gamma2, mass2, xmass); 756 BrWig(gamma2, mass2, xmass);
775 return result; 757 return result;
776 } 758 }
777 759
778 G4double G4KineticTrack::IntegrandFunction2(G4 760 G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
779 { 761 {
780 G4double mass = theDefinition->GetPDGMass(); 762 G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
781 G4double mass1 = theDaughterMass[0]; 763 G4double mass1 = theDaughterMass[0];
782 G4double mass2 = theDaughterMass[1]; 764 G4double mass2 = theDaughterMass[1];
783 G4double gamma2 = theDaughterWidth[1]; 765 G4double gamma2 = theDaughterWidth[1];
784 G4double result = (1. / (2 * mass)) * 766 G4double result = (1. / (2 * mass)) *
785 std::sqrt(std::max((((mass * mass) - (mass 767 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
786 ((mass * mass) - (mass1 - xmass) * (mas 768 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
787 BrWig(gamma2, mass2, xmass); 769 BrWig(gamma2, mass2, xmass);
788 return result; 770 return result;
789 } 771 }
790 772
791 G4double G4KineticTrack::IntegrandFunction3(G4 773 G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
792 { 774 {
793 const G4double mass = G4KineticTrack_Gmass; 775 const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
794 // const G4double mass1 = theDaughterMass[0]; 776 // const G4double mass1 = theDaughterMass[0];
795 const G4double mass2 = theDaughterMass[1]; 777 const G4double mass2 = theDaughterMass[1];
796 const G4double gamma2 = theDaughterWidth[1]; 778 const G4double gamma2 = theDaughterWidth[1];
797 779
798 const G4double result = (1. / (2 * mass)) * 780 const G4double result = (1. / (2 * mass)) *
799 std::sqrt(((mass * mass) - (G4KineticTrack 781 std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
800 ((mass * mass) - (G4KineticTrack_xmass1 - x 782 ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
801 BrWig(gamma2, mass2, xmass); 783 BrWig(gamma2, mass2, xmass);
802 return result; 784 return result;
803 } 785 }
804 786
805 G4double G4KineticTrack::IntegrandFunction4(G4 787 G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
806 { 788 {
807 const G4double mass = G4KineticTrack_Gmass; 789 const G4double mass = G4KineticTrack_Gmass;
808 const G4double mass1 = theDaughterMass[0]; 790 const G4double mass1 = theDaughterMass[0];
809 const G4double gamma1 = theDaughterWidth[0]; 791 const G4double gamma1 = theDaughterWidth[0];
810 // const G4double mass2 = theDaughterMass[1]; 792 // const G4double mass2 = theDaughterMass[1];
811 793
812 G4KineticTrack_xmass1 = xmass; 794 G4KineticTrack_xmass1 = xmass;
813 795
814 const G4double theLowerLimit = 0.0; 796 const G4double theLowerLimit = 0.0;
815 const G4double theUpperLimit = mass - xmass; 797 const G4double theUpperLimit = mass - xmass;
816 const G4int nIterations = 100; 798 const G4int nIterations = 100;
817 799
818 G4Integrator<const G4KineticTrack, G4double( 800 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
819 G4double result = BrWig(gamma1, mass1, xmass 801 G4double result = BrWig(gamma1, mass1, xmass)*
820 integral.Simpson(this, &G4KineticTrack::In 802 integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
821 803
822 return result; 804 return result;
823 } 805 }
824 806
825 G4double G4KineticTrack::IntegrateCMMomentum(c 807 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
826 { 808 {
827 const G4double theUpperLimit = theActualMass 809 const G4double theUpperLimit = theActualMass - theDaughterMass[0];
828 const G4int nIterations = 100; 810 const G4int nIterations = 100;
829 811
830 if (theLowerLimit>=theUpperLimit) return 0.0 812 if (theLowerLimit>=theUpperLimit) return 0.0;
831 813
832 G4Integrator<const G4KineticTrack, G4double( 814 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
833 G4double theIntegralOverMass2 = integral.Sim 815 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
834 theLowerLimit, theUpperLimit, n 816 theLowerLimit, theUpperLimit, nIterations);
835 return theIntegralOverMass2; 817 return theIntegralOverMass2;
836 } 818 }
837 819
838 G4double G4KineticTrack::IntegrateCMMomentum(c 820 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
839 { 821 {
840 const G4double theUpperLimit = poleMass - th 822 const G4double theUpperLimit = poleMass - theDaughterMass[0];
841 const G4int nIterations = 100; 823 const G4int nIterations = 100;
842 824
843 if (theLowerLimit>=theUpperLimit) return 0.0 825 if (theLowerLimit>=theUpperLimit) return 0.0;
844 826
845 G4Integrator<const G4KineticTrack, G4double( 827 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
846 const G4double theIntegralOverMass2 = integr 828 const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
847 theLowerLimit, theUpperLimit, 829 theLowerLimit, theUpperLimit, nIterations);
848 return theIntegralOverMass2; 830 return theIntegralOverMass2;
849 } 831 }
850 832
851 833
852 G4double G4KineticTrack::IntegrateCMMomentum2( 834 G4double G4KineticTrack::IntegrateCMMomentum2() const
853 { 835 {
854 const G4double theLowerLimit = 0.0; 836 const G4double theLowerLimit = 0.0;
855 const G4double theUpperLimit = theActualMass 837 const G4double theUpperLimit = theActualMass;
856 const G4int nIterations = 100; 838 const G4int nIterations = 100;
857 839
858 if (theLowerLimit>=theUpperLimit) return 0.0 840 if (theLowerLimit>=theUpperLimit) return 0.0;
859 841
860 G4Integrator<const G4KineticTrack, G4double( 842 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
861 G4double theIntegralOverMass2 = integral.Sim 843 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
862 theLowerLimit, theUpperLimit, n 844 theLowerLimit, theUpperLimit, nIterations);
863 return theIntegralOverMass2; 845 return theIntegralOverMass2;
864 } 846 }
>> 847
>> 848
>> 849
>> 850
>> 851
>> 852
>> 853
>> 854
865 855