EVGEN/AliGenParam.cxx

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