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