]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenHijing.cxx
Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric...
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.cxx
index 75e0b45dee5e52d623201ce4d5a5f14737c64292..9bd097b9daf447a3d6a98261d599f92bc5e6dd68 100644 (file)
 
 /*
 $Log$
+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)
+
+Revision 1.20  2001/02/14 15:50:40  hristov
+The last particle in event marked using SetHighWaterMark
+
+Revision 1.19  2000/12/21 16:24:06  morsch
+Coding convention clean-up
+
+Revision 1.18  2000/12/06 17:46:30  morsch
+Avoid random numbers 1 and 0.
+
+Revision 1.17  2000/12/04 11:22:03  morsch
+Init of sRandom as in 1.15
+
+Revision 1.16  2000/12/02 11:41:39  morsch
+Use SetRandom() to initialize random number generator in constructor.
+
+Revision 1.15  2000/11/30 20:29:02  morsch
+Initialise static variable sRandom in constructor: sRandom = fRandom;
+
+Revision 1.14  2000/11/30 07:12:50  alibrary
+Introducing new Rndm and QA classes
+
+Revision 1.13  2000/11/09 17:40:27  morsch
+Possibility to select/unselect spectator protons and neutrons.
+Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
+
+Revision 1.12  2000/10/20 13:38:38  morsch
+Debug printouts commented.
+
+Revision 1.11  2000/10/20 13:22:26  morsch
+- skip particle type 92 (string)
+- Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
+  mesons.
+
+Revision 1.10  2000/10/17 15:10:20  morsch
+Write first all the parent particles to the stack and then the final state particles.
+
+Revision 1.9  2000/10/17 13:38:59  morsch
+Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
+
 Revision 1.8  2000/10/17 12:46:31  morsch
 Protect EvaluateCrossSections() against division by zero.
 
@@ -43,25 +118,36 @@ AliGenerator interface class to HIJING using THijing (test version)
 
 */
 
+
+
+// Generator using HIJING as an external generator
+// The main HIJING options are accessable for the user through this interface.
+// Uses the THijing implementation of TGenerator.
+//
+// andreas.morsch@cern.ch
+
 #include "AliGenHijing.h"
 #include "AliGenHijingEventHeader.h"
 #include "AliRun.h"
-#include "AliMC.h"
+#include "AliPDG.h"
 
 #include <TArrayI.h>
 #include <TParticle.h>
 #include <THijing.h>
+#include <TGraph.h>
+#include <TLorentzVector.h>
+
 
  ClassImp(AliGenHijing)
 
 AliGenHijing::AliGenHijing()
-                 :AliGenerator()
+                 :AliGenMC()
 {
 // Constructor
 }
 
 AliGenHijing::AliGenHijing(Int_t npart)
-                 :AliGenerator(npart)
+    :AliGenMC(npart)
 {
 // Default PbPb collisions at 5. 5 TeV
 //
@@ -69,14 +155,26 @@ AliGenHijing::AliGenHijing(Int_t npart)
     SetImpactParameterRange();
     SetTarget();
     SetProjectile();
-    fKeep=0;
-    fQuench=1;
-    fShadowing=1;
-    fTrigger=0;
-    fDecaysOff=1;
-    fEvaluate=0;
-    fSelectAll=0;
-    fFlavor=0;
+    SetBoostLHC();
+    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  =  1;
+    fEventVertex.Set(3);
+//
+// Set random number generator   
+    sRandom = fRandom;
+    fHijing = 0;
+
 }
 
 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
@@ -88,6 +186,8 @@ AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
 AliGenHijing::~AliGenHijing()
 {
 // Destructor
+    if ( fDsigmaDb) delete  fDsigmaDb;  
+    if ( fDnDb)     delete  fDnDb;  
 }
 
 void AliGenHijing::Init()
@@ -102,184 +202,208 @@ 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->Rluset(50,0);
+    fHijing->SetHIPR1(10, fPtMinJet);  
     fHijing->Initialize();
 
     
 //
     if (fEvaluate) EvaluateCrossSections();
