Pileup generation (C.Cheshkov)

[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 8c4d5cd145cec9054aad6db2a22cf001c925148c..90f714450ce72c695c8cebeea3e119829da9cd64 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -55,6 +55,9 @@ AliGenPythia::AliGenPythia()
   fReadFromFile = 0;
   fEventTime = 0.;
   fInteractionRate = 0.;
+  fTimeWindow = 0.;
+  fEventsTime = 0;
+  fCurSubEvent = 0;
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
@@ -83,6 +86,9 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fReadFromFile = 0;
     fEventTime = 0.;
     fInteractionRate = 0.;
+    fTimeWindow = 0.;
+    fEventsTime = 0;
+    fCurSubEvent = 0;
     SetProcess();
     SetStrucFunc();
     SetForceDecay();
@@ -129,6 +135,57 @@ AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 AliGenPythia::~AliGenPythia()
 {
 // Destructor
+  if(fEventsTime) delete fEventsTime;
+}
+
+void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
+{
+// Generate pileup using user specified rate
+    fInteractionRate = rate;
+    fTimeWindow = timewindow;
+    GeneratePileup();
+}
+
+void AliGenPythia::GeneratePileup()
+{
+// Generate sub events time for pileup
+    fEventsTime = 0;
+    if(fInteractionRate == 0.) {
+      Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
+      return;
+    }
+
+    Int_t npart = NumberParticles();
+    if(npart < 0) {
+      Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
+      return;
+    }
+
+    if(fEventsTime) delete fEventsTime;
+    fEventsTime = new TArrayF(npart);
+    TArrayF &array = *fEventsTime;
+    for(Int_t ipart = 0; ipart < npart; ipart++)
+      array[ipart] = 0.;
+
+    Float_t eventtime = 0.;
+    while(1)
+      {
+       eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(eventtime > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    eventtime = 0.;
+    while(1)
+      {
+       eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(TMath::Abs(eventtime) > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    SetNumberParticles(fEventsTime->GetSize());
 }
 
 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
@@ -1020,8 +1077,13 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 void AliGenPythia::GetSubEventTime()
 {
   // Calculates time of the next subevent
-  if (fInteractionRate != 0.)
-    fEventTime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+  fEventTime = 0.;
+  if (fEventsTime) {
+    TArrayF &array = *fEventsTime;
+    fEventTime = array[fCurSubEvent++];
+  }
+  //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
+  return;
 }
 
 #ifdef never