]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - THijing/AliGenHijing.cxx
Kalman filter vertex in Psi2s task
[u/mrichter/AliRoot.git] / THijing / AliGenHijing.cxx
index 70b924dbf78f986d66c250bf33a68072c40dc61c..43cbbd8fd793212b2ee754ec5ea5cdf518e84e7c 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.1  2003/03/15 14:45:57  morsch
-Classes imported from EVGEN
-
-Revision 1.47  2003/01/14 10:50:18  alibrary
-Cleanup of STEER coding conventions
-
-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)
-
-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.
-
-Revision 1.7  2000/10/02 21:28:06  fca
-Removal of useless dependecies via forward declarations
-
-Revision 1.6  2000/09/11 13:23:37  morsch
-Write last seed to file (fortran lun 50) and reed back from same lun using calls to
-luget_hijing and luset_hijing.
-
-Revision 1.5  2000/09/07 16:55:40  morsch
-fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
-
-Revision 1.4  2000/07/11 18:24:56  fca
-Coding convention corrections + few minor bug fixes
-
-Revision 1.3  2000/06/30 12:08:36  morsch
-In member data: char* replaced by TString, Init takes care of resizing the strings to
-8 characters required by Hijing.
-
-Revision 1.2  2000/06/15 14:15:05  morsch
-Add possibility for heavy flavor selection: charm and beauty.
-
-Revision 1.1  2000/06/09 20:47:27  morsch
-AliGenerator interface class to HIJING using THijing (test version)
-
-*/
-
-
+/* $Id$ */
 
 // 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.
+// Author:
+// Andreas Morsch    (andreas.morsch@cern.ch)
 //
-// andreas.morsch@cern.ch
 
-#include <TArrayI.h>
+#include <TClonesArray.h>
 #include <TGraph.h>
 #include <THijing.h>
 #include <TLorentzVector.h>
@@ -194,73 +31,116 @@ AliGenerator interface class to HIJING using THijing (test version)
 
 #include "AliGenHijing.h"
 #include "AliGenHijingEventHeader.h"
+#include "AliHijingRndm.h"
+#include "AliLog.h"
 #include "AliRun.h"
 
-
-  ClassImp(AliGenHijing)
+ClassImp(AliGenHijing)
 
 AliGenHijing::AliGenHijing()
-                 :AliGenMC()
+    :AliGenMC(),
+     fFrame("CMS"),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fKeep(0),
+     fQuench(1),
+     fShadowing(1),
+     fDecaysOff(3),
+     fTrigger(0),     
+     fEvaluate(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fKineBias(0.),
+     fTrials(0),
+     fXsection(0.),
+     fHijing(0),
+     fPtHardMin(2.0),
+     fPtHardMax(-1),
+     fSpectators(1),
+     fDsigmaDb(0),
+     fDnDb(0),
+     fPtMinJet(-2.5),
+     fEtaMinJet(-20.),
+     fEtaMaxJet(+20.),
+     fPhiMinJet(0.),
+     fPhiMaxJet(2. * TMath::Pi()),
+     fRadiation(3),
+     fSimpleJet(kFALSE),
+     fNoGammas(kFALSE),
+     fProjectileSpecn(0),
+     fProjectileSpecp(0),
+     fTargetSpecn(0),
+     fTargetSpecp(0),
+     fLHC(kFALSE),
+     fRandomPz(kFALSE),
+     fNoHeavyQuarks(kFALSE),
+     fHeader(AliGenHijingEventHeader("Hijing")),
+     fSigmaNN(-1),
+     fNoElas(0)
 {
-// Constructor
-    fParticles = 0;
-    fHijing    = 0;
-    fDsigmaDb  = 0;
-    fDnDb      = 0;
+  // Constructor
+  fEnergyCMS = 5500.;
+  AliHijingRndm::SetHijingRandom(GetRandom());
 }
 
 AliGenHijing::AliGenHijing(Int_t npart)
-    :AliGenMC(npart)
+    :AliGenMC(npart),
+     fFrame("CMS"),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fKeep(0),
+     fQuench(1),
+     fShadowing(1),
+     fDecaysOff(3),
+     fTrigger(0),     
+     fEvaluate(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fKineBias(0.),
+     fTrials(0),
+     fXsection(0.),
+     fHijing(0),
+     fPtHardMin(2.0),
+     fPtHardMax(-1),
+     fSpectators(1),
+     fDsigmaDb(0),
+     fDnDb(0),
+     fPtMinJet(-2.5),
+     fEtaMinJet(-20.),
+     fEtaMaxJet(+20.),
+     fPhiMinJet(0.),
+     fPhiMaxJet(2. * TMath::Pi()),
+     fRadiation(3),
+     fSimpleJet(kFALSE),
+     fNoGammas(kFALSE),
+     fProjectileSpecn(0),
+     fProjectileSpecp(0),
+     fTargetSpecn(0),
+     fTargetSpecp(0),
+     fLHC(kFALSE),
+     fRandomPz(kFALSE),
+     fNoHeavyQuarks(kFALSE),
+     fHeader(AliGenHijingEventHeader("Hijing")),
+     fSigmaNN(-1),
+     fNoElas(0)
 {
 // Default PbPb collisions at 5. 5 TeV
 //
+    fEnergyCMS = 5500.;
     fName = "Hijing";
     fTitle= "Particle Generator using HIJING";
-
-    SetEnergyCMS();
-    SetImpactParameterRange();
-    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)
-{
-// copy constructor
+    AliHijingRndm::SetHijingRandom(GetRandom());
+    
 }
 
-
 AliGenHijing::~AliGenHijing()
 {
 // Destructor
     if ( fDsigmaDb) delete  fDsigmaDb;  
     if ( fDnDb)     delete  fDnDb;  
-    delete fParticles;
 }
 
 void AliGenHijing::Init()
@@ -274,14 +154,22 @@ void AliGenHijing::Init()
                      fAProjectile, fZProjectile, fATarget, fZTarget, 
                      fMinImpactParam, fMaxImpactParam));
 
