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