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