Possibility to simulate pile-up event with correct time distribution.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 May 2005 04:18:19 +0000 (04:18 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 May 2005 04:18:19 +0000 (04:18 +0000)
PYTHIA6/AliGenPythia.cxx
PYTHIA6/AliGenPythia.h

index 44084fb..7b65e67 100644 (file)
@@ -53,6 +53,7 @@ AliGenPythia::AliGenPythia()
   fPythia    = 0;
   fHeader = 0;
   fReadFromFile = 0;
+  fEventTime = 0.;
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
@@ -79,6 +80,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fTitle= "Particle Generator using PYTHIA";
     fXsection  = 0.;
     fReadFromFile = 0;
+    fEventTime = 0.;
     SetProcess();
     SetStrucFunc();
     SetForceDecay();
@@ -324,6 +326,8 @@ void AliGenPythia::Generate()
     Int_t jev=0;
     Int_t j, kf;
     fTrials=0;
+    fEventTime = 0.;
+    
     
 
     //  Set collision vertex position 
@@ -414,6 +418,7 @@ void AliGenPythia::Generate()
                kf = TMath::Abs(kf);
                Int_t kfl = kf;
                // Resonance
+
                if (kfl > 100000) kfl %= 100000;
                if (kfl > 10000)  kfl %= 10000;
                // meson ?
@@ -422,11 +427,29 @@ void AliGenPythia::Generate()
                if (kfl > 10) kfl/=10;
                Int_t ipa = iparticle->GetFirstMother()-1;
                Int_t kfMo = 0;
+//
+// Establish mother daughter relation between heavy quarks and mesons
+//
+               if (kf >= fFlavorSelect && kf <= 6) {
+                   Int_t idau = iparticle->GetFirstDaughter() - 1;
+                   if (idau > -1) {
+                       TParticle* daughter = (TParticle *) fParticles->At(idau);
+                       Int_t pdgD = daughter->GetPdgCode();
+                       if (pdgD == 91 || pdgD == 92) {
+                           Int_t jmin = daughter->GetFirstDaughter() - 1;
+                           Int_t jmax = daughter->GetLastDaughter()  - 1;                          
+                           for (Int_t j = jmin; j <= jmax; j++)
+                               ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+                       } // is string or cluster
+                   } // has daughter
+               } // heavy quark
                
+
                if (ipa > -1) {
                    TParticle *  mother = (TParticle *) fParticles->At(ipa);
                    kfMo = TMath::Abs(mother->GetPdgCode());
                }
+               
                // What to keep in Stack?
                Bool_t flavorOK = kFALSE;
                Bool_t selectOK = kFALSE;
@@ -545,6 +568,8 @@ void AliGenPythia::Generate()
        } else {
            nc = GenerateMB();
        } // mb ?
+       
+       GetSubEventTime();
 
        if (pParent)   delete[] pParent;
        if (pSelected) delete[] pSelected;
@@ -612,7 +637,7 @@ Int_t  AliGenPythia::GenerateMB()
        if (!CheckTrigger(jet1, jet2)) return 0;
     }
     
-    for (i = 0; i<np; i++) {
+    for (i = 0; i < np; i++) {
        Int_t trackIt = 0;
        TParticle *  iparticle = (TParticle *) fParticles->At(i);
        kf = CheckPDGCode(iparticle->GetPdgCode());
@@ -639,7 +664,7 @@ Int_t  AliGenPythia::GenerateMB()
            origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
            origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
            
-           Float_t tof=kconv*iparticle->T();
+           Float_t tof = fEventTime + kconv * iparticle->T();
 
            PushTrack(fTrackIt*trackIt, iparent, kf, 
                      p[0], p[1], p[2], p[3], 
@@ -662,13 +687,15 @@ Int_t  AliGenPythia::GenerateMB()
            
            KeepTrack(nt);
            pParent[i] = nt;
+           SetHighWaterMark(nt);
+           
        } // select particle
     } // particle loop 
 
     if (pParent) delete[] pParent;
     
     printf("\n I've put %i particles on the stack \n",nc);
-    return nc;
+    return 1;
 }
 
 
@@ -715,7 +742,7 @@ void AliGenPythia::MakeHeader()
 {
   if (gAlice) {
     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
-       gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
+       gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
   }
 
 // Builds the event header, to be called after each event
@@ -732,6 +759,7 @@ void AliGenPythia::MakeHeader()
     fHeader->SetPrimaryVertex(fVertex);
 //
 // Jets that have triggered
+
     if (fProcess == kPyJets)
     {
        Int_t ntrig, njet;
@@ -987,6 +1015,11 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
     }
 }
 
+void AliGenPythia::GetSubEventTime()
+{
+  // Calculates time of the next subevent
+  fEventTime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+}
 
 #ifdef never
 void AliGenPythia::Streamer(TBuffer &R__b)
index c91fea1..9aa53b7 100644 (file)
@@ -91,7 +91,11 @@ class AliGenPythia : public AliGenMC
     virtual void SetQuench(Int_t flag = 0) {fQuench = flag;}
     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 cross section of process
     virtual Float_t GetXsection() const {return fXsection;}
     // get triggered jets
@@ -126,10 +130,10 @@ class AliGenPythia : public AliGenMC
     AliGenPythia & operator=(const AliGenPythia & rhs);
  protected:
     // adjust the weight from kinematic cuts
-    void   AdjustWeights();
-    Int_t  GenerateMB();
-    void   MakeHeader();    
-
+    void     AdjustWeights();
+    Int_t    GenerateMB();
+    void     MakeHeader();    
+    void     GetSubEventTime();
     
     Process_t   fProcess;           //Process type
     StrucFunc_t fStrucFunc;         //Structure Function
@@ -140,6 +144,8 @@ class AliGenPythia : public AliGenMC
     Float_t     fQ;                 //Mean Q
     Float_t     fX1;                //Mean x1
     Float_t     fX2;                //Mean x2
+    Float_t     fEventTime;         //Time of the subevent
+    Float_t     fInteractionRate;   //Interaction rate (set by user)
     Int_t       fNev;               //Number of events 
     Int_t       fFlavorSelect;      //Heavy Flavor Selection
     Float_t     fXsection;          //Cross-section