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

Last change on this file since 819 was 819, checked in by garnier, 16 years ago

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