-    fHijing=(THijing*) fgMCEvGen;
+    fHijing=(THijing*) fMCEvGen;
     fHijing->SetIHPR2(2,  fRadiation);
     fHijing->SetIHPR2(3,  fTrigger);
     fHijing->SetIHPR2(6,  fShadowing);
     fHijing->SetIHPR2(12, fDecaysOff);    
     fHijing->SetIHPR2(21, fKeep);
-    fHijing->SetHIPR1(10, fPtMinJet);  
+    fHijing->SetHIPR1(8,  fPtHardMin);         
+    fHijing->SetHIPR1(9,  fPtHardMax);         
+    fHijing->SetHIPR1(10, fPtMinJet); 
+    if (fSigmaNN>0)
+      fHijing->SetHIPR1(31, fSigmaNN/2.);      
     fHijing->SetHIPR1(50, fSimpleJet);
+    //
+    // Switching off elastic scattering
+    if (fNoElas)
+      fHijing->SetIHPR2(14, 0);
 //
 //  Quenching
 //
@@ -314,7 +202,17 @@ void AliGenHijing::Init()
        fHijing->SetHIPR1(11, 2.5);
     }
     
+//
+// Heavy quarks
+//    
+    if (fNoHeavyQuarks) {
+       fHijing->SetIHPR2(49, 1);
+    } else {
+       fHijing->SetIHPR2(49, 0);
+    }
+    
     
+    AliGenMC::Init();
     
 //
 //  Initialize Hijing  
@@ -332,58 +230,54 @@ void AliGenHijing::Generate()
   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 time0 = 0.;
+  Float_t p[3];
   Float_t tof;
 
 //  converts from mm/c to s
-  const Float_t kconv = 0.001/2.999792458e8;
+  const Float_t kconv = 0.001/2.99792458e8;
 //
   Int_t nt  = 0;
   Int_t jev = 0;
-  Int_t j, kf, ks, imo;
+  Int_t j, kf, ks, ksp, imo;
   kf = 0;
     
 
     
   fTrials = 0;
+  
   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
+  time0 = fTimeOrigin;
+
   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);
-      for (j = 0; j < 3; j++) {
-//           fHijing->SetPARP(151+j, fOsigma[j]*10.);
-      }
-  }
+      Vertex();
+      for (j=0; j < 3; j++) origin0[j] = fVertex[j];
+      time0 = fTime;
+  } 
+
+
+  Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
+
   while(1)
   {
 //    Generate one event
 // --------------------------------------------------------------------------
-      fSpecn   = 0;  
-      fSpecp   = 0;
+      fProjectileSpecn    = 0;  
+      fProjectileSpecp    = 0;
+      fTargetSpecn        = 0;  
+      fTargetSpecp        = 0;
 // --------------------------------------------------------------------------
       fHijing->GenerateEvent();
       fTrials++;
-      fHijing->ImportParticles(fParticles,"All");
+      fNprimaries = 0;
+      fHijing->ImportParticles(&fParticles,"All");
       if (fTrigger != kNoTrigger) {
          if (!CheckTrigger()) continue;
       }
-      Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
-                                         (Double_t(fZTarget)    * Double_t(fAProjectile)));
-      if (fLHC) Boost(dy);
+      if (fLHC) Boost();
       
       
-      Int_t np = fParticles->GetEntriesFast();
-      printf("\n **************************************************%d\n",np);
+      Int_t np = fParticles.GetEntriesFast();
       Int_t nc = 0;
       if (np == 0 ) continue;
       Int_t i;
