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