From 9d7108a742034479e0f375d007337835dd2f1ccc Mon Sep 17 00:00:00 2001
From: hristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Fri, 10 Jun 2005 06:52:23 +0000
Subject: [PATCH] Pileup generation (C.Cheshkov)

---
 PYTHIA6/AliGenPythia.cxx | 66 ++++++++++++++++++++++++++++++++++++++--
 PYTHIA6/AliGenPythia.h   | 14 ++++++---
 2 files changed, 74 insertions(+), 6 deletions(-)

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 8c4d5cd145c..90f714450ce 100644
--- 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
diff --git a/PYTHIA6/AliGenPythia.h b/PYTHIA6/AliGenPythia.h
index 9aa53b7bab9..c3381a518e2 100644
--- a/PYTHIA6/AliGenPythia.h
+++ b/PYTHIA6/AliGenPythia.h
@@ -92,9 +92,9 @@ class AliGenPythia : public AliGenMC
     virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
     virtual void SetReadFromFile(const Text_t *filname) {fFileName = filname;  fReadFromFile = 1;}    
     
-    // 
-    virtual void     SetInteractionRate(Float_t rate) {fInteractionRate = rate;}
-    virtual Float_t  GetInteractionRate()             {return fInteractionRate;}
+    // Get interaction rate for pileup studies
+    virtual void    SetInteractionRate(Float_t rate,Float_t timewindow = 90.e-6);
+    virtual Float_t GetInteractionRate() const {return fInteractionRate;}
     
     // get cross section of process
     virtual Float_t GetXsection() const {return fXsection;}
@@ -128,12 +128,15 @@ class AliGenPythia : public AliGenMC
     
     // Assignment Operator
     AliGenPythia & operator=(const AliGenPythia & rhs);
+
+    void     GetSubEventTime();
+
  protected:
     // adjust the weight from kinematic cuts
     void     AdjustWeights();
     Int_t    GenerateMB();
     void     MakeHeader();    
-    void     GetSubEventTime();
+    void     GeneratePileup();
     
     Process_t   fProcess;           //Process type
     StrucFunc_t fStrucFunc;         //Structure Function
@@ -146,6 +149,9 @@ class AliGenPythia : public AliGenMC
     Float_t     fX2;                //Mean x2
     Float_t     fEventTime;         //Time of the subevent
     Float_t     fInteractionRate;   //Interaction rate (set by user)
+    Float_t     fTimeWindow;        //Time window for pileup events (set by user)
+    Int_t       fCurSubEvent;       //Index of the current sub-event
+    TArrayF     *fEventsTime;       //Subevents time for pileup
     Int_t       fNev;               //Number of events 
     Int_t       fFlavorSelect;      //Heavy Flavor Selection
     Float_t     fXsection;          //Cross-section
-- 
2.43.0