Modifications needed to use EXODUS Decayer
authormorsch <andreas.morsch@cern.ch>
Sun, 4 May 2014 06:48:09 +0000 (08:48 +0200)
committermorsch <andreas.morsch@cern.ch>
Sun, 4 May 2014 06:48:09 +0000 (08:48 +0200)
Irem Erdemir

PYTHIA6/AliDecayerPythia.cxx
PYTHIA6/AliDecayerPythia.h
PYTHIA6/AliGenPythia.cxx
PYTHIA6/AliGenPythia.h
PYTHIA6/AliPythia.cxx
PYTHIA6/AliPythia.h

index 6b9a023..0b9be55 100644 (file)
@@ -142,19 +142,58 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p)
 {
 //  Decay a particle
 // 
-    Float_t energy = p->Energy();
-    Float_t theta  = p->Theta();
-    Float_t phi    = p->Phi();
-    if (!fPatchOmegaDalitz) {
-       Lu1Ent(0, idpart, energy, theta, phi);
-    } else {
-       fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
-       Lu1Ent(0, idpart, energy, theta, phi);
-       fPythia->DalitzDecays();
-       fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
-       fPythia->Pyexec();
+  Float_t energy = p->Energy();
+  Float_t theta  = p->Theta();
+  Float_t phi    = p->Phi();
+
+  if(!fDecayerExodus) {
+    Lu1Ent(0, idpart, energy, theta, phi);
+  } else {
+
+  // EXODUS decayer
+    if(idpart == 111){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->PizeroDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 221){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->EtaDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 113){
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->RhoDirect();
     }
-    fPythia->GetPrimaries();
+    else if(idpart == 223){
+      fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->OmegaDirect();
+      fPythia->OmegaDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+    }
+    else if(idpart == 331){
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->EtaprimeDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+    }
+    else if(idpart == 333){
+      fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->PhiDirect();
+      fPythia->PhiDalitz();
+      fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+    }
+    else if(idpart == 443){
+      Lu1Ent(0, idpart, energy, theta, phi);
+      fPythia->JPsiDirect();
+    }  
+  }
+
+  fPythia->GetPrimaries();
 }
 
 Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
index 5545a54..7ea934a 100644 (file)
@@ -29,6 +29,7 @@ public AliDecayer
       {SetForceDecay((Decay_t) decay);}
     virtual void    ForceDecay();
     virtual void    SetPatchOmegaDalitz() {fPatchOmegaDalitz = 1;}
+    virtual void    SetDecayerExodus()    {fDecayerExodus = 1;}
     virtual void    HeavyFlavourOff() {fHeavyFlavour = kFALSE;}
     virtual void    DecayLongLivedParticles()  {fLongLived    = kTRUE;}
     virtual Float_t GetPartialBranchingRatio(Int_t ipart);
@@ -64,10 +65,11 @@ public AliDecayer
     Bool_t      fHeavyFlavour;    //! Flag for heavy flavors
     Bool_t      fLongLived;       //! Flag for long lived particle decay
     Bool_t      fPatchOmegaDalitz;//! Flag to patch the omega Dalitz decays 
+    Bool_t      fDecayerExodus;  //! Flag for EXODUS decayer
     Bool_t      fPi0;             //! Flag for pi0 decay 
     static Bool_t fgInit;         //! initialization flag 
     
-    ClassDef(AliDecayerPythia, 4) // AliDecayer implementation using Pythia  
+    ClassDef(AliDecayerPythia, 5) // AliDecayer implementation using Pythia  
 };
 #endif
 
index 3de870a..eccf0f7 100644 (file)
@@ -79,6 +79,7 @@ AliGenPythia::AliGenPythia():
     fCRoff(0),
     fHadronisation(1),
     fPatchOmegaDalitz(0), 
+    fDecayerExodus(0),
     fNpartons(0),
     fReadFromFile(0),
     fReadLHEF(0),
@@ -195,7 +196,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fGfinal(kTRUE),
      fCRoff(kFALSE),
      fHadronisation(kTRUE),
-     fPatchOmegaDalitz(0), 
+     fPatchOmegaDalitz(0),
+     fDecayerExodus(0), 
      fNpartons(0),
      fReadFromFile(kFALSE),
      fReadLHEF(0),
