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