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 <TParticle.h>
  11
  12 ClassImp(AliGenParam)
  13
  14 //------------------------------------------------------------
  15
  16   //Begin_Html
  17   /*
  18     <img src="picts/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="picts/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
 139   Float_t polar[3]= {0,0,0};
 140   //
 141   Float_t origin[3], origin0[3];
 142   Float_t pt, pl, ptot;
 143   Float_t phi, theta;
 144   Float_t p[3];
 145   Float_t ty, xmt;
 146   Int_t i, nt, j;
 147   Float_t  wgtp, wgtch;
 148   Double_t dummy;
 149   //
 150   //
 151   Float_t random[6];
 152   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 153   if(fVertexSmear==perEvent) {
 154       gMC->Rndm(random,6);
 155       for (j=0;j<3;j++) {
 156           origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 157               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 158       }
 159   }
 160   for(i=0;i<fNpart;i++) {
 161       while(1) {
 162 //
 163 // particle type
 164           Int_t Ipart = fIpParaFunc();
 165           fChildWeight=(fPythia->GetBraPart(Ipart))*fParentWeight;
 166           Float_t am=fPythia->GetPMAS(fPythia->LuComp(Ipart),1);
 167           gMC->Rndm(random,2);
 168 //
 169 // phi
 170           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
 171 //
 172 // y
 173           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
 174 //
 175 // pT
 176           if (fAnalog) {
 177               pt=fPtPara->GetRandom();
 178               wgtp=fParentWeight;
 179               wgtch=fChildWeight;
 180           } else {
 181               pt=fPtMin+random[1]*(fPtMax-fPtMin);
 182               Double_t ptd=pt;
 183               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
 184               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
 185           }
 186           xmt=sqrt(pt*pt+am*am);
 187           pl=xmt*ty/sqrt(1.-ty*ty);
 188           theta=TMath::ATan2(pt,pl);
 189           if(theta<fThetaMin || theta>fThetaMax) continue;
 190           ptot=TMath::Sqrt(pt*pt+pl*pl);
 191           if(ptot<fPMin || ptot>fPMax) continue;
 192           p[0]=pt*TMath::Cos(phi);
 193           p[1]=pt*TMath::Sin(phi);
 194           p[2]=pl;
 195           if(fVertexSmear==perTrack) {
 196               gMC->Rndm(random,6);
 197               for (j=0;j<3;j++) {
 198                   origin0[j]=
 199                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 200                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 201               }
 202           }
 203
 204 //
 205 // parent
 206
 207           gAlice->
 208               SetTrack(0,-1,Ipart,p,origin,polar,0,"Primary",nt,wgtp);
 209           Int_t iparent=nt;
 210 //
 211 // use lujet to decay particle
 212
 213           Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 214           fPythia->DecayParticle(Ipart,energy,theta,phi);
 215 //        fPythia->LuList(1);
 216
 217 //
 218 // select muons
 219           TObjArray* particles = fPythia->ImportParticles() ;
 220           Int_t np = particles->GetEntriesFast();
 221           for (Int_t i = 0; i<np; i++) {
 222               TParticle *  iparticle = (TParticle *) particles->At(i);
 223               Int_t kf = iparticle->GetPdgCode();
 224 //
 225 // children
 226               if (ChildSelected(TMath::Abs(kf)))
 227               {
 228                   p[0]=iparticle->Px();
 229                   p[1]=iparticle->Py();
 230                   p[2]=iparticle->Pz();
 231                   origin[0]=origin0[0]+iparticle->Vx()/10;
 232                   origin[1]=origin0[1]+iparticle->Vy()/10;
 233                   origin[2]=origin0[2]+iparticle->Vz()/10;
 234                   gAlice->SetTrack(fTrackIt,iparent,kf,
 235                                    p,origin,polar,
 236                                    0,"Decay",nt,wgtch);
 237                   gAlice->KeepTrack(nt);
 238               } // select muon
 239           } // decay particle loop
 240           break;
 241       } // kinematic selection
 242   } // event loop
 243 }
 244
 245 Bool_t AliGenParam::ChildSelected(Int_t ip)
 246 {
 247     for (Int_t i=0; i<5; i++)
 248     {
 249         if (fChildSelect[i]==ip) return kTRUE;
 250     }
 251     return kFALSE;
 252 }
 253
 254 Bool_t AliGenParam::KinematicSelection(TParticle *particle)
 255 {
 256     Float_t px=particle->Px();
 257     Float_t py=particle->Py();
 258     Float_t pz=particle->Pz();
 259 //
 260 // momentum cut
 261     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
 262     if (p > fPMax || p < fPMin)
 263     {
 264 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 265         return kFALSE;
 266     }
 267     Float_t pt=TMath::Sqrt(px*px+py*py);
 268
 269 //
 270 // theta cut
 271     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
 272     if (theta > fThetaMax || theta < fThetaMin)
 273     {
 274 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 275         return kFALSE;
 276     }
 277
 278     return kTRUE;
 279 }
 280
 281
 282
 283