@@ -704,8 +706,29 @@ void AliGenPythia::Generate()
              fPythia->DalitzDecays();
              fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
            } 
-           fPythia->Pyexec();
-       }
+         
+      else  if (fDecayerExodus) {
+
+        fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+        fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+        fPythia->Pyexec();
+        fPythia->OmegaDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+        fPythia->PizeroDalitz();
+        fPythia->PhiDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+        fPythia->EtaDalitz();
+        fPythia->EtaprimeDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+        fPythia->RhoDirect();
+        fPythia->OmegaDirect();
+        fPythia->PhiDirect();
+        fPythia->JPsiDirect();
+      }
+    
+      fPythia->Pyexec();
+  }
        fTrials++;
        fPythia->ImportParticles(&fParticles,"All");
        
index e888bba..99fbf85 100644 (file)
@@ -200,6 +200,7 @@ class AliGenPythia : public AliGenMC
     {fpyquenT = t0; fpyquenTau = tau0; fpyquenNf=nf;fpyquenEloss=iengl;fpyquenAngle=iangl;}
     virtual void SetHadronisation(Int_t flag = 1) {fHadronisation = flag;}
     virtual void SetPatchOmegaDalitz(Int_t flag = 1) {fPatchOmegaDalitz = flag;}
+    virtual void SetDecayerExodus(Int_t flag = 1) {fDecayerExodus = flag;}
     virtual void SetReadFromFile(const Text_t *filname) {fkFileName = filname;  fReadFromFile = 1;}    
     virtual void SetReadLHEF(const Text_t *filename) {fkNameLHEF = filename; fReadLHEF = 1;}    
 
@@ -284,6 +285,7 @@ class AliGenPythia : public AliGenMC
     Int_t       fCRoff;             //color reconnection off in the pythia6 annealying model
     Int_t       fHadronisation;     //hadronisation
     Bool_t      fPatchOmegaDalitz;  //flag for omega dalitz decay patch
+    Bool_t      fDecayerExodus;     //flag for exodus decayer
     Int_t       fNpartons;          //Number of partons before hadronisation
     Int_t       fReadFromFile;      //read partons from file
     Int_t       fReadLHEF;          //read lhef file
index 360cd8e..d8a2679 100644 (file)
@@ -20,6 +20,7 @@
 #include "AliFastGlauber.h"
 #include "AliQuenchingWeights.h"
 #include "AliOmegaDalitz.h"
+#include "AliDecayerExodus.h"
 #include "AliLog.h"
 #include "TVector3.h"
 #include "TLorentzVector.h"
@@ -79,7 +80,8 @@ AliPythia::AliPythia():
     fGlauber(0),
     fQuenchingWeights(0),
     fItune(-1),
-    fOmegaDalitz()
+    fOmegaDalitz(),
+    fExodus()
 {
 // Default Constructor
 //
@@ -109,7 +111,8 @@ AliPythia::AliPythia(const AliPythia& pythia):
     fGlauber(0),
     fQuenchingWeights(0),
     fItune(-1),
-    fOmegaDalitz()
+    fOmegaDalitz(),
+    fExodus()
 {
     // Copy Constructor
     Int_t i;
@@ -674,6 +677,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       Initialize("CMS",fProjectile,fTarget,fEcms);
     }
     fOmegaDalitz.Init();
+    fExodus.Init();
 }
 
 Int_t AliPythia::CheckedLuComp(Int_t kf)
@@ -1548,3 +1552,235 @@ void AliPythia::DalitzDecays()
     }
   }
 }