@@ -397,17 +291,16 @@ void AliGenHijing::Generate()
       
 //      Get event vertex
 //
-      TParticle *  iparticle = (TParticle *) fParticles->At(0);
-      fEventVertex[0] = origin0[0];
-      fEventVertex[1] = origin0[1];    
-      fEventVertex[2] = origin0[2];
-      
+      fVertex[0] = origin0[0];
+      fVertex[1] = origin0[1]; 
+      fVertex[2] = origin0[2];
+      fTime = time0;
 //
 //      First select parent particles
 //
-
+      TParticle *  iparticle = 0;
       for (i = 0; i < np; i++) {
-         iparticle = (TParticle *) fParticles->At(i);
+         iparticle = (TParticle *) fParticles.At(i);
 
 // Is this a parent particle ?
          if (Stable(iparticle)) continue;
@@ -437,26 +330,31 @@ void AliGenHijing::Generate()
 //
 
       for (i = 0; i<np; i++) {
-         TParticle *  iparticle = (TParticle *) fParticles->At(i);
+         iparticle = (TParticle *) fParticles.At(i);
 // Is this a final state particle ?
          if (!Stable(iparticle)) continue;
       
          Bool_t  selected             =  kTRUE;
          kf        = iparticle->GetPdgCode();
          ks        = iparticle->GetStatusCode();
+         ksp       = iparticle->GetUniqueID();
          
 // --------------------------------------------------------------------------
 // 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(ksp == 0 || ksp == 1){
+             if(kf == kNeutron) fProjectileSpecn += 1;
+             if(kf == kProton)  fProjectileSpecp += 1;
+         }
+         else if(ksp == 10 || ksp == 11){
+             if(kf == kNeutron) fTargetSpecn += 1;
+             if(kf == kProton)  fTargetSpecp += 1;
          }
 // --------------------------------------------------------------------------
 //         
          if (!fSelectAll) {
              selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
-             if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
-                                                       && ks != 11);
+             if (!fSpectators && selected) selected = (ksp != 0 && ksp != 1 && ksp != 10
+                                                       && ksp != 11);
          }
 //
 // Put particle on the stack if selected
@@ -466,34 +364,35 @@ void AliGenHijing::Generate()
              pSelected[i] = 1;
          } // selected
       } // particle loop final state
+
 //
 // Write particles to stack
-//
+
       for (i = 0; i<np; i++) {
-         TParticle *  iparticle = (TParticle *) fParticles->At(i);
+         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();
+             p[2] = iparticle->Pz() * sign;
              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();
+             tof = time0+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;
+                 mother = (TParticle *) fParticles.At(imo);
+                 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
              } // if has mother   
              Bool_t tFlag = (fTrackIt && !hasDaughter);
-             SetTrack(tFlag,imo,kf,p,origin,polar,
-                      tof,kPNoProcess,nt, 1., ks);
+             PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
+             fNprimaries++;
              KeepTrack(nt);
              newPos[i] = nt;
          } // if selected
@@ -501,16 +400,17 @@ void AliGenHijing::Generate()
       delete[] newPos;
       delete[] pSelected;
       
-      printf("\n I've put %i particles on the stack \n",nc);
+      AliInfo(Form("\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);
+             AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
              break;
          }
       }
   } // event loop
+
   MakeHeader();
   SetHighWaterMark(nt);
 }
@@ -537,13 +437,19 @@ void AliGenHijing::EvaluateCrossSections()
 
     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
+    printf("\n Inelastic and total cross section (mb) %f %f \n",fHijing->GetHINT1(12), fHijing->GetHINT1(13));    
     Int_t i;
     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++){
+      b[i] = 0.;
+      si1[i] = 0.;
+      si2[i] = 0.;
+    }
+
     for (i = 0; i < kMax; i++)
     {
        Float_t xb  = bMin+i*kdib;
@@ -586,7 +492,7 @@ void AliGenHijing::EvaluateCrossSections()
     fDnDb      = new TGraph(i, b, si2);
 }
 
-Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
+Bool_t AliGenHijing::DaughtersSelection(const TParticle* iparticle)
 {
 //
 // Looks recursively if one of the daughters has been selected
@@ -601,7 +507,7 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
        imin = iparticle->GetFirstDaughter();
        imax = iparticle->GetLastDaughter();       
        for (i = imin; i <= imax; i++){
-           TParticle *  jparticle = (TParticle *) fParticles->At(i);   
+           TParticle *  jparticle = (TParticle *) fParticles.At(i);    
            Int_t ip = jparticle->GetPdgCode();
            if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
@@ -638,7 +544,7 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid)
     return res;
 }
 
-Bool_t AliGenHijing::Stable(TParticle*  particle)
+Bool_t AliGenHijing::Stable(const TParticle*  particle) const
 {
 // Return true for a stable particle
 //
@@ -656,17 +562,19 @@ Bool_t AliGenHijing::Stable(TParticle*  particle)
 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());
-    ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
+    fHeader.SetNProduced(fNprimaries);
+    fHeader.SetImpactParameter(fHijing->GetHINT1(19));
+    fHeader.SetTotalEnergy(fHijing->GetEATT());
+    fHeader.SetHardScatters(fHijing->GetJATT());
+    fHeader.SetParticipants(fHijing->GetNP(), fHijing->GetNT());
+    fHeader.SetCollisions(fHijing->GetN0(),
+                         fHijing->GetN01(),
+                         fHijing->GetN10(),
+                         fHijing->GetN11());
+    fHeader.SetSpectators(fProjectileSpecn, fProjectileSpecp,
+                         fTargetSpecn,fTargetSpecp);
+    fHeader.SetReactionPlaneAngle(fHijing->GetHINT1(20));
+    fHeader.SetTrueNPart(fHijing->GetNPART());
 
 // 4-momentum vectors of the triggered jets.
 //
@@ -690,15 +598,34 @@ void AliGenHijing::MakeHeader()
                                              fHijing->GetHINT1(37),
                                              fHijing->GetHINT1(38),
                                              fHijing->GetHINT1(39));
-    ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
+    fHeader.SetJets(jet1, jet2, jet3, jet4);
 // Bookkeeping for kinematic bias
-    ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
+    fHeader.SetTrials(fTrials);
 // Event Vertex
-    header->SetPrimaryVertex(fEventVertex);
-    gAlice->SetGenEventHeader(header);   
-    fCollisionGeometry = (AliGenHijingEventHeader*)  header;
+    fHeader.SetPrimaryVertex(fVertex);
+    fHeader.SetInteractionTime(fTime);
+
+    Int_t nsd1 = 0,nsd2 = 0,ndd = 0;
+    Int_t nT = fHijing->GetNT();
+    Int_t nP = fHijing->GetNP();
+    for (Int_t i = 1; i <= nP; ++i) {
+      for (Int_t j = 1; j <= nT; ++j) {
+      Int_t tp = fHijing->GetNFP(i, 5);
+      Int_t tt = fHijing->GetNFT(j, 5);
+      if (tp == 2)
+        nsd1++;
+      if (tt == 2)
+        nsd2++;
+      if (tp == 2 && tt == 2)
+        ndd++;
+      }
+    }
+    fHeader.SetNDiffractive(nsd1, nsd2, ndd);
+    AddHeader(&fHeader);
+    fCollisionGeometry = &fHeader;
 }
 
+
 Bool_t AliGenHijing::CheckTrigger()
 {
 // Check the kinematic trigger condition
@@ -735,12 +662,12 @@ Bool_t AliGenHijing::CheckTrigger()
     } else if (fTrigger == 2) {
 //  Gamma Jet
 //
-       Int_t np = fParticles->GetEntriesFast();
+       Int_t np = fParticles.GetEntriesFast();
        for (Int_t i = 0; i < np; i++) {
-           TParticle* part = (TParticle*) fParticles->At(i);
+           TParticle* part = (TParticle*) fParticles.At(i);
            Int_t kf = part->GetPdgCode();
-           Int_t ks = part->GetStatusCode();
-           if (kf == 22 && ks == 40) {
+           Int_t ksp = part->GetUniqueID();
+           if (kf == 22 && ksp == 40) {
                Float_t phi = part->Phi();
                Float_t eta = part->Eta();
                if  (eta < fEtaMaxJet && 
@@ -755,40 +682,3 @@ Bool_t AliGenHijing::CheckTrigger()
     } // fTrigger == 2
     return triggered;
 }
-
-
-
-
-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;
-  }
-}