From d7d0c637440195733f0c3e9a9a45b570e4c612da Mon Sep 17 00:00:00 2001 From: morsch Date: Sun, 4 May 2014 08:48:09 +0200 Subject: [PATCH] Modifications needed to use EXODUS Decayer Irem Erdemir --- PYTHIA6/AliDecayerPythia.cxx | 63 +++++++-- PYTHIA6/AliDecayerPythia.h | 4 +- PYTHIA6/AliGenPythia.cxx | 29 ++++- PYTHIA6/AliGenPythia.h | 2 + PYTHIA6/AliPythia.cxx | 240 ++++++++++++++++++++++++++++++++++- PYTHIA6/AliPythia.h | 12 ++ 6 files changed, 332 insertions(+), 18 deletions(-) diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx index 6b9a0237866..0b9be557115 100644 --- a/PYTHIA6/AliDecayerPythia.cxx +++ b/PYTHIA6/AliDecayerPythia.cxx @@ -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) diff --git a/PYTHIA6/AliDecayerPythia.h b/PYTHIA6/AliDecayerPythia.h index 5545a54cf2e..7ea934a7c5b 100644 --- a/PYTHIA6/AliDecayerPythia.h +++ b/PYTHIA6/AliDecayerPythia.h @@ -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 diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx index 3de870a0465..eccf0f72495 100644 --- a/PYTHIA6/AliGenPythia.cxx +++ b/PYTHIA6/AliGenPythia.cxx @@ -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"); diff --git a/PYTHIA6/AliGenPythia.h b/PYTHIA6/AliGenPythia.h index e888bba64b6..99fbf85afef 100644 --- a/PYTHIA6/AliGenPythia.h +++ b/PYTHIA6/AliGenPythia.h @@ -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 diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx index 360cd8e5994..d8a26794082 100644 --- a/PYTHIA6/AliPythia.cxx +++ b/PYTHIA6/AliPythia.cxx @@ -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]; + } + } + } +} + + diff --git a/PYTHIA6/AliPythia.h b/PYTHIA6/AliPythia.h index b294d6d67e9..bc153c72fbd 100644 --- a/PYTHIA6/AliPythia.h +++ b/PYTHIA6/AliPythia.h @@ -10,6 +10,7 @@ #include #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(); -- 2.43.0