+
+
+//
+// Replace all dalitz(pi0,eta,omega,eta',phi) and resonance(rho,omega,phi,jpsi) decays with the correct matrix element decays
+// for di-electron cocktail calculations
+//
+
+
+void AliPythia::PizeroDalitz()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 111)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 2)                                continue;
+  if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11) )
+  continue;
+  TLorentzVector pizero(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &pizero);
+  for (Int_t j = 0; j < 3; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_pion())[2-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+      }
+  }
+}
+
+
+void AliPythia::EtaDalitz()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 221)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 2)                                continue;
+  if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+  continue;
+  TLorentzVector eta(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &eta);
+  for (Int_t j = 0; j < 3; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_eta())[2-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::RhoDirect()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 113)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 1)                                continue;
+  if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+  continue;
+  TLorentzVector rho(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &rho);
+  for (Int_t j = 0; j < 2; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_rho())[1-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::OmegaDalitz()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 223)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 2)                                continue;
+  if ((fPyjets->K[1][fd] != 111) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+  continue;
+  TLorentzVector omegadalitz(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &omegadalitz);
+  for (Int_t j = 0; j < 3; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_omega_dalitz())[2-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::OmegaDirect()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 223)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 1)                                continue;
+  if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+  continue;
+  TLorentzVector omegadirect(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &omegadirect);
+  for (Int_t j = 0; j < 2; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_omega())[1-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::EtaprimeDalitz()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 331)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 2)                                continue;
+  if ((fPyjets->K[1][fd] != 22) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+  continue;
+  TLorentzVector etaprime(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &etaprime);
+  for (Int_t j = 0; j < 3; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_etaprime())[2-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::PhiDalitz()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 333)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 2)                                continue;
+  if ((fPyjets->K[1][fd] != 221) || (TMath::Abs(fPyjets->K[1][fd+1]) != 11))
+  continue;
+  TLorentzVector phidalitz(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &phidalitz);
+  for (Int_t j = 0; j < 3; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_phi_dalitz())[2-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
+void AliPythia::PhiDirect()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 333)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 1)                                continue;
+  if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+  continue;
+  TLorentzVector phi(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &phi);
+  for (Int_t j = 0; j < 2; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_phi())[1-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+void AliPythia::JPsiDirect()
+{
+
+  Int_t nt = fPyjets->N;
+  for (Int_t i = 0; i < nt; i++) {
+  if (fPyjets->K[1][i] != 443)                       continue;
+  Int_t fd = fPyjets->K[3][i] - 1;
+  Int_t ld = fPyjets->K[4][i] - 1;
+  if (fd < 0)                                        continue;
+  if ((ld - fd) != 1)                                continue;
+  if ((TMath::Abs(fPyjets->K[1][fd]) != 11))
+  continue;
+  TLorentzVector jpsi(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+  Int_t pdg = TMath::Abs(fPyjets->K[1][i]);
+  fExodus.Decay(pdg, &jpsi);
+  for (Int_t j = 0; j < 2; j++) {
+  for (Int_t k = 0; k < 4; k++) {
+  TLorentzVector vec = (fExodus.Products_jpsi())[1-j];
+  fPyjets->P[k][fd+j] = vec[k];
+        }
+     }
+  }
+}
+
+
index b294d6d..bc153c7 100644 (file)
@@ -10,6 +10,7 @@
 #include <AliStructFuncType.h>
 #include "PythiaProcesses.h"
 #include "AliOmegaDalitz.h"
+#include "AliDecayerExodus.h"
 class AliFastGlauber;
 class AliQuenchingWeights;
 
@@ -48,6 +49,16 @@ class AliPythia : public TPythia6, public AliRndm
     static  AliPythia* Instance();
     virtual void Quench();
     void DalitzDecays();
+    // Dalitz and resonance decays from EXODUS
+    void PizeroDalitz();
+    void EtaDalitz();
+    void RhoDirect();
+    void OmegaDalitz();
+    void OmegaDirect();
+    void EtaprimeDalitz();
+    void PhiDalitz();
+    void PhiDirect();
+    void JPsiDirect();
     // Assignment Operator
     AliPythia & operator=(const AliPythia & rhs);
     void Copy(TObject&) const;
@@ -68,6 +79,7 @@ class AliPythia : public TPythia6, public AliRndm
     AliQuenchingWeights*  fQuenchingWeights;  //  ! The Quenching Weights model
     Int_t                 fItune;             //  ! Pythia tune 
     AliOmegaDalitz        fOmegaDalitz;       //  ! omega dalitz decayer
+    AliDecayerExodus      fExodus;            // ! EXODUS decayer
     static AliPythia*     fgAliPythia;        // Pointer to single instance
  private: 
     AliPythia();