| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class implementation file
29 //
30 // History: first implementation, A. Feli
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_
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 construc
93 G4cerr << " ===============================
94 */
95 }
96
97
98
99 //
100 // Copy constructor
101 //
102
103 G4KineticTrack::G4KineticTrack(const G4Kinetic
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.theProjectileP
123 theCreatorModel = right.GetCreatorModelID();
124 theParentResonanceDef = right.GetParentResona
125 theParentResonanceID = right.GetParentResonan
126
127 ////////////////
128 // DEBUG //
129 ////////////////
130
131 /*
132 G4cerr << G4endl << G4endl << G4endl;
133 G4cerr << " G4KineticTrack copy constructor
134 G4cerr << " ===============================
135 */
136 }
137
138
139 //
140 // By argument constructor
141 //
142
143 G4KineticTrack::G4KineticTrack(const G4Particl
144 G4double aForma
145 const G4ThreeVe
146 const G4Lorentz
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() == theDefin
162 {
163 if(G4UniformRand()<0.5)
164 {
165 theDefinition = G4KaonZeroShort::KaonZer
166 }
167 else
168 {
169 theDefinition = G4KaonZeroLong::KaonZero
170 }
171 }
172
173 //
174 // Get the number of decay channels
175 //
176
177 G4DecayTable* theDecayTable = theDefinition->
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 parti
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 G4doub
205
206 // cout << " ****CONSTR*** ActualMass *****
207 G4int index;
208 for (index = nChannels - 1; index >= 0; --in
209 {
210 G4VDecayChannel* theChannel = theDecayTa
211 G4int nDaughters = theChannel->GetNumber
212 G4double theMotherWidth;
213 if (nDaughters == 2 || nDaughters == 3)
214 {
215 G4double thePoleMass = theDefinitio
216 theMotherWidth = theDefinition->GetP
217 G4double thePoleWidth = theChannel->
218 const G4ParticleDefinition* aDaughte
219 theDaughterMass = new G4double[nDaug
220 theDaughterWidth = new G4double[nDau
221 theDaughterIsShortLived = new G4bool[nDaug
222 for (G4int n = 0; n < nDaughters; ++
223 {
224 aDaughter = theChannel->GetDaugh
225 theDaughterMass[n] = aDaughter->
226 theDaughterWidth[n] = aDaughter-
227 theDaughterIsShortLived[n] = aDaughter
228 }
229
230 //
231 // Check whether both the decay prod
232 //
233
234 G4double theActualMom = 0.0;
235 G4double thePoleMom = 0.0;
236 G4SampleResonance aSampler;
237 if (nDaughters==2)
238 {
239 if ( !theDaughterIsShortLived[0] && !t
240 {
241
242 // G4cout << G4endl << "Bot
243 // " decay
244 // cout << " LB: Both de
245 // cout << " parent:
246 // cout << " particle1:
247 // cout << " particle2:
248
249 theActualMom = EvaluateCMMomentum(theAct
250 theDaughterMass);
251 thePoleMom = EvaluateCMMomentum(thePoleM
252 theDaughterMass);
253 // cout << G4endl;
254 // cout << " LB: ActualMass/Daugh
255 // cout << " LB: ActualMom " << t
256 // cout << " LB: PoleMom " << t
257 // cout << G4endl;
258 }
259 else if ( !theDaughterIsShortLived[0]
260 {
261
262 // G4cout << G4endl << "Onl
263 // cout << " LB: only the first
264 // cout << " parent: " << th
265 // cout << " particle1: " << th
266 // cout << " particle2: " << th
267
268 // global variable definition
269 G4double lowerLimit = aSampler.GetMinimu
270 theActualMom = IntegrateCMMomentum(lower
271 thePoleMom = IntegrateCMMomentum(lowerLi
272 // cout << " LB Parent Mass = " <
273 // cout << " LB Actual Mass = " <
274 // cout << " LB Daughter1 Mass =
275 // cout << " LB Daughter2 Mass =
276 // cout << " The Actual Momentum
277 // cout << " The Pole Momentum
278 // cout << G4endl;
279
280 }
281 else if ( theDaughterIsShortLived[0] &
282 {
283
284 // G4cout << G4endl << "Onl
285 // " decay
286 // cout << " LB: only th
287 // cout << " parent: " << th
288 // cout << " particle1: " << th
289 // cout << " particle2: " << th
290
291 //
292 // Swap the content of the theDa
293 //
294
295 G4SwapObj(theDaughterMass, theDaughterMa
296 G4SwapObj(theDaughterWidth, theDaughterW
297
298 // global variable definition
299 G4double lowerLimit = aSampler.GetMinimu
300 theActualMom = IntegrateCMMomentum(lower
301 thePoleMom = IntegrateCMMomentum(lowerLi
302 // cout << " LB Parent Mass = " <
303 // cout << " LB Actual Mass = " <
304 // cout << " LB Daughter1 Mass =
305 // cout << " LB Daughter2 Mass =
306 // cout << " The Actual Momentum
307 // cout << " The Pole Momentum
308 // cout << G4endl;
309
310 }
311 else if ( theDaughterIsShortLived[0] &
312 {
313
314 // G4cout << G4endl << "Both the
315 // " decay produc
316 // cout << " LB: both decay prod
317 // cout << " parent: " << th
318 // cout << " particle1: " << th
319 // cout << " particle2: " << th
320
321 // global variable definition
322 G4KineticTrack_Gmass = theActualMass;
323 theActualMom = IntegrateCMMomentum2();
324 G4KineticTrack_Gmass = thePoleMass;
325 thePoleMom = IntegrateCMMomentum2();
326 // cout << " LB Parent Mass = " <
327 // cout << " LB Daughter1 Mass =
328 // cout << " LB Daughter2 Mass =
329 // cout << " The Actual Mom
330 // cout << " The Pole Momen
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, theDaughterMa
347 G4SwapObj(theDaughterWidth, theDaughterW
348 }
349 if ( theDaughterIsShortLived[2] )
350 {
351 ++nShortLived;
352 G4SwapObj(theDaughterMass, theDaughterMa
353 G4SwapObj(theDaughterWidth, theDaughterW
354 }
355 if ( nShortLived == 0 )
356 {
357 theDaughterMass[1]+=theDaughterMass[2];
358 theActualMom = EvaluateCMMomentum(theAct
359 theDaughterMass);
360 thePoleMom = EvaluateCMMomentum(thePoleM
361 theDaughterMass);
362 }
363 // else if ( nShortLived == 1 )
364 else if ( nShortLived >= 1 )
365 {
366 // need the shortlived particle in slot
367 G4SwapObj(theDaughterMass, theDaughterMa
368 G4SwapObj(theDaughterWidth, theDaughterW
369 theDaughterMass[0] += theDaughterMass[2]
370 theActualMom = IntegrateCMMomentum(0.0);
371 thePoleMom = IntegrateCMMomentum(0.0, th
372 }
373 // else
374 // {
375 // throw G4HadronicException(__FILE__, __
376 // }
377
378 }
379
380 //if(nDaughters<3) theChannel->GetAngularM
381 G4double theMassRatio = thePoleMass / theA
382 G4double theMomRatio = theActualMom
383 // VI 11.06.2015: for l=0 one not need use
384 //G4double l=0;
385 //theActualWidth[index] = thePoleWid
386 // std::pow(t
387 // (1.2 / (1+
388 theActualWidth[index] = thePoleWidth
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
399 {
400 theMotherWidth = theDefinition->GetPDGWidt
401 theActualWidth[index] = theChannel->GetBR(
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 cons
418 G4cerr << " ===============================
419 */
420
421 }
422
423 G4KineticTrack::G4KineticTrack(G4Nucleon * nuc
424 const G4ThreeVector& aPosition,
425 const G4Lorentz
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()->GetP
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 [] theActualW
451 if (theDaughterMass != 0) delete [] theDaught
452 if (theDaughterWidth != 0) delete [] theDaugh
453 }
454
455
456
457 G4KineticTrack& G4KineticTrack::operator=(cons
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.theFermi3Moment
466 theTotal4Momentum = right.theTotal4Moment
467 theNucleon = right.theNucleon;
468 theStateToNucleus = right.theStateToNucle
469 delete [] theActualWidth;
470 nChannels = right.GetnChannels();
471 theActualWidth = new G4double[nChannels];
472 for (G4int i = 0; i < nChannels; ++i) the
473 theCreatorModel = right.GetCreatorModelID
474 theParentResonanceDef = right.GetParentRe
475 theParentResonanceID = right.GetParentRes
476 }
477 return *this;
478 }
479
480
481
482 G4bool G4KineticTrack::operator==(const G4Kine
483 {
484 return (this == & right);
485 }
486
487
488
489 G4bool G4KineticTrack::operator!=(const G4Kine
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;
504 G4cout << "DECAY Actual Width IND/Actual
505 G4cout << "DECAY Actual Mass " << theAct
506 */
507 const G4ParticleDefinition* thisDefinition =
508 if(!thisDefinition)
509 {
510 G4cerr << "Error condition encountered in
511 G4cerr << " track has no particle definit
512 return nullptr;
513 }
514 G4DecayTable* theDecayTable = thisDefinition
515 if(!theDecayTable)
516 {
517 G4cerr << "Error condition encountered in
518 G4cerr << " particle definition has no de
519 G4cerr << " particle was "<<thisDefinitio
520 return nullptr;
521 }
522
523 G4int chargeBalance = G4lrint(theDefinition->
524 G4int baryonBalance = G4lrint(theDefinition->
525 G4LorentzVector energyMomentumBalance(Get4Mom
526 G4double theTotalActualWidth = this->Evaluate
527 if (theTotalActualWidth !=0)
528 {
529
530 //AR-16Aug2016 : Repeat the sampling of t
531 // kinematically above thre
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 G4dou
551 for (G4int index = nChannels - 1; index
552 {
553 theSumActualWidth += theActualWidth
554 theCumActualWidth[index] = theSumAc
555 // cout << "DECAY Cum. Width " << index
556 }
557 // cout << "DECAY Total Width " << the
558 // cout << "DECAY Total Width " << the
559 G4double r = theTotalActualWidth * G4Un
560 G4VDecayChannel* theDecayChannel(nullpt
561 chosench=-1;
562 for (G4int index = nChannels - 1; index
563 {
564 if (r < theCumActualWidth[index])
565 {
566 theDecayChannel = theDecayTable
567 // cout << "DECAY SELECTED CHANN
568 chosench=index;
569 break;
570 }
571 }
572
573 delete [] theCumActualWidth;
574
575 if(theDecayChannel == nullptr)
576 {
577 G4cerr << "Error condition encountere
578 G4cerr << " decay channel has 0x0 ch
579 G4cerr << " particle was "<<thisDefi
580 G4cerr << " channel index "<< chosen
581 return nullptr;
582 }
583 theParentName = theDecayChannel->GetPar
584 theParentMass = this->GetActualMass();
585 theBR = theActualWidth[chosench];
586 // cout << "**BR*** DECAYNEW " <<
587 theNumberOfDaughters = theDecayChannel-
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 < theNumberOfDaught
598 {
599 const G4ParticleDefinition* aDaughte
600 masses[aD] = aDaughter->GetPDGMass()
601 if ( aDaughter->IsShortLived() )
602 {
603 shortlivedDaughters[numberOfShortlived
604 ++numberOfShortliveds;
605 } else {
606 SumLongLivedMass += aDaughter->GetPDGM
607 }
608
609 }
610 switch (theNumberOfDaughters)
611 {
612 case 0:
613 break;
614 case 1:
615 theDaughtersName1 = theDecayChan
616 theDaughtersName2 = "";
617 theDaughtersName3 = "";
618 theDaughtersName4 = "";
619 break;
620 case 2:
621 theDaughtersName1 = theDecayChan
622 theDaughtersName2 = theDecayChan
623 theDaughtersName3 = "";
624 theDaughtersName4 = "";
625 if ( numberOfShortliveds == 1)
626 { G4SampleResonance aSampler;
627 G4double massmax=theParentMa
628 const G4ParticleDefinition * aDaughter=t
629 masses[shortlivedDaughters[0]]= aS
630 } else if ( numberOfShortliveds == 2)
631 // choose masses one after the oth
632 G4int zero= (G4UniformRand() > 0.5
633 G4int one = 1-zero;
634 G4SampleResonance aSampler;
635 G4double massmax=theParentMass - aSample
636 const G4ParticleDefinition * aDaughter=t
637 masses[shortlivedDaughters[zero]]=aSampl
638 massmax=theParentMass - masses[shortlive
639 aDaughter=theDecayChannel->GetDaughter(s
640 masses[shortlivedDaughters[one]]=aSample
641 }
642 break;
643 case 3:
644 theDaughtersName1 = theDecayChan
645 theDaughtersName2 = theDecayChan
646 theDaughtersName3 = theDecayChan
647 theDaughtersName4 = "";
648 if ( numberOfShortliveds == 1)
649 { G4SampleResonance aSampler;
650 G4double massmax=theParentMa
651 const G4ParticleDefinition * aDaught
652 masses[shortlivedDaughters[0]]= aS
653 }
654 break;
655 default:
656 theDaughtersName1 = theDecayChan
657 theDaughtersName2 = theDecayChan
658 theDaughtersName3 = theDecayChan
659 theDaughtersName4 = theDecayChan
660 if ( numberOfShortliveds == 1)
661 { G4SampleResonance aSampler;
662 G4double massmax=theParentMa
663 const G4ParticleDefinition * aDaught
664 masses[shortlivedDaughters[0]]= aS
665 }
666 if ( theNumberOfDaughters > 4 )
667 G4ExceptionDescription ed;
668 ed << "More than 4 decay daugh
669 G4Exception( "G4KineticTrack::
670 }
671 break;
672 }
673
674 //AR-16Aug2016 : Check whether the s
675 // If this is still no
676 // then the subsequent
677 G4double sumDaughterMasses = 0.0;
678 for (G4int i=0; i < 4; ++i) sumDaugh
679 if ( theParentMass - sumDaughterMass
680
681 } while ( isChannelBelowThreshold && ++lo
682
683 //
684 // Get the decay products List
685 //
686
687 G4GeneralPhaseSpaceDecay thePhaseSpaceDec
688
689
690
691
692 th
693 th
694 th
695 masses);
696 G4DecayProducts* theDecayProducts = thePh
697 if(theDecayProducts == nullptr)
698 {
699 G4ExceptionDescription ed;
700 ed << "Error condition encountered: pha
701 << "\t the decaying particle is: " <
702 << "\t the channel index is: "<< cho
703 << "\t " << theNumberOfDaughters <<
704 << theDaughtersName1 << " " << theDa
705 << theDaughtersName4 << G4endl;
706 G4Exception( "G4KineticTrack::Decay ",
707 return nullptr;
708 }
709
710 //
711 // Create the kinetic track List associat
712 //
713 // For the decay products of hadronic res
714 // the same as their parent
715 G4LorentzRotation toMoving(Get4Momentum()
716 G4DynamicParticle* theDynamicParticle;
717 G4double formationTime = 0.0;
718 G4ThreeVector position = this->GetPositio
719 G4LorentzVector momentum;
720 G4LorentzVector momentumBalanceCMS(0);
721 G4KineticTrackVector* theDecayProductList
722 G4int dEntries = theDecayProducts->entrie
723 const G4ParticleDefinition * aProduct = n
724 // Use the integer round mass in keV to g
725 G4int uniqueID = static_cast< G4int >( ro
726 for (G4int i=dEntries; i > 0; --i)
727 {
728 theDynamicParticle = theDecayProducts->PopP
729 aProduct = theDynamicParticle->GetDef
730 chargeBalance -= G4lrint(aProduct->Ge
731 baryonBalance -= G4lrint(aProduct->Ge
732 momentumBalanceCMS += theDynamicParticle->G
733 momentum = toMoving*theDynamicParticl
734 energyMomentumBalance -= momentum;
735 G4KineticTrack* aDaughter = new G4Kin
736
737
738
739 if (aDaughter != nullptr)
740 {
741 aDaughter->SetCreatorModelID(GetC
742 aDaughter->SetParentResonanceDef(
743 aDaughter->SetParentResonanceID(u
744 }
745 theDecayProductList->push_back(aDaugh
746 delete theDynamicParticle;
747 }
748 delete theDecayProducts;
749 if(std::getenv("DecayEnergyBalanceCheck")
750 std::cout << "DEBUGGING energy 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(G4
765 {
766 G4double mass = theActualMass; /* the actu
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) - (mass
773 ((mass * mass) - (mass1 - xmass) * (mas
774 BrWig(gamma2, mass2, xmass);
775 return result;
776 }
777
778 G4double G4KineticTrack::IntegrandFunction2(G4
779 {
780 G4double mass = theDefinition->GetPDGMass();
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) - (mass
786 ((mass * mass) - (mass1 - xmass) * (mas
787 BrWig(gamma2, mass2, xmass);
788 return result;
789 }
790
791 G4double G4KineticTrack::IntegrandFunction3(G4
792 {
793 const G4double mass = G4KineticTrack_Gmass;
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
800 ((mass * mass) - (G4KineticTrack_xmass1 - x
801 BrWig(gamma2, mass2, xmass);
802 return result;
803 }
804
805 G4double G4KineticTrack::IntegrandFunction4(G4
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(
819 G4double result = BrWig(gamma1, mass1, xmass
820 integral.Simpson(this, &G4KineticTrack::In
821
822 return result;
823 }
824
825 G4double G4KineticTrack::IntegrateCMMomentum(c
826 {
827 const G4double theUpperLimit = theActualMass
828 const G4int nIterations = 100;
829
830 if (theLowerLimit>=theUpperLimit) return 0.0
831
832 G4Integrator<const G4KineticTrack, G4double(
833 G4double theIntegralOverMass2 = integral.Sim
834 theLowerLimit, theUpperLimit, n
835 return theIntegralOverMass2;
836 }
837
838 G4double G4KineticTrack::IntegrateCMMomentum(c
839 {
840 const G4double theUpperLimit = poleMass - th
841 const G4int nIterations = 100;
842
843 if (theLowerLimit>=theUpperLimit) return 0.0
844
845 G4Integrator<const G4KineticTrack, G4double(
846 const G4double theIntegralOverMass2 = integr
847 theLowerLimit, theUpperLimit,
848 return theIntegralOverMass2;
849 }
850
851
852 G4double G4KineticTrack::IntegrateCMMomentum2(
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(
861 G4double theIntegralOverMass2 = integral.Sim
862 theLowerLimit, theUpperLimit, n
863 return theIntegralOverMass2;
864 }
865