source: trunk/source/processes/hadronic/models/util/src/G4KineticTrack.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 28.2 KB
Line 
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
61static G4double  G4KineticTrack_Gmass, G4KineticTrack_xmass1;
62
63//
64//   Default constructor
65//
66
67G4KineticTrack::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
100G4KineticTrack::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
138G4KineticTrack::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
413G4KineticTrack::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
435G4KineticTrack::~G4KineticTrack()
436{
437 if (theActualWidth != 0) delete [] theActualWidth;
438 if (theDaughterMass != 0) delete [] theDaughterMass;
439 if (theDaughterWidth != 0) delete [] theDaughterWidth;
440}
441
442
443
444const 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
470G4int G4KineticTrack::operator==(const G4KineticTrack& right) const
471{
472 return (this == & right);
473}
474
475
476
477G4int G4KineticTrack::operator!=(const G4KineticTrack& right) const
478{
479 return (this != & right);
480}
481
482
483
484G4KineticTrackVector* G4KineticTrack::Decay()
485{
486//
487//   Select a possible decay channel
488//
489/*
490    G4int index1;
491    for (index1 = nChannels - 1; index1 >= 0; index1--)
492      G4cout << "DECAY Actual Width IND/ActualW " << index1 << "  " << theActualWidth[index1] << G4endl;
493      G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
494*/
495  G4ParticleDefinition* thisDefinition = this->GetDefinition();
496  if(!thisDefinition)
497  {
498    G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
499    G4cerr << "  track has no particle definition associated."<<G4endl;
500    return 0;
501  }
502  G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
503  if(!theDecayTable)
504  {
505    G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
506    G4cerr << "  particle definiton has no decay table associated."<<G4endl;
507    G4cerr << "  particle was "<<thisDefinition->GetParticleName()<<G4endl;
508    return 0;
509  }
510 
511 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );     
512 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
513 G4LorentzVector energyMomentumBalance(Get4Momentum());
514 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
515 if (theTotalActualWidth !=0)
516    {
517     G4int index;
518     G4double theSumActualWidth = 0.0;
519     G4double* theCumActualWidth = new G4double[nChannels];
520     for (index = nChannels - 1; index >= 0; index--)
521        {
522         theSumActualWidth += theActualWidth[index];
523         theCumActualWidth[index] = theSumActualWidth;
524         //      cout << "DECAY Cum. Width " << index << "  " << theCumActualWidth[index] << G4endl;
525        }
526     //  cout << "DECAY Total Width " << theSumActualWidth << G4endl;
527     //  cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
528     G4double r = theTotalActualWidth * G4UniformRand();
529     G4VDecayChannel* theDecayChannel(0);
530     for (index = nChannels - 1; index >= 0; index--)
531        {
532         if (r < theCumActualWidth[index])
533            {
534             theDecayChannel = theDecayTable->GetDecayChannel(index);
535             //      cout << "DECAY SELECTED CHANNEL" << index << G4endl;
536             chosench=index;
537             break; 
538            }
539        }
540
541     delete [] theCumActualWidth;
542   
543     if(!theDecayChannel)
544     {
545       G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
546       G4cerr << "  decay channel has 0x0 channel associated."<<G4endl;
547       G4cerr << "  particle was "<<thisDefinition->GetParticleName()<<G4endl;
548       G4cerr << "  channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
549       return 0;
550     }
551     G4String theParentName = theDecayChannel->GetParentName();
552     G4double theParentMass = this->GetActualMass();
553     G4double theBR = theActualWidth[index];
554     //     cout << "**BR*** DECAYNEW  " << theBR << G4endl;
555     G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
556     G4String theDaughtersName1 = "";
557     G4String theDaughtersName2 = "";
558     G4String theDaughtersName3 = "";
559     
560     G4double masses[3]={0.,0.,0.};
561     G4int shortlivedDaughters[3];
562     G4int numberOfShortliveds(0);
563     G4double SumLongLivedMass(0);
564     for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
565     {
566        G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
567        masses[aD] = aDaughter->GetPDGMass();
568        if ( aDaughter->IsShortLived() ) 
569        {
570            shortlivedDaughters[numberOfShortliveds]=aD;
571            numberOfShortliveds++;
572        } else {
573            SumLongLivedMass += aDaughter->GetPDGMass();
574        }
575               
576     }   
577     switch (theNumberOfDaughters)
578        {
579         case 0:
580            break;
581         case 1:
582            theDaughtersName1 = theDecayChannel->GetDaughterName(0);
583            theDaughtersName2 = "";
584            theDaughtersName3 = "";
585            break;
586         case 2:   
587            theDaughtersName1 = theDecayChannel->GetDaughterName(0);       
588            theDaughtersName2 = theDecayChannel->GetDaughterName(1);
589            theDaughtersName3 = "";
590            if (  numberOfShortliveds == 1) 
591            {   G4SampleResonance aSampler;
592                G4double massmax=theParentMass - SumLongLivedMass;
593                G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);     
594                masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
595            } else if (  numberOfShortliveds == 2) {
596                // choose masses one after the other, start with randomly choosen
597                G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
598                G4int one = 1-zero;
599                G4SampleResonance aSampler;
600                G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
601                G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);           
602                masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
603                massmax=theParentMass - masses[shortlivedDaughters[zero]];
604                aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
605                masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
606            }
607            break;             
608         default:   
609            theDaughtersName1 = theDecayChannel->GetDaughterName(0);
610            theDaughtersName2 = theDecayChannel->GetDaughterName(1);
611            theDaughtersName3 = theDecayChannel->GetDaughterName(2);
612            if (  numberOfShortliveds == 1) 
613            {   G4SampleResonance aSampler;
614                G4double massmax=theParentMass - SumLongLivedMass;
615                G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);     
616                masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
617            }
618            break;
619        }
620
621//     
622//      Get the decay products List
623//
624     
625     G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
626                                                        theParentMass,
627                                                        theBR,
628                                                        theNumberOfDaughters,
629                                                        theDaughtersName1,                 
630                                                        theDaughtersName2,
631                                                        theDaughtersName3,
632                                                        masses);
633     G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
634     if(!theDecayProducts)
635     {
636       G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
637       G4cerr << "  phase-space decay failed."<<G4endl;
638       G4cerr << "  particle was "<<thisDefinition->GetParticleName()<<G4endl;
639       G4cerr << "  channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
640       G4cerr << "  "<<theNumberOfDaughters<< " Daughter particles: "
641              << theDaughtersName1<<" "<<theDaughtersName2<<" "<<theDaughtersName3<<G4endl;
642       return 0;
643     }
644                                                       
645//
646//      Create the kinetic track List associated to the decay products
647//
648     G4LorentzRotation toMoving(Get4Momentum().boostVector());
649     G4DynamicParticle* theDynamicParticle;
650     G4double formationTime = 0.0;
651     G4ThreeVector position = this->GetPosition();
652     G4LorentzVector momentum;
653     G4LorentzVector momentumBalanceCMS(0);
654     G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
655     G4int dEntries = theDecayProducts->entries();
656     G4ParticleDefinition * aProduct = 0;
657     for (G4int i=dEntries; i > 0; i--)
658        {
659         theDynamicParticle = theDecayProducts->PopProducts();
660         aProduct = theDynamicParticle->GetDefinition();
661         chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
662         baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
663         momentumBalanceCMS += theDynamicParticle->Get4Momentum();
664         momentum = toMoving*theDynamicParticle->Get4Momentum();
665         energyMomentumBalance -= momentum;
666         theDecayProductList->push_back(new G4KineticTrack (aProduct,
667                                                         formationTime,
668                                                         position,
669                                                         momentum));
670         delete theDynamicParticle;
671        }
672     delete theDecayProducts;
673     if(getenv("DecayEnergyBalanceCheck"))
674       std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
675                 << momentumBalanceCMS << " " 
676                 <<energyMomentumBalance << " " 
677                 <<chargeBalance<<" "
678                 <<baryonBalance<<" "
679                 <<G4endl;
680     return theDecayProductList;
681    }
682 else
683    {
684     return 0;
685    }
686}
687
688G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const 
689{
690  G4double mass = theActualMass;   /* the actual mass value */
691  G4double mass1 = theDaughterMass[0];
692  G4double mass2 = theDaughterMass[1];
693  G4double gamma2 = theDaughterWidth[1];
694 
695  G4double result = (1. / (2 * mass)) *
696    std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
697             ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
698    BrWig(gamma2, mass2, xmass);
699  return result;
700}
701
702G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
703{
704  G4double mass = theDefinition->GetPDGMass();   /* the pole mass value */
705  G4double mass1 = theDaughterMass[0];
706  G4double mass2 = theDaughterMass[1];
707  G4double gamma2 = theDaughterWidth[1];
708  G4double result = (1. / (2 * mass)) *
709    std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
710             ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
711    BrWig(gamma2, mass2, xmass);
712 return result;
713}
714
715G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
716{
717  const G4double mass =  G4KineticTrack_Gmass;   /* the actual mass value */
718//  const G4double mass1 = theDaughterMass[0];
719  const G4double mass2 = theDaughterMass[1];
720  const G4double gamma2 = theDaughterWidth[1];
721
722  const G4double result = (1. / (2 * mass)) *
723    std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
724         ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
725    BrWig(gamma2, mass2, xmass);
726  return result;
727}
728
729G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
730{
731  const G4double mass =  G4KineticTrack_Gmass;
732  const G4double mass1 = theDaughterMass[0];
733  const G4double gamma1 = theDaughterWidth[0];
734//  const G4double mass2 = theDaughterMass[1];
735 
736  G4KineticTrack_xmass1 = xmass;
737 
738  const G4double theLowerLimit = 0.0;
739  const G4double theUpperLimit = mass - xmass;
740  const G4int nIterations = 100;
741 
742  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
743  G4double result = BrWig(gamma1, mass1, xmass)*
744    integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
745
746  return result;
747}
748
749G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
750{
751  const G4double theUpperLimit = theActualMass - theDaughterMass[0];
752  const G4int nIterations = 100;
753 
754  if (theLowerLimit>=theUpperLimit) return 0.0;
755
756  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
757  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1, 
758                                                   theLowerLimit, theUpperLimit, nIterations);
759  return theIntegralOverMass2;
760}
761
762G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
763{
764  const G4double theUpperLimit = poleMass - theDaughterMass[0];
765  const G4int nIterations = 100;
766 
767  if (theLowerLimit>=theUpperLimit) return 0.0;
768
769  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
770  const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
771                                                    theLowerLimit, theUpperLimit, nIterations);
772  return theIntegralOverMass2;
773}
774
775
776G4double G4KineticTrack::IntegrateCMMomentum2() const
777{
778  const G4double theLowerLimit = 0.0;
779  const G4double theUpperLimit = theActualMass;
780  const G4int nIterations = 100;
781 
782  if (theLowerLimit>=theUpperLimit) return 0.0;
783 
784  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
785  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
786                                                   theLowerLimit, theUpperLimit, nIterations);
787  return theIntegralOverMass2;
788}
789
790
791
792
793
794
795
796
Note: See TracBrowser for help on using the repository browser.