1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // ----------------------------------------------------------------------------- |
---|
28 | // GEANT 4 class implementation file |
---|
29 | // |
---|
30 | // History: first implementation, A. Feliciello, 20th May 1998 |
---|
31 | // ----------------------------------------------------------------------------- |
---|
32 | |
---|
33 | #include "globals.hh" |
---|
34 | #include "G4ios.hh" |
---|
35 | #include <cmath> |
---|
36 | |
---|
37 | #include "Randomize.hh" |
---|
38 | #include "G4SimpleIntegration.hh" |
---|
39 | #include "G4ThreeVector.hh" |
---|
40 | #include "G4LorentzVector.hh" |
---|
41 | #include "G4KineticTrack.hh" |
---|
42 | #include "G4KineticTrackVector.hh" |
---|
43 | #include "G4ParticleDefinition.hh" |
---|
44 | #include "G4DecayTable.hh" |
---|
45 | #include "G4GeneralPhaseSpaceDecay.hh" |
---|
46 | #include "G4DecayProducts.hh" |
---|
47 | #include "G4LorentzRotation.hh" |
---|
48 | #include "G4SampleResonance.hh" |
---|
49 | #include "G4Integrator.hh" |
---|
50 | #include "G4KaonZero.hh" |
---|
51 | #include "G4KaonZeroShort.hh" |
---|
52 | #include "G4KaonZeroLong.hh" |
---|
53 | #include "G4AntiKaonZero.hh" |
---|
54 | |
---|
55 | #include "G4HadTmpUtil.hh" |
---|
56 | |
---|
57 | // |
---|
58 | // Some static clobal for integration |
---|
59 | // |
---|
60 | |
---|
61 | static G4double G4KineticTrack_Gmass, G4KineticTrack_xmass1; |
---|
62 | |
---|
63 | // |
---|
64 | // Default constructor |
---|
65 | // |
---|
66 | |
---|
67 | G4KineticTrack::G4KineticTrack() : |
---|
68 | theDefinition(0), |
---|
69 | theFormationTime(0), |
---|
70 | thePosition(0), |
---|
71 | the4Momentum(0), |
---|
72 | theFermi3Momentum(0), |
---|
73 | theTotal4Momentum(0), |
---|
74 | theNucleon(0), |
---|
75 | nChannels(0), |
---|
76 | theActualMass(0), |
---|
77 | theActualWidth(0), |
---|
78 | theDaughterMass(0), |
---|
79 | theDaughterWidth(0), |
---|
80 | theStateToNucleus(undefined), |
---|
81 | theProjectilePotential(0) |
---|
82 | { |
---|
83 | //////////////// |
---|
84 | // DEBUG // |
---|
85 | //////////////// |
---|
86 | |
---|
87 | /* |
---|
88 | G4cerr << G4endl << G4endl << G4endl; |
---|
89 | G4cerr << " G4KineticTrack default constructor invoked! \n"; |
---|
90 | G4cerr << " =========================================== \n" << G4endl; |
---|
91 | */ |
---|
92 | } |
---|
93 | |
---|
94 | |
---|
95 | |
---|
96 | // |
---|
97 | // Copy constructor |
---|
98 | // |
---|
99 | |
---|
100 | G4KineticTrack::G4KineticTrack(const G4KineticTrack &right) : G4VKineticNucleon() |
---|
101 | { |
---|
102 | G4int i; |
---|
103 | theDefinition = right.GetDefinition(); |
---|
104 | theFormationTime = right.GetFormationTime(); |
---|
105 | thePosition = right.GetPosition(); |
---|
106 | the4Momentum = right.GetTrackingMomentum(); |
---|
107 | theFermi3Momentum = right.theFermi3Momentum; |
---|
108 | theTotal4Momentum = right.theTotal4Momentum; |
---|
109 | theNucleon=right.theNucleon; |
---|
110 | nChannels = right.GetnChannels(); |
---|
111 | theActualMass = right.GetActualMass(); |
---|
112 | theActualWidth = new G4double[nChannels]; |
---|
113 | for (i = 0; i < nChannels; i++) |
---|
114 | { |
---|
115 | theActualWidth[i] = right.theActualWidth[i]; |
---|
116 | } |
---|
117 | theDaughterMass = 0; |
---|
118 | theDaughterWidth = 0; |
---|
119 | theStateToNucleus=right.theStateToNucleus; |
---|
120 | theProjectilePotential=right.theProjectilePotential; |
---|
121 | |
---|
122 | //////////////// |
---|
123 | // DEBUG // |
---|
124 | //////////////// |
---|
125 | |
---|
126 | /* |
---|
127 | G4cerr << G4endl << G4endl << G4endl; |
---|
128 | G4cerr << " G4KineticTrack copy constructor invoked! \n"; |
---|
129 | G4cerr << " ======================================== \n" <<G4endl; |
---|
130 | */ |
---|
131 | } |
---|
132 | |
---|
133 | |
---|
134 | // |
---|
135 | // By argument constructor |
---|
136 | // |
---|
137 | |
---|
138 | G4KineticTrack::G4KineticTrack(G4ParticleDefinition* aDefinition, |
---|
139 | G4double aFormationTime, |
---|
140 | G4ThreeVector aPosition, |
---|
141 | G4LorentzVector& a4Momentum) : |
---|
142 | theDefinition(aDefinition), |
---|
143 | theFormationTime(aFormationTime), |
---|
144 | thePosition(aPosition), |
---|
145 | the4Momentum(a4Momentum), |
---|
146 | theFermi3Momentum(0), |
---|
147 | theTotal4Momentum(a4Momentum), |
---|
148 | theNucleon(0), |
---|
149 | theStateToNucleus(undefined), |
---|
150 | theProjectilePotential(0) |
---|
151 | { |
---|
152 | if(G4KaonZero::KaonZero() == theDefinition || |
---|
153 | G4AntiKaonZero::AntiKaonZero() == theDefinition) |
---|
154 | { |
---|
155 | if(G4UniformRand()<0.5) |
---|
156 | { |
---|
157 | theDefinition = G4KaonZeroShort::KaonZeroShort(); |
---|
158 | } |
---|
159 | else |
---|
160 | { |
---|
161 | theDefinition = G4KaonZeroLong::KaonZeroLong(); |
---|
162 | } |
---|
163 | } |
---|
164 | |
---|
165 | // |
---|
166 | // Get the number of decay channels |
---|
167 | // |
---|
168 | |
---|
169 | G4DecayTable* theDecayTable = theDefinition->GetDecayTable(); |
---|
170 | if (theDecayTable != 0) |
---|
171 | { |
---|
172 | nChannels = theDecayTable->entries(); |
---|
173 | |
---|
174 | } |
---|
175 | else |
---|
176 | { |
---|
177 | nChannels = 0; |
---|
178 | } |
---|
179 | |
---|
180 | // |
---|
181 | // Get the actual mass value |
---|
182 | // |
---|
183 | |
---|
184 | theActualMass = GetActualMass(); |
---|
185 | |
---|
186 | // |
---|
187 | // Create an array to Store the actual partial widths |
---|
188 | // of the decay channels |
---|
189 | // |
---|
190 | |
---|
191 | theDaughterMass = 0; |
---|
192 | theDaughterWidth = 0; |
---|
193 | theActualWidth = 0; |
---|
194 | G4bool * theDaughterIsShortLived = 0; |
---|
195 | |
---|
196 | if(nChannels!=0) theActualWidth = new G4double[nChannels]; |
---|
197 | |
---|
198 | // cout << " ****CONSTR*** ActualMass ******* " << theActualMass << G4endl; |
---|
199 | G4int index; |
---|
200 | for (index = nChannels - 1; index >= 0; index--) |
---|
201 | { |
---|
202 | G4VDecayChannel* theChannel = theDecayTable->GetDecayChannel(index); |
---|
203 | G4int nDaughters = theChannel->GetNumberOfDaughters(); |
---|
204 | G4double theMotherWidth; |
---|
205 | if (nDaughters == 2 || nDaughters == 3) |
---|
206 | { |
---|
207 | G4double thePoleMass = theDefinition->GetPDGMass(); |
---|
208 | theMotherWidth = theDefinition->GetPDGWidth(); |
---|
209 | G4double thePoleWidth = theChannel->GetBR()*theMotherWidth; |
---|
210 | G4ParticleDefinition* aDaughter; |
---|
211 | theDaughterMass = new G4double[nDaughters]; |
---|
212 | theDaughterWidth = new G4double[nDaughters]; |
---|
213 | theDaughterIsShortLived = new G4bool[nDaughters]; |
---|
214 | G4int n; |
---|
215 | for (n = 0; n < nDaughters; n++) |
---|
216 | { |
---|
217 | aDaughter = theChannel->GetDaughter(n); |
---|
218 | theDaughterMass[n] = aDaughter->GetPDGMass(); |
---|
219 | theDaughterWidth[n] = aDaughter->GetPDGWidth(); |
---|
220 | theDaughterIsShortLived[n] = aDaughter->IsShortLived(); |
---|
221 | } |
---|
222 | |
---|
223 | // |
---|
224 | // Check whether both the decay products are stable |
---|
225 | // |
---|
226 | |
---|
227 | G4double theActualMom = 0.0; |
---|
228 | G4double thePoleMom = 0.0; |
---|
229 | G4SampleResonance aSampler; |
---|
230 | if (nDaughters==2) |
---|
231 | { |
---|
232 | if ( !theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] ) |
---|
233 | { |
---|
234 | |
---|
235 | // G4cout << G4endl << "Both the " << nDaughters << |
---|
236 | // " decay products are stable!"; |
---|
237 | // cout << " LB: Both decay products STABLE !" << G4endl; |
---|
238 | // cout << " parent: " << theChannel->GetParentName() << G4endl; |
---|
239 | // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl; |
---|
240 | // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl; |
---|
241 | |
---|
242 | theActualMom = EvaluateCMMomentum(theActualMass, |
---|
243 | theDaughterMass); |
---|
244 | thePoleMom = EvaluateCMMomentum(thePoleMass, |
---|
245 | theDaughterMass); |
---|
246 | // cout << G4endl; |
---|
247 | // cout << " LB: ActualMass/DaughterMass " << theActualMass << " " << theDaughterMass << G4endl; |
---|
248 | // cout << " LB: ActualMom " << theActualMom << G4endl; |
---|
249 | // cout << " LB: PoleMom " << thePoleMom << G4endl; |
---|
250 | // cout << G4endl; |
---|
251 | } |
---|
252 | else if ( !theDaughterIsShortLived[0] && theDaughterIsShortLived[1] ) |
---|
253 | { |
---|
254 | |
---|
255 | // G4cout << G4endl << "Only the first of the " << nDaughters <<" decay products is stable!"; |
---|
256 | // cout << " LB: only the first decay product is STABLE !" << G4endl; |
---|
257 | // cout << " parent: " << theChannel->GetParentName() << G4endl; |
---|
258 | // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl; |
---|
259 | // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl; |
---|
260 | |
---|
261 | // global variable definition |
---|
262 | G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(1)); |
---|
263 | theActualMom = IntegrateCMMomentum(lowerLimit); |
---|
264 | thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass); |
---|
265 | // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl; |
---|
266 | // cout << " LB Actual Mass = " << theActualMass << G4endl; |
---|
267 | // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl; |
---|
268 | // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl; |
---|
269 | // cout << " The Actual Momentum = " << theActualMom << G4endl; |
---|
270 | // cout << " The Pole Momentum = " << thePoleMom << G4endl; |
---|
271 | // cout << G4endl; |
---|
272 | |
---|
273 | } |
---|
274 | else if ( theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] ) |
---|
275 | { |
---|
276 | |
---|
277 | // G4cout << G4endl << "Only the second of the " << nDaughters << |
---|
278 | // " decay products is stable!"; |
---|
279 | // cout << " LB: only the second decay product is STABLE !" << G4endl; |
---|
280 | // cout << " parent: " << theChannel->GetParentName() << G4endl; |
---|
281 | // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl; |
---|
282 | // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl; |
---|
283 | |
---|
284 | // |
---|
285 | // Swap the content of the theDaughterMass and theDaughterWidth arrays!!! |
---|
286 | // |
---|
287 | |
---|
288 | G4SwapObj(theDaughterMass, theDaughterMass + 1); |
---|
289 | G4SwapObj(theDaughterWidth, theDaughterWidth + 1); |
---|
290 | |
---|
291 | // global variable definition |
---|
292 | G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(0)); |
---|
293 | theActualMom = IntegrateCMMomentum(lowerLimit); |
---|
294 | thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass); |
---|
295 | // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl; |
---|
296 | // cout << " LB Actual Mass = " << theActualMass << G4endl; |
---|
297 | // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl; |
---|
298 | // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl; |
---|
299 | // cout << " The Actual Momentum = " << theActualMom << G4endl; |
---|
300 | // cout << " The Pole Momentum = " << thePoleMom << G4endl; |
---|
301 | // cout << G4endl; |
---|
302 | |
---|
303 | } |
---|
304 | else if ( theDaughterIsShortLived[0] && theDaughterIsShortLived[1] ) |
---|
305 | { |
---|
306 | |
---|
307 | // G4cout << G4endl << "Both the " << nDaughters << |
---|
308 | // " decay products are resonances!"; |
---|
309 | // cout << " LB: both decay products are RESONANCES !" << G4endl; |
---|
310 | // cout << " parent: " << theChannel->GetParentName() << G4endl; |
---|
311 | // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl; |
---|
312 | // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl; |
---|
313 | |
---|
314 | // global variable definition |
---|
315 | G4KineticTrack_Gmass = theActualMass; |
---|
316 | theActualMom = IntegrateCMMomentum2(); |
---|
317 | G4KineticTrack_Gmass = thePoleMass; |
---|
318 | thePoleMom = IntegrateCMMomentum2(); |
---|
319 | // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl; |
---|
320 | // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl; |
---|
321 | // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl; |
---|
322 | // cout << " The Actual Momentum = " << theActualMom << G4endl; |
---|
323 | // cout << " The Pole Momentum = " << thePoleMom << G4endl; |
---|
324 | // cout << G4endl; |
---|
325 | |
---|
326 | } |
---|
327 | } |
---|
328 | else // (nDaughter==3) |
---|
329 | { |
---|
330 | |
---|
331 | int nShortLived = 0; |
---|
332 | if ( theDaughterIsShortLived[0] ) |
---|
333 | { |
---|
334 | nShortLived++; |
---|
335 | } |
---|
336 | if ( theDaughterIsShortLived[1] ) |
---|
337 | { |
---|
338 | nShortLived++; |
---|
339 | G4SwapObj(theDaughterMass, theDaughterMass + 1); |
---|
340 | G4SwapObj(theDaughterWidth, theDaughterWidth + 1); |
---|
341 | } |
---|
342 | if ( theDaughterIsShortLived[2] ) |
---|
343 | { |
---|
344 | nShortLived++; |
---|
345 | G4SwapObj(theDaughterMass, theDaughterMass + 2); |
---|
346 | G4SwapObj(theDaughterWidth, theDaughterWidth + 2); |
---|
347 | } |
---|
348 | if ( nShortLived == 0 ) |
---|
349 | { |
---|
350 | theDaughterMass[1]+=theDaughterMass[2]; |
---|
351 | theActualMom = EvaluateCMMomentum(theActualMass, |
---|
352 | theDaughterMass); |
---|
353 | thePoleMom = EvaluateCMMomentum(thePoleMass, |
---|
354 | theDaughterMass); |
---|
355 | } |
---|
356 | // else if ( nShortLived == 1 ) |
---|
357 | else if ( nShortLived >= 1 ) |
---|
358 | { |
---|
359 | // need the shortlived particle in slot 1! (very bad style...) |
---|
360 | G4SwapObj(theDaughterMass, theDaughterMass + 1); |
---|
361 | G4SwapObj(theDaughterWidth, theDaughterWidth + 1); |
---|
362 | theDaughterMass[0] += theDaughterMass[2]; |
---|
363 | theActualMom = IntegrateCMMomentum(0.0); |
---|
364 | thePoleMom = IntegrateCMMomentum(0.0, thePoleMass); |
---|
365 | } |
---|
366 | // else |
---|
367 | // { |
---|
368 | // throw G4HadronicException(__FILE__, __LINE__, ("can't handle more than one shortlived in 3 particle output channel"); |
---|
369 | // } |
---|
370 | |
---|
371 | } |
---|
372 | |
---|
373 | G4double l=0; |
---|
374 | //if(nDaughters<3) theChannel->GetAngularMomentum(); |
---|
375 | G4double theMassRatio = thePoleMass / theActualMass; |
---|
376 | G4double theMomRatio = theActualMom / thePoleMom; |
---|
377 | theActualWidth[index] = thePoleWidth * theMassRatio * |
---|
378 | std::pow(theMomRatio, (2 * l + 1)) * |
---|
379 | (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l)))); |
---|
380 | delete [] theDaughterMass; |
---|
381 | theDaughterMass = 0; |
---|
382 | delete [] theDaughterWidth; |
---|
383 | theDaughterWidth = 0; |
---|
384 | delete [] theDaughterIsShortLived; |
---|
385 | theDaughterIsShortLived = 0; |
---|
386 | } |
---|
387 | |
---|
388 | else // nDaughter = 1 ( e.g. K0 decays 50% to Kshort, 50% Klong |
---|
389 | { |
---|
390 | theMotherWidth = theDefinition->GetPDGWidth(); |
---|
391 | theActualWidth[index] = theChannel->GetBR()*theMotherWidth; |
---|
392 | } |
---|
393 | } |
---|
394 | |
---|
395 | //////////////// |
---|
396 | // DEBUG // |
---|
397 | //////////////// |
---|
398 | |
---|
399 | // for (G4int y = nChannels - 1; y >= 0; y--) |
---|
400 | // { |
---|
401 | // G4cout << G4endl << theActualWidth[y]; |
---|
402 | // } |
---|
403 | // G4cout << G4endl << G4endl << G4endl; |
---|
404 | |
---|
405 | /* |
---|
406 | G4cerr << G4endl << G4endl << G4endl; |
---|
407 | G4cerr << " G4KineticTrack by argument constructor invoked! \n"; |
---|
408 | G4cerr << " =============================================== \n" << G4endl; |
---|
409 | */ |
---|
410 | |
---|
411 | } |
---|
412 | |
---|
413 | G4KineticTrack::G4KineticTrack(G4Nucleon * nucleon, |
---|
414 | G4ThreeVector aPosition, |
---|
415 | G4LorentzVector& a4Momentum) |
---|
416 | : theDefinition(nucleon->GetDefinition()), |
---|
417 | theFormationTime(0), |
---|
418 | thePosition(aPosition), |
---|
419 | the4Momentum(a4Momentum), |
---|
420 | theFermi3Momentum(nucleon->GetMomentum()), |
---|
421 | theNucleon(nucleon), |
---|
422 | nChannels(0), |
---|
423 | theActualMass(nucleon->GetDefinition()->GetPDGMass()), |
---|
424 | theActualWidth(0), |
---|
425 | theDaughterMass(0), |
---|
426 | theDaughterWidth(0), |
---|
427 | theStateToNucleus(undefined), |
---|
428 | theProjectilePotential(0) |
---|
429 | { |
---|
430 | theFermi3Momentum.setE(0); |
---|
431 | Set4Momentum(a4Momentum); |
---|
432 | } |
---|
433 | |
---|
434 | |
---|
435 | G4KineticTrack::~G4KineticTrack() |
---|
436 | { |
---|
437 | if (theActualWidth != 0) delete [] theActualWidth; |
---|
438 | if (theDaughterMass != 0) delete [] theDaughterMass; |
---|
439 | if (theDaughterWidth != 0) delete [] theDaughterWidth; |
---|
440 | } |
---|
441 | |
---|
442 | |
---|
443 | |
---|
444 | const G4KineticTrack& G4KineticTrack::operator=(const G4KineticTrack& right) |
---|
445 | { |
---|
446 | G4int i; |
---|
447 | if (this != &right) |
---|
448 | { |
---|
449 | theDefinition = right.GetDefinition(); |
---|
450 | theFormationTime = right.GetFormationTime(); |
---|
451 | the4Momentum = right.the4Momentum; |
---|
452 | the4Momentum = right.GetTrackingMomentum(); |
---|
453 | theFermi3Momentum = right.theFermi3Momentum; |
---|
454 | theTotal4Momentum = right.theTotal4Momentum; |
---|
455 | theNucleon=right.theNucleon; |
---|
456 | theStateToNucleus=right.theStateToNucleus; |
---|
457 | if (theActualWidth != 0) delete [] theActualWidth; |
---|
458 | nChannels = right.GetnChannels(); |
---|
459 | theActualWidth = new G4double[nChannels]; |
---|
460 | for (i = 0; i < nChannels; i++) |
---|
461 | { |
---|
462 | theActualWidth[i] = right.theActualWidth[i]; |
---|
463 | } |
---|
464 | } |
---|
465 | return *this; |
---|
466 | } |
---|
467 | |
---|
468 | |
---|
469 | |
---|
470 | G4int G4KineticTrack::operator==(const G4KineticTrack& right) const |
---|
471 | { |
---|
472 | return (this == & right); |
---|
473 | } |
---|
474 | |
---|
475 | |
---|
476 | |
---|
477 | G4int G4KineticTrack::operator!=(const G4KineticTrack& right) const |
---|
478 | { |
---|
479 | return (this != & right); |
---|
480 | } |
---|
481 | |
---|
482 | |
---|
483 | |
---|
484 | G4KineticTrackVector* G4KineticTrack::Decay() |
---|
485 | { |
---|
486 | // |
---|
487 | // Select a possible decay channel |
---|
488 | // |
---|
489 | // G4int index1; |
---|
490 | // for (index1 = nChannels - 1; index1 >= 0; index1--) |
---|
491 | // cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl; |
---|
492 | // cout << "DECAY Actual Mass " << theActualMass << G4endl; |
---|
493 | |
---|
494 | G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() ); |
---|
495 | G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() ); |
---|
496 | G4LorentzVector energyMomentumBalance(Get4Momentum()); |
---|
497 | G4double theTotalActualWidth = this->EvaluateTotalActualWidth(); |
---|
498 | if (theTotalActualWidth !=0) |
---|
499 | { |
---|
500 | G4int index; |
---|
501 | G4double theSumActualWidth = 0.0; |
---|
502 | G4double* theCumActualWidth = new G4double[nChannels]; |
---|
503 | for (index = nChannels - 1; index >= 0; index--) |
---|
504 | { |
---|
505 | theSumActualWidth += theActualWidth[index]; |
---|
506 | theCumActualWidth[index] = theSumActualWidth; |
---|
507 | // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl; |
---|
508 | } |
---|
509 | // cout << "DECAY Total Width " << theSumActualWidth << G4endl; |
---|
510 | // cout << "DECAY Total Width " << theTotalActualWidth << G4endl; |
---|
511 | G4double r = theTotalActualWidth * G4UniformRand(); |
---|
512 | G4ParticleDefinition* thisDefinition = this->GetDefinition(); |
---|
513 | if(!thisDefinition) |
---|
514 | { |
---|
515 | G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl; |
---|
516 | G4cerr << " track has no particle definition associated."<<G4endl; |
---|
517 | return 0; |
---|
518 | } |
---|
519 | G4DecayTable* theDecayTable = thisDefinition->GetDecayTable(); |
---|
520 | if(!theDecayTable) |
---|
521 | { |
---|
522 | G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl; |
---|
523 | G4cerr << " particle definiton has no decay table associated."<<G4endl; |
---|
524 | G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl; |
---|
525 | return 0; |
---|
526 | } |
---|
527 | G4VDecayChannel* theDecayChannel=NULL; |
---|
528 | for (index = nChannels - 1; index >= 0; index--) |
---|
529 | { |
---|
530 | if (r < theCumActualWidth[index]) |
---|
531 | { |
---|
532 | theDecayChannel = theDecayTable->GetDecayChannel(index); |
---|
533 | // cout << "DECAY SELECTED CHANNEL" << index << G4endl; |
---|
534 | chosench=index; |
---|
535 | break; |
---|
536 | } |
---|
537 | } |
---|
538 | |
---|
539 | if(!theDecayChannel) |
---|
540 | { |
---|
541 | G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl; |
---|
542 | G4cerr << " decay channel has 0x0 channel associated."<<G4endl; |
---|
543 | G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl; |
---|
544 | G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl; |
---|
545 | return 0; |
---|
546 | } |
---|
547 | G4String theParentName = theDecayChannel->GetParentName(); |
---|
548 | G4double theParentMass = this->GetActualMass(); |
---|
549 | G4double theBR = theActualWidth[index]; |
---|
550 | // cout << "**BR*** DECAYNEW " << theBR << G4endl; |
---|
551 | G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters(); |
---|
552 | G4String theDaughtersName1 = ""; |
---|
553 | G4String theDaughtersName2 = ""; |
---|
554 | G4String theDaughtersName3 = ""; |
---|
555 | |
---|
556 | G4double masses[3]={0.,0.,0.}; |
---|
557 | G4int shortlivedDaughters[3]; |
---|
558 | G4int numberOfShortliveds(0); |
---|
559 | G4double SumLongLivedMass(0); |
---|
560 | for (G4int aD=0; aD < theNumberOfDaughters ; aD++) |
---|
561 | { |
---|
562 | G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD); |
---|
563 | masses[aD] = aDaughter->GetPDGMass(); |
---|
564 | if ( aDaughter->IsShortLived() ) |
---|
565 | { |
---|
566 | shortlivedDaughters[numberOfShortliveds]=aD; |
---|
567 | numberOfShortliveds++; |
---|
568 | } else { |
---|
569 | SumLongLivedMass += aDaughter->GetPDGMass(); |
---|
570 | } |
---|
571 | |
---|
572 | } |
---|
573 | switch (theNumberOfDaughters) |
---|
574 | { |
---|
575 | case 0: |
---|
576 | break; |
---|
577 | case 1: |
---|
578 | theDaughtersName1 = theDecayChannel->GetDaughterName(0); |
---|
579 | theDaughtersName2 = ""; |
---|
580 | theDaughtersName3 = ""; |
---|
581 | break; |
---|
582 | case 2: |
---|
583 | theDaughtersName1 = theDecayChannel->GetDaughterName(0); |
---|
584 | theDaughtersName2 = theDecayChannel->GetDaughterName(1); |
---|
585 | theDaughtersName3 = ""; |
---|
586 | if ( numberOfShortliveds == 1) |
---|
587 | { G4SampleResonance aSampler; |
---|
588 | G4double massmax=theParentMass - SumLongLivedMass; |
---|
589 | G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]); |
---|
590 | masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax); |
---|
591 | } else if ( numberOfShortliveds == 2) { |
---|
592 | // choose masses one after the other, start with randomly choosen |
---|
593 | G4int zero= (G4UniformRand() > 0.5) ? 0 : 1; |
---|
594 | G4int one = 1-zero; |
---|
595 | G4SampleResonance aSampler; |
---|
596 | G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one])); |
---|
597 | G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]); |
---|
598 | masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax); |
---|
599 | massmax=theParentMass - masses[shortlivedDaughters[zero]]; |
---|
600 | aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]); |
---|
601 | masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax); |
---|
602 | } |
---|
603 | break; |
---|
604 | default: |
---|
605 | theDaughtersName1 = theDecayChannel->GetDaughterName(0); |
---|
606 | theDaughtersName2 = theDecayChannel->GetDaughterName(1); |
---|
607 | theDaughtersName3 = theDecayChannel->GetDaughterName(2); |
---|
608 | if ( numberOfShortliveds == 1) |
---|
609 | { G4SampleResonance aSampler; |
---|
610 | G4double massmax=theParentMass - SumLongLivedMass; |
---|
611 | G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]); |
---|
612 | masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax); |
---|
613 | } |
---|
614 | break; |
---|
615 | } |
---|
616 | |
---|
617 | // |
---|
618 | // Get the decay products List |
---|
619 | // |
---|
620 | |
---|
621 | G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName, |
---|
622 | theParentMass, |
---|
623 | theBR, |
---|
624 | theNumberOfDaughters, |
---|
625 | theDaughtersName1, |
---|
626 | theDaughtersName2, |
---|
627 | theDaughtersName3, |
---|
628 | masses); |
---|
629 | G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt(); |
---|
630 | if(!theDecayProducts) |
---|
631 | { |
---|
632 | G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl; |
---|
633 | G4cerr << " phase-space decay failed."<<G4endl; |
---|
634 | G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl; |
---|
635 | G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl; |
---|
636 | G4cerr << " "<<theNumberOfDaughters<< " Daughter particles: " |
---|
637 | << theDaughtersName1<<" "<<theDaughtersName2<<" "<<theDaughtersName3<<G4endl; |
---|
638 | return 0; |
---|
639 | } |
---|
640 | |
---|
641 | // |
---|
642 | // Create the kinetic track List associated to the decay products |
---|
643 | // |
---|
644 | G4LorentzRotation toMoving(Get4Momentum().boostVector()); |
---|
645 | G4DynamicParticle* theDynamicParticle; |
---|
646 | G4double formationTime = 0.0; |
---|
647 | G4ThreeVector position = this->GetPosition(); |
---|
648 | G4LorentzVector momentum; |
---|
649 | G4LorentzVector momentumBalanceCMS(0); |
---|
650 | G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector; |
---|
651 | G4int dEntries = theDecayProducts->entries(); |
---|
652 | G4ParticleDefinition * aProduct = 0; |
---|
653 | for (G4int i=dEntries; i > 0; i--) |
---|
654 | { |
---|
655 | theDynamicParticle = theDecayProducts->PopProducts(); |
---|
656 | aProduct = theDynamicParticle->GetDefinition(); |
---|
657 | chargeBalance -= G4lrint(aProduct->GetPDGCharge() ); |
---|
658 | baryonBalance -= G4lrint(aProduct->GetBaryonNumber() ); |
---|
659 | momentumBalanceCMS += theDynamicParticle->Get4Momentum(); |
---|
660 | momentum = toMoving*theDynamicParticle->Get4Momentum(); |
---|
661 | energyMomentumBalance -= momentum; |
---|
662 | theDecayProductList->push_back(new G4KineticTrack (aProduct, |
---|
663 | formationTime, |
---|
664 | position, |
---|
665 | momentum)); |
---|
666 | delete theDynamicParticle; |
---|
667 | } |
---|
668 | delete theDecayProducts; |
---|
669 | delete [] theCumActualWidth; |
---|
670 | if(getenv("DecayEnergyBalanceCheck")) |
---|
671 | std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : " |
---|
672 | << momentumBalanceCMS << " " |
---|
673 | <<energyMomentumBalance << " " |
---|
674 | <<chargeBalance<<" " |
---|
675 | <<baryonBalance<<" " |
---|
676 | <<G4endl; |
---|
677 | return theDecayProductList; |
---|
678 | } |
---|
679 | else |
---|
680 | { |
---|
681 | return 0; |
---|
682 | } |
---|
683 | } |
---|
684 | |
---|
685 | G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const |
---|
686 | { |
---|
687 | G4double mass = theActualMass; /* the actual mass value */ |
---|
688 | G4double mass1 = theDaughterMass[0]; |
---|
689 | G4double mass2 = theDaughterMass[1]; |
---|
690 | G4double gamma2 = theDaughterWidth[1]; |
---|
691 | |
---|
692 | G4double result = (1. / (2 * mass)) * |
---|
693 | std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) * |
---|
694 | ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) * |
---|
695 | BrWig(gamma2, mass2, xmass); |
---|
696 | return result; |
---|
697 | } |
---|
698 | |
---|
699 | G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const |
---|
700 | { |
---|
701 | G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */ |
---|
702 | G4double mass1 = theDaughterMass[0]; |
---|
703 | G4double mass2 = theDaughterMass[1]; |
---|
704 | G4double gamma2 = theDaughterWidth[1]; |
---|
705 | G4double result = (1. / (2 * mass)) * |
---|
706 | std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) * |
---|
707 | ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) * |
---|
708 | BrWig(gamma2, mass2, xmass); |
---|
709 | return result; |
---|
710 | } |
---|
711 | |
---|
712 | G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const |
---|
713 | { |
---|
714 | const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */ |
---|
715 | // const G4double mass1 = theDaughterMass[0]; |
---|
716 | const G4double mass2 = theDaughterMass[1]; |
---|
717 | const G4double gamma2 = theDaughterWidth[1]; |
---|
718 | |
---|
719 | const G4double result = (1. / (2 * mass)) * |
---|
720 | std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) * |
---|
721 | ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) * |
---|
722 | BrWig(gamma2, mass2, xmass); |
---|
723 | return result; |
---|
724 | } |
---|
725 | |
---|
726 | G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const |
---|
727 | { |
---|
728 | const G4double mass = G4KineticTrack_Gmass; |
---|
729 | const G4double mass1 = theDaughterMass[0]; |
---|
730 | const G4double gamma1 = theDaughterWidth[0]; |
---|
731 | // const G4double mass2 = theDaughterMass[1]; |
---|
732 | |
---|
733 | G4KineticTrack_xmass1 = xmass; |
---|
734 | |
---|
735 | const G4double theLowerLimit = 0.0; |
---|
736 | const G4double theUpperLimit = mass - xmass; |
---|
737 | const G4int nIterations = 100; |
---|
738 | |
---|
739 | G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral; |
---|
740 | G4double result = BrWig(gamma1, mass1, xmass)* |
---|
741 | integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations); |
---|
742 | |
---|
743 | return result; |
---|
744 | } |
---|
745 | |
---|
746 | G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const |
---|
747 | { |
---|
748 | const G4double theUpperLimit = theActualMass - theDaughterMass[0]; |
---|
749 | const G4int nIterations = 100; |
---|
750 | |
---|
751 | if (theLowerLimit>=theUpperLimit) return 0.0; |
---|
752 | |
---|
753 | G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral; |
---|
754 | G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1, |
---|
755 | theLowerLimit, theUpperLimit, nIterations); |
---|
756 | return theIntegralOverMass2; |
---|
757 | } |
---|
758 | |
---|
759 | G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const |
---|
760 | { |
---|
761 | const G4double theUpperLimit = poleMass - theDaughterMass[0]; |
---|
762 | const G4int nIterations = 100; |
---|
763 | |
---|
764 | if (theLowerLimit>=theUpperLimit) return 0.0; |
---|
765 | |
---|
766 | G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral; |
---|
767 | const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2, |
---|
768 | theLowerLimit, theUpperLimit, nIterations); |
---|
769 | return theIntegralOverMass2; |
---|
770 | } |
---|
771 | |
---|
772 | |
---|
773 | G4double G4KineticTrack::IntegrateCMMomentum2() const |
---|
774 | { |
---|
775 | const G4double theLowerLimit = 0.0; |
---|
776 | const G4double theUpperLimit = theActualMass; |
---|
777 | const G4int nIterations = 100; |
---|
778 | |
---|
779 | if (theLowerLimit>=theUpperLimit) return 0.0; |
---|
780 | |
---|
781 | G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral; |
---|
782 | G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4, |
---|
783 | theLowerLimit, theUpperLimit, nIterations); |
---|
784 | return theIntegralOverMass2; |
---|
785 | } |
---|
786 | |
---|
787 | |
---|
788 | |
---|
789 | |
---|
790 | |
---|
791 | |
---|
792 | |
---|
793 | |
---|