Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KineticTrack.cc
Go to the documentation of this file.
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"
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 
62 
63 //
64 // Default constructor
65 //
66 
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 {
84 // DEBUG //
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 
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 
123 // DEBUG //
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 
139  G4double aFormationTime,
140  const G4ThreeVector& aPosition,
141  const 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  const 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  G4int 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  //if(nDaughters<3) theChannel->GetAngularMomentum();
374  G4double theMassRatio = thePoleMass / theActualMass;
375  G4double theMomRatio = theActualMom / thePoleMom;
376  // VI 11.06.2015: for l=0 one not need use pow
377  //G4double l=0;
378  //theActualWidth[index] = thePoleWidth * theMassRatio *
379  // std::pow(theMomRatio, (2 * l + 1)) *
380  // (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l))));
381  theActualWidth[index] = thePoleWidth * theMassRatio *
382  theMomRatio;
383  delete [] theDaughterMass;
384  theDaughterMass = 0;
385  delete [] theDaughterWidth;
386  theDaughterWidth = 0;
387  delete [] theDaughterIsShortLived;
388  theDaughterIsShortLived = 0;
389  }
390 
391  else // nDaughter = 1 ( e.g. K0 decays 50% to Kshort, 50% Klong
392  {
393  theMotherWidth = theDefinition->GetPDGWidth();
394  theActualWidth[index] = theChannel->GetBR()*theMotherWidth;
395  }
396  }
397 
399 // DEBUG //
401 
402 // for (G4int y = nChannels - 1; y >= 0; y--)
403 // {
404 // G4cout << G4endl << theActualWidth[y];
405 // }
406 // G4cout << G4endl << G4endl << G4endl;
407 
408  /*
409  G4cerr << G4endl << G4endl << G4endl;
410  G4cerr << " G4KineticTrack by argument constructor invoked! \n";
411  G4cerr << " =============================================== \n" << G4endl;
412  */
413 
414 }
415 
417  const G4ThreeVector& aPosition,
418  const G4LorentzVector& a4Momentum)
419  : theDefinition(nucleon->GetDefinition()),
420  theFormationTime(0),
421  thePosition(aPosition),
422  the4Momentum(a4Momentum),
423  theFermi3Momentum(nucleon->GetMomentum()),
424  theNucleon(nucleon),
425  nChannels(0),
426  theActualMass(nucleon->GetDefinition()->GetPDGMass()),
427  theActualWidth(0),
428  theDaughterMass(0),
429  theDaughterWidth(0),
430  theStateToNucleus(undefined),
431  theProjectilePotential(0)
432 {
433  theFermi3Momentum.setE(0);
434  Set4Momentum(a4Momentum);
435 }
436 
437 
439 {
440  if (theActualWidth != 0) delete [] theActualWidth;
441  if (theDaughterMass != 0) delete [] theDaughterMass;
442  if (theDaughterWidth != 0) delete [] theDaughterWidth;
443 }
444 
445 
446 
448 {
449  if (this != &right)
450  {
451  theDefinition = right.GetDefinition();
452  theFormationTime = right.GetFormationTime();
453  the4Momentum = right.the4Momentum;
454  the4Momentum = right.GetTrackingMomentum();
455  theFermi3Momentum = right.theFermi3Momentum;
456  theTotal4Momentum = right.theTotal4Momentum;
457  theNucleon=right.theNucleon;
458  theStateToNucleus=right.theStateToNucleus;
459  if (theActualWidth != 0) delete [] theActualWidth;
460  nChannels = right.GetnChannels();
461  theActualWidth = new G4double[nChannels];
462  for ( G4int i = 0; i < nChannels; i++)
463  {
464  theActualWidth[i] = right.theActualWidth[i];
465  }
466  }
467  return *this;
468 }
469 
470 
471 
473 {
474  return (this == & right);
475 }
476 
477 
478 
480 {
481  return (this != & right);
482 }
483 
484 
485 
487 {
488 //
489 // Select a possible decay channel
490 //
491 /*
492  G4int index1;
493  for (index1 = nChannels - 1; index1 >= 0; index1--)
494  G4cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl;
495  G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
496 */
497  const G4ParticleDefinition* thisDefinition = this->GetDefinition();
498  if(!thisDefinition)
499  {
500  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
501  G4cerr << " track has no particle definition associated."<<G4endl;
502  return 0;
503  }
504  G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
505  if(!theDecayTable)
506  {
507  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
508  G4cerr << " particle definiton has no decay table associated."<<G4endl;
509  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
510  return 0;
511  }
512 
513  G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
514  G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
515  G4LorentzVector energyMomentumBalance(Get4Momentum());
516  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
517  if (theTotalActualWidth !=0)
518  {
519 
520  //AR-16Aug2016 : Repeat the sampling of the decay channel until is
521  // kinematically above threshold or a max number of attempts is reached
522  G4bool isChannelBelowThreshold = true;
523  const G4int maxNumberOfLoops = 10000;
524  G4int loopCounter = 0;
525 
526  G4int chosench;
527  G4String theParentName;
528  G4double theParentMass;
529  G4double theBR;
530  G4int theNumberOfDaughters;
531  G4String theDaughtersName1;
532  G4String theDaughtersName2;
533  G4String theDaughtersName3;
534  G4String theDaughtersName4;
535  G4double masses[4]={0.,0.,0.,0.};
536 
537  do {
538 
539  G4int index;
540  G4double theSumActualWidth = 0.0;
541  G4double* theCumActualWidth = new G4double[nChannels];
542  for (index = nChannels - 1; index >= 0; index--)
543  {
544  theSumActualWidth += theActualWidth[index];
545  theCumActualWidth[index] = theSumActualWidth;
546  // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
547  }
548  // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
549  // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
550  G4double r = theTotalActualWidth * G4UniformRand();
551  G4VDecayChannel* theDecayChannel(0);
552  chosench=-1;
553  for (index = nChannels - 1; index >= 0; index--)
554  {
555  if (r < theCumActualWidth[index])
556  {
557  theDecayChannel = theDecayTable->GetDecayChannel(index);
558  // cout << "DECAY SELECTED CHANNEL" << index << G4endl;
559  chosench=index;
560  break;
561  }
562  }
563 
564  delete [] theCumActualWidth;
565 
566  if(!theDecayChannel)
567  {
568  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
569  G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
570  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
571  G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
572  return 0;
573  }
574  theParentName = theDecayChannel->GetParentName();
575  theParentMass = this->GetActualMass();
576  theBR = theActualWidth[index];
577  // cout << "**BR*** DECAYNEW " << theBR << G4endl;
578  theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
579  theDaughtersName1 = "";
580  theDaughtersName2 = "";
581  theDaughtersName3 = "";
582  theDaughtersName4 = "";
583 
584  for (G4int i=0; i < 4; i++) masses[i]=0.;
585  G4int shortlivedDaughters[4];
586  G4int numberOfShortliveds(0);
587  G4double SumLongLivedMass(0);
588  for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
589  {
590  const G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
591  masses[aD] = aDaughter->GetPDGMass();
592  if ( aDaughter->IsShortLived() )
593  {
594  shortlivedDaughters[numberOfShortliveds]=aD;
595  numberOfShortliveds++;
596  } else {
597  SumLongLivedMass += aDaughter->GetPDGMass();
598  }
599 
600  }
601  switch (theNumberOfDaughters)
602  {
603  case 0:
604  break;
605  case 1:
606  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
607  theDaughtersName2 = "";
608  theDaughtersName3 = "";
609  theDaughtersName4 = "";
610  break;
611  case 2:
612  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
613  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
614  theDaughtersName3 = "";
615  theDaughtersName4 = "";
616  if ( numberOfShortliveds == 1)
617  { G4SampleResonance aSampler;
618  G4double massmax=theParentMass - SumLongLivedMass;
619  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
620  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
621  } else if ( numberOfShortliveds == 2) {
622  // choose masses one after the other, start with randomly choosen
623  G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
624  G4int one = 1-zero;
625  G4SampleResonance aSampler;
626  G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
627  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
628  masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
629  massmax=theParentMass - masses[shortlivedDaughters[zero]];
630  aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
631  masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
632  }
633  break;
634  case 3:
635  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
636  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
637  theDaughtersName3 = theDecayChannel->GetDaughterName(2);
638  theDaughtersName4 = "";
639  if ( numberOfShortliveds == 1)
640  { G4SampleResonance aSampler;
641  G4double massmax=theParentMass - SumLongLivedMass;
642  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
643  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
644  }
645  break;
646  default:
647  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
648  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
649  theDaughtersName3 = theDecayChannel->GetDaughterName(2);
650  theDaughtersName4 = theDecayChannel->GetDaughterName(3);
651  if ( numberOfShortliveds == 1)
652  { G4SampleResonance aSampler;
653  G4double massmax=theParentMass - SumLongLivedMass;
654  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
655  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
656  }
657  if ( theNumberOfDaughters > 4 ) {
659  ed << "More than 4 decay daughters: kept only the first 4" << G4endl;
660  G4Exception( "G4KineticTrack::Decay()", "KINTRK5", JustWarning, ed );
661  }
662  break;
663  }
664 
665  //AR-16Aug2016 : Check whether the sum of the masses of the daughters is smaller than the parent mass.
666  // If this is still not the case, but the max number of attempts has been reached,
667  // then the subsequent call thePhaseSpaceDecayChannel.DecayIt() will throw an exception.
668  G4double sumDaughterMasses = 0.0;
669  for (G4int i=0; i < 4; i++) sumDaughterMasses += masses[i];
670  if ( theParentMass - sumDaughterMasses > 0.0 ) isChannelBelowThreshold = false;
671 
672  } while ( isChannelBelowThreshold && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 16.08.2016, A.Ribon */
673 
674 //
675 // Get the decay products List
676 //
677 
678  G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
679  theParentMass,
680  theBR,
681  theNumberOfDaughters,
682  theDaughtersName1,
683  theDaughtersName2,
684  theDaughtersName3,
685  theDaughtersName4,
686  masses);
687  G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
688  if(!theDecayProducts)
689  {
691  ed << "Error condition encountered: phase-space decay failed." << G4endl
692  << "\t the decaying particle is: " << thisDefinition->GetParticleName() << G4endl
693  << "\t the channel index is: "<< chosench << " of "<< nChannels << "channels" << G4endl
694  << "\t " << theNumberOfDaughters << " daughter particles: "
695  << theDaughtersName1 << " " << theDaughtersName2 << " " << theDaughtersName3 << " "
696  << theDaughtersName4 << G4endl;
697  G4Exception( "G4KineticTrack::Decay ", "HAD_KINTRACK_001", JustWarning, ed );
698  return 0;
699  }
700 
701 //
702 // Create the kinetic track List associated to the decay products
703 //
704  G4LorentzRotation toMoving(Get4Momentum().boostVector());
705  G4DynamicParticle* theDynamicParticle;
706  G4double formationTime = 0.0;
708  G4LorentzVector momentum;
709  G4LorentzVector momentumBalanceCMS(0);
710  G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
711  G4int dEntries = theDecayProducts->entries();
712  const G4ParticleDefinition * aProduct = 0;
713  for (G4int i=dEntries; i > 0; i--)
714  {
715  theDynamicParticle = theDecayProducts->PopProducts();
716  aProduct = theDynamicParticle->GetDefinition();
717  chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
718  baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
719  momentumBalanceCMS += theDynamicParticle->Get4Momentum();
720  momentum = toMoving*theDynamicParticle->Get4Momentum();
721  energyMomentumBalance -= momentum;
722  theDecayProductList->push_back(new G4KineticTrack (aProduct,
723  formationTime,
724  position,
725  momentum));
726  delete theDynamicParticle;
727  }
728  delete theDecayProducts;
729  if(getenv("DecayEnergyBalanceCheck"))
730  std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
731  << momentumBalanceCMS << " "
732  <<energyMomentumBalance << " "
733  <<chargeBalance<<" "
734  <<baryonBalance<<" "
735  <<G4endl;
736  return theDecayProductList;
737  }
738  else
739  {
740  return 0;
741  }
742 }
743 
744 G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const
745 {
746  G4double mass = theActualMass; /* the actual mass value */
747  G4double mass1 = theDaughterMass[0];
748  G4double mass2 = theDaughterMass[1];
749  G4double gamma2 = theDaughterWidth[1];
750 
751  G4double result = (1. / (2 * mass)) *
752  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
753  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
754  BrWig(gamma2, mass2, xmass);
755  return result;
756 }
757 
758 G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
759 {
760  G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
761  G4double mass1 = theDaughterMass[0];
762  G4double mass2 = theDaughterMass[1];
763  G4double gamma2 = theDaughterWidth[1];
764  G4double result = (1. / (2 * mass)) *
765  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
766  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
767  BrWig(gamma2, mass2, xmass);
768  return result;
769 }
770 
771 G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
772 {
773  const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
774 // const G4double mass1 = theDaughterMass[0];
775  const G4double mass2 = theDaughterMass[1];
776  const G4double gamma2 = theDaughterWidth[1];
777 
778  const G4double result = (1. / (2 * mass)) *
779  std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
780  ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
781  BrWig(gamma2, mass2, xmass);
782  return result;
783 }
784 
785 G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
786 {
787  const G4double mass = G4KineticTrack_Gmass;
788  const G4double mass1 = theDaughterMass[0];
789  const G4double gamma1 = theDaughterWidth[0];
790 // const G4double mass2 = theDaughterMass[1];
791 
792  G4KineticTrack_xmass1 = xmass;
793 
794  const G4double theLowerLimit = 0.0;
795  const G4double theUpperLimit = mass - xmass;
796  const G4int nIterations = 100;
797 
799  G4double result = BrWig(gamma1, mass1, xmass)*
800  integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
801 
802  return result;
803 }
804 
805 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
806 {
807  const G4double theUpperLimit = theActualMass - theDaughterMass[0];
808  const G4int nIterations = 100;
809 
810  if (theLowerLimit>=theUpperLimit) return 0.0;
811 
813  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
814  theLowerLimit, theUpperLimit, nIterations);
815  return theIntegralOverMass2;
816 }
817 
818 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
819 {
820  const G4double theUpperLimit = poleMass - theDaughterMass[0];
821  const G4int nIterations = 100;
822 
823  if (theLowerLimit>=theUpperLimit) return 0.0;
824 
826  const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
827  theLowerLimit, theUpperLimit, nIterations);
828  return theIntegralOverMass2;
829 }
830 
831 
832 G4double G4KineticTrack::IntegrateCMMomentum2() const
833 {
834  const G4double theLowerLimit = 0.0;
835  const G4double theUpperLimit = theActualMass;
836  const G4int nIterations = 100;
837 
838  if (theLowerLimit>=theUpperLimit) return 0.0;
839 
841  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
842  theLowerLimit, theUpperLimit, nIterations);
843  return theIntegralOverMass2;
844 }
845 
846 
847 
848 
849 
850 
851 
852 
G4double G4ParticleHPJENDLHEData::G4double result
static G4KaonZero * KaonZero()
Definition: G4KaonZero.cc:104
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
G4double GetBR() const
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4int operator==(const G4KineticTrack &right) const
G4int GetnChannels() const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4String & GetParentName() const
G4int GetNumberOfDaughters() const
const G4ThreeVector & GetPosition() const
G4KineticTrackVector * Decay()
G4ParticleDefinition * GetDaughter(G4int anIndex)
static G4ThreadLocal G4double G4KineticTrack_Gmass
static G4KaonZeroLong * KaonZeroLong()
G4VDecayChannel * GetDecayChannel(G4int index) const
G4double GetActualMass() const
G4ParticleDefinition * GetDefinition() const
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4int entries() const
G4DecayTable * GetDecayTable() const
G4int operator!=(const G4KineticTrack &right) const
G4double GetPDGWidth() const
#define G4UniformRand()
Definition: Randomize.hh:97
static G4KaonZeroShort * KaonZeroShort()
G4bool nucleon(G4int ityp)
G4double GetFormationTime() const
bool G4bool
Definition: G4Types.hh:79
G4KineticTrack & operator=(const G4KineticTrack &right)
static G4ThreadLocal G4double G4KineticTrack_xmass1
const G4int n
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4LorentzVector Get4Momentum() const
const G4String & GetDaughterName(G4int anIndex) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetMinimumMass(const G4ParticleDefinition *p) const
virtual G4DecayProducts * DecayIt(G4double mass=0.0)
static G4AntiKaonZero * AntiKaonZero()
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4LorentzVector & GetTrackingMomentum() const
G4DynamicParticle * PopProducts()
#define G4endl
Definition: G4ios.hh:61
G4double Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
G4int entries() const
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
const G4ParticleDefinition * GetDefinition() const
G4double BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
G4GLOB_DLL std::ostream G4cerr