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