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