From d971512bde5ec8c6bbebc7ecf25b3afe5e86af98 Mon Sep 17 00:00:00 2001 From: ffionda Date: Thu, 19 Oct 2023 09:06:02 +0200 Subject: [PATCH] update prompt / non-prompt charmonia simulation for introducing trigger gap --- .../generator/GeneratorPromptCharmonia.C | 460 ++++++++++++++++++ ...pythia8_NonPromptSignals_gaptriggered_dq.C | 276 +++++++++++ ...ithInjectedPromptSignals_gaptriggered_dq.C | 92 ++++ ...neratorHF_bbbarToBplus_midy_triggerGap.ini | 20 + ...torHF_bbbar_PsiAndJpsi_midy_triggerGap.ini | 20 + ...InjectedPromptCharmoniaMidy_TriggerGap.ini | 7 + ...GeneratorHF_bbbarToBplus_midy_triggerGap.C | 104 ++++ ...ratorHF_bbbar_PsiAndJpsi_midy_triggerGap.C | 92 ++++ ...r_InjectedPromptCharmoniaMidy_TriggerGap.C | 84 ++++ .../generator/pythia8_inel_triggerGap.cfg | 11 + ...unBeautyToPsiAndJpsi_midy_pp_triggerGap.sh | 21 + .../runBplusToJpsi_midy_pp_triggerGap.sh | 21 + .../runPromptCharmonia_midy_pp_triggerGap.sh | 21 + 13 files changed, 1229 insertions(+) create mode 100644 MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C create mode 100644 MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C create mode 100644 MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C create mode 100644 MC/config/PWGDQ/ini/GeneratorHF_bbbarToBplus_midy_triggerGap.ini create mode 100644 MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.ini create mode 100755 MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaMidy_TriggerGap.ini create mode 100644 MC/config/PWGDQ/ini/tests/GeneratorHF_bbbarToBplus_midy_triggerGap.C create mode 100644 MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.C create mode 100644 MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaMidy_TriggerGap.C create mode 100644 MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg create mode 100755 MC/run/PWGDQ/runBeautyToPsiAndJpsi_midy_pp_triggerGap.sh create mode 100755 MC/run/PWGDQ/runBplusToJpsi_midy_pp_triggerGap.sh create mode 100755 MC/run/PWGDQ/runPromptCharmonia_midy_pp_triggerGap.sh diff --git a/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C new file mode 100644 index 000000000..c92f803dd --- /dev/null +++ b/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C @@ -0,0 +1,460 @@ +// +// generators for prompt charmonia considering different cases (prompt jpsi, prompt psi2S, prompt jpsi+psi2S) at midrapidity and forward rapidity +// +// usage: +// Jpsi+Psi2S midy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp13TeV()" +// Jpsi midy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorParamPromptJpsiToElectronEvtGen_pp13TeV()" +// Psi2S midy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorParamPromptPSiToElectronEvtGen_pp13TeV()" +// Jpsi+Psi2S fwdy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV()" +// Jpsi fwdy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV()" +// Psi2S fwdy: o2-sim -j 4 -n 10 -g external -o sgn --configKeyValues "GeneratorExternal.fileName=$O2DPG_ROOT/MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C;GeneratorExternal.funcName=GeneratorParamPromptPSiToMuonEvtGen_pp13TeV()" +// + +R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen) +#include "GeneratorCocktail.C" +#include "GeneratorEvtGen.C" + +namespace o2 +{ +namespace eventgen +{ + +class O2_GeneratorParamJpsiMidY : public GeneratorTGenerator +{ + + public: + O2_GeneratorParamJpsiMidY() : GeneratorTGenerator("ParamJpsiMidY") + { + paramJpsi = new GeneratorParam(1, -1, PtJPsipp13TeV, YJPsipp13TeV, V2JPsipp13TeV, IpJPsipp13TeV); + paramJpsi->SetMomentumRange(0., 1.e6); + paramJpsi->SetPtRange(0., 1000.); + paramJpsi->SetYRange(-1.0, 1.0); + paramJpsi->SetPhiRange(0., 360.); + paramJpsi->SetDecayer(new TPythia6Decayer()); // Pythia + paramJpsi->SetForceDecay(kNoDecay); // particle left undecayed + setTGenerator(paramJpsi); + }; + + ~O2_GeneratorParamJpsiMidY() + { + delete paramJpsi; + }; + + Bool_t Init() override + { + GeneratorTGenerator::Init(); + paramJpsi->Init(); + return true; + } + + void SetNSignalPerEvent(Int_t nsig) { paramJpsi->SetNumberParticles(nsig); } + + //-------------------------------------------------------------------------// + static Double_t PtJPsipp13TeV(const Double_t* px, const Double_t* /*dummy*/) + { + // prompt J/Psi pT + // pp, 13TeV (tuned on pp 13 TeV, 2016-2018) + // + const Double_t kC = 2.28550e+00; + const Double_t kpt0 = 3.73619e+00; + const Double_t kn = 2.81708e+00; + Double_t pt = px[0]; + + return kC * pt / TMath::Power((1. + (pt / kpt0) * (pt / kpt0)), kn); + } + + //-------------------------------------------------------------------------// + static Double_t YJPsipp13TeV(const Double_t* py, const Double_t* /*dummy*/) + { + // jpsi y in pp at 13 TeV, tuned on data, prompt jpsi ALICE+LHCb, 13 TeV + Double_t y = *py; + Float_t p0, p1, p2; + p0 = 7.79382e+00; + p1 = 2.87827e-06; + p2 = 4.41847e+00; + return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2)); + } + + //-------------------------------------------------------------------------// + static Double_t V2JPsipp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/) + { + // jpsi v2 + return 0.; + } + + //-------------------------------------------------------------------------// + static Int_t IpJPsipp13TeV(TRandom*) + { + return 443; + } + + private: + GeneratorParam* paramJpsi = nullptr; +}; + +class O2_GeneratorParamPsiMidY : public GeneratorTGenerator +{ + + public: + O2_GeneratorParamPsiMidY() : GeneratorTGenerator("ParamPsi") + { + paramPsi = new GeneratorParam(1, -1, PtPsipp13TeV, YPsipp13TeV, V2Psipp13TeV, IpPsipp13TeV); + paramPsi->SetMomentumRange(0., 1.e6); // Momentum range added from me + paramPsi->SetPtRange(0., 1000.); // transverse of momentum range + paramPsi->SetYRange(-1.0, 1.0); // rapidity range + paramPsi->SetPhiRange(0., 360.); // phi range + paramPsi->SetDecayer(new TPythia6Decayer()); // Pythia decayer + paramPsi->SetForceDecay(kNoDecay); // particle left undecayed + setTGenerator(paramPsi); // Setting parameters to ParamPsi for Psi(2S) + }; + + ~O2_GeneratorParamPsiMidY() + { + delete paramPsi; + }; + + Bool_t Init() override + { + GeneratorTGenerator::Init(); + paramPsi->Init(); + return true; + } + void SetNSignalPerEvent(Int_t nsig) { paramPsi->SetNumberParticles(nsig); } + + //-------------------------------------------------------------------------// + static Double_t PtPsipp13TeV(const Double_t* px, const Double_t* /*dummy*/) + { + // prompt J/Psi pT + // pp, 13TeV (tuned on pp 13 TeV, 2016-2018) + // + const Double_t kC = 2.28550e+00; + const Double_t kpt0 = 3.73619e+00; + const Double_t kn = 2.81708e+00; + Double_t pt = px[0]; + + return kC * pt / TMath::Power((1. + (pt / kpt0) * (pt / kpt0)), kn); + } + + //-------------------------------------------------------------------------// + static Double_t YPsipp13TeV(const Double_t* py, const Double_t* /*dummy*/) + { + // jpsi y in pp at 13 TeV, tuned on data, prompt jpsi ALICE+LHCb, 13 TeV + Double_t y = *py; + Float_t p0, p1, p2; + p0 = 7.79382e+00; + p1 = 2.87827e-06; + p2 = 4.41847e+00; + return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2)); + } + + //-------------------------------------------------------------------------// + static Double_t V2Psipp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/) + { + // jpsi v2 + return 0.; + } + + //-------------------------------------------------------------------------// + static Int_t IpPsipp13TeV(TRandom*) + { + return 100443; + } + + private: + GeneratorParam* paramPsi = nullptr; +}; + +class O2_GeneratorParamJpsiFwdY : public GeneratorTGenerator +{ + + public: + O2_GeneratorParamJpsiFwdY() : GeneratorTGenerator("ParamJpsi") + { + paramJpsi = new GeneratorParam(1, -1, PtJPsipp13TeV, YJPsipp13TeV, V2JPsipp13TeV, IpJPsipp13TeV); + paramJpsi->SetMomentumRange(0., 1.e6); + paramJpsi->SetPtRange(0, 999.); + paramJpsi->SetYRange(-4.2, -2.3); + paramJpsi->SetPhiRange(0., 360.); + paramJpsi->SetDecayer(new TPythia6Decayer()); + paramJpsi->SetForceDecay(kNoDecay); // particle left undecayed + // - - paramJpsi->SetTrackingFlag(1); // (from AliGenParam) -> check this + setTGenerator(paramJpsi); + }; + + ~O2_GeneratorParamJpsiFwdY() + { + delete paramJpsi; + }; + + Bool_t Init() override + { + GeneratorTGenerator::Init(); + paramJpsi->Init(); + return true; + } + + void SetNSignalPerEvent(Int_t nsig) { paramJpsi->SetNumberParticles(nsig); } + + //-------------------------------------------------------------------------// + static Double_t PtJPsipp13TeV(const Double_t* px, const Double_t* /*dummy*/) + { + // jpsi pT in pp at 13 TeV, tuned on data (2015) + Double_t x = *px; + Float_t p0, p1, p2, p3; + p0 = 1; + p1 = 4.75208; + p2 = 1.69247; + p3 = 4.49224; + return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3); + } + + //-------------------------------------------------------------------------// + static Double_t YJPsipp13TeV(const Double_t* py, const Double_t* /*dummy*/) + { + // jpsi y in pp at 13 TeV, tuned on data (2015) + Double_t y = *py; + Float_t p0, p1, p2; + p0 = 1; + p1 = 0; + p2 = 2.98887; + return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2)); + } + + //-------------------------------------------------------------------------// + static Double_t V2JPsipp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/) + { + // jpsi v2 + return 0.; + } + + //-------------------------------------------------------------------------// + static Int_t IpJPsipp13TeV(TRandom*) + { + return 443; + } + + private: + GeneratorParam* paramJpsi = nullptr; +}; + +class O2_GeneratorParamPsiFwdY : public GeneratorTGenerator +{ + + public: + O2_GeneratorParamPsiFwdY() : GeneratorTGenerator("ParamPsi") + { + paramPsi = new GeneratorParam(1, -1, PtPsipp13TeV, YPsipp13TeV, V2Psipp13TeV, IpPsipp13TeV); + paramPsi->SetMomentumRange(0., 1.e6); + paramPsi->SetPtRange(0, 999.); + paramPsi->SetYRange(-4.2, -2.3); + paramPsi->SetPhiRange(0., 360.); + paramPsi->SetDecayer(new TPythia6Decayer()); + paramPsi->SetForceDecay(kNoDecay); // particle left undecayed + // - - paramJpsi->SetTrackingFlag(1); // check this + setTGenerator(paramPsi); + }; + + ~O2_GeneratorParamPsiFwdY() + { + delete paramPsi; + }; + + Bool_t Init() override + { + GeneratorTGenerator::Init(); + paramPsi->Init(); + return true; + } + + void SetNSignalPerEvent(Int_t nsig) { paramPsi->SetNumberParticles(nsig); } + + //-------------------------------------------------------------------------// + static Double_t PtPsipp13TeV(const Double_t* px, const Double_t* /*dummy*/) + { + // jpsi pT in pp at 13 TeV, tuned on data (2015) + Double_t x = *px; + Float_t p0, p1, p2, p3; + p0 = 1; + p1 = 4.75208; + p2 = 1.69247; + p3 = 4.49224; + return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3); + } + + //-------------------------------------------------------------------------// + static Double_t YPsipp13TeV(const Double_t* py, const Double_t* /*dummy*/) + { + // jpsi y in pp at 13 TeV, tuned on data (2015) + Double_t y = *py; + Float_t p0, p1, p2; + p0 = 1; + p1 = 0; + p2 = 2.98887; + return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2)); + } + + //-------------------------------------------------------------------------// + static Double_t V2Psipp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/) + { + // jpsi v2 + return 0.; + } + + //-------------------------------------------------------------------------// + static Int_t IpPsipp13TeV(TRandom*) + { + return 100443; + } + + private: + GeneratorParam* paramPsi = nullptr; +}; + + +} // namespace eventgen +} // namespace o2 + +FairGenerator* + GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp13TeV() +{ + auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen(); + + auto genJpsi = new o2::eventgen::O2_GeneratorParamJpsiMidY; + genJpsi->SetNSignalPerEvent(1); // signal per event for J/Psi + auto genPsi = new o2::eventgen::O2_GeneratorParamPsiMidY; + genPsi->SetNSignalPerEvent(1); // signal per event for Psi(2s) + genCocktailEvtGen->AddGenerator(genJpsi, 1); // add cocktail --> J/Psi + genCocktailEvtGen->AddGenerator(genPsi, 1); // add cocktail --> Psi(2s) + + TString pdgs = "443;100443"; + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + genCocktailEvtGen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + genCocktailEvtGen->SetForceDecay(kEvtDiElectron); + + // print debug + genCocktailEvtGen->PrintDebug(); + + return genCocktailEvtGen; +} + +FairGenerator* + GeneratorParamPromptJpsiToElectronEvtGen_pp13TeV(TString pdgs = "443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + gen->SetForceDecay(kEvtDiElectron); + + // print debug + gen->PrintDebug(); + + return gen; +} + +FairGenerator* + GeneratorParamPromptPsiToElectronEvtGen_pp13TeV(TString pdgs = "100443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + gen->SetForceDecay(kEvtDiElectron); + + // print debug + gen->PrintDebug(); + + return gen; +} + + +FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV() +{ + + auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen(); + + auto genJpsi = new o2::eventgen::O2_GeneratorParamJpsiFwdY; + genJpsi->SetNSignalPerEvent(4); // 4 J/psi generated per event by GeneratorParam + auto genPsi = new o2::eventgen::O2_GeneratorParamPsiFwdY; + genPsi->SetNSignalPerEvent(2); // 2 Psi(2S) generated per event by GeneratorParam + genCocktailEvtGen->AddGenerator(genJpsi, 1); // 2/3 J/psi + genCocktailEvtGen->AddGenerator(genPsi, 1); // 1/3 Psi(2S) + + TString pdgs = "443;100443"; + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + genCocktailEvtGen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + genCocktailEvtGen->SetForceDecay(kEvtDiMuon); + + return genCocktailEvtGen; +} + +FairGenerator* + GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV(TString pdgs = "443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + gen->SetForceDecay(kEvtDiMuon); + + // print debug + gen->PrintDebug(); + + return gen; +} + +FairGenerator* + GeneratorParamPromptPsiToMuonEvtGen_pp13TeV(TString pdgs = "100443") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->SetNSignalPerEvent(1); // number of jpsis per event + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + } + gen->SetForceDecay(kEvtDiMuon); + + // print debug + gen->PrintDebug(); + + return gen; +} + diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C b/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C new file mode 100644 index 000000000..8c5859997 --- /dev/null +++ b/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C @@ -0,0 +1,276 @@ +#include "FairGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "Pythia8/Pythia.h" +#include "TRandom.h" + +R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen) +#include "GeneratorEvtGen.C" + +#include + +using namespace o2::eventgen; +using namespace Pythia8; + +namespace o2 +{ +namespace eventgen +{ + +class GeneratorPythia8NonPromptInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 { +public: + + /// constructor + GeneratorPythia8NonPromptInjectedGapTriggeredDQ(int inputTriggerRatio = 5) { + + mGeneratedEvents = 0; + mInverseTriggerRatio = inputTriggerRatio; + // define minimum bias event generator + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + TString pathconfigMB = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg"); + pythiaMBgen.readFile(pathconfigMB.Data()); + pythiaMBgen.readString("Random:setSeed on"); + pythiaMBgen.readString("Random:seed " + std::to_string(seed)); + mConfigMBdecays = ""; + mPDG = 5; + mRapidityMin = -1.; + mRapidityMax = 1.; + mHadronMultiplicity = -1; + mHadronRapidityMin = -1.; + mHadronRapidityMax = 1.; + mVerbose = false; + } + + /// Destructor + ~GeneratorPythia8NonPromptInjectedGapTriggeredDQ() = default; + + void setPDG(int val) { mPDG = val; }; + void addHadronPDGs(int pdg) { mHadronsPDGs.push_back(pdg); }; + void setHadronMultiplicity(int val) { mHadronMultiplicity = val; }; + void setRapidity(double valMin, double valMax) + { + mRapidityMin = valMin; + mRapidityMax = valMax; + }; + + void setRapidityHadron(double valMin, double valMax) + { + mHadronRapidityMin = valMin; + mHadronRapidityMax = valMax; + }; + + void setConfigMBdecays(TString val){mConfigMBdecays = val;} + + void setVerbose(bool val) { mVerbose = val; }; + +protected: + +Bool_t generateEvent() override + { + // reset event + bool genOk = false; + if (mGeneratedEvents % mInverseTriggerRatio == 0){ + bool ancestor = false; + while (! (genOk && ancestor) ){ + /// reset event + mPythia.event.reset(); + genOk = GeneratorPythia8::generateEvent(); + // find the q-qbar ancestor + ancestor = findHeavyQuarkPair(mPythia.event); + } + } else { + /// reset event + pythiaMBgen.event.reset(); + while (!genOk) { + genOk = pythiaMBgen.next(); + } + mPythia.event = pythiaMBgen.event; + } + mGeneratedEvents++; + if (mVerbose) mOutputEvent.list(); + return true; + } + +Bool_t Init() override + { + if(mConfigMBdecays.Contains("cfg")) pythiaMBgen.readFile(mConfigMBdecays.Data()); + GeneratorPythia8::Init(); + pythiaMBgen.init(); + return true; + } + + // search for q-qbar mother with at least one q in a selected rapidity window + bool findHeavyQuarkPair(Pythia8::Event& event) + { + int countH[mHadronsPDGs.size()]; for(int ipdg=0; ipdg < mHadronsPDGs.size(); ipdg++) countH[ipdg]=0; + bool hasq = false, hasqbar = false, atSelectedY = false, isOkAtPartonicLevel = false; + for (int ipa = 0; ipa < event.size(); ++ipa) { + + if(!isOkAtPartonicLevel){ + auto daughterList = event[ipa].daughterList(); + hasq = false; hasqbar = false; atSelectedY = false; + for (auto ida : daughterList) { + if (event[ida].id() == mPDG) + hasq = true; + if (event[ida].id() == -mPDG) + hasqbar = true; + if ((event[ida].y() > mRapidityMin) && (event[ida].y() < mRapidityMax)) + atSelectedY = true; + } + if (hasq && hasqbar && atSelectedY) isOkAtPartonicLevel = true; + } + + if( (mHadronMultiplicity <= 0) && isOkAtPartonicLevel) return true; // no selection at hadron level + + /// check at hadron level if needed + int ipdg=0; + for (auto& pdgVal : mHadronsPDGs){ + if ( (TMath::Abs(event[ipa].id()) == pdgVal) && (event[ipa].y() > mHadronRapidityMin) && (event[ipa].y() < mHadronRapidityMax) ) { countH[ipdg]++; printf("ndaughters %d \n",event[ipa].daughterList().size()); } + if(isOkAtPartonicLevel && countH[ipdg] >= mHadronMultiplicity) return true; + ipdg++; + } + } + return false; + }; + + +private: +// Interface to override import particles +Pythia8::Event mOutputEvent; + + // Control gap-triggering + unsigned long long mGeneratedEvents; + int mInverseTriggerRatio; + Pythia pythiaMBgen; // minimum bias event + TString mConfigMBdecays; + int mPDG; + std::vector mHadronsPDGs; + int mHadronMultiplicity; + double mRapidityMin; + double mRapidityMax; + double mHadronRapidityMin; + double mHadronRapidityMax; + bool mVerbose; + }; + +} + +} + +// Predefined generators: +FairGenerator* + GeneratorBeautyToJpsi_EvtGenMidY(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "511;521;531;541;5112;5122;5232;5132;5332") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->setRapidity(rapidityMin, rapidityMax); + gen->setPDG(5); + gen->setRapidityHadron(-1.5,1.5); + gen->setHadronMultiplicity(1); + TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBhadrons.cfg"); + gen->readFile(pathO2table.Data()); + gen->setConfigMBdecays(pathO2table); + gen->setVerbose(verbose); + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + gen->addHadronPDGs(std::stoi(spdg)); + printf("PDG %d \n", std::stoi(spdg)); + } + gen->SetForceDecay(kEvtBJpsiDiElectron); + + // set random seed + gen->readString("Random:setSeed on"); + uint random_seed; + unsigned long long int random_value = 0; + ifstream urandom("/dev/urandom", ios::in|ios::binary); + urandom.read(reinterpret_cast(&random_value), sizeof(random_seed)); + gen->readString(Form("Random:seed = %d", random_value % 900000001)); + + // print debug + // gen->PrintDebug(); + + return gen; +} + +FairGenerator* + GeneratorBeautyToPsiAndJpsi_EvtGenMidY(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "511;521;531;541;5112;5122;5232;5132;5332") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->setRapidity(rapidityMin, rapidityMax); + gen->setPDG(5); + gen->setRapidityHadron(rapidityMin,rapidityMax); + gen->setHadronMultiplicity(1); + TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBhadrons.cfg"); + gen->readFile(pathO2table.Data()); + gen->setConfigMBdecays(pathO2table); + gen->setVerbose(verbose); + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + printf("PDG %d \n", std::stoi(spdg)); + gen->addHadronPDGs(std::stoi(spdg)); + } + gen->SetForceDecay(kEvtBPsiAndJpsiDiElectron); + + // set random seed + gen->readString("Random:setSeed on"); + uint random_seed; + unsigned long long int random_value = 0; + ifstream urandom("/dev/urandom", ios::in|ios::binary); + urandom.read(reinterpret_cast(&random_value), sizeof(random_seed)); + gen->readString(Form("Random:seed = %d", random_value % 900000001)); + + // print debug + // gen->PrintDebug(); + + return gen; +} + +FairGenerator* + GeneratorBplusToJpsiKaon_EvtGen(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "521") +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->setRapidity(rapidityMin, rapidityMax); + gen->setPDG(5); + gen->addHadronPDGs(521); + gen->setRapidityHadron(rapidityMin,rapidityMax); + gen->setHadronMultiplicity(2); + TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBplus.cfg"); + gen->readFile(pathO2table.Data()); + gen->setConfigMBdecays(pathO2table); + gen->setVerbose(verbose); + + std::string spdg; + TObjArray* obj = pdgs.Tokenize(";"); + gen->SetSizePdg(obj->GetEntriesFast()); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + spdg = obj->At(i)->GetName(); + gen->AddPdg(std::stoi(spdg), i); + gen->addHadronPDGs(std::stoi(spdg)); + printf("PDG %d \n", std::stoi(spdg)); + } + + TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen"); + gen->SetDecayTable(Form("%s/BPLUSTOKAONJPSITOELE.DEC", pathO2.Data())); + // print debug + // gen->PrintDebug(); + // set random seed + gen->readString("Random:setSeed on"); + uint random_seed; + unsigned long long int random_value = 0; + ifstream urandom("/dev/urandom", ios::in|ios::binary); + urandom.read(reinterpret_cast(&random_value), sizeof(random_seed)); + gen->readString(Form("Random:seed = %d", random_value % 900000001)); + + return gen; +} + + + diff --git a/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C new file mode 100644 index 000000000..5901729f4 --- /dev/null +++ b/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C @@ -0,0 +1,92 @@ +#include "FairGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "Pythia8/Pythia.h" +#include "TRandom.h" +#include "GeneratorPromptCharmonia.C" + +#include + +using namespace o2::eventgen; +using namespace Pythia8; + +class GeneratorPythia8PromptInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 { +public: + + /// default constructor + GeneratorPythia8PromptInjectedGapTriggeredDQ() = default; + + /// constructor + GeneratorPythia8PromptInjectedGapTriggeredDQ(int inputTriggerRatio = 5, int gentype = 0) { + + mGeneratedEvents = 0; + mGeneratorParam = 0x0; + mInverseTriggerRatio = inputTriggerRatio; + switch (gentype) { + case 0: // generate prompt jpsi + psi2s cocktail at midrapidity + mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp13TeV(); + break; + case 1: // generate prompt jpsi at midrapidity + mGeneratorParam = (Generator*)GeneratorParamPromptJpsiToElectronEvtGen_pp13TeV("443"); + break; + case 2: // generate prompt psi2S at midrapidity + mGeneratorParam = (Generator*)GeneratorParamPromptPsiToElectronEvtGen_pp13TeV("100443"); + break; + case 3: // generate prompt jpsi + psi2s cocktail at forward rapidity + mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV(); + break; + case 4: // generate prompt jpsi at forward rapidity + mGeneratorParam = (Generator*)GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV("443"); + break; + case 5: // generate prompt psi2S at forward rapidity + mGeneratorParam = (Generator*)GeneratorParamPromptPsiToMuonEvtGen_pp13TeV("100443"); + break; + } + mGeneratorParam->Init(); + } + + /// Destructor + ~GeneratorPythia8PromptInjectedGapTriggeredDQ() = default; + +protected: + +Bool_t importParticles() override + { + GeneratorPythia8::importParticles(); + bool genOk = false; + if (mGeneratedEvents % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack + bool genOk = false; + while (!genOk){ + genOk = (mGeneratorParam->generateEvent() && mGeneratorParam->importParticles()) ? true : false ; + } + int originalSize = mParticles.size(); + for(int ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++){ + TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart)); + if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize); + if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize); + if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize); + mParticles.push_back(part); + // encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C + } + mGeneratorParam->clearParticles(); + } + + mGeneratedEvents++; + return true; + } + + +private: + Generator* mGeneratorParam; + // Control gap-triggering + unsigned long long mGeneratedEvents; + int mInverseTriggerRatio; +}; + +// Predefined generators: +FairGenerator *GeneratorPythia8InjectedPromptCharmoniaGapTriggered(int inputTriggerRatio, int gentype) { + auto myGen = new GeneratorPythia8PromptInjectedGapTriggeredDQ(inputTriggerRatio,gentype); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + return myGen; +} diff --git a/MC/config/PWGDQ/ini/GeneratorHF_bbbarToBplus_midy_triggerGap.ini b/MC/config/PWGDQ/ini/GeneratorHF_bbbarToBplus_midy_triggerGap.ini new file mode 100644 index 000000000..69f0d4a2b --- /dev/null +++ b/MC/config/PWGDQ/ini/GeneratorHF_bbbarToBplus_midy_triggerGap.ini @@ -0,0 +1,20 @@ +### The setup uses an external event generator +### This part sets the path of the file and the function call to retrieve it + +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C +funcName = GeneratorBplusToJpsiKaon_EvtGen() + +### The external generator derives from GeneratorPythia8. +### This part configures the bits of the interface: configuration and user hooks + +[GeneratorPythia8] +config = ${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hf.cfg +hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C +hooksFuncName = pythia8_userhooks_bbbar(-1.5,1.5) + +### The setup inhibits transport of primary particles which are produce at forward rapidity. +### The settings below only transports particles in the barrel, which is currently defined as |eta| < 2 + +[Stack] +transportPrimary = barrel diff --git a/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.ini b/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.ini new file mode 100644 index 000000000..f014af9c8 --- /dev/null +++ b/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.ini @@ -0,0 +1,20 @@ +### The setup uses an external event generator +### This part sets the path of the file and the function call to retrieve it + +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C +funcName = GeneratorBeautyToPsiAndJpsi_EvtGenMidY() + +### The external generator derives from GeneratorPythia8. +### This part configures the bits of the interface: configuration and user hooks + +[GeneratorPythia8] +config = ${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hf.cfg +hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C +hooksFuncName = pythia8_userhooks_bbbar(-1.5,1.5) + +### The setup inhibits transport of primary particles which are produce at forward rapidity. +### The settings below only transports particles in the barrel, which is currently defined as |eta| < 2 + +[Stack] +transportPrimary = barrel diff --git a/MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaMidy_TriggerGap.ini b/MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaMidy_TriggerGap.ini new file mode 100755 index 000000000..f2ed68c40 --- /dev/null +++ b/MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaMidy_TriggerGap.ini @@ -0,0 +1,7 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C +funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,0) + +[GeneratorPythia8] +config=${O2DPG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg diff --git a/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbarToBplus_midy_triggerGap.C b/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbarToBplus_midy_triggerGap.C new file mode 100644 index 000000000..fe282e809 --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbarToBplus_midy_triggerGap.C @@ -0,0 +1,104 @@ +int External() +{ + int checkPdgSignal[] = {443}; + int checkPdgDecay = 11; + int kaonPdg = 321; + std::string path{"o2sim_Kine.root"}; + std::cout << "Check for\nsignal PDG " << checkPdgSignal[0] << "\ndecay PDG " << checkPdgDecay << "\n"; + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree*)file.Get("o2sim"); + std::vector* tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nLeptons{}; + int nAntileptons{}; + int nLeptonPairs{}; + int nLeptonPairsToBeDone{}; + int nSignalJpsi{}; + int nSignalKaons{}; + int nSignalPsi2S{}; + int nSignalJpsiWithinAcc{}; + int nSignalKaonsWithinAcc{}; + auto nEvents = tree->GetEntries(); + o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine); + Int_t bpdgs[] = {521}; + Int_t sizePdg = sizeof(bpdgs)/sizeof(Int_t); + Bool_t hasBeautyMoth = kFALSE; + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + for (auto& track : *tracks) { + auto pdg = track.GetPdgCode(); + auto rapidity = track.GetRapidity(); + auto idMoth = track.getMotherTrackId(); + if (pdg == checkPdgDecay) { + // count leptons + nLeptons++; + } else if(pdg == -checkPdgDecay) { + // count anti-leptons + nAntileptons++; + } else if (pdg == checkPdgSignal[0]) { + // check if mothers are beauty hadrons + hasBeautyMoth = kFALSE; + if(idMoth){ // check beauty mother + auto tdM = mcreader.getTrack(i, idMoth); + for(int i=0; iGetPdgCode()) == bpdgs[i] ) hasBeautyMoth = kTRUE; } + // check that it has 2 daughters, Jpsi + K + auto child0b = o2::mcutils::MCTrackNavigator::getDaughter0(*tdM, *tracks); + auto child1b = o2::mcutils::MCTrackNavigator::getDaughter1(*tdM, *tracks); + if (child0b != nullptr && child1b != nullptr) { + auto pdg0b = child0b->GetPdgCode(); + auto pdg1b = child1b->GetPdgCode(); + std::cout << "First and last children of parent B+ " << tdM->GetPdgCode() << " are PDG0: " << pdg0b << " PDG1: " << pdg1b << "\n"; + if(TMath::Abs(pdg0b) == kaonPdg ) { nSignalKaons++; if( std::abs(track.GetRapidity()) < 1.5) nSignalKaonsWithinAcc++; } + if(TMath::Abs(pdg1b) == kaonPdg ) { nSignalKaons++; if( std::abs(track.GetRapidity()) < 1.5) nSignalKaonsWithinAcc++; } + } + } + if(hasBeautyMoth){ + // count signal PDG + pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++; + // count signal PDG within acceptance + if( (std::abs(rapidity) < 1.5) && pdg == checkPdgSignal[0] ) nSignalJpsiWithinAcc++; + } + auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks); + auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks); + if (child0 != nullptr && child1 != nullptr) { + // check for parent-child relations + auto pdg0 = child0->GetPdgCode(); + auto pdg1 = child1->GetPdgCode(); + std::cout << "First and last children of parent " << checkPdgSignal[0] << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) { + nLeptonPairs++; + if (child0->getToBeDone() && child1->getToBeDone()) { + nLeptonPairsToBeDone++; + } + } + } + } + } + } + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#signal (jpsi <- B+): " << nSignalJpsi << "; within acceptance (|y| < 1.5): " << nSignalJpsiWithinAcc << "\n" + << "#signal (K+ <- B+): " << nSignalKaons << "; within acceptance (|y| < 1.5): " << nSignalKaonsWithinAcc << "\n" + << "#lepton pairs: " << nLeptonPairs << "\n" + << "#lepton pairs to be done: " << nLeptonPairs << "\n"; + + + if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) { + std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n"; + return 1; + } + if (nLeptonPairs != nLeptonPairsToBeDone) { + std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.C b/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.C new file mode 100644 index 000000000..8a0964910 --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.C @@ -0,0 +1,92 @@ +int External() +{ + int checkPdgSignal[] = {443,100443}; + int checkPdgDecay = 11; + std::string path{"o2sim_Kine.root"}; + std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n"; + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree*)file.Get("o2sim"); + std::vector* tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nLeptons{}; + int nAntileptons{}; + int nLeptonPairs{}; + int nLeptonPairsToBeDone{}; + int nSignalJpsi{}; + int nSignalPsi2S{}; + int nSignalJpsiWithinAcc{}; + int nSignalPsi2SWithinAcc{}; + auto nEvents = tree->GetEntries(); + o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine); + Int_t bpdgs[] = {511, 521, 531, 5112, 5122, 5232, 5132}; + Int_t sizePdg = sizeof(bpdgs)/sizeof(Int_t); + Bool_t hasBeautyMoth = kFALSE; + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + for (auto& track : *tracks) { + auto pdg = track.GetPdgCode(); + auto rapidity = track.GetRapidity(); + auto idMoth = track.getMotherTrackId(); + if (pdg == checkPdgDecay) { + // count leptons + nLeptons++; + } else if(pdg == -checkPdgDecay) { + // count anti-leptons + nAntileptons++; + } else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) { + // check if mothers are beauty hadrons + hasBeautyMoth = kFALSE; + if(idMoth){ // check beauty mother + auto tdM = mcreader.getTrack(i, idMoth); + for(int i=0; iGetPdgCode()) == bpdgs[i] ) hasBeautyMoth = kTRUE; } + } + if(hasBeautyMoth){ + // count signal PDG + pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++; + // count signal PDG within acceptance + if(std::abs(rapidity) < 1.0) { pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;} + } + auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks); + auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks); + if (child0 != nullptr && child1 != nullptr) { + // check for parent-child relations + auto pdg0 = child0->GetPdgCode(); + auto pdg1 = child1->GetPdgCode(); + std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) { + nLeptonPairs++; + if (child0->getToBeDone() && child1->getToBeDone()) { + nLeptonPairsToBeDone++; + } + } + } + } + } + } + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#signal (jpsi <- b): " << nSignalJpsi << "; within acceptance (|y| < 1): " << nSignalJpsiWithinAcc << "\n" + << "#signal (psi2S <- b): " << nSignalPsi2S << "; within acceptance (|y| < 1): " << nSignalPsi2SWithinAcc << "\n" + << "#lepton pairs: " << nLeptonPairs << "\n" + << "#lepton pairs to be done: " << nLeptonPairs << "\n"; + + + if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) { + std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n"; + return 1; + } + if (nLeptonPairs != nLeptonPairsToBeDone) { + std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaMidy_TriggerGap.C b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaMidy_TriggerGap.C new file mode 100644 index 000000000..b155ec9ac --- /dev/null +++ b/MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaMidy_TriggerGap.C @@ -0,0 +1,84 @@ +int External() +{ + int checkPdgSignal[] = {443,100443}; + int checkPdgDecay = 11; + std::string path{"o2sim_Kine.root"}; + std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n"; + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree*)file.Get("o2sim"); + std::vector* tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nLeptons{}; + int nAntileptons{}; + int nLeptonPairs{}; + int nLeptonPairsToBeDone{}; + int nSignalJpsi{}; + int nSignalPsi2S{}; + int nSignalJpsiWithinAcc{}; + int nSignalPsi2SWithinAcc{}; + auto nEvents = tree->GetEntries(); + o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine); + Bool_t isInjected = kFALSE; + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + for (auto& track : *tracks) { + auto pdg = track.GetPdgCode(); + auto rapidity = track.GetRapidity(); + auto idMoth = track.getMotherTrackId(); + if (pdg == checkPdgDecay) { + // count leptons + nLeptons++; + } else if(pdg == -checkPdgDecay) { + // count anti-leptons + nAntileptons++; + } else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) { + if(idMoth < 0){ + // count signal PDG + pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++; + // count signal PDG within acceptance + if(std::abs(rapidity) < 1.0) { pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;} + } + auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks); + auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks); + if (child0 != nullptr && child1 != nullptr) { + // check for parent-child relations + auto pdg0 = child0->GetPdgCode(); + auto pdg1 = child1->GetPdgCode(); + std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; + if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) { + nLeptonPairs++; + if (child0->getToBeDone() && child1->getToBeDone()) { + nLeptonPairsToBeDone++; + } + } + } + } + } + } + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance (|y| < 1): " << nSignalJpsiWithinAcc << "\n" + << "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance (|y| < 1): " << nSignalPsi2SWithinAcc << "\n" + << "#lepton pairs: " << nLeptonPairs << "\n" + << "#lepton pairs to be done: " << nLeptonPairs << "\n"; + + + if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) { + std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n"; + return 1; + } + if (nLeptonPairs != nLeptonPairsToBeDone) { + std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg b/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg new file mode 100644 index 000000000..a4c535826 --- /dev/null +++ b/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg @@ -0,0 +1,11 @@ +### beams +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +SoftQCD:inelastic on # all inelastic processes + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. diff --git a/MC/run/PWGDQ/runBeautyToPsiAndJpsi_midy_pp_triggerGap.sh b/MC/run/PWGDQ/runBeautyToPsiAndJpsi_midy_pp_triggerGap.sh new file mode 100755 index 000000000..145c29d6f --- /dev/null +++ b/MC/run/PWGDQ/runBeautyToPsiAndJpsi_midy_pp_triggerGap.sh @@ -0,0 +1,21 @@ +#!/usr/bin/env bash + +# make sure O2DPG + O2 is loaded +[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1 +[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1 + + +# ----------- LOAD UTILITY FUNCTIONS -------------------------- +. ${O2_ROOT}/share/scripts/jobutils.sh + +RNDSEED=${RNDSEED:-0} +NSIGEVENTS=${NSIGEVENTS:-1} +NBKGEVENTS=${NBKGEVENTS:-1} +NWORKERS=${NWORKERS:-8} +NTIMEFRAMES=${NTIMEFRAMES:-1} + +${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM 900 -gen external -j ${NWORKERS} -ns ${NSIGEVENTS} -tf ${NTIMEFRAMES} -e TGeant4 -mod "--skipModules ZDC" \ + -ini $O2DPG_ROOT/MC/config/PWGDQ/ini/GeneratorHF_bbbar_PsiAndJpsi_midy_triggerGap.ini --mft-reco-full + +# run workflow +${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json diff --git a/MC/run/PWGDQ/runBplusToJpsi_midy_pp_triggerGap.sh b/MC/run/PWGDQ/runBplusToJpsi_midy_pp_triggerGap.sh new file mode 100755 index 000000000..0e74920d8 --- /dev/null +++ b/MC/run/PWGDQ/runBplusToJpsi_midy_pp_triggerGap.sh @@ -0,0 +1,21 @@ +#!/usr/bin/env bash + +# make sure O2DPG + O2 is loaded +[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1 +[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1 + + +# ----------- LOAD UTILITY FUNCTIONS -------------------------- +. ${O2_ROOT}/share/scripts/jobutils.sh + +RNDSEED=${RNDSEED:-0} +NSIGEVENTS=${NSIGEVENTS:-1} +NBKGEVENTS=${NBKGEVENTS:-1} +NWORKERS=${NWORKERS:-8} +NTIMEFRAMES=${NTIMEFRAMES:-1} + +${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM 13600 -gen external -j ${NWORKERS} -ns ${NSIGEVENTS} -tf ${NTIMEFRAMES} -e TGeant4 -mod "--skipModules ZDC" \ + -ini $O2DPG_ROOT/MC/config/PWGDQ/ini/GeneratorHF_bbbarToBplus_midy_triggerGap.ini --mft-reco-full + +# run workflow +${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json diff --git a/MC/run/PWGDQ/runPromptCharmonia_midy_pp_triggerGap.sh b/MC/run/PWGDQ/runPromptCharmonia_midy_pp_triggerGap.sh new file mode 100755 index 000000000..0636e7eb6 --- /dev/null +++ b/MC/run/PWGDQ/runPromptCharmonia_midy_pp_triggerGap.sh @@ -0,0 +1,21 @@ +#!/usr/bin/env bash + +# make sure O2DPG + O2 is loaded +[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1 +[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1 + + +# ----------- LOAD UTILITY FUNCTIONS -------------------------- +. ${O2_ROOT}/share/scripts/jobutils.sh + +RNDSEED=${RNDSEED:-0} +NSIGEVENTS=${NSIGEVENTS:-1} +NBKGEVENTS=${NBKGEVENTS:-1} +NWORKERS=${NWORKERS:-8} +NTIMEFRAMES=${NTIMEFRAMES:-1} + +${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM 13600 -gen external -j ${NWORKERS} -ns ${NSIGEVENTS} -tf ${NTIMEFRAMES} -e TGeant4 -mod "--skipModules ZDC" \ + -ini $O2DPG_ROOT/MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaMidy_TriggerGap.ini --mft-reco-full + +# run workflow +${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json