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