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