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