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