ROOT logo
// @(#)root/physics:$Id: TGenPhaseSpace.cxx 20882 2007-11-19 11:31:26Z rdm $
// Author: Rene Brun , Valerio Filippini  06/09/2000 

//_____________________________________________________________________________________
//
//  Utility class to generate n-body event,
//  with constant cross-section (default)
//  or with Fermi energy dependence (opt="Fermi").
//  The event is generated in the center-of-mass frame, 
//  but the decay products are finally boosted
//  using the betas of the original particle.
//
//  The code is based on the GENBOD function (W515 from CERNLIB)
//  using the Raubold and Lynch method
//      F. James, Monte Carlo Phase Space, CERN 68-15 (1968)
//
// see example of use in $ROOTSYS/tutorials/physics/PhaseSpace.C
//
// Note that Momentum, Energy units are Gev/C, GeV

#include "TGenPhaseSpace.h"
#include "TRandom.h"
#include "TMath.h"
#include <stdlib.h>

const Int_t kMAXP = 18;

ClassImp(TGenPhaseSpace)

//_____________________________________________________________________________________
Double_t TGenPhaseSpace::PDK(Double_t a, Double_t b, Double_t c) 
{
   //the PDK function
   Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c);
   x = TMath::Sqrt(x)/(2*a);
   return x;
}

//_____________________________________________________________________________________
Int_t DoubleMax(const void *a, const void *b) 
{
   //special max function
   Double_t aa = * ((Double_t *) a); 
   Double_t bb = * ((Double_t *) b); 
   if (aa > bb) return  1;
   if (aa < bb) return -1;
   return 0;

}

//__________________________________________________________________________________________________
TGenPhaseSpace::TGenPhaseSpace(const TGenPhaseSpace &gen) : TObject(gen)
{
   //copy constructor
   fNt      = gen.fNt;
   fWtMax   = gen.fWtMax;
   fTeCmTm  = gen.fTeCmTm;
   fBeta[0] = gen.fBeta[0];
   fBeta[1] = gen.fBeta[1];
   fBeta[2] = gen.fBeta[2];
   for (Int_t i=0;i<fNt;i++) {
      fMass[i]   = gen.fMass[i];
      fDecPro[i] = gen.fDecPro[i];
   }
}

//__________________________________________________________________________________________________
Double_t TGenPhaseSpace::Generate() 
{
   //  Generate a random final state.
   //  The function returns the weigth of the current event.
   //  The TLorentzVector of each decay product can be obtained using GetDecay(n).
   //
   // Note that Momentum, Energy units are Gev/C, GeV

   Double_t rno[kMAXP];
   rno[0] = 0;
   Int_t n;
   if (fNt>2) {
      for (n=1; n<fNt-1; n++)  rno[n]=gRandom->Rndm();   // fNt-2 random numbers
      qsort(rno+1 ,fNt-2 ,sizeof(Double_t) ,DoubleMax);  // sort them
   }
   rno[fNt-1] = 1;

   Double_t invMas[kMAXP], sum=0;
   for (n=0; n<fNt; n++) {
      sum      += fMass[n];
      invMas[n] = rno[n]*fTeCmTm + sum;
   }

   //
   //-----> compute the weight of the current event
   //
   Double_t wt=fWtMax;
   Double_t pd[kMAXP];
   for (n=0; n<fNt-1; n++) {
      pd[n] = PDK(invMas[n+1],invMas[n],fMass[n+1]);
      wt *= pd[n];
   }

   //
   //-----> complete specification of event (Raubold-Lynch method)
   //
   fDecPro[0].SetPxPyPzE(0, pd[0], 0 , TMath::Sqrt(pd[0]*pd[0]+fMass[0]*fMass[0]) );

   Int_t i=1;
   Int_t j;
   while (1) {
      fDecPro[i].SetPxPyPzE(0, -pd[i-1], 0 , TMath::Sqrt(pd[i-1]*pd[i-1]+fMass[i]*fMass[i]) );

      Double_t cZ   = 2*gRandom->Rndm() - 1;
      Double_t sZ   = TMath::Sqrt(1-cZ*cZ);
      Double_t angY = 2*TMath::Pi() * gRandom->Rndm();
      Double_t cY   = TMath::Cos(angY);
      Double_t sY   = TMath::Sin(angY);
      for (j=0; j<=i; j++) {
         TLorentzVector *v = fDecPro+j;
         Double_t x = v->Px();
         Double_t y = v->Py();
         v->SetPx( cZ*x - sZ*y );
         v->SetPy( sZ*x + cZ*y );   // rotation around Z
         x = v->Px();
         Double_t z = v->Pz();
         v->SetPx( cY*x - sY*z );
         v->SetPz( sY*x + cY*z );   // rotation around Y
      }

      if (i == (fNt-1)) break;

      Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]);
      for (j=0; j<=i; j++) fDecPro[j].Boost(0,beta,0);
      i++;
   }
 
   //
   //---> final boost of all particles  
   //
   for (n=0;n<fNt;n++) fDecPro[n].Boost(fBeta[0],fBeta[1],fBeta[2]);

   //
   //---> return the weigth of event
   //
   return wt;
}

