]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - THijing/AliGenHijing.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / THijing / AliGenHijing.cxx
index ac40008c32885244d54e7286348c37ad7a913f96..43cbbd8fd793212b2ee754ec5ea5cdf518e84e7c 100644 (file)
@@ -31,8 +31,9 @@
 
 #include "AliGenHijing.h"
 #include "AliGenHijingEventHeader.h"
-#include "AliRun.h"
 #include "AliHijingRndm.h"
+#include "AliLog.h"
+#include "AliRun.h"
 
 ClassImp(AliGenHijing)
 
@@ -44,18 +45,17 @@ AliGenHijing::AliGenHijing()
      fKeep(0),
      fQuench(1),
      fShadowing(1),
-     fDecaysOff(1),
+     fDecaysOff(3),
      fTrigger(0),     
      fEvaluate(0),
      fSelectAll(0),
      fFlavor(0),
-     fEnergyCMS(5500.),
      fKineBias(0.),
      fTrials(0),
      fXsection(0.),
      fHijing(0),
-     fPtHardMin(0.),
-     fPtHardMax(1.e4),
+     fPtHardMin(2.0),
+     fPtHardMax(-1),
      fSpectators(1),
      fDsigmaDb(0),
      fDnDb(0),
@@ -73,10 +73,14 @@ AliGenHijing::AliGenHijing()
      fTargetSpecp(0),
      fLHC(kFALSE),
      fRandomPz(kFALSE),
-     fNoHeavyQuarks(kFALSE)
+     fNoHeavyQuarks(kFALSE),
+     fHeader(AliGenHijingEventHeader("Hijing")),
+     fSigmaNN(-1),
+     fNoElas(0)
 {
-// Constructor
-    AliHijingRndm::SetHijingRandom(GetRandom());
+  // Constructor
+  fEnergyCMS = 5500.;
+  AliHijingRndm::SetHijingRandom(GetRandom());
 }
 
 AliGenHijing::AliGenHijing(Int_t npart)
@@ -87,18 +91,17 @@ AliGenHijing::AliGenHijing(Int_t npart)
      fKeep(0),
      fQuench(1),
      fShadowing(1),
-     fDecaysOff(1),
+     fDecaysOff(3),
      fTrigger(0),     
      fEvaluate(0),
      fSelectAll(0),
      fFlavor(0),
-     fEnergyCMS(5500.),
      fKineBias(0.),
      fTrials(0),
      fXsection(0.),
      fHijing(0),
-     fPtHardMin(0.),
-     fPtHardMax(1.e4),
+     fPtHardMin(2.0),
+     fPtHardMax(-1),
      fSpectators(1),
      fDsigmaDb(0),
      fDnDb(0),
@@ -116,16 +119,21 @@ AliGenHijing::AliGenHijing(Int_t npart)
      fTargetSpecp(0),
      fLHC(kFALSE),
      fRandomPz(kFALSE),
-     fNoHeavyQuarks(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";
 //
 //
 // Set random number generator   
     AliHijingRndm::SetHijingRandom(GetRandom());
+    
 }
 
 AliGenHijing::~AliGenHijing()
@@ -152,8 +160,16 @@ void AliGenHijing::Init()
     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
 //
@@ -214,6 +230,7 @@ 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 time0 = 0.;
   Float_t p[3];
   Float_t tof;
 
@@ -230,13 +247,17 @@ void AliGenHijing::Generate()
   fTrials = 0;
   
   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
+  time0 = fTimeOrigin;
+
   if(fVertexSmear == kPerEvent) {
       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
@@ -257,7 +278,6 @@ void AliGenHijing::Generate()
       
       
       Int_t np = fParticles.GetEntriesFast();
-      printf("\n **************************************************%d\n",np);
       Int_t nc = 0;
       if (np == 0 ) continue;
       Int_t i;
@@ -271,15 +291,14 @@ void AliGenHijing::Generate()
       
 //      Get event vertex
 //
-      TParticle *  iparticle = (TParticle *) fParticles.At(0);
       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);
 
@@ -311,7 +330,7 @@ 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;
       
@@ -346,16 +365,11 @@ void AliGenHijing::Generate()
          } // selected
       } // particle loop final state
 
-//
-//    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]) {
@@ -367,9 +381,8 @@ 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() + sign * origin0[2] / 3.e10;
-             if (fPileUpTimeWindow > 0.) tof += tInt;
-             
+             tof = time0+kconv * iparticle->T();
+
              imo = -1;
              TParticle* mother = 0;
              if (hasMother) {
@@ -387,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);
 }
@@ -423,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;
@@ -472,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
@@ -524,7 +544,7 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid)
     return res;
 }
 
-Bool_t AliGenHijing::Stable(TParticle*  particle) const
+Bool_t AliGenHijing::Stable(const TParticle*  particle) const
 {
 // Return true for a stable particle
 //
@@ -542,21 +562,19 @@ Bool_t AliGenHijing::Stable(TParticle*  particle) const
 void AliGenHijing::MakeHeader()
 {
 // Builds the event header, to be called after each event
-    AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
-    ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
-    ((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(fProjectileSpecn, fProjectileSpecp,
-                                                      fTargetSpecn,fTargetSpecp);
-    ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(fHijing->GetHINT1(20));
-    
-
+    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.
 //
@@ -580,13 +598,31 @@ 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(fVertex);
-    AddHeader(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;
 }