EVGEN/AliGenHijing.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.1  2000/06/09 20:47:27  morsch
  19 AliGenerator interface class to HIJING using THijing (test version)
  20
  21 */
  22
  23 #include "AliGenHijing.h"
  24 #include "AliRun.h"
  25
  26 #include <TArrayI.h>
  27 #include <TParticle.h>
  28 #include <THijing.h>
  29
  30  ClassImp(AliGenHijing)
  31
  32 AliGenHijing::AliGenHijing()
  33                  :AliGenerator()
  34 {
  35 // Constructor
  36 }
  37
  38 AliGenHijing::AliGenHijing(Int_t npart)
  39                  :AliGenerator(npart)
  40 {
  41 // Default PbPb collisions at 5. 5 TeV
  42 //
  43     SetEnergyCMS();
  44     SetImpactParameterRange();
  45     SetTarget();
  46     SetProjectile();
  47     fKeep=0;
  48     fQuench=1;
  49     fShadowing=1;
  50     fTrigger=0;
  51     fDecaysOff=1;
  52     fEvaluate=0;
  53     fSelectAll=0;
  54     fFlavor=0;
  55 }
  56
  57 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
  58 {
  59 // copy constructor
  60 }
  61
  62
  63 AliGenHijing::~AliGenHijing()
  64 {
  65 // Destructor
  66 }
  67
  68 void AliGenHijing::Init()
  69 {
  70 // Initialisation
  71     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
  72                       fAProjectile, fZProjectile, fATarget, fZTarget,
  73                       fMinImpactParam, fMaxImpactParam));
  74
  75     fHijing=(THijing*) fgMCEvGen;
  76     fHijing->Initialize();
  77     fHijing->SetIHPR2(3,  fTrigger);
  78     fHijing->SetIHPR2(4,  fQuench);
  79     fHijing->SetIHPR2(6,  fShadowing);
  80     fHijing->SetIHPR2(12, fDecaysOff);
  81     fHijing->SetIHPR2(21, fKeep);
  82 //
  83     if (fEvaluate) EvaluateCrossSections();
  84 }
  85
  86 void AliGenHijing::Generate()
  87 {
  88 // Generate one event
  89
  90     Float_t polar[3] =   {0,0,0};
  91     Float_t origin[3]=   {0,0,0};
  92     Float_t origin0[3]=  {0,0,0};
  93     Float_t p[3], random[6];
  94     Float_t tof;
  95
  96     static TClonesArray *particles;
  97 //  converts from mm/c to s
  98     const Float_t kconv=0.001/2.999792458e8;
  99 //
 100     Int_t nt=0;
 101     Int_t jev=0;
 102     Int_t j, kf, ks, imo;
 103     kf=0;
 104
 105     if(!particles) particles=new TClonesArray("TParticle",10000);
 106
 107     fTrials=0;
 108     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 109     if(fVertexSmear==perEvent) {
 110         gMC->Rndm(random,6);
 111         for (j=0;j<3;j++) {
 112             origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 113                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 114 //          fHijing->SetMSTP(151,0);
 115         }
 116     } else if (fVertexSmear==perTrack) {
 117 //      fHijing->SetMSTP(151,0);
 118         for (j=0;j<3;j++) {
 119 //          fHijing->SetPARP(151+j, fOsigma[j]*10.);
 120         }
 121     }
 122     while(1)
 123     {
 124
 125         fHijing->GenerateEvent();
 126         fTrials++;
 127         fHijing->ImportParticles(particles,"All");
 128         Int_t np = particles->GetEntriesFast();
 129         printf("\n **************************************************%d\n",np);
 130         Int_t nc=0;
 131         if (np == 0 ) continue;
 132         Int_t i;
 133         Int_t * newPos = new Int_t[np];
 134         for (i = 0; i<np; i++) *(newPos+i)=i;
 135
 136         for (i = 0; i<np; i++) {
 137             TParticle *  iparticle       = (TParticle *) particles->At(i);
 138
 139             Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
 140             Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
 141             Bool_t  selected             =  kTRUE;
 142             Bool_t  hasSelectedDaughters =  kFALSE;
 143
 144
 145             kf        = iparticle->GetPdgCode();
 146             if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
 147             if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
 148 //
 149 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
 150 //
 151
 152             if (selected || hasSelectedDaughters) {
 153                 nc++;
 154                 ks        = iparticle->GetStatusCode();
 155                 p[0]=iparticle->Px();
 156                 p[1]=iparticle->Py();
 157                 p[2]=iparticle->Pz();
 158                 origin[0]=origin0[0]+iparticle->Vx()/10;
 159                 origin[1]=origin0[1]+iparticle->Vy()/10;
 160                 origin[2]=origin0[2]+iparticle->Vz()/10;
 161                 tof=kconv*iparticle->T();
 162                 imo=-1;
 163                 if (hasMother) {
 164                     imo=iparticle->GetFirstMother();
 165                     imo=*(newPos+imo);
 166                 }
 167
 168 //              printf("\n selected iparent %d %d %d \n",i, kf, imo);
 169                 if (hasDaughter) {
 170                     gAlice->SetTrack(0,imo,kf,p,origin,polar,
 171                                      tof,"Primary",nt);
 172                 } else {
 173                     gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
 174                                      tof,"Secondary",nt);
 175                 }
 176                 *(newPos+i)=nt;
 177             } // selected
 178         } // particle loop
 179         delete newPos;
 180         printf("\n I've put %i particles on the stack \n",nc);
 181         if (nc > 0) {
 182             jev+=nc;
 183             if (jev >= fNpart || fNpart == -1) {
 184                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 185                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 186                 break;
 187             }
 188         }
 189     } // event loop
 190 }
 191
 192 Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
 193 {
 194 // Perform kinematic selection
 195     Float_t px=particle->Px();
 196     Float_t py=particle->Py();
 197     Float_t pz=particle->Pz();
 198     Float_t  e=particle->Energy();
 199
 200 //
 201 //  transverse momentum cut
 202     Float_t pt=TMath::Sqrt(px*px+py*py);
 203     if (pt > fPtMax || pt < fPtMin)
 204     {
 205 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
 206         return kFALSE;
 207     }
 208 //
 209 // momentum cut
 210     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
 211     if (p > fPMax || p < fPMin)
 212     {
 213 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 214         return kFALSE;
 215     }
 216
 217 //
 218 // theta cut
 219     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
 220     if (theta > fThetaMax || theta < fThetaMin)
 221     {
 222
 223 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 224         return kFALSE;
 225     }
 226
 227 //
 228 // rapidity cut
 229     Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
 230     if (y > fYMax || y < fYMin)
 231     {
 232 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
 233         return kFALSE;
 234     }
 235
 236 //
 237 // phi cut
 238     Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
 239     if (phi > fPhiMax || phi < fPhiMin)
 240     {
 241 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
 242         return kFALSE;
 243     }
 244
 245     return kTRUE;
 246 }
 247
 248 void AliGenHijing::KeepFullEvent()
 249 {
 250     fKeep=1;
 251 }
 252
 253 void AliGenHijing::EvaluateCrossSections()
 254 {
 255 //     Glauber Calculation of geometrical x-section
 256 //
 257     Float_t xTot=0.;          // barn
 258     Float_t xTotHard=0.;      // barn
 259     Float_t xPart=0.;         // barn
 260     Float_t xPartHard=0.;     // barn
 261     Float_t sigmaHard=0.1;    // mbarn
 262     Float_t bMin=0.;
 263     Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
 264     const Float_t kdib=0.2;
 265     Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;
 266
 267
 268     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
 269     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));
 270     Int_t i;
 271     Float_t oldvalue=0.;
 272
 273     for (i=0; i<kMax; i++)
 274     {
 275         Float_t xb=bMin+i*kdib;
 276         Float_t ov;
 277         ov=fHijing->Profile(xb);
 278         Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
 279         Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
 280         xTot+=gb;
 281         xTotHard+=gbh;
 282         if (xb > fMinImpactParam && xb < fMaxImpactParam)
 283         {
 284             xPart+=gb;
 285             xPartHard+=gbh;
 286         }
 287
 288         if ((xTot-oldvalue)/oldvalue<0.0001) break;
 289         oldvalue=xTot;
 290         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
 291         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
 292     }
 293     printf("\n Total cross section (barn): %f \n",xTot);
 294     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
 295     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
 296     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
 297 }
 298
 299 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
 300 {
 301 //
 302 // Looks recursively if one of the daughters has been selected
 303 //
 304 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 305     Int_t imin=-1;
 306     Int_t imax=-1;
 307     Int_t i;
 308     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
 309     Bool_t selected=kFALSE;
 310     if (hasDaughters) {
 311         imin=iparticle->GetFirstDaughter();
 312         imax=iparticle->GetLastDaughter();
 313         for (i=imin; i<= imax; i++){
 314             TParticle *  jparticle       = (TParticle *) particles->At(i);
 315             Int_t ip=jparticle->GetPdgCode();
 316             if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
 317             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
 318         }
 319     } else {
 320         return kFALSE;
 321     }
 322
 323     return selected;
 324 }
 325
 326
 327 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
 328 {
 329 // Select flavor of particle
 330 // 0: all
 331 // 4: charm and beauty
 332 // 5: beauty
 333     if (fFlavor == 0) return kTRUE;
 334
 335     Int_t ifl=TMath::Abs(pid/100);
 336     if (ifl > 10) ifl/=10;
 337     return ((fFlavor==4 && (ifl==4 || ifl==5))  ||
 338             (fFlavor==5 &&  ifl==5));
 339
 340 }
 341
 342 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
 343 {
 344 // Assignment operator
 345     return *this;
 346 }