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