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