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