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