+//
+//
+//  Initialize random generator
 }
 
 void AliGenHijing::Generate()
 {
 // Generate one event
 
-    Float_t polar[3] =   {0,0,0};
-    Float_t origin[3]=   {0,0,0};
-    Float_t origin0[3]=  {0,0,0};
-    Float_t p[3], random[6];
-    Float_t tof;
+  Float_t polar[3]    =   {0,0,0};
+  Float_t origin[3]   =   {0,0,0};
+  Float_t origin0[3]  =   {0,0,0};
+  Float_t p[3], random[6];
+  Float_t tof;
 
-    static TClonesArray *particles;
+  static TClonesArray *particles;
 //  converts from mm/c to s
-    const Float_t kconv=0.001/2.999792458e8;
+  const Float_t kconv = 0.001/2.999792458e8;
 //
-    Int_t nt=0;
-    Int_t jev=0;
-    Int_t j, kf, ks, imo;
-    kf=0;
+  Int_t nt  = 0;
+  Int_t jev = 0;
+  Int_t j, kf, ks, imo;
+  kf = 0;
     
-    if(!particles) particles=new TClonesArray("TParticle",10000);
+  if(!particles) particles = new TClonesArray("TParticle",10000);
     
-    fTrials=0;
-    for (j=0;j<3;j++) origin0[j]=fOrigin[j];
-    if(fVertexSmear==kPerEvent) {
-       gMC->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]));
+  fTrials = 0;
+  for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
+  if(fVertexSmear == kPerEvent) {
+      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);
-       }
-    } else if (fVertexSmear==kPerTrack) {
-//     fHijing->SetMSTP(151,0);
-       for (j=0;j<3;j++) {
-//         fHijing->SetPARP(151+j, fOsigma[j]*10.);
-       }
-    }
-    while(1)
-    {
-
-       fHijing->GenerateEvent();
-       fTrials++;
-       fHijing->ImportParticles(particles,"All");
-       Int_t np = particles->GetEntriesFast();
-       printf("\n **************************************************%d\n",np);
-       Int_t nc=0;
-       if (np == 0 ) continue;
-       Int_t i;
-       Int_t * newPos = new Int_t[np];
-       for (i = 0; i<np; i++) *(newPos+i)=i;
-       
-       for (i = 0; i<np; i++) {
-           TParticle *  iparticle       = (TParticle *) particles->At(i);
-
-           Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
-           Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
-           Bool_t  selected             =  kTRUE;
-           Bool_t  hasSelectedDaughters =  kFALSE;
-
+      for (j = 0; j < 3; j++) {
+//           fHijing->SetPARP(151+j, fOsigma[j]*10.);
+      }
+  }
+  while(1)
+  {
+//    Generate one event
+// --------------------------------------------------------------------------
+      fSpecn   = 0;  
+      fSpecp   = 0;
+// --------------------------------------------------------------------------
+      fHijing->GenerateEvent();
+      fTrials++;
+      fHijing->ImportParticles(particles,"All");
+      if (fLHC) Boost(particles);
+      
+      Int_t np = particles->GetEntriesFast();
+      printf("\n **************************************************%d\n",np);
+      Int_t nc = 0;
+      if (np == 0 ) continue;
+      Int_t i;
+      Int_t * newPos = new Int_t[np];
+
+      for (i = 0; i < np; i++) *(newPos+i) = i;
+//      Get event vertex
+//
+      TParticle *  iparticle = (TParticle *) particles->At(0);
+      fEventVertex[0] = origin0[0];
+      fEventVertex[1] = origin0[1];    
+      fEventVertex[2] = origin0[2];
+      
+//
+//      First write parent particles
+//
 
-           kf        = iparticle->GetPdgCode();
-           if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
-           if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
+      for (i = 0; i < np; i++) {
+         iparticle       = (TParticle *) particles->At(i);
+// Is this a parent particle ?
+         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;
+           
+        if (!fSelectAll) selected = KinematicSelection(iparticle, 0)&&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++;
-               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;
-               if (hasMother) {
-                   imo=iparticle->GetFirstMother();
-                   imo=*(newPos+imo);
-               }
+        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;
+           TParticle* mother = 0;
+           if (hasMother) {
+               imo = iparticle->GetFirstMother();
+               mother = (TParticle *) particles->At(imo);
+               imo = (mother->GetPdgCode() != 92) ? imo =* (newPos+imo) : -1;
                
-//             printf("\n selected iparent %d %d %d \n",i, kf, imo);
-               if (hasDaughter) {
-                   gAlice->SetTrack(0,imo,kf,p,origin,polar,
-                                    tof,"Primary",nt);
-               } else {
-                   gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
-                                    tof,"Secondary",nt);
-               }
-               *(newPos+i)=nt;
-           } // selected
-       } // particle loop 
-       delete newPos;
-       printf("\n I've put %i particles on the stack \n",nc);
-       if (nc > 0) {
-           jev+=nc;
-           if (jev >= fNpart || fNpart == -1) {
-               fKineBias=Float_t(fNpart)/Float_t(fTrials);
-               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
-               break;
            }
-       }
-    } // event loop
-    fHijing->Rluget(50,-1);
-}
-
-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();
-
+// Put particle on the stack ... 
+//             printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
+           
+           SetTrack(0,imo,kf,p,origin,polar, tof,kPPrimary,nt);
+// ... and keep it there
+           KeepTrack(nt);
 //
