diff --git a/MC/config/PWGEM/external/generator/Generator_pythia8_HFLepton_pp.C b/MC/config/PWGEM/external/generator/Generator_pythia8_HFLepton_pp.C new file mode 100644 index 000000000..854d6ff42 --- /dev/null +++ b/MC/config/PWGEM/external/generator/Generator_pythia8_HFLepton_pp.C @@ -0,0 +1,202 @@ +#include "Pythia8/Pythia.h" +#include "Pythia8/HeavyIons.h" +#include "FairGenerator.h" +#include "FairPrimaryGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "TRandom3.h" +#include "TParticlePDG.h" +#include "TDatabasePDG.h" + +#include +#include + +using namespace Pythia8; + +class GeneratorPythia8HFLeptonpp : public o2::eventgen::GeneratorPythia8 +{ +public: + /// default constructor + GeneratorPythia8HFLeptonpp() = default; + + /// constructor + GeneratorPythia8HFLeptonpp(TString configsignal, int quarkPdg = 4, int lInputExternalID = 0) + { + + lGeneratedEvents = 0; + lExternalID = lInputExternalID; + mQuarkPdg = quarkPdg; + + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + + int offset = (int)(gRandom->Uniform(1)); // create offset to mitigate edge effects due to small number of events per job + lGeneratedEvents += offset; + + cout << "Initalizing PYTHIA object used to generate signal events..." << endl; + TString pathconfigSignal = gSystem->ExpandPathName(configsignal.Data()); + pythiaObjectSignal.readFile(pathconfigSignal.Data()); + pythiaObjectSignal.readString("Random:setSeed on"); + pythiaObjectSignal.readString("Random:seed " + std::to_string(seed)); + pythiaObjectSignal.readString("Beams:eCM = 5360.0"); + pythiaObjectSignal.init(); + cout << "Initalization of signal event is complete" << endl; + + // flag the generators using type + // addCocktailConstituent(type, "interesting"); + // addCocktailConstitent(0, "minbias"); + // Add Sub generators + addSubGenerator(1, "charm lepton"); + addSubGenerator(2, "beauty forced decay"); + addSubGenerator(3, "beauty no foced decay"); + } + + /// Destructor + ~GeneratorPythia8HFLeptonpp() = default; + + void addTriggerOnDaughter(int nb, int pdg) + { + mNbDaughter = nb; + mPdgDaughter = pdg; + }; + void setQuarkRapidity(float yMin, float yMax) + { + mQuarkRapidityMin = yMin; + mQuarkRapidityMax = yMax; + }; + void setDaughterRapidity(float yMin, float yMax) + { + mDaughterRapidityMin = yMin; + mDaughterRapidityMax = yMax; + }; + +protected: + //__________________________________________________________________ + Bool_t generateEvent() override + { + /// reset event + mPythia.event.reset(); + + // Generate event of interest + Bool_t lGenerationOK = kFALSE; + while (!lGenerationOK) { + if (pythiaObjectSignal.next()) { + lGenerationOK = selectEvent(pythiaObjectSignal.event); + } + } + mPythia.event = pythiaObjectSignal.event; + notifySubGenerator(lExternalID); + + lGeneratedEvents++; + // mPythia.next(); + + return true; + } + + bool selectEvent(const Pythia8::Event& event) + { + bool isGoodAtPartonLevel = false, isGoodAtDaughterLevel = (mPdgDaughter != 0) ? false : true; + int nbDaughter = 0; + for (auto iPart{0}; iPart < event.size(); ++iPart) { + // search for Q-Qbar mother with at least one Q in rapidity window + if (!isGoodAtPartonLevel) { + auto daughterList = event[iPart].daughterList(); + bool hasQ = false, hasQbar = false, atSelectedY = false; + for (auto iDau : daughterList) { + if (event[iDau].id() == mQuarkPdg) { + hasQ = true; + } + if (event[iDau].id() == -mQuarkPdg) { + hasQbar = true; + } + if ((std::abs(event[iDau].id()) == mQuarkPdg) && (event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax)) + atSelectedY = true; + } + if (hasQ && hasQbar && atSelectedY) { + isGoodAtPartonLevel = true; + } + } + // search for mNbDaughter daughters of type mPdgDaughter in rapidity window + if (!isGoodAtDaughterLevel) { + int id = std::abs(event[iPart].id()); + float rap = event[iPart].y(); + if (id == mPdgDaughter) { + int motherindexa = event[iPart].mother1(); + if (motherindexa > 0) { + int idmother = std::abs(event[motherindexa].id()); + if (int(std::abs(idmother) / 100.) == 4 || int(std::abs(idmother) / 1000.) == 4 || int(std::abs(idmother) / 100.) == 5 || int(std::abs(idmother) / 1000.) == 5) { + if (rap > mDaughterRapidityMin && rap < mDaughterRapidityMax) { + nbDaughter++; + if (nbDaughter >= mNbDaughter) isGoodAtDaughterLevel = true; + } + } + } + } + } + // we send the trigger + if (isGoodAtPartonLevel && isGoodAtDaughterLevel) { + return true; + } + } + return false; + }; + +private: + // Interface to override import particles + Pythia8::Event mOutputEvent; + + // Properties of selection + int mQuarkPdg; + float mQuarkRapidityMin; + float mQuarkRapidityMax; + int mPdgDaughter; + int mNbDaughter; + float mDaughterRapidityMin; + float mDaughterRapidityMax; + + Long64_t lGeneratedEvents; + // ID for different generators + int lExternalID; + + // Base event generators + Pythia8::Pythia pythiaObjectSignal; ///Signal collision generator +}; + +// Predefined generators: + +// Charm-enriched forced decay +FairGenerator* GeneratorPythia8CharmLepton(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_pp_cr2_forceddecayscharm.cfg", 4, inputExternalID); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + myGen->setQuarkRapidity(yMinQ, yMaxQ); + myGen->addTriggerOnDaughter(2, pdgLepton); + myGen->setDaughterRapidity(yMinL, yMaxL); + return myGen; +} + +// Beauty-enriched forced decay +FairGenerator* GeneratorPythia8BeautyForcedDecays(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar_forceddecayscharmbeauty.cfg", 5, inputExternalID); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + myGen->setQuarkRapidity(yMinQ, yMaxQ); + myGen->addTriggerOnDaughter(2, pdgLepton); + myGen->setDaughterRapidity(yMinL, yMaxL); + return myGen; +} + +// Beauty-enriched no forced decay +FairGenerator* GeneratorPythia8BeautyNoForcedDecays(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar.cfg", 5, inputExternalID); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + myGen->setQuarkRapidity(yMinQ, yMaxQ); + myGen->addTriggerOnDaughter(2, pdgLepton); + myGen->setDaughterRapidity(yMinL, yMaxL); + return myGen; +} diff --git a/MC/config/PWGEM/external/generator/Generator_pythia8_embed_HFLepton.C b/MC/config/PWGEM/external/generator/Generator_pythia8_embed_HFLepton.C new file mode 100644 index 000000000..fe55587aa --- /dev/null +++ b/MC/config/PWGEM/external/generator/Generator_pythia8_embed_HFLepton.C @@ -0,0 +1,396 @@ +/////////////////////////////////////////////////////////////////////////////// +/// /// +/// HF MC generator for Pb-Pb /// +/// Option 1: generate N PYTHIA events triggered on ccbar and/or bbbar /// +/// to be embedded with a underlying Pb-Pb event /// +/// /// +/////////////////////////////////////////////////////////////////////////////// + +#include "Generator_pythia8_HFLepton_pp.C" + +using namespace Pythia8; + +#include +namespace hf_generators +{ + enum GenType : int { + Charm = 0, // --> GeneratorPythia8CharmLepton: charm enriched + Beauty, // --> GeneratorPythia8BeautyForcedDecays: beauty enriched with forced decays + BeautyNoForcedDecays, // --> GeneratorPythia8BeautyNoForcedDecays: beauty enriched no forced decays + NGenType + }; +} + +class GeneratorPythia8EmbedHFLepton : public o2::eventgen::GeneratorPythia8 +{ +public: + + /// default constructor + GeneratorPythia8EmbedHFLepton() = default; + + /// Destructor + ~GeneratorPythia8EmbedHFLepton() { + // Clean up the internally created HF generator if any + if (mGeneratorEvHFLepton) { + delete mGeneratorEvHFLepton; + mGeneratorEvHFLepton = nullptr; + } + } + + /// Init + bool Init() override + { + return o2::eventgen::GeneratorPythia8::Init(); + } + + /// @brief setup the event generator for HF signals + /// \param gentype generator type (only ccbar, only bbbar, both) + /// \param usePtHardBins flag to enable/disable pt-hard bins + /// \param yQuarkMin minimum quark rapidity + /// \param yQuarkMax maximum quark rapidity + /// \param yHadronMin minimum hadron rapidity + /// \param yHadronMax maximum hadron rapidity + /// \param hadronPdgList list of PDG codes for hadrons to be used in trigger + /// \param quarkPdgList list of PDG codes for quarks to be enriched in the trigger + void setupGeneratorEvHFLepton(int genType, int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) { + mGeneratorEvHFLepton = nullptr; + switch (genType) + { + case hf_generators::Charm: + LOG(info) << "********** [GeneratorPythia8EmbedHFLepton] configuring GeneratorPythia8CharmLepton **********"; + LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs; + mGeneratorEvHFLepton = dynamic_cast(GeneratorPythia8CharmLepton(inputExternalID, pdgLepton, yMinQ, yMaxQ, yMinL, yMaxL)); + break; + case hf_generators::Beauty: + LOG(info) << "********** [GeneratorPythia8EmbedHFLepton] configuring GeneratorPythia8BeautyForcedDecays **********"; + LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs; + mGeneratorEvHFLepton = dynamic_cast(GeneratorPythia8BeautyForcedDecays(inputExternalID, pdgLepton, yMinQ, yMaxQ, yMinL, yMaxL)); + break; + case hf_generators::BeautyNoForcedDecays: + LOG(info) << "********** [GeneratorPythia8EmbedHFLepton] configuring GeneratorPythia8BeautyNoForcedDecays **********"; + LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs; + mGeneratorEvHFLepton = dynamic_cast(GeneratorPythia8BeautyNoForcedDecays(inputExternalID, pdgLepton, yMinQ, yMaxQ, yMinL, yMaxL)); + break; + default: + LOG(fatal) << "********** [GeneratorPythia8EmbedHFLepton] bad configuration, fix it! **********"; + break; + } + mGeneratorEvHFLepton->Init(); + } + + // This function is called by the primary generator + // for each event in case we are in embedding mode. + // We use it to setup the number of signal events + // to be generated and to be embedded on the background. + void notifyEmbedding(const o2::dataformats::MCEventHeader* bkgHeader) override + { + LOG(info) << "[notifyEmbedding] ----- Function called"; + + /// Impact parameter between the two nuclei + const float x = bkgHeader->GetB(); + LOG(info) << "[notifyEmbedding] ----- Collision impact parameter: " << x; + + /// number of events to be embedded in a background event + // gen->setFormula("max(1.,120.*(x<5.)+80.*(1.-x/20.)*(x>5.)*(x<11.)+240.*(1.-x/13.)*(x>11.))"); + //mNumSigEvs = static_cast(std::lround(5.0 + 0.886202881 * std::pow(std::max(0.0f, 17.5f - x), 1.7))); + mNumSigEvs = static_cast(std::max(1.,120.*(x<5.)+80.*(1.-x/20.)*(x>5.)*(x<11.)+240.*(1.-x/13.)*(x>11.))); + //mNumSigEvs = 1; + LOG(info) << "[notifyEmbedding] ----- generating " << mNumSigEvs << " signal events " << std::endl; + }; + +protected: + +/// @brief Main function for event generation +bool generateEvent() override +{ + /// Overriding that from GeneratorPythia8, to avoid the simulation of an untriggered event as first + return true; +} + +/// @brief Main function to find out whether the particle comes charm or beauty quark +/// @param partId is the index of the particle under study +/// @param particles are the particles of the full event +bool isFromCharmOrBeauty(const int partId, std::vector const& particles) { + + // Let's check wheter this is already a c or b quark? + const TParticle& part = particles.at(partId); + const int pdgAbs = std::abs(part.GetPdgCode()); + if(pdgAbs == 4 || pdgAbs == 5) { + return true; + } + + // Let's check the mother particles of the hadron at all stages + // and look for the charm or beauty quark + std::vector> arrayIds{}; + std::vector initVec{partId}; + arrayIds.push_back(initVec); // the first vector contains the index of the original particle + int stage = 0; + while(arrayIds[-stage].size() > 0) { + + //LOG(info) << "### stage " << stage << ", arrayIds[-stage].size() = " << arrayIds[-stage].size(); + + std::vector arrayIdsStage{}; + + for (auto& iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage + const TParticle& partStage = particles.at(iPart); + + // check the first mother + const int firstMotherId = partStage.GetFirstMother(); + if( firstMotherId >= 0) { + const TParticle& firstMother = particles.at(firstMotherId); + const int pdgAbsFirstMother = std::abs(firstMother.GetPdgCode()); + if(pdgAbsFirstMother == 4 || pdgAbsFirstMother == 5) { + return true; + } + // the first mother is not a charm or beauty quark + arrayIdsStage.push_back(firstMotherId); + } + + // let's check all other mothers, if any + const int lastMotherId = partStage.GetSecondMother(); + if(lastMotherId >=0 && lastMotherId != firstMotherId) { + for(int motherId = firstMotherId+1 /*first mother already considered*/; motherId <= lastMotherId; motherId++) { + const TParticle& mother = particles.at(motherId); + const int pdgAbsMother = std::abs(mother.GetPdgCode()); + if(pdgAbsMother == 4 || pdgAbsMother == 5) { + return true; + } + // this mother is not a charm or beauty quark + arrayIdsStage.push_back(motherId); + } + } + + } + + /* + All light-flavour mothers are not considered with this approach + eg: D+ coming from c and uBar --> uBar lost + --> TODO: check if the current particle has a charm or beauty hadron as daughter. If yes, keep it + >>> we can ignore it! This might be useful only for jet analyses, however this approach of embedding N pp events into a Pb-Pb one might be not ideal for them + */ + + // none of the particle mothers is a charm or beauty quark, let's consider their indices for the next stage + arrayIds.push_back(arrayIdsStage); + stage--; // ready to go to next stage + + } /// end while(arrayIds[-stage].size() > 0) + + return false; +} + + +void printParticleVector(std::vector v) { + for(int id=0; id idFirstMother=" << idFirstMother << ", idLastMother=" << idLastMother << ", idFirstDaughter=" << idFirstDaughter << ", idLastDaughter=" << idLastDaughter; + } +} + + +int findKey(const std::map& m, int value) { + for(std::pair p : m) { + if(p.second /*index*/ == value) { + return p.first; // key --> it becomes the new index + } + } + return -1; +} + + +/// @brief Main function to copy the generated particles in mPythia.event into the stack (this.mParticles) +Bool_t importParticles() override +{ + /// Import particles from generated event + /// This should not do anything now, since we override generateEvent + GeneratorPythia8::importParticles(); + + LOG(info) << ""; + LOG(info) << "*************************************************************"; + LOG(info) << "************** New signal event considered **************"; + LOG(info) << "*************************************************************"; + LOG(info) << ""; + + /// Generate mNumSigEvs HF events to be merged in one + int nEvsHF = 0; + while(nEvsHF < mNumSigEvs) { + + /// generate the HF event + bool genOk = false; + while(!genOk) { + genOk = (mGeneratorEvHFLepton->generateEvent() && mGeneratorEvHFLepton->importParticles() /*copy particles from mGeneratorEvHF.mPythia.event to mGeneratorEvHF.mParticles*/ ); + } + + int originalSize = mParticles.size(); // stack of this event generator + + // for debug + // LOG(info) << ""; + // LOG(info) << "============ Before HF event " << nEvsHF; + // LOG(info) << "Full stack (size " << originalSize << "):"; + // printParticleVector(mParticles); + + /// copy the particles from the HF event in the particle stack + auto particlesHfEvent = mGeneratorEvHFLepton->getParticles(); + std::map mapHfParticles = {}; + int counterHfParticles = 0; + + // for debug + // LOG(info) << "-----------------------------------------------"; + // LOG(info) << ">>> HF event " << nEvsHF; + // LOG(info) << " HF event stack:"; + // printParticleVector(particlesHfEvent); + + for(int iPart=0; iPart>>"; + // LOG(info) << " >>> printing mapHfParticles:"; + // for(auto& p : mapHfParticles) { + // const int pdgCodeFromMap = particlesHfEvent.at(p.second).GetPdgCode(); + // LOG(info) << " >>> entry " << p.first << ", original id = " << p.second << ", pdgCode=" << pdgCodeFromMap << " --> firstMotherId=" << particlesHfEvent.at(p.second).GetFirstMother() << ", lastMotherId=" << particlesHfEvent.at(p.second).GetSecondMother() << ", firstDaughterId=" << particlesHfEvent.at(p.second).GetFirstDaughter() << ", lastDaughterId=" << particlesHfEvent.at(p.second).GetLastDaughter(); + // } + + + // In the map we have only the particles from charm or beauty + // Let's readapt the mother/daughter indices accordingly + int offset = originalSize; + for(int iHfPart=0; iHfPart= 0) { + idFirstMother = findKey(mapHfParticles, idFirstMother); + /// If idFirstMother>=0, the 1st mother is from charm or beauty, i.e. is not a light-flavoured parton + /// Instead, if idFirstMother==-1 from findKey this means that the first mother was a light-flavoured parton --> not stored in the map + if(idFirstMother >=0) { + /// the 1st mother is from charm or beauty, i.e. is not a light-flavoured parton + if(idLastMother != idFirstMotherOrig) { + if(idLastMother != -1) { + /// idLastMother is >= 0 + } + } else { + /// idLastMother is equal to idFirstMother + idLastMother = idFirstMother; + } + isFirstMotherOk = true; + } + } + if(!isFirstMotherOk) { + /// - If we are here, it means that the 1st mother was not from charm or beauty + /// - No need to check whether idLastMother>=0, + /// because this would mean that none of the mother is from charm or beauty and this was checked already in isFromCharmOrBeauty + /// - Need to loop between 1st and last mother, to treat cases like these + /// [11:52:13][INFO] id = 565, pdgCode = -2 --> idFirstMother=519, idLastMother=519 + /// [11:52:13][INFO] id = 566, pdgCode = -4 --> idFirstMother=520, idLastMother=520 + /// [11:52:13][INFO] id = 567, pdgCode = -1 --> idFirstMother=518, idLastMother=518 + /// [11:52:13][INFO] id = 568, pdgCode = -311 --> idFirstMother=565, idLastMother=567 + /// [11:52:13][INFO] id = 569, pdgCode = -4212 --> idFirstMother=565, idLastMother=567 + /// --> w/o loop between 1st and last mother, the mother Ids assigned to this Sc+ (4212) by findKey are -1, both first and last + bool foundAnyMother = false; + for(int idMotherOrig=(idFirstMotherOrig+1); idMotherOrig<=idLastMother; idMotherOrig++) { + const int idMother = findKey(mapHfParticles, idMotherOrig); + if(idMother >= 0) { + /// this should mean that the mother is from HF, i.e. that we found the correct one + idFirstMother = idMother; + idLastMother = idFirstMother; + foundAnyMother = true; + break; + } + } + // set last mother to -1 if no mother has been found so far + if (!foundAnyMother) { + idLastMother = -1; + } + } + + /// fix daughter indices + idFirstDaughter = findKey(mapHfParticles, idFirstDaughter); + idLastDaughter = findKey(mapHfParticles, idLastDaughter); + + /// adjust the particle mother and daughter indices + particle.SetFirstMother((idFirstMother >= 0) ? idFirstMother + offset : idFirstMother); + particle.SetLastMother((idLastMother >= 0) ? idLastMother + offset : idLastMother); + particle.SetFirstDaughter((idFirstDaughter >= 0) ? idFirstDaughter + offset : idFirstDaughter); + particle.SetLastDaughter((idLastDaughter >= 0) ? idLastDaughter + offset : idLastDaughter); + + /// copy inside this.mParticles from mGeneratorEvHF.mParticles, i.e. the particles generated in mGeneratorEvHF + mParticles.push_back(particle); + + } + + // for debug + // LOG(info) << "-----------------------------------------------"; + // LOG(info) << "============ After HF event " << nEvsHF; + // LOG(info) << "Full stack:"; + // printParticleVector(mParticles); + + /// one more event generated, let's update the counter and clear it, to allow the next generation + nEvsHF++; + //mGeneratedEvents++; + mGeneratorEvHFLepton->clearParticles(); + } + + return true; +} + +private: + + Generator* mGeneratorEvHFLepton; // to generate HF signal events + + int mNumSigEvs{1}; // number of HF signal events to be merged in one Pythia event + //unsigned long long mGeneratedEvents; + +}; + +// Charm enriched +FairGenerator * GeneratorPythia8EmbedHFLeptonCharm(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8EmbedHFLepton(); + + /// setup the internal generator for HF events + myGen->setupGeneratorEvHFLepton(hf_generators::Charm, inputExternalID, pdgLepton, yMinQ = -1.5, yMaxQ = 1.5, yMinL = -1, yMaxL = 1); + + return myGen; +} + +// Beauty enriched +FairGenerator * GeneratorPythia8EmbedHFLeptonBeauty(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8EmbedHFLepton(); + + /// setup the internal generator for HF events + myGen->setupGeneratorEvHFLepton(hf_generators::Beauty, inputExternalID, pdgLepton, yMinQ = -1.5, yMaxQ = 1.5, yMinL = -1, yMaxL = 1); + + return myGen; +} + +// Charm and beauty enriched (with same ratio) +FairGenerator * GeneratorPythia8EmbedHFLeptonBeautyNoForcedDecays(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1) +{ + auto myGen = new GeneratorPythia8EmbedHFLepton(); + + /// setup the internal generator for HF events + myGen->setupGeneratorEvHFLepton(hf_generators::BeautyNoForcedDecays, inputExternalID, pdgLepton, yMinQ = -1.5, yMaxQ = 1.5, yMinL = -1, yMaxL = 1); + + return myGen; +} diff --git a/MC/config/PWGEM/ini/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.ini b/MC/config/PWGEM/ini/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.ini new file mode 100644 index 000000000..f24c66d3e --- /dev/null +++ b/MC/config/PWGEM/ini/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.ini @@ -0,0 +1,9 @@ +#NEV_TEST> 10 +### The external generator derives from GeneratorPythia8 +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_embed_HFLepton.C +funcName=GeneratorPythia8EmbedHFLeptonBeautyNoForcedDecays(3, 11) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg +includePartonEvent=true diff --git a/MC/config/PWGEM/ini/GeneratorHF_Charm_PbPb_electron.ini b/MC/config/PWGEM/ini/GeneratorHF_Charm_PbPb_electron.ini new file mode 100644 index 000000000..b72cfb523 --- /dev/null +++ b/MC/config/PWGEM/ini/GeneratorHF_Charm_PbPb_electron.ini @@ -0,0 +1,9 @@ +#NEV_TEST> 10 +### The external generator derives from GeneratorPythia8 +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_embed_HFLepton.C +funcName=GeneratorPythia8EmbedHFLeptonCharm(1, 11) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg +includePartonEvent=true diff --git a/MC/config/PWGEM/ini/tests/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.C b/MC/config/PWGEM/ini/tests/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.C new file mode 100644 index 000000000..375e79b75 --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorHF_BeautyNoForcedDecay_PbPb_electron.C @@ -0,0 +1,111 @@ +int External() +{ + int checkPdgDecay = 11; + std::string path{"o2sim_Kine.root"}; + 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 nLeptonsInAcceptance{}; + int nLeptons{}; + int nAntileptons{}; + int nLeptonsToBeDone{}; + int nAntileptonsToBeDone{}; + int nSignalPairs{}; + int nLeptonPairs{}; + int nLeptonPairsToBeDone{}; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + int nElectrons = 0; + int nPositrons = 0; + int nElectronsToBeDone = 0; + int nPositronsToBeDone = 0; + int nOpenBeautyPos = 0; + int nOpenBeautyNeg = 0; + int nPositronsElectronsInAcceptance = 0; + for (auto& track : *tracks) { + auto pdg = track.GetPdgCode(); + auto y = track.GetRapidity(); + if (pdg == checkPdgDecay) { + int igmother = track.getMotherTrackId(); + if (igmother > 0) { + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 || int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) { + nLeptons++; + nElectrons++; + if (-1 < y && y < 1) nPositronsElectronsInAcceptance++; + if (track.getToBeDone()) { + nLeptonsToBeDone++; + nElectronsToBeDone++; + } + } + } + } else if (pdg == -checkPdgDecay) { + int igmother = track.getMotherTrackId(); + if (igmother > 0) { + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5) { + nAntileptons++; + nPositrons++; + if (-1 < y && y < 1) nPositronsElectronsInAcceptance++; + if (track.getToBeDone()) { + nAntileptonsToBeDone++; + nPositronsToBeDone++; + } + } + } + } else if (pdg == 511 || pdg == 521 || pdg == 531 || pdg == 5122 || pdg == 5132 || pdg == 5232 || pdg == 5332) { + nOpenBeautyPos++; + } else if (pdg == -511 || pdg == -521 || pdg == -531 || pdg == -5122 || pdg == -5132 || pdg == -5232 || pdg == -5332) { + nOpenBeautyNeg++; + } + } + if (nOpenBeautyPos > 0 && nOpenBeautyNeg > 0) { + nSignalPairs++; + } + if (nPositronsElectronsInAcceptance > 1) { + nLeptonsInAcceptance++; + } + if (nElectrons > 0 && nPositrons > 0) { + nLeptonPairs++; + } + if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++; + } + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#leptons to be done: " << nLeptonsToBeDone << "\n" + << "#antileptons to be done: " << nAntileptonsToBeDone << "\n" + << "#Open-beauty hadron pairs: " << nSignalPairs << "\n" + << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n" + << "#Electron-positron pairs: " << nLeptonPairs << "\n" + << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n"; + if (nLeptons == 0 && nAntileptons == 0) { + std::cerr << "Number of leptons, number of anti-leptons 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; + } + if (nLeptons != nLeptonsToBeDone) { + std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n"; + return 1; + } + if (nLeptonsInAcceptance != nEvents) { + std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorHF_Charm_PbPb_electron.C b/MC/config/PWGEM/ini/tests/GeneratorHF_Charm_PbPb_electron.C new file mode 100644 index 000000000..d55c16c11 --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorHF_Charm_PbPb_electron.C @@ -0,0 +1,116 @@ +int External() +{ + + int checkPdgDecay = 11; + std::string path{"o2sim_Kine.root"}; + 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 nLeptonsInAcceptance{}; + int nLeptons{}; + int nAntileptons{}; + int nLeptonsToBeDone{}; + int nAntileptonsToBeDone{}; + int nSignalPairs{}; + int nLeptonPairs{}; + int nLeptonPairsToBeDone{}; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + int nElectrons = 0; + int nPositrons = 0; + int nElectronsToBeDone = 0; + int nPositronsToBeDone = 0; + int nOpenCharmPos = 0; + int nOpenCharmNeg = 0; + int nPositronsElectronsInAcceptance = 0; + for (auto& track : *tracks) { + auto pdg = track.GetPdgCode(); + auto y = track.GetRapidity(); + if (pdg == checkPdgDecay) { + int igmother = track.getMotherTrackId(); + if (igmother > 0) { + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4) { + nLeptons++; + nElectrons++; + if (-1 < y && y < 1) nPositronsElectronsInAcceptance++; + if (track.getToBeDone()) { + nLeptonsToBeDone++; + nElectronsToBeDone++; + } + } + } + } else if (pdg == -checkPdgDecay) { + int igmother = track.getMotherTrackId(); + if (igmother > 0) { + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(TMath::Abs(gmpdg)/100.) == 4 || int(TMath::Abs(gmpdg)/1000.) == 4) { + nAntileptons++; + nPositrons++; + if (-1 < y && y < 1) nPositronsElectronsInAcceptance++; + if (track.getToBeDone()) { + nAntileptonsToBeDone++; + nPositronsToBeDone++; + } + } + } + } else if (pdg == 411 || pdg == 421 || pdg == 431 || pdg == 4122 || pdg == 4132 || pdg == 4232 || pdg == 4332) { + nOpenCharmPos++; + } else if (pdg == -411 || pdg == -421 || pdg == -431 || pdg == -4122 || pdg == -4132 || pdg == -4232 || pdg == -4332) { + nOpenCharmNeg++; + } + } + if (nOpenCharmPos > 0 && nOpenCharmNeg > 0) { + nSignalPairs++; + } + if (nPositronsElectronsInAcceptance > 1) { + nLeptonsInAcceptance++; + } + if (nElectrons > 0 && nPositrons > 0) { + nLeptonPairs++; + } + if (nElectronsToBeDone > 0 && nPositronsToBeDone > 0) nLeptonPairsToBeDone++; + } + std::cout << "#events: " << nEvents << "\n" + << "#leptons: " << nLeptons << "\n" + << "#antileptons: " << nAntileptons << "\n" + << "#leptons to be done: " << nLeptonsToBeDone << "\n" + << "#antileptons to be done: " << nAntileptonsToBeDone << "\n" + << "#open-charm hadron pairs: " << nSignalPairs << "\n" + << "#leptons in acceptance: " << nLeptonsInAcceptance << "\n" + << "#Electron-positron pairs: " << nLeptonPairs << "\n" + << "#Electron-positron pairs to be done: " << nLeptonPairsToBeDone << "\n"; + if (nLeptons == 0 && nAntileptons == 0) { + std::cerr << "Number of leptons, number of anti-leptons should all be greater than 1.\n"; + return 1; + } + if (nLeptonPairs != nLeptonPairsToBeDone) { + std::cerr << "The number of electron-positron pairs should be the same as the number of electron-positron pairs which should be transported.\n"; + return 1; + } + if (nLeptons != nLeptonsToBeDone) { + std::cerr << "The number of leptons should be the same as the number of leptons which should be transported.\n"; + return 1; + } + if (nLeptonsInAcceptance != nEvents) { + std::cerr << "The number of leptons in acceptance should be at least equaled to the number of events.\n"; + return 1; + } + if (nLeptonPairs < nLeptonsInAcceptance) { + std::cerr << "The number of positron-electron pairs should be at least equaled to the number of leptons in acceptance.\n"; + return 1; + } + + return 0; +} diff --git a/MC/run/PWGEM/runHFToDielectrons_PbPb.sh b/MC/run/PWGEM/runHFToDielectrons_PbPb.sh new file mode 100644 index 000000000..4193f6751 --- /dev/null +++ b/MC/run/PWGEM/runHFToDielectrons_PbPb.sh @@ -0,0 +1,51 @@ +#!/usr/bin/env bash + +RNDSEED=${RNDSEED:-0} +NWORKERS=${NWORKERS:-8} + + +export ALIEN_JDL_LPMPRODUCTIONTYPE=MC +export ALIEN_JDL_CPULIMIT=20 +export ALIEN_JDL_LPMANCHORPASSNAME=apass4 +export ALIEN_JDL_MCANCHOR=apass4 +export ALIEN_JDL_COLLISIONSYSTEM=PbPb +export ALIEN_JDL_LPMPASSNAME=apass4 +export ALIEN_JDL_LPMRUNNUMBER=544474 +export ALIEN_JDL_LPMANCHORRUN=544474 + +export ALIEN_JDL_LPMINTERACTIONTYPE=PbPb +export ALIEN_JDL_LPMPRODUCTIONTAG="TestProd17032026" +export ALIEN_JDL_LPMANCHORPRODUCTION="LHC26ac" +export ALIEN_JDL_LPMANCHORYEAR=2023 + +export ALIEN_JDL_O2DPGWORKFLOWTARGET="aod" + +export PRODSPLIT=${ALIEN_O2DPG_GRIDSUBMIT_PRODSPLIT:-10} +export SPLITID=${ALIEN_O2DPG_GRIDSUBMIT_SUBJOBID:-10} +export CYCLE=0 +export NTIMEFRAMES=8 +export NSIGEVENTS=20 +export NBKGEVENTS=10 + +# define the generator via ini file +# use 20/80 sampling for different generators +# generate random number +RNDSIG=$(($RANDOM % 100)) + + +if [[ $RNDSIG -ge 0 && $RNDSIG -lt 20 ]]; +then + CONFIGNAME="GeneratorHF_Charm_PbPb_electron.ini" +elif [[ $RNDSIG -ge 20 && $RNDSIG -lt 100 ]]; +then + CONFIGNAME="GeneratorHF_BeautyNoForcedDecay_PbPb_electron.ini" +fi + + +# generator and other sim configuration; parallel world for ITS +SIM_OPTIONS="-eCM 5360 -gen external -j ${NWORKERS} -tf ${NTIMEFRAMES} -e TGeant4 -seed 0 -ini ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/ini/$CONFIGNAME -genBkg pythia8 -procBkg \"heavy_ion\" -colBkg PbPb --embedding -nb $NBKGEVENTS" + +export ALIEN_JDL_ANCHOR_SIM_OPTIONS="${SIM_OPTIONS}" + +${O2DPG_ROOT}/MC/run/ANCHOR/anchorMC.sh +