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