From 173264b0f4f94e25ac1fa1c6a08153cee847786c Mon Sep 17 00:00:00 2001 From: DanielSamitz Date: Wed, 18 Dec 2024 00:17:31 +0100 Subject: [PATCH] add psi(2S) and upsilon --- GeneratorParam/ExodusDecayer.cxx | 61 ++++++-- GeneratorParam/ExodusDecayer.h | 8 + GeneratorParam/GeneratorParamEMlibV2.cxx | 182 +++++++++++++++++------ GeneratorParam/GeneratorParamEMlibV2.h | 36 +++-- GeneratorParam/PythiaDecayerConfig.cxx | 60 ++++++++ GeneratorParam/PythiaDecayerConfig.h | 2 + 6 files changed, 282 insertions(+), 67 deletions(-) diff --git a/GeneratorParam/ExodusDecayer.cxx b/GeneratorParam/ExodusDecayer.cxx index dfeedca..b5de3b9 100644 --- a/GeneratorParam/ExodusDecayer.cxx +++ b/GeneratorParam/ExodusDecayer.cxx @@ -29,6 +29,8 @@ ExodusDecayer::ExodusDecayer(): fEPMassPhiDalitz(0), fEPMassPhiDalitz_toPi0(0), fEPMassJPsi(0), + fEPMassPsi2S(0), + fEPMassUpsilon(0), fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)), /* Polarization Function for resonances */ fInit(0), fDecayToDimuon(0) @@ -51,6 +53,8 @@ ExodusDecayer::~ExodusDecayer() delete fEPMassPhiDalitz; delete fEPMassPhiDalitz_toPi0; delete fEPMassJPsi; + delete fEPMassPsi2S; + delete fEPMassUpsilon; delete fPol; } @@ -74,8 +78,8 @@ void ExodusDecayer::Init() Float_t binwidth; Float_t mass_bin, mass_min, mass_max; - Double_t vmass_eta, vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_eta, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi; - Double_t weight_eta, weight_rho, weight_omega, weight_phi, weight_jpsi; + Double_t vmass_eta, vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vmass_psi2s, vmass_upsilon, vwidth_eta, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi, vwidth_psi2s, vwidth_upsilon; + Double_t weight_eta, weight_rho, weight_omega, weight_phi, weight_jpsi, weight_psi2s, weight_upsilon; //================================================================================// // Create electron pair mass histograms from dalitz decays // @@ -83,7 +87,7 @@ void ExodusDecayer::Init() // Get the particle masses // parent - nbins = 1000; + nbins = 2000; mass_min = 0.; mass_max = 0.; pionmass = (TDatabasePDG::Instance()->GetParticle(111))->Mass(); @@ -315,6 +319,8 @@ void ExodusDecayer::Init() vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass(); vmass_phi = (TDatabasePDG::Instance()->GetParticle(333))->Mass(); vmass_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Mass(); + vmass_psi2s = (TDatabasePDG::Instance()->GetParticle(100443))->Mass(); + vmass_upsilon = (TDatabasePDG::Instance()->GetParticle(553))->Mass(); // Get the particle widths // parent vwidth_eta = (TDatabasePDG::Instance()->GetParticle(221))->Width(); @@ -322,22 +328,32 @@ void ExodusDecayer::Init() vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width(); vwidth_phi = (TDatabasePDG::Instance()->GetParticle(333))->Width(); vwidth_jpsi = (TDatabasePDG::Instance()->GetParticle(443))->Width(); + vwidth_psi2s = (TDatabasePDG::Instance()->GetParticle(100443))->Width(); + vwidth_upsilon = (TDatabasePDG::Instance()->GetParticle(553))->Width(); if( vwidth_jpsi < 1e-6 ){ printf(Form("ExodusDecayer: Warning: Width JPsi = %f (Mass = %f) < 1e-6, set to 1e-6\n",vwidth_jpsi,vmass_jpsi)); vwidth_jpsi = 1e-6; } + if( vwidth_psi2s < 1e-6 ){ + printf(Form("ExodusDecayer: Warning: Width Psi2S = %f (Mass = %f) < 1e-6, set to 1e-6\n",vwidth_psi2s,vmass_psi2s)); + vwidth_psi2s = 1e-6; + } + if( vwidth_upsilon < 1e-6 ){ + printf(Form("ExodusDecayer: Warning: Width Upsilon = %f (Mass = %f) < 1e-6, set to 1e-6\n",vwidth_upsilon,vmass_upsilon)); + vwidth_upsilon = 1e-6; + } if ( mass_min == 0. && mass_max == 0. ) { mass_min = 2.*pimass; - mass_max = 5.; + mass_max = 10.; } binwidth = (mass_max-mass_min)/(Double_t)nbins; - // create pair mass histograms for resonances of rho, omega, phi and jpsi + // create pair mass histograms for resonances of rho, omega, phi, jpsi, psi2s and upsilon fEPMassEta = new TH1F("fEPMassEta","mass eta",nbins,mass_min,mass_max); fEPMassEta->SetDirectory(0); fEPMassRho = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max); @@ -348,6 +364,10 @@ void ExodusDecayer::Init() fEPMassPhi->SetDirectory(0); fEPMassJPsi = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max); fEPMassJPsi->SetDirectory(0); + fEPMassPsi2S = new TH1F("fEPMassPsi2S","mass psi2s",nbins,mass_min,mass_max); + fEPMassPsi2S->SetDirectory(0); + fEPMassUpsilon = new TH1F("fEPMassUpsilon","mass upsilon",nbins,mass_min,mass_max); + fEPMassUpsilon->SetDirectory(0); for (ibin=1; ibin<=nbins; ibin++ ) { @@ -360,6 +380,8 @@ void ExodusDecayer::Init() weight_omega = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass); weight_phi = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass); weight_jpsi = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi); + weight_psi2s = (Float_t)Lorentz(mass_bin,vmass_psi2s,vwidth_psi2s); + weight_upsilon = (Float_t)Lorentz(mass_bin,vmass_upsilon,vwidth_upsilon); // Fill histograms of electron pair masses from resonance decays fEPMassEta ->AddBinContent(ibin,weight_eta); @@ -367,6 +389,8 @@ void ExodusDecayer::Init() fEPMassOmega->AddBinContent(ibin,weight_omega); fEPMassPhi ->AddBinContent(ibin,weight_phi); fEPMassJPsi ->AddBinContent(ibin,weight_jpsi); + fEPMassPsi2S ->AddBinContent(ibin,weight_psi2s); + fEPMassUpsilon ->AddBinContent(ibin,weight_upsilon); } fInit=1; } @@ -478,7 +502,11 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) // -------- get id of the partner from idpart ------- // Int_t idpartner=idpart/1000; - idpart=idpart%1000; + if (idpart!=100443){ //protect psi(2S) pdg code (100443) + idpart=idpart%1000; + } else { + idpartner=0; + } //local variables for dalitz/2-body decay: Double_t pmass, epmass, realp_mass, e1, p1, e3, p3; @@ -492,6 +520,8 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) Int_t idPi0=111; Int_t idEta=221; Int_t idEtaPrime=331; + Int_t idPsi2S=100443; + Int_t idUpsilon=553; // Get the particle masses of daughters Double_t emass, proton_mass, omass_pion, omass_eta, omass_gamma, omass_omega; @@ -651,11 +681,9 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) //-----------------------------------------------------------------------------// -// Generate 2-body resonance decays: Rho/Omega/Phi/JPsi // +// Generate 2-body resonance decays: Rho/Omega/Phi/JPsi/Psi2S/Upsilon // //-----------------------------------------------------------------------------// - - if(((idpart==idEta&&fDecayToDimuon==1)||idpart==idRho||idpart==idOmega||idpart==idPhi||idpart==idJPsi)&&(idpartner==0)){ - + if(((idpart==idEta&&fDecayToDimuon==1)||idpart==idRho||idpart==idOmega||idpart==idPhi||idpart==idJPsi||idpart==idPsi2S||idpart==idUpsilon)&&(idpartner==0)){ //get the parent mass mp_res = pparent->M(); @@ -666,7 +694,6 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) printf("ExodusDecayer: res into ee Decay kinematically impossible! \n"); return; } - // Sample the electron pair mass from a histogram and set Polarization for( ;; ) { if(idpart==idEta&&fDecayToDimuon==1){ @@ -687,6 +714,12 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) }else if(idpart==idJPsi){ epmass_res = fEPMassJPsi->GetRandom(); PolPar=0.; + }else if(idpart==idPsi2S){ + epmass_res = fEPMassPsi2S->GetRandom(); + PolPar=0.; + } else if(idpart==idUpsilon){ + epmass_res = fEPMassUpsilon->GetRandom(); + PolPar=0.; }else{ printf("ExodusDecayer: ERROR: Resonance mass G-S parametrization not found \n"); return; } @@ -757,6 +790,12 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent) }else if(idpart==idJPsi){ fProducts_jpsi[0]=fProducts_res[0]; fProducts_jpsi[1]=fProducts_res[1]; + }else if(idpart==idPsi2S){ + fProducts_psi2s[0]=fProducts_res[0]; + fProducts_psi2s[1]=fProducts_res[1]; + }else if(idpart==idUpsilon){ + fProducts_upsilon[0]=fProducts_res[0]; + fProducts_upsilon[1]=fProducts_res[1]; } } diff --git a/GeneratorParam/ExodusDecayer.h b/GeneratorParam/ExodusDecayer.h index 6329ab5..ac4ee66 100644 --- a/GeneratorParam/ExodusDecayer.h +++ b/GeneratorParam/ExodusDecayer.h @@ -63,6 +63,8 @@ class ExodusDecayer : public TVirtualMCDecayer virtual TH1F* ElectronPairMassHistoPhiDalitz() {return fEPMassPhiDalitz;} virtual TH1F* ElectronPairMassHistoPhiDalitz_toPi0() {return fEPMassPhiDalitz_toPi0;} virtual TH1F* ElectronPairMassHistoJPsi() {return fEPMassJPsi;} + virtual TH1F* ElectronPairMassHistoPsi2S() {return fEPMassPsi2S;} + virtual TH1F* ElectronPairMassHistoUpsilon() {return fEPMassUpsilon;} virtual const TLorentzVector* Products_pion() const {return fProducts_pion;} virtual const TLorentzVector* Products_eta() const {return fProducts_eta;} @@ -76,6 +78,8 @@ class ExodusDecayer : public TVirtualMCDecayer virtual const TLorentzVector* Products_phi_dalitz() const {return fProducts_phi_dalitz;} virtual const TLorentzVector* Products_phi_dalitz_toPi0() const {return fProducts_phi_dalitz_toPi0;} virtual const TLorentzVector* Products_jpsi() const {return fProducts_jpsi;} + virtual const TLorentzVector* Products_psi2s() const {return fProducts_psi2s;} + virtual const TLorentzVector* Products_upsilon() const {return fProducts_upsilon;} protected: // Histograms for electron pair mass @@ -91,6 +95,8 @@ class ExodusDecayer : public TVirtualMCDecayer TH1F* fEPMassPhiDalitz; TH1F* fEPMassPhiDalitz_toPi0; TH1F* fEPMassJPsi; + TH1F* fEPMassPsi2S; + TH1F* fEPMassUpsilon; TF1* fPol; @@ -107,6 +113,8 @@ class ExodusDecayer : public TVirtualMCDecayer TLorentzVector fProducts_phi_dalitz[3]; TLorentzVector fProducts_phi_dalitz_toPi0[3]; TLorentzVector fProducts_jpsi[2]; + TLorentzVector fProducts_psi2s[2]; + TLorentzVector fProducts_upsilon[2]; Bool_t fInit; diff --git a/GeneratorParam/GeneratorParamEMlibV2.cxx b/GeneratorParam/GeneratorParamEMlibV2.cxx index 588aa8a..19a019f 100644 --- a/GeneratorParam/GeneratorParamEMlibV2.cxx +++ b/GeneratorParam/GeneratorParamEMlibV2.cxx @@ -108,10 +108,10 @@ const Double_t GeneratorParamEMlibV2::fgkThermPtParam[kCentralities][2] = { ,{ 8.088291e-01, 2.013231e+00 } // 20-50 }; -// MASS 0=>PIZERO, 1=>ETA, 2=>RHO0, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI, 7=>SIGMA, 8=>K0s, 9=>DELTA++, 10=>DELTA+, 11=>DELTA-, 12=>DELTA0, 13=>Rho+, 14=>Rho-, 15=>K0*, 16=>K0l, 17=>Lambda, 18=>K+, 19=>K-, 20=>Omega+, 21=>Omega-, 22=>Xi+, 23=>Xi-, 24=>Sigma+, 25=>Sigma- -const Double_t GeneratorParamEMlibV2::fgkHM[26] = {0.1349766, 0.547853, 0.77549, 0.78265, 0.95778, 1.019455, 3.096916, 1.192642, 0.497614, 1.2311, 1.2349, 1.2349, 1.23340, 0.77549, 0.77549, 0.896, 0.497614, 1.115683, 0.493677, 0.493677, 1.67245, 1.67245, 1.32171, 1.32171, 1.3828, 1.3872}; +// MASS PIZERO, ETA, RHO0, OMEGA, ETAPRIME, PHI, JPSI, PSI(2S) UPSILON SIGMA, K0s, DELTA++, DELTA+, DELTA-, DELTA0, Rho+, Rho-, K0*, K0l, Lambda, K+, K-, Omega+, Omega-, Xi+, Xi-, Sigma+, Sigma- +const Double_t GeneratorParamEMlibV2::fgkHM[kNHadrons] = {0.1349768, 0.547862, 0.77526, 0.78266, 0.95778, 1.019461, 3.096900, 3.686097, 9.46040, 1.192642, 0.497614, 1.2311, 1.2349, 1.2349, 1.23340, 0.77549, 0.77549, 0.896, 0.497614, 1.115683, 0.493677, 0.493677, 1.67245, 1.67245, 1.32171, 1.32171, 1.3828, 1.3872}; -const Double_t GeneratorParamEMlibV2::fgkMtFactor[3][26] = { +const Double_t GeneratorParamEMlibV2::fgkMtFactor[3][kNHadrons] = { // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054}, // factor for pp from arXiv:1110.3929 // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004} // factor for PbPb from arXiv:1110.3929 //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values) @@ -125,9 +125,9 @@ const Double_t GeneratorParamEMlibV2::fgkMtFactor[3][26] = { /*best guess: - pp values for eta/pi0 [arXiv:1205.5724], omega/pi0 [arXiv:1210.5749], phi/(pi+/-) [arXiv:1208.5717], K+-/pi+- [arXiv:1504.00024v2] from measured 7 Tev data */ - {1., 0.476, 1.0, 0.85, 0.4, 0.13, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, //pp - {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, //pPb - {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0} //PbPb + {1., 0.476, 1.0, 0.85, 0.4, 0.13, 1.,1., 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, //pp + {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1.,1., 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}, //pPb + {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1.,1., 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0, 0.575, 0.18, 0.41, 0.41, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0} //PbPb }; // Exponential @@ -529,6 +529,72 @@ Double_t GeneratorParamEMlibV2::V2Jpsi( const Double_t *px, const Double_t */*du return 0; } +//-------------------------------------------------------------------------- +// +// Psi2S +// +//-------------------------------------------------------------------------- +Int_t GeneratorParamEMlibV2::IpPsi2S(TRandom *) +{ + // Return Psi2S pdg code + return 100443; +} + +Double_t GeneratorParamEMlibV2::PtPsi2S( const Double_t *px, const Double_t */*dummy*/ ) +{ + const double &pt=px[0]; + return fPtParametrization[kPsi2S]->Eval(pt); +} + +Double_t GeneratorParamEMlibV2::YPsi2S( const Double_t *py, const Double_t */*dummy*/ ) +{ + return YFlat(*py); +} + +Double_t GeneratorParamEMlibV2::V2Psi2S( const Double_t *px, const Double_t */*dummy*/ ) +{ + //If there are parameterizations read from file, use them + if(fV2Parametrization[kPsi2S]) + return fV2Parametrization[kPsi2S]->Eval(EtScalingV2(px[0], kPsi2S,fV2RefParameterization[kPsi2S])) ; + + //else use build-in parameterizations + return KEtScal(*px,kPsi2S); +} + + + +//-------------------------------------------------------------------------- +// +// Upsilon +// +//-------------------------------------------------------------------------- +Int_t GeneratorParamEMlibV2::IpUpsilon(TRandom *) +{ + // Return Upsilon pdg code + return 553; +} + +Double_t GeneratorParamEMlibV2::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ ) +{ + const double &pt=px[0]; + return fPtParametrization[kUpsilon]->Eval(pt); +} + +Double_t GeneratorParamEMlibV2::YUpsilon( const Double_t *py, const Double_t */*dummy*/ ) +{ + return YFlat(*py); +} + +Double_t GeneratorParamEMlibV2::V2Upsilon( const Double_t *px, const Double_t */*dummy*/ ) +{ + //If there are parameterizations read from file, use them + if(fV2Parametrization[kUpsilon]) + return fV2Parametrization[kUpsilon]->Eval(EtScalingV2(px[0], kUpsilon,fV2RefParameterization[kUpsilon])) ; + + //else use build-in parameterizations + return KEtScal(*px,kUpsilon); +} + //-------------------------------------------------------------------------- // @@ -1352,7 +1418,7 @@ Bool_t GeneratorParamEMlibV2::SetPtParametrizations(TString fileName, TString di TRandom* rndm=NULL; // get parametrizations from file - for (Int_t i=1; i<26; i++) { + for (Int_t i=1; iGet(Form("%d_pt", ip)); @@ -1364,7 +1430,7 @@ Bool_t GeneratorParamEMlibV2::SetPtParametrizations(TString fileName, TString di } fPtParametrization[i]->SetName(Form("%d_pt", ip)); } else { - if (i==7 || i==9 || i==10 || i==11 || i==12 || i==17 || (i>=20 && i<=25)) + if (i==kSigma0 || i==kDeltaPlPl || i==kDeltaPl || i==kDeltaZero || i==kDeltaMi || i==kLambda || i==kOmegaPl || i==kOmegaMi || i==kXiPl || i==kXiMi || i==kSigmaPl || i==kSigmaMi) fPtParametrization[i] = (TF1*)MtScal(i, Form("%d_pt_mtScaled", ip), 0); else fPtParametrization[i] = (TF1*)MtScal(i, Form("%d_pt_mtScaled", ip), 1); @@ -1416,7 +1482,7 @@ Bool_t GeneratorParamEMlibV2::SetPtParametrizations(TString fileName, TString di TRandom* rndm=NULL; // get parametrizations from file - for (Int_t i=1; i<26; i++) { + for (Int_t i=1; i=20 && i<=25)) + if (i==kSigma0 || i==kDeltaPlPl || i==kDeltaPl || i==kDeltaZero || i==kDeltaMi || i==kLambda || i==kOmegaPl || i==kOmegaMi || i==kXiPl || i==kXiMi || i==kSigmaPl || i==kSigmaMi) fPtParametrization[i] = (TF1*)MtScal(i, Form("%d_pt_mtScaled", ip), 0); else fPtParametrization[i] = (TF1*)MtScal(i, Form("%d_pt_mtScaled", ip), 1); @@ -1566,9 +1632,9 @@ Bool_t GeneratorParamEMlibV2::SetFlowParametrizations(TString fileName, TString // //-------------------------------------------------------------------------- TF1* GeneratorParamEMlibV2::GetPtParametrization(Int_t np) { - if (np<26) + if (npGetYaxis()->SetTitle("mt scaling factor"); - fMtFactorHisto->GetXaxis()->SetBinLabel(1,"111"); - fMtFactorHisto->GetXaxis()->SetBinLabel(2,"221"); - fMtFactorHisto->GetXaxis()->SetBinLabel(3,"113"); - fMtFactorHisto->GetXaxis()->SetBinLabel(4,"223"); - fMtFactorHisto->GetXaxis()->SetBinLabel(5,"331"); - fMtFactorHisto->GetXaxis()->SetBinLabel(6,"333"); - fMtFactorHisto->GetXaxis()->SetBinLabel(7,"443"); - fMtFactorHisto->GetXaxis()->SetBinLabel(8,"3212"); - fMtFactorHisto->GetXaxis()->SetBinLabel(9,"310"); - fMtFactorHisto->GetXaxis()->SetBinLabel(10,"2224"); - fMtFactorHisto->GetXaxis()->SetBinLabel(11,"2214"); - fMtFactorHisto->GetXaxis()->SetBinLabel(12,"1114"); - fMtFactorHisto->GetXaxis()->SetBinLabel(13,"2114"); - fMtFactorHisto->GetXaxis()->SetBinLabel(14,"213"); - fMtFactorHisto->GetXaxis()->SetBinLabel(15,"-213"); - fMtFactorHisto->GetXaxis()->SetBinLabel(16,"313"); - fMtFactorHisto->GetXaxis()->SetBinLabel(17,"130"); - fMtFactorHisto->GetXaxis()->SetBinLabel(18,"3122"); - fMtFactorHisto->GetXaxis()->SetBinLabel(19,"321"); - fMtFactorHisto->GetXaxis()->SetBinLabel(20,"-321"); - fMtFactorHisto->GetXaxis()->SetBinLabel(21,"-3334"); - fMtFactorHisto->GetXaxis()->SetBinLabel(22,"3334"); - fMtFactorHisto->GetXaxis()->SetBinLabel(23,"-3312"); - fMtFactorHisto->GetXaxis()->SetBinLabel(24,"3312"); - fMtFactorHisto->GetXaxis()->SetBinLabel(25,"3224"); - fMtFactorHisto->GetXaxis()->SetBinLabel(26,"3114"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kPizero+1,"111"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kEta+1,"221"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kRho0+1,"113"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kOmega+1,"223"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kEtaprime+1,"331"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kPhi+1,"333"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kJpsi+1,"443"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kPsi2S+1,"100443"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kUpsilon+1,"553"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kSigma0+1,"3212"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kK0s+1,"310"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kDeltaPlPl+1,"2224"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kDeltaPl+1,"2214"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kDeltaMi+1,"1114"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kDeltaZero+1,"2114"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kRhoPl+1,"213"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kRhoMi+1,"-213"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kK0star+1,"313"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kK0l+1,"130"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kLambda+1,"3122"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kKPl+1,"321"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kKMi+1,"-321"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kOmegaPl+1,"-3334"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kOmegaMi+1,"3334"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kXiPl+1,"-3312"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kXiMi+1,"3312"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kSigmaPl+1,"3224"); + fMtFactorHisto->GetXaxis()->SetBinLabel(kSigmaMi+1,"3114"); fMtFactorHisto->SetDirectory(0); if (fileName.EndsWith(".root")){ //read the old ROOT files @@ -1648,7 +1716,7 @@ void GeneratorParamEMlibV2::SetMtScalingFactors(TString fileName, TString dirNam if (fMtFactorHistoTemp) { GeneratorParamEMlibV2 lib; TRandom* rndm=NULL; - for (Int_t i=0; i<26; i++) { + for (Int_t i=0; iGetNbinsX()+1; j++) { @@ -1671,7 +1739,7 @@ void GeneratorParamEMlibV2::SetMtScalingFactors(TString fileName, TString dirNam delete fMtFactorFile; } else{ // read the JSON file - for (Int_t i=0; i<26; i++) { + for (Int_t i=0; iSetBinContent(i+1, fgkMtFactor[selectedCol][i]); //first set all factors to hard coded value } @@ -1683,7 +1751,7 @@ void GeneratorParamEMlibV2::SetMtScalingFactors(TString fileName, TString dirNam if (dir.contains("histoMtScaleFactor")){ GeneratorParamEMlibV2 lib; TRandom* rndm=NULL; - for (Int_t i=0; i<26; i++) { + for (Int_t i=0; iGet(Form("%d_pt_y", ip)); if (ptYTemp) { @@ -1771,7 +1839,7 @@ Bool_t GeneratorParamEMlibV2::SetPtYDistributions(TString fileName, TString dirN // check for pt-y parametrizations GeneratorParamEMlibV2 lib; TRandom* rndm=NULL; - for (Int_t i=0; i<26; i++) { + for (Int_t i=0; iPy1ent(0, idpart, energy, theta, phi); JPsiDirect(); } + else if(idpart == 100443){ + fPythia->Py1ent(0, idpart, energy, theta, phi); + Psi2SDirect(); + } + else if(idpart == 553){ + fPythia->Py1ent(0, idpart, energy, theta, phi); + UpsilonDirect(); + } } fPythia->SetMSTU(10,2); fPythia->GetPrimaries(); @@ -1234,6 +1244,56 @@ void PythiaDecayerConfig::JPsiDirect(){ } } +void PythiaDecayerConfig::Psi2SDirect(){ + Int_t nt = fPythia->GetN(); + for (Int_t i = 0; i < nt; i++) { + if (fPythia->GetK(i+1,2) != 100443) continue; + Int_t fd = fPythia->GetK(i+1,4) - 1; + Int_t ld = fPythia->GetK(i+1,5) - 1; + if (fd < 0) continue; + if ((ld - fd) != 1) continue; + if (fDecayToDimuon == 0) { + if ((TMath::Abs(fPythia->GetK(fd+1,2)) != 11)) continue; + } else { + if ((TMath::Abs(fPythia->GetK(fd+1,2)) != 13)) continue; + } + TLorentzVector psi2s(fPythia->GetP(i+1,1), fPythia->GetP(i+1,2), fPythia->GetP(i+1,3), fPythia->GetP(i+1,4)); + Int_t pdg = TMath::Abs(fPythia->GetK(i+1,2)); + fDecayerExodus->Decay(pdg, &psi2s); + for (Int_t j = 0; j < 2; j++) { + for (Int_t k = 0; k < 4; k++) { + TLorentzVector vec = (fDecayerExodus->Products_psi2s())[1-j]; + fPythia->SetP(fd+j+1,k+1,vec[k]); + } + } + } +} + +void PythiaDecayerConfig::UpsilonDirect(){ + Int_t nt = fPythia->GetN(); + for (Int_t i = 0; i < nt; i++) { + if (fPythia->GetK(i+1,2) != 553) continue; + Int_t fd = fPythia->GetK(i+1,4) - 1; + Int_t ld = fPythia->GetK(i+1,5) - 1; + if (fd < 0) continue; + if ((ld - fd) != 1) continue; + if (fDecayToDimuon == 0) { + if ((TMath::Abs(fPythia->GetK(fd+1,2)) != 11)) continue; + } else { + if ((TMath::Abs(fPythia->GetK(fd+1,2)) != 13)) continue; + } + TLorentzVector upsilon(fPythia->GetP(i+1,1), fPythia->GetP(i+1,2), fPythia->GetP(i+1,3), fPythia->GetP(i+1,4)); + Int_t pdg = TMath::Abs(fPythia->GetK(i+1,2)); + fDecayerExodus->Decay(pdg, &upsilon); + for (Int_t j = 0; j < 2; j++) { + for (Int_t k = 0; k < 4; k++) { + TLorentzVector vec = (fDecayerExodus->Products_upsilon())[1-j]; + fPythia->SetP(fd+j+1,k+1,vec[k]); + } + } + } +} + void PythiaDecayerConfig::Copy(TObject &) const { // // Copy *this onto PythiaDecayerConfig -- not implemented diff --git a/GeneratorParam/PythiaDecayerConfig.h b/GeneratorParam/PythiaDecayerConfig.h index 726602c..7aa1809 100644 --- a/GeneratorParam/PythiaDecayerConfig.h +++ b/GeneratorParam/PythiaDecayerConfig.h @@ -118,6 +118,8 @@ class PythiaDecayerConfig : public TVirtualMCDecayer { void PhiDalitz(); void PhiDirect(); void JPsiDirect(); + void Psi2SDirect(); + void UpsilonDirect(); private: Int_t CountProducts(Int_t channel, Int_t particle);