EVGEN/AliGenPythia.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.11  1999/09/29 09:24:14  fca
  19 Introduction of the Copyright and cvs Log
  20 */
  21
  22 #include "AliGenerator.h"
  23 #include "AliGenPythia.h"
  24 #include "AliRun.h"
  25 #include "AliPythia.h"
  26 #include <TDirectory.h>
  27 #include <TFile.h>
  28 #include <TTree.h>
  29 #include <stdlib.h>
  30 #include <AliPythia.h>
  31 #include <TParticle.h>
  32 //#include <GParticle.h>
  33  ClassImp(AliGenPythia)
  34
  35 AliGenPythia::AliGenPythia()
  36                  :AliGenerator()
  37 {
  38 }
  39
  40 AliGenPythia::AliGenPythia(Int_t npart)
  41                  :AliGenerator(npart)
  42 {
  43 // default charm production at 5. 5 TeV
  44 // semimuonic decay
  45 // structure function GRVHO
  46 //
  47     fXsection  = 0.;
  48     fParentSelect.Set(5);
  49     fChildSelect.Set(5);
  50     for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
  51     SetProcess();
  52     SetStrucFunc();
  53     SetForceDecay();
  54     SetPtHard();
  55     SetEnergyCMS();
  56 }
  57
  58 AliGenPythia::~AliGenPythia()
  59 {
  60 }
  61
  62 void AliGenPythia::Init()
  63 {
  64     SetMC(new AliPythia());
  65     fPythia=(AliPythia*) fgMCEvGen;
  66 //
  67     fParentWeight=1./Float_t(fNpart);
  68 //
  69 //  Forward Paramters to the AliPythia object
  70     fPythia->DefineParticles();
  71     fPythia->SetCKIN(3,fPtHardMin);
  72     fPythia->SetCKIN(4,fPtHardMax);
  73     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
  74     fPythia->ForceDecay(fForceDecay);
  75     fPythia->LuList(0);
  76     fPythia->PyStat(2);
  77 //  Parent and Children Selection
  78     switch (fProcess)
  79     {
  80     case charm:
  81
  82         fParentSelect[0]=411;
  83         fParentSelect[1]=421;
  84         fParentSelect[2]=431;
  85         fParentSelect[3]=4122;
  86         break;
  87     case charm_unforced:
  88
  89         fParentSelect[0]=411;
  90         fParentSelect[1]=421;
  91         fParentSelect[2]=431;
  92         fParentSelect[3]=4122;
  93         break;
  94     case beauty:
  95         fParentSelect[0]=511;
  96         fParentSelect[1]=521;
  97         fParentSelect[2]=531;
  98         fParentSelect[3]=5122;
  99         break;
 100     case beauty_unforced:
 101         fParentSelect[0]=511;
 102         fParentSelect[1]=521;
 103         fParentSelect[2]=531;
 104         fParentSelect[3]=5122;
 105         break;
 106     case jpsi_chi:
 107     case jpsi:
 108         fParentSelect[0]=443;
 109         break;
 110     case mb:
 111         break;
 112     }
 113
 114     switch (fForceDecay)
 115     {
 116     case semielectronic:
 117     case dielectron:
 118     case b_jpsi_dielectron:
 119     case b_psip_dielectron:
 120         fChildSelect[0]=11;
 121         break;
 122     case semimuonic:
 123     case dimuon:
 124     case b_jpsi_dimuon:
 125     case b_psip_dimuon:
 126     case pitomu:
 127     case katomu:
 128         fChildSelect[0]=13;
 129         break;
 130     case all:
 131     case nodecay:
 132         break;
 133     }
 134 }
 135
 136 void AliGenPythia::Generate()
 137 {
 138
 139     Float_t polar[3] =   {0,0,0};
 140     Float_t origin[3]=   {0,0,0};
 141     Float_t origin_p[3]= {0,0,0};
 142     Float_t origin0[3]=  {0,0,0};
 143     Float_t p[3], p_p[4], random[6];
 144     static TClonesArray *particles;
 145 //  converts from mm/c to s
 146     const Float_t kconv=0.001/2.999792458e8;
 147
 148
 149 //
 150     Int_t nt=0;
 151     Int_t nt_p=0;
 152     Int_t jev=0;
 153     Int_t j, kf;
 154
 155     if(!particles) particles=new TClonesArray("TParticle",1000);
 156
 157     fTrials=0;
 158     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 159     if(fVertexSmear==perEvent) {
 160         gMC->Rndm(random,6);
 161         for (j=0;j<3;j++) {
 162             origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 163                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 164             fPythia->SetMSTP(151,0);
 165         }
 166     } else if (fVertexSmear==perTrack) {
 167         fPythia->SetMSTP(151,0);
 168         for (j=0;j<3;j++) {
 169             fPythia->SetPARP(151+j, fOsigma[j]*10.);
 170         }
 171     }
 172     while(1)
 173     {
 174         fPythia->PyEvnt();
 175         fTrials++;
 176         fPythia->ImportParticles(particles,"All");
 177         Int_t np = particles->GetEntriesFast();
 178         printf("\n **************************************************%d\n",np);
 179         Int_t nc=0;
 180         if (np == 0 ) continue;
 181         if (fProcess != mb) {
 182             for (Int_t i = 0; i<np; i++) {
 183                 TParticle *  iparticle = (TParticle *) particles->At(i);
 184                 kf = CheckPDGCode(iparticle->GetPdgCode());
 185                 fChildWeight=(fPythia->GetBraPart(kf))*fParentWeight;
 186 //
 187 // Parent
 188                 if (ParentSelected(TMath::Abs(kf))) {
 189                     if (KinematicSelection(iparticle)) {
 190                         if (nc==0) {
 191 //
 192 // Store information concerning the hard scattering process
 193 //
 194                             Float_t mass_p = fPythia->GetPARI(13);
 195                             Float_t   pt_p = fPythia->GetPARI(17);
 196                             Float_t    y_p = fPythia->GetPARI(37);
 197                             Float_t  xmt_p = sqrt(pt_p*pt_p+mass_p*mass_p);
 198                             Float_t     ty = Float_t(TMath::TanH(y_p));
 199                             p_p[0] = pt_p;
 200                             p_p[1] = 0;
 201                             p_p[2] = xmt_p*ty/sqrt(1.-ty*ty);
 202                             p_p[3] = mass_p;
 203                             gAlice->SetTrack(0,-1,-1,
 204                                              p_p,origin_p,polar,
 205                                              0,"Hard Scat.",nt_p,fParentWeight);
 206                             gAlice->KeepTrack(nt_p);
 207                         }
 208                         nc++;
 209 //
 210 // store parent track information
 211                         p[0]=iparticle->Px();
 212                         p[1]=iparticle->Py();
 213                         p[2]=iparticle->Pz();
 214                         origin[0]=origin0[0]+iparticle->Vx()/10;
 215                         origin[1]=origin0[1]+iparticle->Vy()/10;
 216                         origin[2]=origin0[2]+iparticle->Vz()/10;
 217
 218                         Int_t ifch=iparticle->GetFirstDaughter();
 219                         Int_t ilch=iparticle->GetLastDaughter();
 220                         if (ifch !=0 && ilch !=0) {
 221                             gAlice->SetTrack(0,nt_p,kf,
 222                                              p,origin,polar,
 223                                              0,"Primary",nt,fParentWeight);
 224                             gAlice->KeepTrack(nt);
 225                             Int_t iparent = nt;
 226 //
 227 // Children
 228
 229                             for (j=ifch; j<=ilch; j++)
 230                             {
 231                                 TParticle *  ichild =
 232                                     (TParticle *) particles->At(j-1);
 233                                 kf = CheckPDGCode(ichild->GetPdgCode());
 234 //
 235 //
 236                                 if (ChildSelected(TMath::Abs(kf))) {
 237                                     origin[0]=ichild->Vx();
 238                                     origin[1]=ichild->Vy();
 239                                     origin[2]=ichild->Vz();
 240                                     p[0]=ichild->Px();
 241                                     p[1]=ichild->Py();
 242                                     p[2]=ichild->Pz();
 243                                     Float_t tof=kconv*ichild->T();
 244                                     gAlice->SetTrack(fTrackIt, iparent, kf,
 245                                                      p,origin,polar,
 246                                                      tof,"Decay",nt,fChildWeight);
 247                                     gAlice->KeepTrack(nt);
 248                                 } // select child
 249                             } // child loop
 250                         }
 251                     } // kinematic selection
 252                 } // select particle
 253             } // particle loop
 254         } else {
 255             for (Int_t i = 0; i<np; i++) {
 256                 TParticle *  iparticle = (TParticle *) particles->At(i);
 257                 kf = CheckPDGCode(iparticle->GetPdgCode());
 258                 Int_t ks = iparticle->GetStatusCode();
 259                 if (ks==1 && kf!=0 && KinematicSelection(iparticle)) {
 260                         nc++;
 261 //
 262 // store track information
 263                         p[0]=iparticle->Px();
 264                         p[1]=iparticle->Py();
 265                         p[2]=iparticle->Pz();
 266                         origin[0]=origin0[0]+iparticle->Vx()/10;
 267                         origin[1]=origin0[1]+iparticle->Vy()/10;
 268                         origin[2]=origin0[2]+iparticle->Vz()/10;
 269                         Float_t tof=kconv*iparticle->T();
 270                         gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
 271                                          tof,"Primary",nt);
 272                         gAlice->KeepTrack(nt);
 273                 } // select particle
 274             } // particle loop
 275             printf("\n I've put %i particles on the stack \n",nc);
 276         } // mb ?
 277         if (nc > 0) {
 278             jev+=nc;
 279             if (jev >= fNpart) {
 280                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 281                 printf("\n Trials: %i\n",fTrials);
 282                 break;
 283             }
 284         }
 285     } // event loop
 286 //  adjust weight due to kinematic selection
 287     AdjustWeights();
 288 //  get cross-section
 289     fXsection=fPythia->GetPARI(1);
 290 }
 291
 292 Bool_t AliGenPythia::ParentSelected(Int_t ip)
 293 {
 294     for (Int_t i=0; i<5; i++)
 295     {
 296         if (fParentSelect[i]==ip) return kTRUE;
 297     }
 298     return kFALSE;
 299 }
 300
 301 Bool_t AliGenPythia::ChildSelected(Int_t ip)
 302 {
 303     for (Int_t i=0; i<5; i++)
 304     {
 305         if (fChildSelect[i]==ip) return kTRUE;
 306     }
 307     return kFALSE;
 308 }
 309
 310 Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
 311 {
 312     Float_t px=particle->Px();
 313     Float_t py=particle->Py();
 314     Float_t pz=particle->Pz();
 315     Float_t  e=particle->Energy();
 316
 317 //
 318 //  transverse momentum cut
 319     Float_t pt=TMath::Sqrt(px*px+py*py);
 320     if (pt > fPtMax || pt < fPtMin)
 321     {
 322 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
 323         return kFALSE;
 324     }
 325 //
 326 // momentum cut
 327     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
 328     if (p > fPMax || p < fPMin)
 329     {
 330 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 331         return kFALSE;
 332     }
 333
 334 //
 335 // theta cut
 336     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
 337     if (theta > fThetaMax || theta < fThetaMin)
 338     {
 339 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 340         return kFALSE;
 341     }
 342
 343 //
 344 // rapidity cut
 345     Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
 346     if (y > fYMax || y < fYMin)
 347     {
 348 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
 349         return kFALSE;
 350     }
 351
 352 //
 353 // phi cut
 354     Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)))+TMath::Pi();
 355     if (phi > fPhiMax || phi < fPhiMin)
 356     {
 357 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
 358         return kFALSE;
 359     }
 360
 361     return kTRUE;
 362 }
 363 void AliGenPythia::AdjustWeights()
 364 {
 365     TClonesArray *PartArray = gAlice->Particles();
 366     TParticle *Part;
 367     Int_t ntrack=gAlice->GetNtrack();
 368     for (Int_t i=0; i<ntrack; i++) {
 369         Part= (TParticle*) PartArray->UncheckedAt(i);
 370         Part->SetWeight(Part->GetWeight()*fKineBias);
 371     }
 372 }
 373
 374 Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
 375 {
 376 //
 377 //  If the particle is in a diffractive state, then take action accordigly
 378   switch (pdgcode) {
 379   case 110:
 380     //rho_diff0 -- difficult to translate, return rho0
 381     return 113;
 382   case 210:
 383     //pi_diffr+ -- change to pi+
 384     return 211;
 385   case 220:
 386     //omega_di0 -- change to omega0
 387     return 223;
 388   case 330:
 389     //phi_diff0 -- return phi0
 390     return 333;
 391   case 440:
 392     //J/psi_di0 -- return J/psi
 393     return 443;
 394   case 2110:
 395     //n_diffr -- return neutron
 396     return 2112;
 397   case 2210:
 398     //p_diffr+ -- return proton
 399     return 2212;
 400   }
 401   //non diffractive state -- return code unchanged
 402   return pdgcode;
 403 }
 404
 405
 406
 407
 408
 409
 410
 411
 412