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