//__________________________________________________________________________________
TLorentzVector *TGenPhaseSpace::GetDecay(Int_t n) 
{ 
   //return Lorentz vector corresponding to decay n
   if (n>fNt) return 0;
   return fDecPro+n;
}


//_____________________________________________________________________________________
Bool_t TGenPhaseSpace::SetDecay(TLorentzVector &P, Int_t nt, 
   Double_t *mass, Option_t *opt) 
{
   // input:
   // TLorentzVector &P:    decay particle (Momentum, Energy units are Gev/C, GeV)
   // Int_t nt:             number of decay products
   // Double_t *mass:       array of decay product masses
   // Option_t *opt:        default -> constant cross section
   //                       "Fermi" -> Fermi energy dependece
   // return value:
   // kTRUE:      the decay is permitted by kinematics
   // kFALSE:     the decay is forbidden by kinematics
   //

   Int_t n;
   fNt = nt;
   if (fNt<2 || fNt>18) return kFALSE;  // no more then 18 particle

   //
   //
   //
   fTeCmTm = P.Mag();           // total energy in C.M. minus the sum of the masses
   for (n=0;n<fNt;n++) {
      fMass[n]  = mass[n];
      fTeCmTm  -= mass[n];
   }

   if (fTeCmTm<=0) return kFALSE;    // not enough energy for this decay

   //
   //------> the max weigth depends on opt:
   //   opt == "Fermi"  --> fermi energy dependence for cross section
   //   else            --> constant cross section as function of TECM (default)
   //
   if (strcasecmp(opt,"fermi")==0) {  
      // ffq[] = pi * (2*pi)**(FNt-2) / (FNt-2)!
      Double_t ffq[] = {0 
                     ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
                     ,256.3704, 268.4705, 240.9780, 189.2637
                     ,132.1308,  83.0202,  47.4210,  24.8295
                     ,12.0006,   5.3858,   2.2560,   0.8859 };
      fWtMax = TMath::Power(fTeCmTm,fNt-2) * ffq[fNt-1] / P.Mag();

   } else {
      Double_t emmax = fTeCmTm + fMass[0];
      Double_t emmin = 0;
      Double_t wtmax = 1;
      for (n=1; n<fNt; n++) {
         emmin += fMass[n-1];
         emmax += fMass[n];
         wtmax *= PDK(emmax, emmin, fMass[n]);
      }
      fWtMax = 1/wtmax;
   }

   //
   //---->  save the betas of the decaying particle
   //
   if (P.Beta()) {
      Double_t w = P.Beta()/P.Rho();
      fBeta[0] = P(0)*w;
      fBeta[1] = P(1)*w;
      fBeta[2] = P(2)*w;
   }
   else fBeta[0]=fBeta[1]=fBeta[2]=0; 

   return kTRUE; 
}
 TGenPhaseSpace.cxx:1
 TGenPhaseSpace.cxx:2
 TGenPhaseSpace.cxx:3
 TGenPhaseSpace.cxx:4
 TGenPhaseSpace.cxx:5
 TGenPhaseSpace.cxx:6
 TGenPhaseSpace.cxx:7
 TGenPhaseSpace.cxx:8
 TGenPhaseSpace.cxx:9
 TGenPhaseSpace.cxx:10
 TGenPhaseSpace.cxx:11
 TGenPhaseSpace.cxx:12
 TGenPhaseSpace.cxx:13
 TGenPhaseSpace.cxx:14
 TGenPhaseSpace.cxx:15
 TGenPhaseSpace.cxx:16
 TGenPhaseSpace.cxx:17
 TGenPhaseSpace.cxx:18
 TGenPhaseSpace.cxx:19
 TGenPhaseSpace.cxx:20
 TGenPhaseSpace.cxx:21
 TGenPhaseSpace.cxx:22
 TGenPhaseSpace.cxx:23
 TGenPhaseSpace.cxx:24
 TGenPhaseSpace.cxx:25
 TGenPhaseSpace.cxx:26
 TGenPhaseSpace.cxx:27
 TGenPhaseSpace.cxx:28
 TGenPhaseSpace.cxx:29
 TGenPhaseSpace.cxx:30
 TGenPhaseSpace.cxx:31
 TGenPhaseSpace.cxx:32
 TGenPhaseSpace.cxx:33
 TGenPhaseSpace.cxx:34
 TGenPhaseSpace.cxx:35
 TGenPhaseSpace.cxx:36
 TGenPhaseSpace.cxx:37
 TGenPhaseSpace.cxx:38
 TGenPhaseSpace.cxx:39
 TGenPhaseSpace.cxx:40
 TGenPhaseSpace.cxx:41
 TGenPhaseSpace.cxx:42
 TGenPhaseSpace.cxx:43
 TGenPhaseSpace.cxx:44
 TGenPhaseSpace.cxx:45
 TGenPhaseSpace.cxx:46
 TGenPhaseSpace.cxx:47
 TGenPhaseSpace.cxx:48
 TGenPhaseSpace.cxx:49
 TGenPhaseSpace.cxx:50
 TGenPhaseSpace.cxx:51
 TGenPhaseSpace.cxx:52
 TGenPhaseSpace.cxx:53
 TGenPhaseSpace.cxx:54
 TGenPhaseSpace.cxx:55
 TGenPhaseSpace.cxx:56
 TGenPhaseSpace.cxx:57
 TGenPhaseSpace.cxx:58
 TGenPhaseSpace.cxx:59
 TGenPhaseSpace.cxx:60
 TGenPhaseSpace.cxx:61
 TGenPhaseSpace.cxx:62
 TGenPhaseSpace.cxx:63
 TGenPhaseSpace.cxx:64
 TGenPhaseSpace.cxx:65
 TGenPhaseSpace.cxx:66
 TGenPhaseSpace.cxx:67
 TGenPhaseSpace.cxx:68
 TGenPhaseSpace.cxx:69
 TGenPhaseSpace.cxx:70
 TGenPhaseSpace.cxx:71
 TGenPhaseSpace.cxx:72
 TGenPhaseSpace.cxx:73
 TGenPhaseSpace.cxx:74
 TGenPhaseSpace.cxx:75
 TGenPhaseSpace.cxx:76
 TGenPhaseSpace.cxx:77
 TGenPhaseSpace.cxx:78
 TGenPhaseSpace.cxx:79
 TGenPhaseSpace.cxx:80
 TGenPhaseSpace.cxx:81
 TGenPhaseSpace.cxx:82
 TGenPhaseSpace.cxx:83
 TGenPhaseSpace.cxx:84
 TGenPhaseSpace.cxx:85
 TGenPhaseSpace.cxx:86
 TGenPhaseSpace.cxx:87
 TGenPhaseSpace.cxx:88
 TGenPhaseSpace.cxx:89
 TGenPhaseSpace.cxx:90
 TGenPhaseSpace.cxx:91
 TGenPhaseSpace.cxx:92
 TGenPhaseSpace.cxx:93
 TGenPhaseSpace.cxx:94
 TGenPhaseSpace.cxx:95
 TGenPhaseSpace.cxx:96
 TGenPhaseSpace.cxx:97
 TGenPhaseSpace.cxx:98
 TGenPhaseSpace.cxx:99
 TGenPhaseSpace.cxx:100
 TGenPhaseSpace.cxx:101
 TGenPhaseSpace.cxx:102
 TGenPhaseSpace.cxx:103
 TGenPhaseSpace.cxx:104
 TGenPhaseSpace.cxx:105
 TGenPhaseSpace.cxx:106
 TGenPhaseSpace.cxx:107
 TGenPhaseSpace.cxx:108
 TGenPhaseSpace.cxx:109
 TGenPhaseSpace.cxx:110
 TGenPhaseSpace.cxx:111
 TGenPhaseSpace.cxx:112
 TGenPhaseSpace.cxx:113
 TGenPhaseSpace.cxx:114
 TGenPhaseSpace.cxx:115
 TGenPhaseSpace.cxx:116
 TGenPhaseSpace.cxx:117
 TGenPhaseSpace.cxx:118
 TGenPhaseSpace.cxx:119
 TGenPhaseSpace.cxx:120
 TGenPhaseSpace.cxx:121
 TGenPhaseSpace.cxx:122
 TGenPhaseSpace.cxx:123
 TGenPhaseSpace.cxx:124
 TGenPhaseSpace.cxx:125
 TGenPhaseSpace.cxx:126
 TGenPhaseSpace.cxx:127
 TGenPhaseSpace.cxx:128
 TGenPhaseSpace.cxx:129
 TGenPhaseSpace.cxx:130
 TGenPhaseSpace.cxx:131
 TGenPhaseSpace.cxx:132
 TGenPhaseSpace.cxx:133
 TGenPhaseSpace.cxx:134
 TGenPhaseSpace.cxx:135
 TGenPhaseSpace.cxx:136
 TGenPhaseSpace.cxx:137
 TGenPhaseSpace.cxx:138
 TGenPhaseSpace.cxx:139
 TGenPhaseSpace.cxx:140
 TGenPhaseSpace.cxx:141
 TGenPhaseSpace.cxx:142
 TGenPhaseSpace.cxx:143
 TGenPhaseSpace.cxx:144
 TGenPhaseSpace.cxx:145
 TGenPhaseSpace.cxx:146
 TGenPhaseSpace.cxx:147
 TGenPhaseSpace.cxx:148
 TGenPhaseSpace.cxx:149
 TGenPhaseSpace.cxx:150
 TGenPhaseSpace.cxx:151
 TGenPhaseSpace.cxx:152
 TGenPhaseSpace.cxx:153
 TGenPhaseSpace.cxx:154
 TGenPhaseSpace.cxx:155
 TGenPhaseSpace.cxx:156
 TGenPhaseSpace.cxx:157
 TGenPhaseSpace.cxx:158
 TGenPhaseSpace.cxx:159
 TGenPhaseSpace.cxx:160
 TGenPhaseSpace.cxx:161
 TGenPhaseSpace.cxx:162
 TGenPhaseSpace.cxx:163
 TGenPhaseSpace.cxx:164
 TGenPhaseSpace.cxx:165
 TGenPhaseSpace.cxx:166
 TGenPhaseSpace.cxx:167
 TGenPhaseSpace.cxx:168
 TGenPhaseSpace.cxx:169
 TGenPhaseSpace.cxx:170
 TGenPhaseSpace.cxx:171
 TGenPhaseSpace.cxx:172
 TGenPhaseSpace.cxx:173
 TGenPhaseSpace.cxx:174
 TGenPhaseSpace.cxx:175
 TGenPhaseSpace.cxx:176
 TGenPhaseSpace.cxx:177
 TGenPhaseSpace.cxx:178
 TGenPhaseSpace.cxx:179
 TGenPhaseSpace.cxx:180
 TGenPhaseSpace.cxx:181
 TGenPhaseSpace.cxx:182
 TGenPhaseSpace.cxx:183
 TGenPhaseSpace.cxx:184
 TGenPhaseSpace.cxx:185
 TGenPhaseSpace.cxx:186
 TGenPhaseSpace.cxx:187
 TGenPhaseSpace.cxx:188
 TGenPhaseSpace.cxx:189
 TGenPhaseSpace.cxx:190
 TGenPhaseSpace.cxx:191
 TGenPhaseSpace.cxx:192
 TGenPhaseSpace.cxx:193
 TGenPhaseSpace.cxx:194
 TGenPhaseSpace.cxx:195
 TGenPhaseSpace.cxx:196
 TGenPhaseSpace.cxx:197
 TGenPhaseSpace.cxx:198
 TGenPhaseSpace.cxx:199
 TGenPhaseSpace.cxx:200
 TGenPhaseSpace.cxx:201
 TGenPhaseSpace.cxx:202
 TGenPhaseSpace.cxx:203
 TGenPhaseSpace.cxx:204
 TGenPhaseSpace.cxx:205
 TGenPhaseSpace.cxx:206
 TGenPhaseSpace.cxx:207
 TGenPhaseSpace.cxx:208
 TGenPhaseSpace.cxx:209
 TGenPhaseSpace.cxx:210
 TGenPhaseSpace.cxx:211
 TGenPhaseSpace.cxx:212
 TGenPhaseSpace.cxx:213
 TGenPhaseSpace.cxx:214
 TGenPhaseSpace.cxx:215
 TGenPhaseSpace.cxx:216
 TGenPhaseSpace.cxx:217
 TGenPhaseSpace.cxx:218
 TGenPhaseSpace.cxx:219
 TGenPhaseSpace.cxx:220
 TGenPhaseSpace.cxx:221
 TGenPhaseSpace.cxx:222
 TGenPhaseSpace.cxx:223