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