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