EVGEN/AliGenParam.cxx

   1 #include "AliGenParam.h"
   2 #include "AliGenMUONlib.h"
   3 #include "AliRun.h"
   4 #include "TGeant3.h"
   5 #include "AliPythia.h"
   6 #include <TDirectory.h>
   7 #include <TFile.h>
   8 #include <TTree.h>
   9 #include <stdlib.h>
  10 #include <TMCParticle.h>
  11
  12 ClassImp(AliGenParam)
  13
  14 //------------------------------------------------------------
  15
  16   //Begin_Html
  17   /*
  18     <img src="gif/AliGenParam.gif">
  19   */
  20   //End_Html
  21
  22 //____________________________________________________________
  23 //____________________________________________________________
  24 AliGenParam::AliGenParam()
  25                  :AliGenerator()
  26 {
  27   fPtPara = 0;
  28   fYPara  = 0;
  29   fIpart  = 0;
  30   fAnalog = 1;
  31 }
  32
  33 //____________________________________________________________
  34 AliGenParam::AliGenParam(Int_t npart, Int_t Ipart)
  35 //                                 Double_t (*PtPara)(Double_t*, Double_t*),
  36 //                                 Double_t (*YPara) (Double_t* ,Double_t*))
  37     :AliGenerator(npart)
  38 {
  39   //
  40   //  fName="HMESONpara";
  41   //  fTitle="Heavy Mesons Parametrisation";
  42   fPtParaFunc = AliGenMUONlib::GetPt(Ipart);
  43   fYParaFunc  = AliGenMUONlib::GetY(Ipart);
  44   fIpParaFunc = AliGenMUONlib::GetIp(Ipart);
  45
  46   fPtPara = 0;
  47   fYPara  = 0;
  48   fIpart  = Ipart;
  49   fAnalog = 1;
  50   fChildSelect.Set(5);
  51   for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
  52   ForceDecay();
  53 }
  54
  55 //____________________________________________________________
  56 AliGenParam::~AliGenParam()
  57 {
  58     delete  fPtPara;
  59     delete  fYPara;
  60 }
  61
  62 //____________________________________________________________
  63 void AliGenParam::Init()
  64 {
  65     SetMC(new AliPythia());
  66     fPythia= (AliPythia*) fgMCEvGen;
  67
  68 //  End of the test !!!
  69   //Begin_Html
  70   /*
  71     <img src="gif/AliGenParam.gif">
  72   */
  73   //End_Html
  74
  75   fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
  76   fYPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
  77 //
  78 // dN/dy| y=0
  79   Double_t y1=0;
  80   Double_t y2=0;
  81
  82   fdNdy0=fYParaFunc(&y1,&y2);
  83 //
  84 // Integral over generation region
  85   Float_t IntYS  = fYPara ->Integral(fYMin, fYMax);
  86   Float_t IntPt0 = fPtPara->Integral(0,15);
  87   Float_t IntPtS = fPtPara->Integral(fPtMin,fPtMax);
  88   Float_t PhiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
  89   if (fAnalog) {
  90      fYWgt  = IntYS/fdNdy0;
  91      fPtWgt = IntPtS/IntPt0;
  92      fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
  93   } else {
  94       fYWgt = IntYS/fdNdy0;
  95       fPtWgt = (fPtMax-fPtMin)/IntPt0;
  96       fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
  97   }
  98 //
  99 // particle decay related initialization
 100   fPythia->DefineParticles();
 101 // semimuonic decays of charm and beauty
 102   fPythia->ForceDecay(fForceDecay);
 103 //
 104     switch (fForceDecay)
 105     {
 106     case semielectronic:
 107     case dielectron:
 108     case b_jpsi_dielectron:
 109     case b_psip_dielectron:
 110         fChildSelect[0]=11;
 111         break;
 112     case semimuonic:
 113     case dimuon:
 114     case b_jpsi_dimuon:
 115     case b_psip_dimuon:
 116         fChildSelect[0]=13;
 117         break;
 118     }
 119
 120 }
 121
 122 //____________________________________________________________
 123 void AliGenParam::Generate()
 124 {
 125 // Generate 'npart' of heavy mesons (J/Psi, upsilon or phi) in the
 126 // the desired theta, phi and momentum windows; Gaussian smearing
 127 // on the vertex is done if selected
 128
 129   AliMC* pMC = AliMC::GetMC();
 130
 131   Float_t polar[3]= {0,0,0};
 132   //
 133   Float_t origin[3], origin0[3];
 134   Float_t pt, pl, ptot;
 135   Float_t phi, theta;
 136   Float_t p[3];
 137   Float_t ty, xmt;
 138   Int_t i, nt, j;
 139   Float_t  wgtp, wgtch;
 140   Double_t dummy;
 141   //
 142   //
 143   Float_t random[6];
 144   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 145   if(fVertexSmear==perEvent) {
 146       pMC->Rndm(random,6);
 147       for (j=0;j<3;j++) {
 148           origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 149               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 150       }
 151   }
 152   for(i=0;i<fNpart;i++) {
 153       while(1) {
 154 //
 155 // particle type
 156           Int_t Ipart = fIpParaFunc();
 157           fChildWeight=(fPythia->GetBraPart(Ipart))*fParentWeight;
 158           Float_t am=fPythia->GetPMAS(fPythia->LuComp(Ipart),1);
 159           pMC->Rndm(random,2);
 160 //
 161 // phi
 162           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
 163 //
 164 // y
 165           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
 166 //
 167 // pT
 168           if (fAnalog) {
 169               pt=fPtPara->GetRandom();
 170               wgtp=fParentWeight;
 171               wgtch=fChildWeight;
 172           } else {
 173               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 174               Double_t ptd=pt;
 175               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 176               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 177           }
 178           xmt=sqrt(pt*pt+am*am);
 179           pl=xmt*ty/sqrt(1.-ty*ty);
 180           theta=TMath::ATan2(pt,pl);
 181           if(theta<fThetaMin || theta>fThetaMax) continue;
 182           ptot=TMath::Sqrt(pt*pt+pl*pl);
 183           if(ptot<fPMin || ptot>fPMax) continue;
 184           p[0]=pt*TMath::Cos(phi);
 185           p[1]=pt*TMath::Sin(phi);
 186           p[2]=pl;
 187           if(fVertexSmear==perTrack) {
 188               pMC->Rndm(random,6);
 189               for (j=0;j<3;j++) {
 190                   origin0[j]=
 191                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 192                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 193               }
 194           }
 195
 196           Int_t kg=fPythia->GetGeantCode(Ipart);
 197 //
 198 // parent
 199
 200           gAlice->
 201               SetTrack(0,-1,kg,p,origin,polar,0,"Primary",nt,wgtp);
 202           Int_t iparent=nt;
 203 //
 204 // use lujet to decay particle
 205
 206           Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 207           fPythia->DecayParticle(Ipart,energy,theta,phi);
 208 //
 209 // select muons
 210           TObjArray* particles = fPythia->GetPrimaries() ;
 211           Int_t np = particles->GetEntriesFast();
 212           for (Int_t i = 0; i<np; i++) {
 213               TMCParticle *  iparticle = (TMCParticle *) particles->At(i);
 214               Int_t kf = iparticle->GetKF();
 215 //
 216 // children
 217               if (ChildSelected(TMath::Abs(kf)))
 218               {
 219                   p[0]=iparticle->GetPx();
 220                   p[1]=iparticle->GetPy();
 221                   p[2]=iparticle->GetPz();
 222                   origin[0]=origin0[0]+iparticle->GetVx()/10;
 223                   origin[1]=origin0[1]+iparticle->GetVy()/10;
 224                   origin[2]=origin0[2]+iparticle->GetVz()/10;
 225                   gAlice->SetTrack(1,iparent,fPythia->GetGeantCode(kf),
 226                                    p,origin,polar,
 227                                    0,"Decay",nt,wgtch);
 228                   gAlice->KeepTrack(nt);
 229               } // select muon
 230           } // decay particle loop
 231           break;
 232       } // kinematic selection
 233   } // event loop
 234 }
 235
 236 Bool_t AliGenParam::ChildSelected(Int_t ip)
 237 {
 238     for (Int_t i=0; i<5; i++)
 239     {
 240         if (fChildSelect[i]==ip) return kTRUE;
 241     }
 242     return kFALSE;
 243 }
 244
 245 Bool_t AliGenParam::KinematicSelection(TMCParticle *particle)
 246 {
 247     Float_t px=particle->GetPx();
 248     Float_t py=particle->GetPy();
 249     Float_t pz=particle->GetPz();
 250 //
 251 // momentum cut
 252     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
 253     if (p > fPMax || p < fPMin)
 254     {
 255 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 256         return kFALSE;
 257     }
 258     Float_t pt=TMath::Sqrt(px*px+py*py);
 259
 260 //
 261 // theta cut
 262     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
 263     if (theta > fThetaMax || theta < fThetaMin)
 264     {
 265 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 266         return kFALSE;
 267     }
 268
 269     return kTRUE;
 270 }
 271
 272
 273
 274