diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.EE.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.EE.DEC new file mode 100644 index 000000000..3bfe3a10b --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.EE.DEC @@ -0,0 +1,5 @@ +Decay J/psi +1.000 e+ e- PHOTOS VLL; +Enddecay +End + diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.MUMU.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.MUMU.DEC new file mode 100644 index 000000000..7e5ab4de6 --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/JPSI.MUMU.DEC @@ -0,0 +1,5 @@ +Decay J/psi +1.000 mu+ mu- PHOTOS VLL; +Enddecay +End + diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OMEGA.3PI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OMEGA.3PI.DEC new file mode 100644 index 000000000..c6641a920 --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OMEGA.3PI.DEC @@ -0,0 +1,11 @@ +Alias mypi0 pi0 + +Decay omega +1.000 pi- pi+ mypi0 OMEGA_DALITZ; +Enddecay + +Decay mypi0 +#1.000 gamma gamma PHSP; +Enddecay + +End diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PHI.KK.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PHI.KK.DEC new file mode 100644 index 000000000..5fa682ee8 --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PHI.KK.DEC @@ -0,0 +1,4 @@ +Decay phi +1.000 K+ K- VSS; +Enddecay +End diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.EEPIPI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.EEPIPI.DEC new file mode 100644 index 000000000..25eea861a --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.EEPIPI.DEC @@ -0,0 +1,11 @@ +Alias myJ/psi J/psi +## +Decay psi(2S) +1.000 myJ/psi pi+ pi- VVPIPI; +Enddecay + +Decay myJ/psi +1.000 e+ e- PHOTOS VLL; +Enddecay +End + diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.MUMUPIPI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.MUMUPIPI.DEC new file mode 100644 index 000000000..d20651ed5 --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/PSI2S.MUMUPIPI.DEC @@ -0,0 +1,11 @@ +Alias myJ/psi J/psi +## +Decay psi(2S) +1.000 myJ/psi pi+ pi- VVPIPI; +Enddecay + +Decay myJ/psi +1.000 mu+ mu- PHOTOS VLL; +Enddecay +End + diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/RHOPRIME.RHOPIPI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/RHOPRIME.RHOPIPI.DEC new file mode 100644 index 000000000..777ef942a --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/RHOPRIME.RHOPIPI.DEC @@ -0,0 +1,15 @@ +Alias myrho0 rho0 +## +Decay rho(2S)0 +1.000 myrho0 pi+ pi- VVPIPI; +Enddecay + +Decay rho(3S)0 +1.000 myrho0 pi+ pi- VVPIPI; +Enddecay + +Decay myrho0 +1.000 pi+ pi- VSS; +Enddecay +End + diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlight.C b/MC/config/PWGUD/external/generator/GeneratorStarlight.C new file mode 100644 index 000000000..21054104e --- /dev/null +++ b/MC/config/PWGUD/external/generator/GeneratorStarlight.C @@ -0,0 +1,283 @@ +R__LOAD_LIBRARY(libStarlib.so) +R__ADD_INCLUDE_PATH($STARlight_ROOT/include) + +#include "randomgenerator.h" +#include "upcXevent.h" +#include "starlight.h" +#include "inputParameters.h" + +// usage: o2-sim -g external --configKeyValues 'GeneratorExternal.fileName=GeneratorStarlight.C;GeneratorExternal.funcName=GeneratorStarlight("kCohJpsiToMu")' + +namespace o2 +{ +namespace eventgen +{ +class GeneratorStarlight_class : public Generator +{ + public: + GeneratorStarlight_class(){}; + ~GeneratorStarlight_class() = default; + void selectConfiguration(std::string val) { mSelectedConfiguration = val; }; + void setCollisionSystem(float energyCM, int beam1Z, int beam1A, int beam2Z, int beam2A) {eCM = energyCM; projZ=beam1Z; projA=beam1A; targZ=beam2Z; targA=beam2A;}; + bool setParameter(std::string line) { + if (not mInputParameters.setParameter(line)){ + std::cout << " --- [ERROR] cannot set parameter: " << line << std::endl; + return false; + } + return true; + } + int getPdgMother(){return mPdgMother;} + + bool Init() override + { + Generator::Init(); + + float beam1energy = TMath::Sqrt(Double_t(projZ)/projA*targA/targZ)*eCM/2; + float beam2energy = TMath::Sqrt(Double_t(projA)/projZ*targZ/targA)*eCM/2; + float gamma1 = beam1energy/0.938272; + float gamma2 = beam2energy/0.938272; + float rapMax = 4.1 + 0.5*(TMath::ACosH(gamma2)-TMath::ACosH(gamma1)); + + const struct SLConfig { + const char* name; + int prod_mode; + int prod_pid; + int nwbins; + float wmin; + float wmax; + float dy; + int pdg_mother; + bool decay_EvtGen; + } slConfig[] = { + {"kTwoGammaToMuLow", 1, 13, 292, 0.4, 15.0, 0.01, -1, 0 }, // from 0.4 to 15 GeV + {"kTwoGammaToElLow", 1, 11, 292, 0.4, 15.0, 0.01, -1, 0 }, // from 0.4 to 15 GeV + {"kTwoGammaToMuMedium", 1, 13, 264, 1.8, 15.0, 0.01, -1, 0 }, // from 1.8 to 15 GeV + {"kTwoGammaToElMedium", 1, 11, 264, 1.8, 15.0, 0.01, -1, 0 }, // from 1.8 to 15 GeV + {"kTwoGammaToMuHigh", 1, 13, 220, 4.0, 15.0, 0.01, -1, 0 }, // from 4.0 to 15 GeV + {"kTwoGammaToElHigh", 1, 11, 220, 4.0, 15.0, 0.01, -1, 0 }, // from 4.0 to 15 GeV + {"kTwoGammaToRhoRho", 1, 33, 20, -1.0, -1.0, 0.01, -1, 0 }, // + {"kTwoGammaToF2", 1, 225, 20, -1.0, -1.0, 0.01, -1, 0 }, // + {"kCohRhoToPi", 3, 113, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kCohRhoToElEl", 3, 113011, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kCohRhoToMuMu", 3, 113013, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kCohRhoToPiWithCont", 3, 913, 1200, -1.0, -1.0, 0.02, -1, 0 }, // + {"kCohRhoToPiFlat", 3, 113, 1, -1.0, 2.5, 0.02, 113, 0 }, // + {"kCohPhiToKa", 2, 333, 20, -1.0, -1.0, 0.01, 333, 0 }, // + {"kDirectPhiToKaKa", 3, 933, 20, -1.0, -1.0, 0.01, 333, 0 }, // + {"kCohOmegaTo2Pi", 2, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // + {"kCohOmegaTo3Pi", 2, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // + {"kCohOmegaToPiPiPi", 2, 223211111, 20, -1.0, -1.0, 0.01, 333, 0 }, // + {"kCohJpsiToMu", 2, 443013, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kCohJpsiToEl", 2, 443011, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kCohJpsiToElRad", 2, 443011, 20, -1.0, -1.0, 0.01, 443, 1 }, // + {"kCohJpsiToProton", 2, 4432212, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kCohPsi2sToMu", 2, 444013, 20, -1.0, -1.0, 0.01, 100443, 0 }, // + {"kCohPsi2sToEl", 2, 444011, 20, -1.0, -1.0, 0.01, 100443, 0 }, // + {"kCohPsi2sToMuPi", 2, 444013, 20, -1.0, -1.0, 0.01, 100443, 1 }, // + {"kCohPsi2sToElPi", 2, 444011, 20, -1.0, -1.0, 0.01, 100443, 1 }, // + {"kCohUpsilonToMu", 2, 553013, 20, -1.0, -1.0, 0.01, 553, 0 }, // + {"kCohUpsilonToEl", 2, 553011, 20, -1.0, -1.0, 0.01, 553, 0 }, // + {"kIncohRhoToPi", 4, 113, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kIncohRhoToElEl", 4, 113011, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kIncohRhoToMuMu", 4, 113013, 1200, -1.0, -1.0, 0.02, 113, 0 }, // + {"kIncohRhoToPiWithCont",4, 913, 1200, -1.0, -1.0, 0.02, -1, 0 }, // + {"kIncohRhoToPiFlat", 4, 113, 1, -1.0, 2.5, 0.02, 113, 0 }, // + {"kIncohPhiToKa", 4, 333, 20, -1.0, -1.0, 0.01, 333, 0 }, // + {"kIncohOmegaTo2Pi", 4, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // + {"kIncohOmegaTo3Pi", 4, 223, 20, -1.0, -1.0, 0.01, 223, 0 }, // + {"kIncohOmegaToPiPiPi", 4, 223211111, 20, -1.0, -1.0, 0.01, 223, 0 }, // + {"kIncohJpsiToMu", 4, 443013, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kIncohJpsiToEl", 4, 443011, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kIncohJpsiToElRad", 4, 443011, 20, -1.0, -1.0, 0.01, 443, 1 }, // + {"kIncohJpsiToProton", 4, 4432212, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kIncohJpsiToLLbar", 4, 4433122, 20, -1.0, -1.0, 0.01, 443, 0 }, // + {"kIncohPsi2sToMu", 4, 444013, 20, -1.0, -1.0, 0.01, 100443, 0 }, // + {"kIncohPsi2sToEl", 4, 444011, 20, -1.0, -1.0, 0.01, 100443, 0 }, // + {"kIncohPsi2sToMuPi", 4, 444013, 20, -1.0, -1.0, 0.01, 100443, 1 }, // + {"kIncohPsi2sToElPi", 4, 444011, 20, -1.0, -1.0, 0.01, 100443, 1 }, // + {"kIncohUpsilonToMu", 4, 553013, 20, -1.0, -1.0, 0.01, 553, 0 }, // + {"kIncohUpsilonToEl", 4, 553011, 20, -1.0, -1.0, 0.01, 553, 0 }, // + }; + + const int nProcess = sizeof(slConfig)/sizeof(SLConfig); + int idx = -1; + for (int i=0; i(&random_value), sizeof(random_seed)); + + setParameter(Form("BEAM_1_Z = %3i #Z of target",targZ)); + setParameter(Form("BEAM_1_A = %3i #A of target",targA)); + setParameter(Form("BEAM_2_Z = %3i #Z of projectile",projZ)); + setParameter(Form("BEAM_2_A = %3i #A of projectile",projA)); + setParameter(Form("BEAM_1_GAMMA = %6.1f #Gamma of the target",gamma1)); + setParameter(Form("BEAM_2_GAMMA = %6.1f #Gamma of the projectile",gamma2)); + setParameter(Form("W_MAX = %.1f #Max value of w",slConfig[idx].wmax)); + setParameter(Form("W_MIN = %.1f #Min value of w",slConfig[idx].wmin)); + setParameter(Form("W_N_BINS = %3i #Bins i w",slConfig[idx].nwbins)); + setParameter(Form("RAP_MAX = %.2f #max y",rapMax)); + setParameter(Form("RAP_N_BINS = %.0f #Bins i y",rapMax*2./slConfig[idx].dy)); + setParameter("CUT_PT = 0 #Cut in pT? 0 = (no, 1 = yes)"); + setParameter("PT_MIN = 0 #Minimum pT in GeV"); + setParameter("PT_MAX = 10 #Maximum pT in GeV"); + setParameter("CUT_ETA = 0 #Cut in pseudorapidity? (0 = no, 1 = yes)"); + setParameter("ETA_MIN = -5 #Minimum pseudorapidity"); + setParameter("ETA_MAX = 5 #Maximum pseudorapidity"); + setParameter(Form("PROD_MODE = %i #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), # 4 = incoherent vector meson, 5 = A+A DPMJet single, 6 = A+A DPMJet double, 7 = p+A DPMJet single, 8 = p+A Pythia single )",slConfig[idx].prod_mode)); + setParameter(Form("PROD_PID = %6i #Channel of interest (not relevant for photonuclear processes)",slConfig[idx].prod_pid)); + setParameter(Form("RND_SEED = %i #Random number seed", random_seed)); + setParameter("BREAKUP_MODE = 5 #Controls the nuclear breakup"); + setParameter("INTERFERENCE = 0 #Interference (0 = off, 1 = on)"); + setParameter("IF_STRENGTH = 1. #% of interfernce (0.0 - 0.1)"); + setParameter("INT_PT_MAX = 0.24 #Maximum pt considered, when interference is turned on"); + setParameter("INT_PT_N_BINS = 120 #Number of pt bins when interference is turned on"); + setParameter("XSEC_METHOD = 1 # Set to 0 to use old method for calculating gamma-gamma luminosity"); //CM + setParameter("BSLOPE_DEFINITION = 0"); // using default slope + setParameter("BSLOPE_VALUE = 4.0"); // default slope value + setParameter("PRINT_VM = 0"); // print cross sections and fluxes vs rapidity in stdout for VM photoproduction processes + + if (not mInputParameters.init()) { + std::cout << "InitStarLight parameter initialization has failed" << std::endl; + return false; + } + + mStarLight = new starlight; + mStarLight->setInputParameters(&mInputParameters); + mRandomGenerator.SetSeed(mInputParameters.randomSeed()); + mStarLight->setRandomGenerator(&mRandomGenerator); + return mStarLight->init(); + + }; + + bool generateEvent() override { + + if (!mStarLight) { + std::cout <<"GenerateEvent: StarLight class/object not properly constructed"<produceEvent(); + // boost event to the experiment CM frame + mEvent.boost(0.5*(TMath::ACosH(mInputParameters.beam1LorentzGamma()) - TMath::ACosH(mInputParameters.beam2LorentzGamma()))); + + return true; + + }; + + // at importParticles we add particles to the output particle vector + // according to the selected configuration + bool importParticles() override + { + int nVtx(0); + float vtx(0), vty(0), vtz(0), vtt(0); + const std::vector* slVtx(mEvent.getVertices()); + if (slVtx == 0) { // not vertex assume 0,0,0,0; + vtx = vty = vtz = vtt = 0.0; + } else { // a vertex exits + slVtx = mEvent.getVertices(); + nVtx = slVtx->size(); + } // end if + + const std::vector* slPartArr(mEvent.getParticles()); + const int npart(mEvent.getParticles()->size()); + + if(mPdgMother != -1){ //Reconstruct mother particle for VM processes + TLorentzVector lmoth; + TLorentzVector ldaug; + for(int ipart=0;ipartat(ipart))); + ldaug.SetPxPyPzE(slPart->GetPx(), slPart->GetPy(), slPart->GetPz(), slPart->GetE()); + lmoth += ldaug; + } + TParticle particle(mPdgMother, + 11, + -1, + -1, + 1, + npart, + lmoth.Px(), + lmoth.Py(), + lmoth.Pz(), + lmoth.E(), + 0,0,0,0); + //particle.Print(); + mParticles.push_back(particle); + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), 11); + } + if(!mDecayEvtGen){ // Don't import daughters in case of external decayer + for(int ipart=0;ipartat(ipart))); + if (nVtx < 1) { // No verticies + vtx = vty = vtz = vtt = 0.0; + } else { + vtx = (slVtx->at((ipart < nVtx ? ipart : 0))).X(); + vty = (slVtx->at((ipart < nVtx ? ipart : 0))).Y(); + vtz = (slVtx->at((ipart < nVtx ? ipart : 0))).Z(); + vtt = 0.0; // no time given. + } // end if + TParticle particle(slPart->getPdgCode(), + 1, + 0, + -1, + slPart->getFirstDaughter(), + slPart->getLastDaughter(), + slPart->GetPx(), + slPart->GetPy(), + slPart->GetPz(), + slPart->GetE(), + vtx,vty,vtz,vtt); + //particle.Print(); + mParticles.push_back(particle); + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), 1); + } + } + return true; + } + + private: + starlight *mStarLight = 0x0; + inputParameters mInputParameters; // simulation input information. + randomGenerator mRandomGenerator; // STARLIGHT's own random generator + upcXEvent mEvent; // object holding STARlight simulated event. + std::string mSelectedConfiguration = ""; + int mPdgMother = -1; + bool mDecayEvtGen = 0; + float eCM = 5020; //CMS energy + int projA=208; //Beam + int targA=208; + int projZ=82; + int targZ=82; + + }; + +} // namespace eventgen +} // namespace o2 + + +FairGenerator* + GeneratorStarlight(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208) +{ + auto gen = new o2::eventgen::GeneratorStarlight_class(); + gen->selectConfiguration(configuration); + gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); + return gen; +} + + + + diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C b/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C new file mode 100644 index 000000000..3fe6e17c0 --- /dev/null +++ b/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C @@ -0,0 +1,33 @@ +// usage (fwdy) : +// o2-sim -j 4 -n 10 -g external -t external -m "PIPE ITS TPC" -o sgn --configFile GeneratorHF_bbbar_fwdy.ini +// usage (midy) : +// o2-sim -j 4 -n 10 -g external -t external -m "PIPE ITS TPC" -o sgn --configFile GeneratorHF_bbbar_midy.ini +// +// +R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen) +R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGUD/external/generator) +#include "GeneratorEvtGen.C" +#include "GeneratorStarlight.C" + +FairGenerator* + GeneratorStarlightToEvtGen(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208) +{ + auto gen = new o2::eventgen::GeneratorEvtGen(); + gen->selectConfiguration(configuration); + gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); + + gen->SetSizePdg(3); + gen->AddPdg(443,0); + gen->AddPdg(100443,1); + gen->AddPdg(223,2); + gen->SetPolarization(1); //Transversal + + TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGUD/external/generator/DecayTablesEvtGen"); + if (configuration.find("Psi2sToMuPi") != std::string::npos) gen->SetDecayTable(Form("%s/PSI2S.MUMUPIPI.DEC",pathO2.Data())); + else if (configuration.find("Psi2sToElPi") != std::string::npos) gen->SetDecayTable(Form("%s/PSI2S.EEPIPI.DEC",pathO2.Data())); + else if (configuration.find("RhoPrime") != std::string::npos) gen->SetDecayTable(Form("%s/RHOPRIME.RHOPIPI.DEC",pathO2.Data())); + else if (configuration.find("OmegaTo3Pi") != std::string::npos) gen->SetDecayTable(Form("%s/OMEGA.3PI.DEC",pathO2.Data())); + else if (configuration.find("JpsiToElRad") != std::string::npos) gen->SetDecayTable(Form("%s/JPSI.EE.DEC",pathO2.Data())); + + return gen; +} \ No newline at end of file diff --git a/MC/config/PWGUD/ini/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.ini b/MC/config/PWGUD/ini/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.ini new file mode 100644 index 000000000..8c38fd160 --- /dev/null +++ b/MC/config/PWGUD/ini/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.ini @@ -0,0 +1,6 @@ +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlight.C +funcName = GeneratorStarlight("kCohJpsiToMu", 5360.000000, 82, 208, 82, 208) +[TriggerExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C +funcName = selectMotherPartInAcc(-4.0,-2.5) diff --git a/MC/config/PWGUD/ini/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.ini b/MC/config/PWGUD/ini/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.ini new file mode 100644 index 000000000..d696abde2 --- /dev/null +++ b/MC/config/PWGUD/ini/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.ini @@ -0,0 +1,6 @@ +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C +funcName = GeneratorStarlightToEvtGen("kCohPsi2sToElPi", 5360.000000, 82, 208, 82, 208) +[TriggerExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C +funcName = selectMotherPartInAcc(-0.9,0.9) diff --git a/MC/config/PWGUD/ini/makeStarlightConfig.py b/MC/config/PWGUD/ini/makeStarlightConfig.py new file mode 100755 index 000000000..6d87e3a9f --- /dev/null +++ b/MC/config/PWGUD/ini/makeStarlightConfig.py @@ -0,0 +1,60 @@ +#! /usr/bin/env python3 + +### @author: Michal Broz +### @email: michal.broz@cern.ch + +import argparse + +parser = argparse.ArgumentParser(description='Make Starlight configuration', + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + +parser.add_argument('--pZ', type=int, default='82', + help='Projectile charge') + +parser.add_argument('--pA', type=int, default='208', + help='Projectile atomic number') + +parser.add_argument('--tZ', type=int, default='82', + help='Target charge') + +parser.add_argument('--tA', type=int, default='208', + help='Target atomic number') + +parser.add_argument('--eCM', type=float, default='5360', + help='Centre-of-mass energy') + +parser.add_argument('--rapidity', default='cent', choices=['cent', 'muon'], + help='Rapidity to select') + +parser.add_argument('--process',default=None, choices=['kTwoGammaToMuLow', 'kTwoGammaToElLow', 'kTwoGammaToMuMedium', 'kTwoGammaToElMedium', 'kTwoGammaToMuHigh', 'kTwoGammaToElHigh', 'kTwoGammaToRhoRho', 'kTwoGammaToF2', 'kCohRhoToPi', 'kCohRhoToElEl', 'kCohRhoToMuMu', 'kCohRhoToPiWithCont', 'kCohRhoToPiFlat', 'kCohPhiToKa', 'kDirectPhiToKaKa','kCohOmegaTo2Pi', 'kCohOmegaTo3Pi', 'kCohOmegaToPiPiPi', 'kCohJpsiToMu', 'kCohJpsiToEl', 'kCohJpsiToElRad', 'kCohJpsiToProton', 'kCohPsi2sToMu','kCohPsi2sToEl', 'kCohPsi2sToMuPi', 'kCohPsi2sToElPi', 'kCohUpsilonToMu', 'kCohUpsilonToEl', 'kIncohRhoToPi', 'kIncohRhoToElEl', 'kIncohRhoToMuMu', 'kIncohRhoToPiWithCont', 'kIncohRhoToPiFlat', 'kIncohPhiToKa', 'kIncohOmegaTo2Pi', 'kIncohOmegaTo3Pi', 'kIncohOmegaToPiPiPi', 'kIncohJpsiToMu', 'kIncohJpsiToEl', 'kIncohJpsiToElRad', 'kIncohJpsiToProton', 'kIncohJpsiToLLbar', 'kIncohPsi2sToMu', 'kIncohPsi2sToEl', 'kIncohPsi2sToMuPi', 'kIncohPsi2sToElPi', 'kIncohUpsilonToMu', 'kIncohUpsilonToEl'], + help='Process to switch on') + + +parser.add_argument('--output', default='GenStarlight.ini', + help='Where to write the configuration') + + +args = parser.parse_args() + +### open output file +fout = open(args.output, 'w') + +### Generator +fout.write('[GeneratorExternal] \n') +if 'Psi2sToMuPi' in args.process or 'Psi2sToElPi' in args.process or 'RhoPrime' in args.process or 'OmegaTo3Pi' in args.process or 'JpsiToElRad' in args.process : + fout.write('fileName = ${O2DPG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C \n') + fout.write('funcName = GeneratorStarlightToEvtGen("%s", %f, %d, %d, %d, %d) \n' % (args.process,args.eCM ,args.pZ ,args.pA,args.tZ,args.tA)) +else: + fout.write('fileName = ${O2DPG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlight.C \n') + fout.write('funcName = GeneratorStarlight("%s", %f, %d, %d, %d, %d) \n' % (args.process,args.eCM ,args.pZ ,args.pA,args.tZ,args.tA)) + +###Trigger +fout.write('[TriggerExternal] \n') +fout.write('fileName = ${O2DPG_ROOT}/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C \n') +if args.rapidity == 'cent': + fout.write('funcName = selectMotherPartInAcc(-0.9,0.9) \n') +if args.rapidity == 'muon': + fout.write('funcName = selectMotherPartInAcc(-4.0,-2.5) \n') + +### close outout file +fout.close() diff --git a/MC/config/PWGUD/ini/tests/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.C b/MC/config/PWGUD/ini/tests/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.C new file mode 100644 index 000000000..2a9424d87 --- /dev/null +++ b/MC/config/PWGUD/ini/tests/GenStarlight_kCohJpsiToMu_PbPb_5360_Muon.C @@ -0,0 +1,37 @@ +int External() +{ + 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"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + bool goodEvent = kFALSE; + auto nEvents = tree->GetEntries(); + for (int i = 0; i < nEvents; i++) + { + goodEvent = kFALSE; + auto check = tree->GetEntry(i); + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 443 && track.getMotherTrackId() == -1){ //Primary J/psi + auto daugh1 = tracks->at(track.getFirstDaughterTrackId()); + auto daugh2 = tracks->at(track.getLastDaughterTrackId()); + if(TMath::Abs(daugh1.GetPdgCode()) == 13 && TMath::Abs(daugh2.GetPdgCode()) == 13) goodEvent = kTRUE; + } + } + if(!goodEvent) return 1; + } + return 0; +} diff --git a/MC/config/PWGUD/ini/tests/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.C b/MC/config/PWGUD/ini/tests/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.C new file mode 100644 index 000000000..74b40771c --- /dev/null +++ b/MC/config/PWGUD/ini/tests/GenStarlight_kCohPsi2sToElPi_PbPb_5360_Cent.C @@ -0,0 +1,51 @@ +int External() +{ + 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"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nPions = 0; + int nElectrons = 0; + int nPhotons = 0; + auto nEvents = tree->GetEntries(); + for (int i = 0; i < nEvents; i++) + { + nPhotons = 0; + nPions = 0; + nElectrons = 0; + auto check = tree->GetEntry(i); + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 100443 && track.getMotherTrackId() == -1){ //Primary psi2s + for(int iDaugh = track.getFirstDaughterTrackId(); iDaugh<=track.getLastDaughterTrackId(); iDaugh++){ + auto daughPsi2s = tracks->at(iDaugh); + if(daughPsi2s.GetPdgCode() == 22) nPhotons++; + if(TMath::Abs(daughPsi2s.GetPdgCode()) == 211) nPions++; + + if(daughPsi2s.GetPdgCode() == 443 ){ + for(int jDaugh = daughPsi2s.getFirstDaughterTrackId(); jDaugh<=daughPsi2s.getLastDaughterTrackId(); jDaugh++){ + auto daughJPsi = tracks->at(jDaugh); + if(daughJPsi.GetPdgCode() == 22) nPhotons++; + if(TMath::Abs(daughJPsi.GetPdgCode()) == 11) nElectrons++; + } + } + } + } + } + if(nElectrons != 2 || nPions != 2)return 1; + } + return 0; +} diff --git a/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C b/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C new file mode 100644 index 000000000..a00bcad6b --- /dev/null +++ b/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C @@ -0,0 +1,31 @@ +R__ADD_INCLUDE_PATH($O2DPG_ROOT) +#include +#include "Generators/Trigger.h" + +/// ================================================================================================================================= +/// Select event with VM or track in a given rapidity or eta window +/// ================================================================================================================================= + +o2::eventgen::Trigger selectMotherPartInAcc(double rapidityMin = -1., double rapidityMax = -1.) +{ + return [rapidityMin, rapidityMax](const std::vector& particles) -> bool { + for (const auto& particle : particles) { + if (particle.GetFirstMother() == -1) + if ((particle.Y() > rapidityMin) && (particle.Y() < rapidityMax)) + return kTRUE; + } + return kFALSE; + }; +} + +o2::eventgen::Trigger selectDaughterPartInAcc(double etaMin = -1., double etaMax = -1.) +{ + return [etaMin, etaMax](const std::vector& particles) -> bool { + for (const auto& particle : particles) { + if (particle.GetFirstMother() != -1) + if ((particle.Eta() < etaMin) || (particle.Eta() > etaMax)) + return kFALSE; + } + return kTRUE; + }; +} \ No newline at end of file