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.45  2002/12/16 09:44:49  morsch
  19 Default for fRadiation is 3.
  20
  21 Revision 1.44  2002/10/14 14:55:35  hristov
  22 Merging the VirtualMC branch to the main development branch (HEAD)
  23
  24 Revision 1.42.4.1  2002/08/28 15:06:50  alibrary
  25 Updating to v3-09-01
  26
  27 Revision 1.43  2002/08/09 12:09:52  morsch
  28 Direct gamma trigger correctly included.
  29
  30 Revision 1.42  2002/03/12 11:07:08  morsch
  31 Add status code of particle to SetTrack call.
  32
  33 Revision 1.41  2002/03/05 11:25:33  morsch
  34 - New quenching options
  35 - Correction in CheckTrigger()
  36
  37 Revision 1.40  2002/02/12 11:05:53  morsch
  38 Get daughter indices right.
  39
  40 Revision 1.39  2002/02/12 09:16:39  morsch
  41 Correction in SelectFlavor()
  42
  43 Revision 1.38  2002/02/12 08:53:21  morsch
  44 SetNoGammas can be used to inhibit writing of gammas and pi0.
  45
  46 Revision 1.37  2002/02/08 16:50:50  morsch
  47 Add name and title in constructor.
  48
  49 Revision 1.36  2002/01/31 20:17:55  morsch
  50 Allow for triggered jets with simplified topology: Exact pT, back-to-back
  51
  52 Revision 1.35  2001/12/13 07:56:25  hristov
  53 Set pointers to zero in the default constructor
  54
  55 Revision 1.34  2001/12/11 16:55:42  morsch
  56 Correct initialization for jet phi-range.
  57
  58 Revision 1.33  2001/12/05 10:18:51  morsch
  59 Possibility of kinematic biasing of jet phi range. (J. Klay)
  60
  61 Revision 1.32  2001/11/28 13:51:11  morsch
  62 Introduce kinematic biasing (etamin, etamax) of jet trigger. Bookkeeping
  63 (number of trials) done in AliGenHijingEventHeader.
  64
  65 Revision 1.31  2001/11/06 12:30:34  morsch
  66 Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric collision systems.
  67
  68 Revision 1.30  2001/10/21 18:35:56  hristov
  69 Several pointers were set to zero in the default constructors to avoid memory management problems
  70
  71 Revision 1.29  2001/10/15 08:12:24  morsch
  72 - Vertex smearing with truncated gaussian.
  73 - Store triggered jet info before and after final state radiation into mc-heade
  74
  75 Revision 1.28  2001/10/08 11:55:25  morsch
  76 Store 4-momenta of trigegred jets in event header.
  77 Possibility to switch of initial and final state radiation.
  78
  79 Revision 1.27  2001/10/08 07:13:14  morsch
  80 Add setter for minimum transverse momentum of triggered jet.
  81
  82 Revision 1.26  2001/10/04 08:12:24  morsch
  83 Redefinition of stable condition.
  84
  85 Revision 1.25  2001/07/27 17:09:36  morsch
  86 Use local SetTrack, KeepTrack and SetHighWaterMark methods
  87 to delegate either to local stack or to stack owned by AliRun.
  88 (Piotr Skowronski, A.M.)
  89
  90 Revision 1.24  2001/07/20 09:34:56  morsch
  91 Count the number of spectator neutrons and protons and add information
  92 to the event header. (Chiara Oppedisano)
  93
  94 Revision 1.23  2001/07/13 17:30:22  morsch
  95 Derive from AliGenMC.
  96
  97 Revision 1.22  2001/06/11 13:09:23  morsch
  98 - Store cross-Section and number of binary collisions as a function of impact parameter
  99 - Pass AliGenHijingEventHeader to gAlice.
 100
 101 Revision 1.21  2001/02/28 17:35:24  morsch
 102 Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
 103
 104 Revision 1.20  2001/02/14 15:50:40  hristov
 105 The last particle in event marked using SetHighWaterMark
 106
 107 Revision 1.19  2000/12/21 16:24:06  morsch
 108 Coding convention clean-up
 109
 110 Revision 1.18  2000/12/06 17:46:30  morsch
 111 Avoid random numbers 1 and 0.
 112
 113 Revision 1.17  2000/12/04 11:22:03  morsch
 114 Init of sRandom as in 1.15
 115
 116 Revision 1.16  2000/12/02 11:41:39  morsch
 117 Use SetRandom() to initialize random number generator in constructor.
 118
 119 Revision 1.15  2000/11/30 20:29:02  morsch
 120 Initialise static variable sRandom in constructor: sRandom = fRandom;
 121
 122 Revision 1.14  2000/11/30 07:12:50  alibrary
 123 Introducing new Rndm and QA classes
 124
 125 Revision 1.13  2000/11/09 17:40:27  morsch
 126 Possibility to select/unselect spectator protons and neutrons.
 127 Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
 128
 129 Revision 1.12  2000/10/20 13:38:38  morsch
 130 Debug printouts commented.
 131
 132 Revision 1.11  2000/10/20 13:22:26  morsch
 133 - skip particle type 92 (string)
 134 - Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
 135   mesons.
 136
 137 Revision 1.10  2000/10/17 15:10:20  morsch
 138 Write first all the parent particles to the stack and then the final state particles.
 139
 140 Revision 1.9  2000/10/17 13:38:59  morsch
 141 Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
 142
 143 Revision 1.8  2000/10/17 12:46:31  morsch
 144 Protect EvaluateCrossSections() against division by zero.
 145
 146 Revision 1.7  2000/10/02 21:28:06  fca
 147 Removal of useless dependecies via forward declarations
 148
 149 Revision 1.6  2000/09/11 13:23:37  morsch
 150 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
 151 luget_hijing and luset_hijing.
 152
 153 Revision 1.5  2000/09/07 16:55:40  morsch
 154 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
 155
 156 Revision 1.4  2000/07/11 18:24:56  fca
 157 Coding convention corrections + few minor bug fixes
 158
 159 Revision 1.3  2000/06/30 12:08:36  morsch
 160 In member data: char* replaced by TString, Init takes care of resizing the strings to
 161 8 characters required by Hijing.
 162
 163 Revision 1.2  2000/06/15 14:15:05  morsch
 164 Add possibility for heavy flavor selection: charm and beauty.
 165
 166 Revision 1.1  2000/06/09 20:47:27  morsch
 167 AliGenerator interface class to HIJING using THijing (test version)
 168
 169 */
 170
 171
 172
 173 // Generator using HIJING as an external generator
 174 // The main HIJING options are accessable for the user through this interface.
 175 // Uses the THijing implementation of TGenerator.
 176 //
 177 // andreas.morsch@cern.ch
 178
 179 #include "AliGenHijing.h"
 180 #include "AliGenHijingEventHeader.h"
 181 #include "AliRun.h"
 182 #include "AliPDG.h"
 183
 184 #include <TArrayI.h>
 185 #include <TParticle.h>
 186 #include <THijing.h>
 187 #include <TGraph.h>
 188 #include <TLorentzVector.h>
 189
 190
 191  ClassImp(AliGenHijing)
 192
 193 AliGenHijing::AliGenHijing()
 194                  :AliGenMC()
 195 {
 196 // Constructor
 197     fParticles = 0;
 198     fHijing    = 0;
 199     fDsigmaDb  = 0;
 200     fDnDb      = 0;
 201 }
 202
 203 AliGenHijing::AliGenHijing(Int_t npart)
 204     :AliGenMC(npart)
 205 {
 206 // Default PbPb collisions at 5. 5 TeV
 207 //
 208     fName = "Hijing";
 209     fTitle= "Particle Generator using HIJING";
 210
 211     SetEnergyCMS();
 212     SetImpactParameterRange();
 213     SetTarget();
 214     SetProjectile();
 215     SetBoostLHC();
 216     SetJetEtaRange();
 217     SetJetPhiRange();
 218
 219     fKeep       =  0;
 220     fQuench     =  1;
 221     fShadowing  =  1;
 222     fTrigger    =  0;
 223     fDecaysOff  =  1;
 224     fEvaluate   =  0;
 225     fSelectAll  =  0;
 226     fFlavor     =  0;
 227     fSpectators =  1;
 228     fDsigmaDb   =  0;
 229     fDnDb       =  0;
 230     fPtMinJet   = -2.5;
 231     fRadiation  =  3;
 232     fEventVertex.Set(3);
 233 //
 234     SetSimpleJets();
 235     SetNoGammas();
 236 //
 237     fParticles = new TClonesArray("TParticle",10000);
 238 //
 239 // Set random number generator
 240     sRandom = fRandom;
 241     fHijing = 0;
 242
 243 }
 244
 245 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
 246 {
 247 // copy constructor
 248 }
 249
 250
 251 AliGenHijing::~AliGenHijing()
 252 {
 253 // Destructor
 254     if ( fDsigmaDb) delete  fDsigmaDb;
 255     if ( fDnDb)     delete  fDnDb;
 256     delete fParticles;
 257 }
 258
 259 void AliGenHijing::Init()
 260 {
 261 // Initialisation
 262     fFrame.Resize(8);
 263     fTarget.Resize(8);
 264     fProjectile.Resize(8);
 265
 266     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
 267                       fAProjectile, fZProjectile, fATarget, fZTarget,
 268                       fMinImpactParam, fMaxImpactParam));
 269
 270     fHijing=(THijing*) fgMCEvGen;
 271     fHijing->SetIHPR2(2,  fRadiation);
 272     fHijing->SetIHPR2(3,  fTrigger);
 273     fHijing->SetIHPR2(6,  fShadowing);
 274     fHijing->SetIHPR2(12, fDecaysOff);
 275     fHijing->SetIHPR2(21, fKeep);
 276     fHijing->SetHIPR1(10, fPtMinJet);
 277     fHijing->SetHIPR1(50, fSimpleJet);
 278 //
 279 //  Quenching
 280 //
 281 //
 282 //  fQuench = 0:  no quenching
 283 //  fQuench = 1:  hijing default
 284 //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
 285 //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
 286 //  fQuench = 4:  new LHC  parameters with log(e) dependence
 287 //  fQuench = 5:  new RHIC parameters with log(e) dependence
 288     fHijing->SetIHPR2(50, 0);
 289     if (fQuench > 0)
 290         fHijing->SetIHPR2(4,  1);
 291     else
 292         fHijing->SetIHPR2(4,  0);
 293 // New LHC parameters from Xin-Nian Wang
 294     if (fQuench == 2) {
 295         fHijing->SetHIPR1(14, 1.1);
 296         fHijing->SetHIPR1(11, 3.7);
 297     } else if (fQuench == 3) {
 298         fHijing->SetHIPR1(14, 0.20);
 299         fHijing->SetHIPR1(11, 2.5);
 300     } else if (fQuench == 4) {
 301         fHijing->SetIHPR2(50, 1);
 302         fHijing->SetHIPR1(14, 4.*0.34);
 303         fHijing->SetHIPR1(11, 3.7);
 304     } else if (fQuench == 5) {
 305         fHijing->SetIHPR2(50, 1);
 306         fHijing->SetHIPR1(14, 0.34);
 307         fHijing->SetHIPR1(11, 2.5);
 308     }
 309
 310
 311
 312 //
 313 //  Initialize Hijing
 314 //
 315     fHijing->Initialize();
 316 //
 317     if (fEvaluate) EvaluateCrossSections();
 318 //
 319 }
 320
 321 void AliGenHijing::Generate()
 322 {
 323 // Generate one event
 324
 325   Float_t polar[3]    =   {0,0,0};
 326   Float_t origin[3]   =   {0,0,0};
 327   Float_t origin0[3]  =   {0,0,0};
 328   Float_t p[3], random[6];
 329   Float_t tof;
 330
 331 //  converts from mm/c to s
 332   const Float_t kconv = 0.001/2.999792458e8;
 333 //
 334   Int_t nt  = 0;
 335   Int_t jev = 0;
 336   Int_t j, kf, ks, imo;
 337   kf = 0;
 338
 339
 340
 341   fTrials = 0;
 342   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
 343   if(fVertexSmear == kPerEvent) {
 344       Float_t dv[3];
 345       dv[2] = 1.e10;
 346       while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
 347           Rndm(random,6);
 348           for (j=0; j < 3; j++) {
 349               dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 350                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 351           }
 352       }
 353       for (j=0; j < 3; j++) origin0[j] += dv[j];
 354   } else if (fVertexSmear == kPerTrack) {
 355 //          fHijing->SetMSTP(151,0);
 356       for (j = 0; j < 3; j++) {
 357 //            fHijing->SetPARP(151+j, fOsigma[j]*10.);
 358       }
 359   }
 360   while(1)
 361   {
 362 //    Generate one event
 363 // --------------------------------------------------------------------------
 364       fSpecn    = 0;
 365       fSpecp    = 0;
 366 // --------------------------------------------------------------------------
 367       fHijing->GenerateEvent();
 368       fTrials++;
 369       fHijing->ImportParticles(fParticles,"All");
 370       if (fTrigger != kNoTrigger) {
 371           if (!CheckTrigger()) continue;
 372       }
 373       if (fLHC) Boost();
 374
 375
 376       Int_t np = fParticles->GetEntriesFast();
 377       printf("\n **************************************************%d\n",np);
 378       Int_t nc = 0;
 379       if (np == 0 ) continue;
 380       Int_t i;
 381       Int_t* newPos     = new Int_t[np];
 382       Int_t* pSelected  = new Int_t[np];
 383
 384       for (i = 0; i < np; i++) {
 385           newPos[i]    = i;
 386           pSelected[i] = 0;
 387       }
 388
 389 //      Get event vertex
 390 //
 391       TParticle *  iparticle = (TParticle *) fParticles->At(0);
 392       fEventVertex[0] = origin0[0];
 393       fEventVertex[1] = origin0[1];
 394       fEventVertex[2] = origin0[2];
 395
 396 //
 397 //      First select parent particles
 398 //
 399
 400       for (i = 0; i < np; i++) {
 401           iparticle = (TParticle *) fParticles->At(i);
 402
 403 // Is this a parent particle ?
 404           if (Stable(iparticle)) continue;
 405 //
 406           Bool_t  selected             =  kTRUE;
 407           Bool_t  hasSelectedDaughters =  kFALSE;
 408
 409
 410           kf        = iparticle->GetPdgCode();
 411           ks        = iparticle->GetStatusCode();
 412           if (kf == 92) continue;
 413
 414           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
 415                                SelectFlavor(kf);
 416           hasSelectedDaughters = DaughtersSelection(iparticle);
 417 //
 418 // Put particle on the stack if it is either selected or
 419 // it is the mother of at least one seleted particle
 420 //
 421           if (selected || hasSelectedDaughters) {
 422               nc++;
 423               pSelected[i] = 1;
 424           } // selected
 425       } // particle loop parents
 426 //
 427 // Now select the final state particles
 428 //
 429
 430       for (i = 0; i<np; i++) {
 431           TParticle *  iparticle = (TParticle *) fParticles->At(i);
 432 // Is this a final state particle ?
 433           if (!Stable(iparticle)) continue;
 434
 435           Bool_t  selected             =  kTRUE;
 436           kf        = iparticle->GetPdgCode();
 437           ks        = iparticle->GetStatusCode();
 438
 439 // --------------------------------------------------------------------------
 440 // Count spectator neutrons and protons
 441           if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
 442               if(kf == kNeutron) fSpecn += 1;
 443               if(kf == kProton)  fSpecp += 1;
 444           }
 445 // --------------------------------------------------------------------------
 446 //
 447           if (!fSelectAll) {
 448               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 449               if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
 450                                                         && ks != 11);
 451           }
 452 //
 453 // Put particle on the stack if selected
 454 //
 455           if (selected) {
 456               nc++;
 457               pSelected[i] = 1;
 458           } // selected
 459       } // particle loop final state
 460 //
 461 // Write particles to stack
 462 //
 463       for (i = 0; i<np; i++) {
 464           TParticle *  iparticle = (TParticle *) fParticles->At(i);
 465           Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
 466           Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
 467
 468           if (pSelected[i]) {
 469               kf   = iparticle->GetPdgCode();
 470               ks   = iparticle->GetStatusCode();
 471               p[0] = iparticle->Px();
 472               p[1] = iparticle->Py();
 473               p[2] = iparticle->Pz();
 474               origin[0] = origin0[0]+iparticle->Vx()/10;
 475               origin[1] = origin0[1]+iparticle->Vy()/10;
 476               origin[2] = origin0[2]+iparticle->Vz()/10;
 477               tof = kconv*iparticle->T();
 478               imo = -1;
 479               TParticle* mother = 0;
 480               if (hasMother) {
 481                   imo = iparticle->GetFirstMother();
 482                   mother = (TParticle *) fParticles->At(imo);
 483                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
 484               } // if has mother
 485               Bool_t tFlag = (fTrackIt && !hasDaughter);
 486               SetTrack(tFlag,imo,kf,p,origin,polar,
 487                        tof,kPNoProcess,nt, 1., ks);
 488               KeepTrack(nt);
 489               newPos[i] = nt;
 490           } // if selected
 491       } // particle loop
 492       delete[] newPos;
 493       delete[] pSelected;
 494
 495       printf("\n I've put %i particles on the stack \n",nc);
 496       if (nc > 0) {
 497           jev += nc;
 498           if (jev >= fNpart || fNpart == -1) {
 499               fKineBias = Float_t(fNpart)/Float_t(fTrials);
 500               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 501               break;
 502           }
 503       }
 504   } // event loop
 505   MakeHeader();
 506   SetHighWaterMark(nt);
 507 }
 508
 509 void AliGenHijing::KeepFullEvent()
 510 {
 511     fKeep=1;
 512 }
 513
 514 void AliGenHijing::EvaluateCrossSections()
 515 {
 516 //     Glauber Calculation of geometrical x-section
 517 //
 518     Float_t xTot       = 0.;          // barn
 519     Float_t xTotHard   = 0.;          // barn
 520     Float_t xPart      = 0.;          // barn
 521     Float_t xPartHard  = 0.;          // barn
 522     Float_t sigmaHard  = 0.1;         // mbarn
 523     Float_t bMin       = 0.;
 524     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
 525     const Float_t kdib = 0.2;
 526     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
 527
 528
 529     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
 530     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));
 531     Int_t i;
 532     Float_t oldvalue= 0.;
 533
 534     Float_t* b   = new Float_t[kMax];
 535     Float_t* si1 = new Float_t[kMax];
 536     Float_t* si2 = new Float_t[kMax];
 537
 538     for (i = 0; i < kMax; i++)
 539     {
 540         Float_t xb  = bMin+i*kdib;
 541         Float_t ov;
 542         ov=fHijing->Profile(xb);
 543         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
 544         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
 545         xTot+=gb;
 546         xTotHard += gbh;
 547         if (xb > fMinImpactParam && xb < fMaxImpactParam)
 548         {
 549             xPart += gb;
 550             xPartHard += gbh;
 551         }
 552
 553         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
 554         oldvalue = xTot;
 555         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
 556         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
 557         if (i>0) {
 558             si1[i] = gb/kdib;
 559             si2[i] = gbh/gb;
 560             b[i]  = xb;
 561         }
 562     }
 563
 564     printf("\n Total cross section (barn): %f \n",xTot);
 565     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
 566     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
 567     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
 568
 569 //  Store result as a graph
 570     b[0] = 0;
 571     si1[0] = 0;
 572     si2[0]=si2[1];
 573
 574     fDsigmaDb  = new TGraph(i, b, si1);
 575     fDnDb      = new TGraph(i, b, si2);
 576 }
 577
 578 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
 579 {
 580 //
 581 // Looks recursively if one of the daughters has been selected
 582 //
 583 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 584     Int_t imin = -1;
 585     Int_t imax = -1;
 586     Int_t i;
 587     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 588     Bool_t selected = kFALSE;
 589     if (hasDaughters) {
 590         imin = iparticle->GetFirstDaughter();
 591         imax = iparticle->GetLastDaughter();
 592         for (i = imin; i <= imax; i++){
 593             TParticle *  jparticle = (TParticle *) fParticles->At(i);
 594             Int_t ip = jparticle->GetPdgCode();
 595             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 596                 selected=kTRUE; break;
 597             }
 598             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 599         }
 600     } else {
 601         return kFALSE;
 602     }
 603     return selected;
 604 }
 605
 606
 607 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
 608 {
 609 // Select flavor of particle
 610 // 0: all
 611 // 4: charm and beauty
 612 // 5: beauty
 613     Bool_t res = 0;
 614
 615     if (fFlavor == 0) {
 616         res = kTRUE;
 617     } else {
 618         Int_t ifl = TMath::Abs(pid/100);
 619         if (ifl > 10) ifl/=10;
 620         res = (fFlavor == ifl);
 621     }
 622 //
 623 //  This part if gamma writing is inhibited
 624     if (fNoGammas)
 625         res = res && (pid != kGamma && pid != kPi0);
 626 //
 627     return res;
 628 }
 629
 630 Bool_t AliGenHijing::Stable(TParticle*  particle)
 631 {
 632 // Return true for a stable particle
 633 //
 634
 635     if (particle->GetFirstDaughter() < 0 )
 636     {
 637         return kTRUE;
 638     } else {
 639         return kFALSE;
 640     }
 641 }
 642
 643
 644 void AliGenHijing::Boost()
 645 {
 646 //
 647 // Boost cms into LHC lab frame
 648 //
 649     Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) /
 650                                       (Double_t(fZTarget)    * Double_t(fAProjectile)));
 651     Double_t beta  = TMath::TanH(dy);
 652     Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
 653     Double_t gb    = gamma * beta;
 654
 655     printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
 656
 657     Int_t i;
 658     Int_t np = fParticles->GetEntriesFast();
 659     for (i = 0; i < np; i++)
 660     {
 661         TParticle* iparticle = (TParticle*) fParticles->At(i);
 662
 663         Double_t e   = iparticle->Energy();
 664         Double_t px  = iparticle->Px();
 665         Double_t py  = iparticle->Py();
 666         Double_t pz  = iparticle->Pz();
 667
 668         Double_t eb  = gamma * e -      gb * pz;
 669         Double_t pzb =   -gb * e +   gamma * pz;
 670
 671         iparticle->SetMomentum(px, py, pzb, eb);
 672     }
 673 }
 674
 675
 676 void AliGenHijing::MakeHeader()
 677 {
 678 // Builds the event header, to be called after each event
 679     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
 680     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
 681     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
 682     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
 683     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
 684     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
 685     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
 686                                                        fHijing->GetN01(),
 687                                                        fHijing->GetN10(),
 688                                                        fHijing->GetN11());
 689     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
 690
 691 // 4-momentum vectors of the triggered jets.
 692 //
 693 // Before final state gluon radiation.
 694     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21),
 695                                               fHijing->GetHINT1(22),
 696                                               fHijing->GetHINT1(23),
 697                                               fHijing->GetHINT1(24));
 698
 699     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31),
 700                                               fHijing->GetHINT1(32),
 701                                               fHijing->GetHINT1(33),
 702                                               fHijing->GetHINT1(34));
 703 // After final state gluon radiation.
 704     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26),
 705                                               fHijing->GetHINT1(27),
 706                                               fHijing->GetHINT1(28),
 707                                               fHijing->GetHINT1(29));
 708
 709     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36),
 710                                               fHijing->GetHINT1(37),
 711                                               fHijing->GetHINT1(38),
 712                                               fHijing->GetHINT1(39));
 713     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
 714 // Bookkeeping for kinematic bias
 715     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
 716 // Event Vertex
 717     header->SetPrimaryVertex(fEventVertex);
 718     gAlice->SetGenEventHeader(header);
 719     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 720 }
 721
 722 Bool_t AliGenHijing::CheckTrigger()
 723 {
 724 // Check the kinematic trigger condition
 725 //
 726     Bool_t   triggered = kFALSE;
 727
 728     if (fTrigger == 1) {
 729 //
 730 //  jet-jet Trigger
 731
 732         TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26),
 733                                                   fHijing->GetHINT1(27),
 734                                                   fHijing->GetHINT1(28),
 735                                                   fHijing->GetHINT1(29));
 736
 737         TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36),
 738                                                   fHijing->GetHINT1(37),
 739                                                   fHijing->GetHINT1(38),
 740                                                   fHijing->GetHINT1(39));
 741         Double_t eta1      = jet1->Eta();
 742         Double_t eta2      = jet2->Eta();
 743         Double_t phi1      = jet1->Phi();
 744         Double_t phi2      = jet2->Phi();
 745 //    printf("\n Trigger: %f %f %f %f",
 746 //         fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
 747         if (
 748             (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
 749              phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
 750             ||
 751             (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
 752              phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
 753             )
 754             triggered = kTRUE;
 755     } else if (fTrigger == 2) {
 756 //  Gamma Jet
 757 //
 758         Int_t np = fParticles->GetEntriesFast();
 759         for (Int_t i = 0; i < np; i++) {
 760             TParticle* part = (TParticle*) fParticles->At(i);
 761             Int_t kf = part->GetPdgCode();
 762             Int_t ks = part->GetStatusCode();
 763             if (kf == 22 && ks == 40) {
 764                 Float_t phi = part->Phi();
 765                 Float_t eta = part->Eta();
 766                 if  (eta < fEtaMaxJet &&
 767                      eta > fEtaMinJet &&
 768                      phi < fPhiMaxJet &&
 769                      phi > fPhiMinJet) {
 770                     triggered = 1;
 771                     break;
 772                 } // check phi,eta within limits
 773             } // direct gamma ?
 774         } // particle loop
 775     } // fTrigger == 2
 776     return triggered;
 777 }
 778
 779
 780
 781
 782 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
 783 {
 784 // Assignment operator
 785     return *this;
 786 }
 787
 788 #ifndef WIN32
 789 # define rluget_hijing rluget_hijing_
 790 # define rluset_hijing rluset_hijing_
 791 # define rlu_hijing    rlu_hijing_
 792 # define type_of_call
 793 #else
 794 # define rluget_hijing RLUGET_HIJING
 795 # define rluset_hijing RLUSET_HIJING
 796 # define rlu_hijing    RLU_HIJING
 797 # define type_of_call _stdcall
 798 #endif
 799
 800
 801 extern "C" {
 802   void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
 803   {printf("Dummy version of rluget_hijing reached\n");}
 804
 805   void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
 806   {printf("Dummy version of rluset_hijing reached\n");}
 807
 808   Double_t type_of_call rlu_hijing(Int_t & /*idum*/)
 809   {
 810       Float_t r;
 811       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
 812       return r;
 813   }
 814 }