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