diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OPENCHARM.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OPENCHARM.DEC new file mode 100644 index 000000000..d35930c4a --- /dev/null +++ b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/OPENCHARM.DEC @@ -0,0 +1,45 @@ +Decay D*+ +1.0 D0 pi+ VSS; +Enddecay + +Decay D*- +1.0 anti-D0 pi- VSS; +Enddecay + +Decay D+ +1.0 K- pi+ pi+ D_DALITZ; #[Reconstructed PDG2011] +Enddecay + +Decay D- +1.0 K+ pi- pi- D_DALITZ; #[Reconstructed PDG2011] +Enddecay + +Decay D0 +1.0 K- pi+ PHSP; #[Reconstructed PDG2011] +Enddecay + +Decay anti-D0 +1.0 K+ pi- PHSP; #[Reconstructed PDG2011] +Enddecay + +Decay Lambda_c+ +1.0 p+ K- pi+ PHSP; #[New mode added] #[Reconstructed PDG2011] +Enddecay + +Decay anti-Lambda_c- +1.0 anti-p- K+ pi- PHSP; #[New mode added] #[Reconstructed PDG2011] +Enddecay + +Decay D_s+ +1.0 phi pi+ SVS; #[Reconstructed PDG2011] +Enddecay + +Decay D_s- +1.0 phi pi- SVS; #[Reconstructed PDG2011] +Enddecay + +Decay phi +1.0 K+ K- VSS; #[Reconstructed PDG2011] +Enddecay + +End diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlight.C b/MC/config/PWGUD/external/generator/GeneratorStarlight.C index 66402130b..9c317ac81 100644 --- a/MC/config/PWGUD/external/generator/GeneratorStarlight.C +++ b/MC/config/PWGUD/external/generator/GeneratorStarlight.C @@ -37,6 +37,27 @@ class GeneratorStarlight_class : public Generator public: GeneratorStarlight_class(){}; ~GeneratorStarlight_class() = default; + + void setupDpmjet(std::string dpmjetconf){ + if(dpmjetconf.size() == 0)return; + //Copy necesary files to the working directory + TString pathDPMJET = gSystem->ExpandPathName("$DPMJET_ROOT/dpmdata"); + system(TString::Format("cp -r %s .",pathDPMJET.Data())); + system(TString::Format("cp %s ./my.input",dpmjetconf.c_str())); + + //Reset four seeds of the DPMJET random generator in the config + std::mt19937 gen(generateRandomSeed()); + std::uniform_int_distribution<> dist(1, 168); + + std::string command = "awk -i inplace -v nums=\""; + for (int i = 0; i < 4; ++i)command += TString::Format("%d.0 ", dist(gen)); + command +=" \" \' "; + command += "BEGIN {split(nums, newvals);}"; + command += "{if ($1 == \"RNDMINIT\") {printf \"%-16s%-9s%-9s%-9s%-9s\\n\", $1, newvals[1], newvals[2], newvals[3], newvals[4];}"; + command += " else {print $0;}}\' \"my.input\" "; + system(command.c_str()); + } + void selectConfiguration(std::string val) { mSelectedConfiguration = val; }; void setExtraParams(std::string val) { mExtraParams = val; }; void setCollisionSystem(float energyCM, int beam1Z, int beam1A, int beam2Z, int beam2A) {eCM = energyCM; projZ=beam1Z; projA=beam1A; targZ=beam2Z; targA=beam2A;}; @@ -48,6 +69,10 @@ class GeneratorStarlight_class : public Generator return true; } int getPdgMother(){return mPdgMother;} + double getPhotonEnergy(){ + //std::cout << mEvent.getGamma().gamma.GetE() << std::endl; + return mEvent.getGamma().gamma.GetE(); + } bool Init() override { @@ -305,7 +330,14 @@ class GeneratorStarlight_class : public Generator } return true; } - + + protected: + float eCM = 5020; //CMS energy + int projA=208; //Beam + int targA=208; + int projZ=82; + int targZ=82; + private: starlight *mStarLight = 0x0; inputParameters mInputParameters; // simulation input information. @@ -316,11 +348,7 @@ class GeneratorStarlight_class : public Generator std::string mExtraParams = ""; int mPdgMother = -1; bool mDecayEvtGen = 0; - float eCM = 5020; //CMS energy - int projA=208; //Beam - int targA=208; - int projZ=82; - int targZ=82; + }; @@ -331,29 +359,12 @@ class GeneratorStarlight_class : public Generator FairGenerator* GeneratorStarlight(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "",std::string dpmjetconf = "") { - if(dpmjetconf.size() != 0){ - //Copy necesary files to the working directory - TString pathDPMJET = gSystem->ExpandPathName("$DPMJET_ROOT/dpmdata"); - system(TString::Format("cp -r %s .",pathDPMJET.Data())); - system(TString::Format("cp %s ./my.input",dpmjetconf.c_str())); - - //Reset four seeds of the DPMJET random generator in the config - std::mt19937 gen(generateRandomSeed()); - std::uniform_int_distribution<> dist(1, 168); - - std::string command = "awk -i inplace -v nums=\""; - for (int i = 0; i < 4; ++i)command += TString::Format("%d.0 ", dist(gen)); - command +=" \" \' "; - command += "BEGIN {split(nums, newvals);}"; - command += "{if ($1 == \"RNDMINIT\") {printf \"%-16s%-9s%-9s%-9s%-9s\\n\", $1, newvals[1], newvals[2], newvals[3], newvals[4];}"; - command += " else {print $0;}}\' \"my.input\" "; - system(command.c_str()); - } auto gen = new o2::eventgen::GeneratorStarlight_class(); gen->selectConfiguration(configuration); gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); gen->setExtraParams(extrapars); + gen->setupDpmjet(dpmjetconf); return gen; } diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C b/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C index f743ff0e7..3bf3afaa6 100644 --- a/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C +++ b/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C @@ -1,29 +1,43 @@ -// 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_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen) R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_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, std::string extrapars = "") + GeneratorStarlightToEvtGen(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "", std::string dpmjetconf = "") { auto gen = new o2::eventgen::GeneratorEvtGen(); gen->selectConfiguration(configuration); gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); gen->setExtraParams(extrapars); - - gen->SetSizePdg(5); - gen->AddPdg(443,0); - gen->AddPdg(100443,1); - gen->AddPdg(223,2); - gen->AddPdg(15,3); - gen->AddPdg(-15,4); - if (configuration.find("kTau") == std::string::npos) gen->SetPolarization(1); //Transversal + gen->setupDpmjet(dpmjetconf); + + if (configuration.find("kTau") != std::string::npos){ + gen->SetSizePdg(2); + gen->AddPdg(15,0); + gen->AddPdg(-15,1); + } + else if(configuration.find("kDpmjet") != std::string::npos){ + gen->SetSizePdg(11); + gen->AddPdg( 411,0); + gen->AddPdg(-411,1); + gen->AddPdg( 421,2); + gen->AddPdg(-421,3); + gen->AddPdg( 413,4); + gen->AddPdg(-413,5); + gen->AddPdg( 431,6); + gen->AddPdg(-431,7); + gen->AddPdg( 4122,8); + gen->AddPdg(-4122,9); + gen->AddPdg( 333,10); + } + else{ + gen->SetPolarization(1); //Transversal + gen->SetSizePdg(3); + gen->AddPdg(443,0); + gen->AddPdg(100443,1); + gen->AddPdg(223,2); + } TString pathO2 = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGUD/external/generator/DecayTablesEvtGen"); if (configuration.find("Psi2sToMuPi") != std::string::npos) gen->SetDecayTable(Form("%s/PSI2S.MUMUPIPI.DEC",pathO2.Data())); @@ -38,6 +52,7 @@ FairGenerator* else if (configuration.find("ToPoPiPi0") != std::string::npos) gen->SetDecayTable(Form("%s/TAUTAU.POPI.DEC",pathO2.Data())); else if (configuration.find("Jpsi4Prong") != std::string::npos) gen->SetDecayTable(Form("%s/JPSI.4PRONG.DEC",pathO2.Data())); else if (configuration.find("Jpsi6Prong") != std::string::npos) gen->SetDecayTable(Form("%s/JPSI.6PRONG.DEC",pathO2.Data())); + else if (configuration.find("Dpmjet") != std::string::npos) gen->SetDecayTable(Form("%s/OPENCHARM.DEC",pathO2.Data())); return gen; } diff --git a/MC/config/PWGUD/external/generator/Generator_nOOn.C b/MC/config/PWGUD/external/generator/Generator_nOOn.C new file mode 100644 index 000000000..dc19f9504 --- /dev/null +++ b/MC/config/PWGUD/external/generator/Generator_nOOn.C @@ -0,0 +1,65 @@ +R__LOAD_LIBRARY(NeutronGenerator_cxx.so) +#include "GeneratorStarlight.C" +#include "NeutronGenerator.h" + +class Generator_nOOn_class : public o2::eventgen::GeneratorStarlight_class +{ +public: + /// constructor + Generator_nOOn_class(){}; + /// Destructor + ~Generator_nOOn_class() = default; + + bool Init() override + { + GeneratorStarlight_class::Init(); + mNeutronGen = new NeutronGenerator(); + mNeutronGen->SetRapidityCut(-6.0,6.0); + + float beam1energy = TMath::Sqrt(Double_t(projZ)/projA*targA/targZ)*eCM/2; + float gamma1 = beam1energy/0.938272; + mNeutronGen->SetRunMode(NeutronGenerator::kInterface); + mNeutronGen->SetBeamParameters(NeutronGenerator::kPb208,gamma1); + mNeutronGen->SetDataPath(gSystem->ExpandPathName("$nOOn_ROOT/include/Data/")); + mNeutronGen->Initialize(); + mNeutronGen->Setup(); + return true; + } + + bool generateEvent() override + { + GeneratorStarlight_class::generateEvent(); + mNeutronGen->GenerateEvent(getPhotonEnergy()); + return true; + } + + bool importParticles() override + { + GeneratorStarlight_class::importParticles(); + + mNeutrons = mNeutronGen->ImportParticles(); + for(Int_t i = 0; iGetEntriesFast(); i++){ + mParticles.push_back(*(TParticle*)(mNeutrons->At(i))); + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), true); + } + mNeutronGen->FinishEvent(); + mNeutrons->Clear("C"); + return true; + } + + +private: + NeutronGenerator *mNeutronGen = 0x0; + TClonesArray *mNeutrons = 0x0; + +}; + +///___________________________________________________________ +FairGenerator *Generator_nOOn(std::string configuration = "empty",float energyCM = 5020, int beam1Z = 82, int beam1A = 208, int beam2Z = 82, int beam2A = 208, std::string extrapars = "") +{ + auto gen = new Generator_nOOn_class(); + gen->selectConfiguration(configuration); + gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A); + gen->setExtraParams(extrapars); + return gen; +} diff --git a/MC/config/PWGUD/ini/makeStarlightConfig.py b/MC/config/PWGUD/ini/makeStarlightConfig.py index 6a6a163df..48feb7548 100755 --- a/MC/config/PWGUD/ini/makeStarlightConfig.py +++ b/MC/config/PWGUD/ini/makeStarlightConfig.py @@ -30,9 +30,15 @@ parser.add_argument('--dpmjetConf', default='', help='DPMJET config file') +parser.add_argument('--nOOn', action='store_true', + help="Enable the neutron production with nOOn") + args = parser.parse_args() +if args.nOOn: + args.extraPars = 'BREAKUP_MODE = 4' + if 'PbPb' in args.collType: pZ = 82 pA = 208 @@ -74,12 +80,16 @@ ### Generator fout.write('[GeneratorExternal] \n') -if 'Psi2sToMuPi' in args.process or 'Psi2sToElPi' in args.process or 'OmegaTo3Pi' in args.process or 'JpsiToElRad' in args.process or 'Jpsi4Prong' in args.process or 'Jpsi6Prong' in args.process or 'kTau' in args.process: - fout.write('fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C \n') - fout.write('funcName = GeneratorStarlightToEvtGen("%s", %f, %d, %d, %d, %d, "%s") \n' % (args.process,args.eCM ,pZ,pA,tZ,tA,args.extraPars)) +if args.nOOn: + fout.write('fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/external/generator/Generator_nOOn.C \n') + fout.write('funcName = Generator_nOOn("%s", %f, %d, %d, %d, %d, "%s") \n' % (args.process,args.eCM ,pZ,pA,tZ,tA,args.extraPars)) else: - fout.write('fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlight.C \n') - fout.write('funcName = GeneratorStarlight("%s", %f, %d, %d, %d, %d, "%s", "%s") \n' % (args.process.split('_')[0],args.eCM ,pZ,pA,tZ,tA,args.extraPars,args.dpmjetConf)) + if 'Psi2sToMuPi' in args.process or 'Psi2sToElPi' in args.process or 'OmegaTo3Pi' in args.process or 'JpsiToElRad' in args.process or 'Jpsi4Prong' in args.process or 'Jpsi6Prong' in args.process or 'kTau' in args.process or 'Dpmjet' in args.process: + fout.write('fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C \n') + fout.write('funcName = GeneratorStarlightToEvtGen("%s", %f, %d, %d, %d, %d, "%s", "%s") \n' % (args.process.split('_')[0],args.eCM ,pZ,pA,tZ,tA,args.extraPars,args.dpmjetConf)) + else: + fout.write('fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGUD/external/generator/GeneratorStarlight.C \n') + fout.write('funcName = GeneratorStarlight("%s", %f, %d, %d, %d, %d, "%s", "%s") \n' % (args.process.split('_')[0],args.eCM ,pZ,pA,tZ,tA,args.extraPars,args.dpmjetConf)) ###Trigger if not 'kDpmjet' in args.process: diff --git a/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C b/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C index ea669c49c..8e74b0811 100644 --- a/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C +++ b/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C @@ -10,6 +10,7 @@ o2::eventgen::Trigger selectMotherPartInAcc(double rapidityMin = -1., double rap { return [rapidityMin, rapidityMax](const std::vector& particles) -> bool { for (const auto& particle : particles) { + if (TMath::Abs(particle.GetPdgCode()) == 2112)continue; if (particle.GetFirstMother() == -1) if ((particle.Y() > rapidityMin) && (particle.Y() < rapidityMax)) return kTRUE; @@ -22,6 +23,7 @@ o2::eventgen::Trigger selectDaughterPartInAcc(double etaMin = -1., double etaMax { return [etaMin, etaMax](const std::vector& particles) -> bool { for (const auto& particle : particles) { + if (TMath::Abs(particle.GetPdgCode()) == 2112)continue; if (particle.GetFirstMother() == -1) if ((particle.Y() < etaMin) || (particle.Y() > etaMax)) return kFALSE; if (particle.GetFirstMother() != -1 && particle.GetFirstDaughter() == -1 && particle.GetPdgCode() != 22 && TMath::Abs(particle.GetPdgCode()) != 12 && TMath::Abs(particle.GetPdgCode()) != 14 && TMath::Abs(particle.GetPdgCode()) != 16) diff --git a/MC/config/PWGUD/trigger/triggerDpmjetParticle.C b/MC/config/PWGUD/trigger/triggerDpmjetParticle.C index 5fe8e2bc8..9d3a99c00 100644 --- a/MC/config/PWGUD/trigger/triggerDpmjetParticle.C +++ b/MC/config/PWGUD/trigger/triggerDpmjetParticle.C @@ -34,7 +34,7 @@ o2::eventgen::Trigger triggerDstar(double rapidityMin = -1., double rapidityMax { return [rapidityMin, rapidityMax](const std::vector& particles) -> bool { for (const auto& particle : particles) { - if ((TMath::Abs(particle.GetPdgCode()) == 413) || (TMath::Abs(particle.GetPdgCode()) == 423)) + if (TMath::Abs(particle.GetPdgCode()) == 413) if ((particle.Y() > rapidityMin) && (particle.Y() < rapidityMax)) return kTRUE; }