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