]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenHijing.cxx
Cleanup of STEER coding conventions
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.cxx
index 21e5974abb72b10c83f3b90bdae231471e45a1c2..156cd52b87cd18be4ed9c9616beebc685d43fe60 100644 (file)
 
 /*
 $Log$
+Revision 1.46  2003/01/07 14:12:33  morsch
+Provides collision geometry.
+
+Revision 1.45  2002/12/16 09:44:49  morsch
+Default for fRadiation is 3.
+
+Revision 1.44  2002/10/14 14:55:35  hristov
+Merging the VirtualMC branch to the main development branch (HEAD)
+
+Revision 1.42.4.1  2002/08/28 15:06:50  alibrary
+Updating to v3-09-01
+
+Revision 1.43  2002/08/09 12:09:52  morsch
+Direct gamma trigger correctly included.
+
+Revision 1.42  2002/03/12 11:07:08  morsch
+Add status code of particle to SetTrack call.
+
+Revision 1.41  2002/03/05 11:25:33  morsch
+- New quenching options
+- Correction in CheckTrigger()
+
+Revision 1.40  2002/02/12 11:05:53  morsch
+Get daughter indices right.
+
+Revision 1.39  2002/02/12 09:16:39  morsch
+Correction in SelectFlavor()
+
+Revision 1.38  2002/02/12 08:53:21  morsch
+SetNoGammas can be used to inhibit writing of gammas and pi0.
+
+Revision 1.37  2002/02/08 16:50:50  morsch
+Add name and title in constructor.
+
+Revision 1.36  2002/01/31 20:17:55  morsch
+Allow for triggered jets with simplified topology: Exact pT, back-to-back
+
+Revision 1.35  2001/12/13 07:56:25  hristov
+Set pointers to zero in the default constructor
+
+Revision 1.34  2001/12/11 16:55:42  morsch
+Correct initialization for jet phi-range.
+
+Revision 1.33  2001/12/05 10:18:51  morsch
+Possibility of kinematic biasing of jet phi range. (J. Klay)
+
+Revision 1.32  2001/11/28 13:51:11  morsch
+Introduce kinematic biasing (etamin, etamax) of jet trigger. Bookkeeping
+(number of trials) done in AliGenHijingEventHeader.
+
+Revision 1.31  2001/11/06 12:30:34  morsch
+Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric collision systems.
+
+Revision 1.30  2001/10/21 18:35:56  hristov
+Several pointers were set to zero in the default constructors to avoid memory management problems
+
+Revision 1.29  2001/10/15 08:12:24  morsch
+- Vertex smearing with truncated gaussian.
+- Store triggered jet info before and after final state radiation into mc-heade
+
+Revision 1.28  2001/10/08 11:55:25  morsch
+Store 4-momenta of trigegred jets in event header.
+Possibility to switch of initial and final state radiation.
+
+Revision 1.27  2001/10/08 07:13:14  morsch
+Add setter for minimum transverse momentum of triggered jet.
+
+Revision 1.26  2001/10/04 08:12:24  morsch
+Redefinition of stable condition.
+
+Revision 1.25  2001/07/27 17:09:36  morsch
+Use local SetTrack, KeepTrack and SetHighWaterMark methods
+to delegate either to local stack or to stack owned by AliRun.
+(Piotr Skowronski, A.M.)
+
+Revision 1.24  2001/07/20 09:34:56  morsch
+Count the number of spectator neutrons and protons and add information
+to the event header. (Chiara Oppedisano)
+
+Revision 1.23  2001/07/13 17:30:22  morsch
+Derive from AliGenMC.
+
+Revision 1.22  2001/06/11 13:09:23  morsch
+- Store cross-Section and number of binary collisions as a function of impact parameter
+- Pass AliGenHijingEventHeader to gAlice.
+
 Revision 1.21  2001/02/28 17:35:24  morsch
 Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
 
@@ -93,47 +179,70 @@ AliGenerator interface class to HIJING using THijing (test version)
 //
 // andreas.morsch@cern.ch
 
+#include <TArrayI.h>
+#include <TGraph.h>
+#include <THijing.h>
+#include <TLorentzVector.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+
 #include "AliGenHijing.h"
 #include "AliGenHijingEventHeader.h"
 #include "AliRun.h"
 
-#include <TArrayI.h>
-#include <TParticle.h>
-#include <THijing.h>
-#include <TGraph.h>
-
 
- ClassImp(AliGenHijing)
 ClassImp(AliGenHijing)
 
 AliGenHijing::AliGenHijing()
-                 :AliGenerator()
+                 :AliGenMC()
 {
 // Constructor
+    fParticles = 0;
+    fHijing    = 0;
+    fDsigmaDb  = 0;
+    fDnDb      = 0;
 }
 
 AliGenHijing::AliGenHijing(Int_t npart)
-    :AliGenerator(npart)
+    :AliGenMC(npart)
 {
 // Default PbPb collisions at 5. 5 TeV
 //
+    fName = "Hijing";
+    fTitle= "Particle Generator using HIJING";
+
     SetEnergyCMS();
     SetImpactParameterRange();
     SetTarget();
     SetProjectile();
-    fKeep       = 0;
-    fQuench     = 1;
-    fShadowing  = 1;
-    fTrigger    = 0;
-    fDecaysOff  = 1;
-    fEvaluate   = 0;
-    fSelectAll  = 0;
-    fFlavor     = 0;
-    fSpectators = 1;
-    fDsigmaDb   = 0;  
-    fDnDb       = 0; 
+    SetBoostLHC();
+    SetJetEtaRange();
+    SetJetPhiRange();
+    
+    fKeep       =  0;
+    fQuench     =  1;
+    fShadowing  =  1;
+    fTrigger    =  0;
+    fDecaysOff  =  1;
+    fEvaluate   =  0;
+    fSelectAll  =  0;
+    fFlavor     =  0;
+    fSpectators =  1;
+    fDsigmaDb   =  0;
+    fDnDb       =  0;
+    fPtMinJet   = -2.5;        
+    fRadiation  =  3;
+    fEventVertex.Set(3);
+//
+    SetSimpleJets();
+    SetNoGammas();
+//
+    fParticles = new TClonesArray("TParticle",10000);    
 //
 // Set random number generator   
     sRandom = fRandom;
+    fHijing = 0;
+
 }
 
 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
@@ -147,6 +256,7 @@ AliGenHijing::~AliGenHijing()
 // Destructor
     if ( fDsigmaDb) delete  fDsigmaDb;  
     if ( fDnDb)     delete  fDnDb;  
+    delete fParticles;
 }
 
 void AliGenHijing::Init()
@@ -161,20 +271,54 @@ void AliGenHijing::Init()
                      fMinImpactParam, fMaxImpactParam));
 
     fHijing=(THijing*) fgMCEvGen;
-
+    fHijing->SetIHPR2(2,  fRadiation);
     fHijing->SetIHPR2(3,  fTrigger);
-    fHijing->SetIHPR2(4,  fQuench);
     fHijing->SetIHPR2(6,  fShadowing);
     fHijing->SetIHPR2(12, fDecaysOff);    
     fHijing->SetIHPR2(21, fKeep);
-    fHijing->Initialize();
-
+    fHijing->SetHIPR1(10, fPtMinJet);  
+    fHijing->SetHIPR1(50, fSimpleJet);
+//
+//  Quenching
+//
+//
+//  fQuench = 0:  no quenching
+//  fQuench = 1:  hijing default
+//  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
+//  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
+//  fQuench = 4:  new LHC  parameters with log(e) dependence
+//  fQuench = 5:  new RHIC parameters with log(e) dependence
+    fHijing->SetIHPR2(50, 0);
+    if (fQuench > 0) 
+       fHijing->SetIHPR2(4,  1);
+    else
+       fHijing->SetIHPR2(4,  0);
+// New LHC parameters from Xin-Nian Wang
+    if (fQuench == 2) {
+       fHijing->SetHIPR1(14, 1.1);
+       fHijing->SetHIPR1(11, 3.7);
+    } else if (fQuench == 3) {
+       fHijing->SetHIPR1(14, 0.20);
+       fHijing->SetHIPR1(11, 2.5);
+    } else if (fQuench == 4) {
+       fHijing->SetIHPR2(50, 1);
+       fHijing->SetHIPR1(14, 4.*0.34);
+       fHijing->SetHIPR1(11, 3.7);
+    } else if (fQuench == 5) {
+       fHijing->SetIHPR2(50, 1);
+       fHijing->SetHIPR1(14, 0.34);
+       fHijing->SetHIPR1(11, 2.5);
+    }
+    
+    
     
 //
-    if (fEvaluate) EvaluateCrossSections();
+//  Initialize Hijing  
+//    
+    fHijing->Initialize();
 //
+    if (fEvaluate) EvaluateCrossSections();
 //
-//  Initialize random generator
 }
 
 void AliGenHijing::Generate()
@@ -187,7 +331,6 @@ void AliGenHijing::Generate()
   Float_t p[3], random[6];
   Float_t tof;
 
-  static TClonesArray *particles;
 //  converts from mm/c to s
   const Float_t kconv = 0.001/2.999792458e8;
 //
@@ -196,16 +339,21 @@ void AliGenHijing::Generate()
   Int_t j, kf, ks, imo;
   kf = 0;
     
-  if(!particles) particles = new TClonesArray("TParticle",10000);
+
     
   fTrials = 0;
   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
   if(fVertexSmear == kPerEvent) {
-       Rndm(random,6);
-       for (j=0; j < 3; j++) {
-           origin0[j] += fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
-               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-       }
+      Float_t dv[3];
+      dv[2] = 1.e10;
+      while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
+         Rndm(random,6);
+         for (j=0; j < 3; j++) {
+             dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
+                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+         }
+      }
+      for (j=0; j < 3; j++) origin0[j] += dv[j];
   } else if (fVertexSmear == kPerTrack) {
 //         fHijing->SetMSTP(151,0);
       for (j = 0; j < 3; j++) {
@@ -215,113 +363,138 @@ void AliGenHijing::Generate()
   while(1)
   {
 //    Generate one event
+// --------------------------------------------------------------------------
+      fSpecn   = 0;  
+      fSpecp   = 0;
+// --------------------------------------------------------------------------
       fHijing->GenerateEvent();
       fTrials++;
-      fHijing->ImportParticles(particles,"All");
-      Int_t np = particles->GetEntriesFast();
+      fHijing->ImportParticles(fParticles,"All");
+      if (fTrigger != kNoTrigger) {
+         if (!CheckTrigger()) continue;
+      }
+      if (fLHC) Boost();
+      
+      
+      Int_t np = fParticles->GetEntriesFast();
       printf("\n **************************************************%d\n",np);
       Int_t nc = 0;
       if (np == 0 ) continue;
       Int_t i;
-      Int_t * newPos = new Int_t[np];
+      Int_t* newPos     = new Int_t[np];
+      Int_t* pSelected  = new Int_t[np];
 
-      for (i = 0; i < np; i++) *(newPos+i) = i;
+      for (i = 0; i < np; i++) {
+         newPos[i]    = i;
+         pSelected[i] = 0;
+      }
+      
+//      Get event vertex
 //
-//      First write parent particles
+      TParticle *  iparticle = (TParticle *) fParticles->At(0);
+      fEventVertex[0] = origin0[0];
+      fEventVertex[1] = origin0[1];    
+      fEventVertex[2] = origin0[2];
+      
+//
+//      First select parent particles
 //
 
       for (i = 0; i < np; i++) {
-         TParticle *  iparticle       = (TParticle *) particles->At(i);
+         iparticle = (TParticle *) fParticles->At(i);
+
 // Is this a parent particle ?
-        if (Stable(iparticle)) continue;
+         if (Stable(iparticle)) continue;
 //
-        Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
-        Bool_t  selected             =  kTRUE;
-        Bool_t  hasSelectedDaughters =  kFALSE;
-           
-           
-        kf        = iparticle->GetPdgCode();
-        ks        = iparticle->GetStatusCode();
-        if (kf == 92) continue;
+         Bool_t  selected             =  kTRUE;
+         Bool_t  hasSelectedDaughters =  kFALSE;
+         
+         
+         kf        = iparticle->GetPdgCode();
+         ks        = iparticle->GetStatusCode();
+         if (kf == 92) continue;
            
-        if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
-        hasSelectedDaughters = DaughtersSelection(iparticle, particles);
-//
-// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
-//
-        if (selected || hasSelectedDaughters) {
-           nc++;
-           p[0] = iparticle->Px();
-           p[1] = iparticle->Py();
-           p[2] = iparticle->Pz();
-           origin[0] = origin0[0]+iparticle->Vx()/10;
-           origin[1] = origin0[1]+iparticle->Vy()/10;
-           origin[2] = origin0[2]+iparticle->Vz()/10;
-           tof = kconv*iparticle->T();
-           imo = -1;
-           if (hasMother) {
-               imo = iparticle->GetFirstMother();
-               TParticle* mother = (TParticle *) particles->At(imo);
-               imo = (mother->GetPdgCode() != 92) ? imo =* (newPos+imo) : -1;
-           }
-// Put particle on the stack ... 
-//             printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
-           
-           gAlice->SetTrack(0,imo,kf,p,origin,polar,
-                            tof,kPPrimary,nt);
-// ... and keep it there
-           gAlice->KeepTrack(nt);
+         if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
+                              SelectFlavor(kf);
+         hasSelectedDaughters = DaughtersSelection(iparticle);
 //
-           *(newPos+i)=nt;
-        } // selected
+// Put particle on the stack if it is either selected or 
+// it is the mother of at least one seleted particle
+//
+         if (selected || hasSelectedDaughters) {
+             nc++;
+             pSelected[i] = 1;
+         } // selected
       } // particle loop parents
 //
-// Now write the final state particles
+// Now select the final state particles
 //
 
       for (i = 0; i<np; i++) {
-        TParticle *  iparticle       = (TParticle *) particles->At(i);
+         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 // Is this a final state particle ?
-        if (!Stable(iparticle)) continue;
+         if (!Stable(iparticle)) continue;
+      
+         Bool_t  selected             =  kTRUE;
+         kf        = iparticle->GetPdgCode();
+         ks        = iparticle->GetStatusCode();
+         
+// --------------------------------------------------------------------------
+// Count spectator neutrons and protons
+         if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
+             if(kf == kNeutron) fSpecn += 1;
+             if(kf == kProton)  fSpecp += 1;
+         }
+// --------------------------------------------------------------------------
 //         
-        Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
-        Bool_t  selected             =  kTRUE;
-        kf        = iparticle->GetPdgCode();
-        ks        = iparticle->GetStatusCode();
-        if (!fSelectAll) {
-           selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
-           if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
-                                                     && ks != 11);
-        }
+         if (!fSelectAll) {
+             selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
+             if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
+                                                       && ks != 11);
+         }
 //
 // Put particle on the stack if selected
 //
-        if (selected) {
-           nc++;
-           p[0] = iparticle->Px();
-           p[1] = iparticle->Py();
-           p[2] = iparticle->Pz();
-           origin[0] = origin0[0]+iparticle->Vx()/10;
-           origin[1] = origin0[1]+iparticle->Vy()/10;
-           origin[2] = origin0[2]+iparticle->Vz()/10;
-           tof = kconv*iparticle->T();
-           imo = -1;
-           
-           if (hasMother) {
-               imo = iparticle->GetFirstMother();
-               TParticle* mother = (TParticle *) particles->At(imo);
-               imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
-           }   
-// Put particle on the stack
-           gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
-                            tof,kPNoProcess,nt);
-           gAlice->KeepTrack(nt);
-           *(newPos+i)=nt;
-        } // selected
+         if (selected) {
+             nc++;
+             pSelected[i] = 1;
+         } // selected
       } // particle loop final state
-      delete newPos;
-
+//
+// Write particles to stack
+//
+      for (i = 0; i<np; i++) {
+         TParticle *  iparticle = (TParticle *) fParticles->At(i);
+         Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
+         Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
+
+         if (pSelected[i]) {
+             kf   = iparticle->GetPdgCode();
+             ks   = iparticle->GetStatusCode();
+             p[0] = iparticle->Px();
+             p[1] = iparticle->Py();
+             p[2] = iparticle->Pz();
+             origin[0] = origin0[0]+iparticle->Vx()/10;
+             origin[1] = origin0[1]+iparticle->Vy()/10;
+             origin[2] = origin0[2]+iparticle->Vz()/10;
+             tof = kconv*iparticle->T();
+             imo = -1;
+             TParticle* mother = 0;
+             if (hasMother) {
+                 imo = iparticle->GetFirstMother();
+                 mother = (TParticle *) fParticles->At(imo);
+                 imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
+             } // if has mother   
+             Bool_t tFlag = (fTrackIt && !hasDaughter);
+             SetTrack(tFlag,imo,kf,p,origin,polar,
+                      tof,kPNoProcess,nt, 1., ks);
+             KeepTrack(nt);
+             newPos[i] = nt;
+         } // if selected
+      } // particle loop
+      delete[] newPos;
+      delete[] pSelected;
+      
       printf("\n I've put %i particles on the stack \n",nc);
       if (nc > 0) {
          jev += nc;
@@ -333,67 +506,7 @@ void AliGenHijing::Generate()
       }
   } // event loop
   MakeHeader();
-  
-  gAlice->SetHighWaterMark(nt);
-}
-
-Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
-{
-// Perform kinematic selection
-    Double_t px = particle->Px();
-    Double_t py = particle->Py();
-    Double_t pz = particle->Pz();
-    Double_t  e = particle->Energy();
-
-//
-//  transverse momentum cut    
-    Double_t pt = TMath::Sqrt(px*px+py*py);
-    if (pt > fPtMax || pt < fPtMin) 
-    {
-//     printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
-       return kFALSE;
-    }
-//
-// momentum cut
-    Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
-    if (p > fPMax || p < fPMin) 
-    {
-//     printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
-       return kFALSE;
-    }
-    
-//
-// theta cut
-    Double_t  theta = Double_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
-    if (theta > fThetaMax || theta < fThetaMin) 
-    {
-       
-//     printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
-       return kFALSE;
-    }
-
-//
-// rapidity cut
-    Double_t y;
-    if(e <= pz) y = 99;
-    else if (e <= -pz)  y = -99;
-    else y = 0.5*TMath::Log((e+pz)/(e-pz));
-    if (y > fYMax || y < fYMin)
-    {
-//     printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
-       return kFALSE;
-    }
-
-//
-// phi cut
-    Double_t phi=Double_t(TMath::ATan2(Double_t(py),Double_t(px)));
-    if (phi > fPhiMax || phi < fPhiMin)
-    {
-//     printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
-       return kFALSE;
-    }
-
-    return kTRUE;
+  SetHighWaterMark(nt);
 }
 
 void AliGenHijing::KeepFullEvent()
@@ -465,7 +578,7 @@ void AliGenHijing::EvaluateCrossSections()
     fDnDb      = new TGraph(i, b, si2);
 }
 
-Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
+Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
 {
 //
 // Looks recursively if one of the daughters has been selected
@@ -480,12 +593,12 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* part
        imin = iparticle->GetFirstDaughter();
        imax = iparticle->GetLastDaughter();       
        for (i = imin; i <= imax; i++){
-           TParticle *  jparticle = (TParticle *) particles->At(i);    
+           TParticle *  jparticle = (TParticle *) fParticles->At(i);   
            Int_t ip = jparticle->GetPdgCode();
-           if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {
+           if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
            }
-           if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
+           if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
        }
     } else {
        return kFALSE;
@@ -500,21 +613,29 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid)
 // 0: all
 // 4: charm and beauty
 // 5: beauty
-    if (fFlavor == 0) return kTRUE;
+    Bool_t res = 0;
     
-    Int_t ifl = TMath::Abs(pid/100);
-    if (ifl > 10) ifl/=10;
-    return (fFlavor == ifl);
+    if (fFlavor == 0) {
+       res = kTRUE;
+    } else {
+       Int_t ifl = TMath::Abs(pid/100);
+       if (ifl > 10) ifl/=10;
+       res = (fFlavor == ifl);
+    }
+//
+//  This part if gamma writing is inhibited
+    if (fNoGammas) 
+       res = res && (pid != kGamma && pid != kPi0);
+//
+    return res;
 }
 
 Bool_t AliGenHijing::Stable(TParticle*  particle)
 {
 // Return true for a stable particle
 //
-    Int_t kf = TMath::Abs(particle->GetPdgCode());
     
-    if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
-        
+    if (particle->GetFirstDaughter() < 0 )
     {
        return kTRUE;
     } else {
@@ -522,23 +643,145 @@ Bool_t AliGenHijing::Stable(TParticle*  particle)
     }
 }
 
+
+void AliGenHijing::Boost()
+{
+//
+// Boost cms into LHC lab frame
+//
+    Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
+                                     (Double_t(fZTarget)    * Double_t(fAProjectile)));
+    Double_t beta  = TMath::TanH(dy);
+    Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
+    Double_t gb    = gamma * beta;
+
+    printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
+    
+    Int_t i;
+    Int_t np = fParticles->GetEntriesFast();
+    for (i = 0; i < np; i++) 
+    {
+       TParticle* iparticle = (TParticle*) fParticles->At(i);
+
+       Double_t e   = iparticle->Energy();
+       Double_t px  = iparticle->Px();
+       Double_t py  = iparticle->Py();
+       Double_t pz  = iparticle->Pz();
+
+       Double_t eb  = gamma * e -      gb * pz;
+       Double_t pzb =   -gb * e +   gamma * pz;
+
+       iparticle->SetMomentum(px, py, pzb, eb);
+    }
+}
+
+
 void AliGenHijing::MakeHeader()
 {
 // Builds the event header, to be called after each event
     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
-   ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
-   ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
-   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
-   ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
-   ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
-   ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
-                         fHijing->GetN01(),
-                         fHijing->GetN10(),
-                         fHijing->GetN11());
-   gAlice->SetGenEventHeader(header);
-   
+    ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
+    ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
+    ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
+    ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
+    ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
+    ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
+                                                      fHijing->GetN01(),
+                                                      fHijing->GetN10(),
+                                                      fHijing->GetN11());
+    ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
+
+// 4-momentum vectors of the triggered jets.
+//
+// Before final state gluon radiation.
+    TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
+                                             fHijing->GetHINT1(22),
+                                             fHijing->GetHINT1(23),
+                                             fHijing->GetHINT1(24));
+
+    TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
+                                             fHijing->GetHINT1(32),
+                                             fHijing->GetHINT1(33),
+                                             fHijing->GetHINT1(34));
+// After final state gluon radiation.
+    TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
+                                             fHijing->GetHINT1(27),
+                                             fHijing->GetHINT1(28),
+                                             fHijing->GetHINT1(29));
+
+    TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
+                                             fHijing->GetHINT1(37),
+                                             fHijing->GetHINT1(38),
+                                             fHijing->GetHINT1(39));
+    ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
+// Bookkeeping for kinematic bias
+    ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
+// Event Vertex
+    header->SetPrimaryVertex(fEventVertex);
+    gAlice->SetGenEventHeader(header);   
+    fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 }
 
+Bool_t AliGenHijing::CheckTrigger()
+{
+// Check the kinematic trigger condition
+//
+    Bool_t   triggered = kFALSE;
+    if (fTrigger == 1) {
+//
+//  jet-jet Trigger    
+       
+       TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), 
+                                                 fHijing->GetHINT1(27),
+                                                 fHijing->GetHINT1(28),
+                                                 fHijing->GetHINT1(29));
+       
+       TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), 
+                                                 fHijing->GetHINT1(37),
+                                                 fHijing->GetHINT1(38),
+                                                 fHijing->GetHINT1(39));
+       Double_t eta1      = jet1->Eta();
+       Double_t eta2      = jet2->Eta();
+       Double_t phi1      = jet1->Phi();
+       Double_t phi2      = jet2->Phi();
+//    printf("\n Trigger: %f %f %f %f",
+//        fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
+       if (
+           (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
+            phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
+           ||
+           (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
+            phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
+           ) 
+           triggered = kTRUE;
+    } else if (fTrigger == 2) {
+//  Gamma Jet
+//
+       Int_t np = fParticles->GetEntriesFast();
+       for (Int_t i = 0; i < np; i++) {
+           TParticle* part = (TParticle*) fParticles->At(i);
+           Int_t kf = part->GetPdgCode();
+           Int_t ks = part->GetStatusCode();
+           if (kf == 22 && ks == 40) {
+               Float_t phi = part->Phi();
+               Float_t eta = part->Eta();
+               if  (eta < fEtaMaxJet && 
+                    eta > fEtaMinJet &&
+                    phi < fPhiMaxJet && 
+                    phi > fPhiMinJet) {
+                   triggered = 1;
+                   break;
+               } // check phi,eta within limits
+           } // direct gamma ? 
+       } // particle loop
+    } // fTrigger == 2
+    return triggered;
+}
+
+
+
+
 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
 {
 // Assignment operator