-//  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;
-    }
+           *(newPos+i)=nt;
+        } // selected
+      } // particle loop parents
 //
-// 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;
-    }
-    
+// Now write the final state particles
 //
-// 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;
-    }
 
+      for (i = 0; i<np; i++) {
+        TParticle *  iparticle       = (TParticle *) particles->At(i);
+// Is this a final state particle ?
+        if (!Stable(iparticle)) continue;
+
+        Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
+        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;
+       }
+// --------------------------------------------------------------------------
+//         
+        if (!fSelectAll) {
+           selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
+           if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
+                                                     && ks != 11);
+        }
 //
-// 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;
-    }
-
+// Put particle on the stack if selected
 //
-// 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;
+        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;
+           TParticle* mother = 0;
+           if (hasMother) {
+               imo = iparticle->GetFirstMother();
+               mother = (TParticle *) particles->At(imo);
+               imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
+           }   
+// Put particle on the stack
+           SetTrack(fTrackIt,imo,kf,p,origin,polar,
+                                                    tof,kPNoProcess,nt);
+           KeepTrack(nt);
+           *(newPos+i)=nt;
+        } // selected
+      } // particle loop final state
+      delete[] newPos;
+
+      printf("\n I've put %i particles on the stack \n",nc);
+      if (nc > 0) {
+         jev += nc;
+         if (jev >= fNpart || fNpart == -1) {
+             fKineBias = Float_t(fNpart)/Float_t(fTrials);
+             printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+             break;
+         }
+      }
+  } // event loop
+  MakeHeader();
+  
+  SetHighWaterMark(nt);
 }
 
 void AliGenHijing::KeepFullEvent()
@@ -291,46 +415,64 @@ void AliGenHijing::EvaluateCrossSections()
 {
 //     Glauber Calculation of geometrical x-section
 //
-    Float_t xTot=0.;          // barn
-    Float_t xTotHard=0.;      // barn 
-    Float_t xPart=0.;         // barn
-    Float_t xPartHard=0.;     // barn 
-    Float_t sigmaHard=0.1;    // mbarn
-    Float_t bMin=0.;
-    Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
-    const Float_t kdib=0.2;
-    Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;
+    Float_t xTot       = 0.;          // barn
+    Float_t xTotHard   = 0.;          // barn 
+    Float_t xPart      = 0.;          // barn
+    Float_t xPartHard  = 0.;          // barn 
+    Float_t sigmaHard  = 0.1;         // mbarn
+    Float_t bMin       = 0.;
+    Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
+    const Float_t kdib = 0.2;
+    Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
 
 
     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
     Int_t i;
-    Float_t oldvalue=0.;
+    Float_t oldvalue= 0.;
+
+    Float_t* b   = new Float_t[kMax];
+    Float_t* si1 = new Float_t[kMax];    
+    Float_t* si2 = new Float_t[kMax];    
     
-    for (i=0; i<kMax; i++)
+    for (i = 0; i < kMax; i++)
     {
-       Float_t xb=bMin+i*kdib;
+       Float_t xb  = bMin+i*kdib;
        Float_t ov;
        ov=fHijing->Profile(xb);
-       Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
-       Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
+       Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
+       Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
        xTot+=gb;
-       xTotHard+=gbh;
+       xTotHard += gbh;
        if (xb > fMinImpactParam && xb < fMaxImpactParam)
        {
-           xPart+=gb;
-           xPartHard+=gbh;
+           xPart += gb;
+           xPartHard += gbh;
        }
        
        if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
-       oldvalue=xTot;
+       oldvalue = xTot;
        printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
        printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
+       if (i>0) {
+           si1[i] = gb/kdib;
+           si2[i] = gbh/gb;
+           b[i]  = xb;
+       }
     }
+
     printf("\n Total cross section (barn): %f \n",xTot);
     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
+
+//  Store result as a graph
+    b[0] = 0;
+    si1[0] = 0;
+    si2[0]=si2[1];
+    
+    fDsigmaDb  = new TGraph(i, b, si1);
+    fDnDb      = new TGraph(i, b, si2);
 }
 
 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
