]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - THijing/AliGenHijing.cxx
Use AliLog for messages
[u/mrichter/AliRoot.git] / THijing / AliGenHijing.cxx
index f456abceef2779c3ac6039b02c726c53f5a3ba3d..645b05da9ed20883250541ff1719bd33d688d7cc 100644 (file)
@@ -22,6 +22,7 @@
 // Andreas Morsch    (andreas.morsch@cern.ch)
 //
 
+#include <TClonesArray.h>
 #include <TGraph.h>
 #include <THijing.h>
 #include <TLorentzVector.h>
 
 #include "AliGenHijing.h"
 #include "AliGenHijingEventHeader.h"
-#include "AliRun.h"
 #include "AliHijingRndm.h"
+#include "AliLog.h"
+#include "AliRun.h"
 
 ClassImp(AliGenHijing)
 
 AliGenHijing::AliGenHijing()
-                 :AliGenMC()
+    :AliGenMC(),
+     fFrame("CMS"),
+     fMinImpactParam(0.),
+     fMaxImpactParam(5.),
+     fKeep(0),
+     fQuench(1),
+     fShadowing(1),
+     fDecaysOff(1),
+     fTrigger(0),     
+     fEvaluate(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fKineBias(0.),
+     fTrials(0),
+     fXsection(0.),
+     fHijing(0),
+     fPtHardMin(0.),
+     fPtHardMax(1.e4),
+     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)
 {
-// Constructor
-    fParticles = 0;
-    fHijing    = 0;
-    fDsigmaDb  = 0;
-    fDnDb      = 0;
-    AliHijingRndm::SetHijingRandom(GetRandom());
+  // 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(1),
+     fTrigger(0),     
+     fEvaluate(0),
+     fSelectAll(0),
+     fFlavor(0),
+     fKineBias(0.),
+     fTrials(0),
+     fXsection(0.),
+     fHijing(0),
+     fPtHardMin(0.),
+     fPtHardMax(1.e4),
+     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)
 {
 // 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;
-    //
-    SetSimpleJets();
-    SetNoGammas();
-    SetRandomPz();
-    SwitchOffHeavyQuarks(kFALSE);
-//
-    fParticles = new TClonesArray("TParticle",10000);    
+//
 //
 // Set random number generator   
     AliHijingRndm::SetHijingRandom(GetRandom());
-    fHijing = 0;
-
 }
 
-AliGenHijing::AliGenHijing(const AliGenHijing & hijing):
-    AliGenMC(hijing)
-{
-// copy constructor
-}
-
-
 AliGenHijing::~AliGenHijing()
 {
 // Destructor
     if ( fDsigmaDb) delete  fDsigmaDb;  
     if ( fDnDb)     delete  fDnDb;  
-    delete fParticles;
 }
 
 void AliGenHijing::Init()
@@ -185,16 +219,17 @@ void AliGenHijing::Generate()
   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];
   if(fVertexSmear == kPerEvent) {
       Vertex();
@@ -214,15 +249,15 @@ void AliGenHijing::Generate()
 // --------------------------------------------------------------------------
       fHijing->GenerateEvent();
       fTrials++;
-      fHijing->ImportParticles(fParticles,"All");
+      fNprimaries = 0;
+      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 np = fParticles.GetEntriesFast();
       Int_t nc = 0;
       if (np == 0 ) continue;
       Int_t i;
@@ -236,7 +271,7 @@ void AliGenHijing::Generate()
       
 //      Get event vertex
 //
-      TParticle *  iparticle = (TParticle *) fParticles->At(0);
+      TParticle *  iparticle = (TParticle *) fParticles.At(0);
       fVertex[0] = origin0[0];
       fVertex[1] = origin0[1]; 
       fVertex[2] = origin0[2];
@@ -246,7 +281,7 @@ void AliGenHijing::Generate()
 //
 
       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;
@@ -276,21 +311,22 @@ 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){
+         if(ksp == 0 || ksp == 1){
              if(kf == kNeutron) fProjectileSpecn += 1;
              if(kf == kProton)  fProjectileSpecp += 1;
          }
-         else if(ks == 10 || ks == 11){
+         else if(ksp == 10 || ksp == 11){
              if(kf == kNeutron) fTargetSpecn += 1;
              if(kf == kProton)  fTargetSpecp += 1;
          }
@@ -298,8 +334,8 @@ void AliGenHijing::Generate()
 //         
          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
@@ -309,14 +345,19 @@ void AliGenHijing::Generate()
              pSelected[i] = 1;
          } // selected
       } // particle loop final state
+
 //
-// Write particles to stack
+//    Time of the interactions
+      Float_t tInt = 0.;
+      if (fPileUpTimeWindow > 0.) tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
+
 //
+// 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();
@@ -326,17 +367,20 @@ void AliGenHijing::Generate()
              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 = kconv * iparticle->T() + sign * origin0[2] / 3.e10;
+             tof = kconv * iparticle->T();
+             if (fPileUpTimeWindow > 0.) tof += tInt;
+             
              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);
-             PushTrack(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
@@ -344,12 +388,12 @@ 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;
          }
       }
@@ -444,7 +488,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;
@@ -500,7 +544,7 @@ 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)->SetNProduced(fNprimaries);
     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
@@ -511,6 +555,9 @@ void AliGenHijing::MakeHeader()
                                                       fHijing->GetN11());
     ((AliGenHijingEventHeader*) header)->SetSpectators(fProjectileSpecn, fProjectileSpecp,
                                                       fTargetSpecn,fTargetSpecp);
+    ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(fHijing->GetHINT1(20));
+    
+
 
 // 4-momentum vectors of the triggered jets.
 //
@@ -539,10 +586,11 @@ void AliGenHijing::MakeHeader()
     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
 // Event Vertex
     header->SetPrimaryVertex(fVertex);
-    if (gAlice) gAlice->SetGenEventHeader(header);   
+    AddHeader(header);
     fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 }
 
+
 Bool_t AliGenHijing::CheckTrigger()
 {
 // Check the kinematic trigger condition
@@ -579,12 +627,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 && 
@@ -599,16 +647,3 @@ Bool_t AliGenHijing::CheckTrigger()
     } // fTrigger == 2
     return triggered;
 }
-
-
-void AliGenHijing::Copy(TObject &) const
-{
-  Fatal("Copy","Not implemented!\n");
-}
-
-AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
-{
-    rhs.Copy(*this); 
-    return (*this);
-}
-