@@ -339,24 +481,25 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* part
 // Looks recursively if one of the daughters has been selected
 //
 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
-    Int_t imin=-1;
-    Int_t imax=-1;
+    Int_t imin = -1;
+    Int_t imax = -1;
     Int_t i;
-    Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
-    Bool_t selected=kFALSE;
+    Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
+    Bool_t selected = kFALSE;
     if (hasDaughters) {
-       imin=iparticle->GetFirstDaughter();
-       imax=iparticle->GetLastDaughter();       
-       for (i=imin; i<= imax; i++){
-           TParticle *  jparticle       = (TParticle *) particles->At(i);      
-           Int_t ip=jparticle->GetPdgCode();
-           if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
+       imin = iparticle->GetFirstDaughter();
+       imax = iparticle->GetLastDaughter();       
+       for (i = imin; i <= imax; i++){
+           TParticle *  jparticle = (TParticle *) particles->At(i);    
+           Int_t ip = jparticle->GetPdgCode();
+           if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
+               selected=kTRUE; break;
+           }
            if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
        }
     } else {
        return kFALSE;
     }
-
     return selected;
 }
 
@@ -369,33 +512,130 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid)
 // 5: beauty
     if (fFlavor == 0) return kTRUE;
     
-    Int_t ifl=TMath::Abs(pid/100);
+    Int_t ifl = TMath::Abs(pid/100);
     if (ifl > 10) ifl/=10;
-    return ((fFlavor==4 && (ifl==4 || ifl==5))  || 
-           (fFlavor==5 &&  ifl==5));
+    return (fFlavor == ifl);
+}
 
+Bool_t AliGenHijing::Stable(TParticle*  particle)
+{
+// Return true for a stable particle
+//
+    if (particle->GetFirstDaughter() < 0 )
+    {
+       return kTRUE;
+    } else {
+       return kFALSE;
+    }
 }
 
+
+void AliGenHijing::Boost(TClonesArray* particles)
+{
+//
+// 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 = particles->GetEntriesFast();
+    for (i = 0; i < np; i++) 
+    {
+       TParticle* iparticle = (TParticle*) particles->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
-    AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
-//    header->SetDate(date);
-//    header->SetRunNumber(run);
-//    header->SetEventNumber(event);
-    header->SetNProduced(fHijing->GetNATT());
-    header->SetImpactParameter(fHijing->GetHINT1(19));
-    header->SetTotalEnergy(fHijing->GetEATT());
-    header->SetHardScatters(fHijing->GetJATT());
-    header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
-    header->SetCollisions(fHijing->GetN0(),
-                         fHijing->GetN01(),
-                         fHijing->GetN10(),
-                         fHijing->GetN11());
+    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());
+    ((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);
+// Event Vertex
+    header->SetPrimaryVertex(fEventVertex);
+    gAlice->SetGenEventHeader(header);    
 }
 
+
 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
 {
 // Assignment operator
     return *this;
 }
+
+#ifndef WIN32
+# define rluget_hijing rluget_hijing_
+# define rluset_hijing rluset_hijing_
+# define rlu_hijing    rlu_hijing_
+# define type_of_call
+#else
+# define rluget_hijing RLUGET_HIJING
+# define rluset_hijing RLUSET_HIJING
+# define rlu_hijing    RLU_HIJING
+# define type_of_call _stdcall
+#endif
+
+
+extern "C" {
+  void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
+  {printf("Dummy version of rluget_hijing reached\n");}
+
+  void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
+  {printf("Dummy version of rluset_hijing reached\n");}
+
+  Double_t type_of_call rlu_hijing(Int_t & /*idum*/) 
+  {
+      Float_t r;
+      do r=sRandom->Rndm(); while(0 >= r || r >= 1);
+      return r;
+  }
+}