Tuning of the dpmjet decays and inclusion of neutron generator

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

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;
 }
 

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<o2::eventgen::GeneratorStarlight_class>();
   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;
 }

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; i<mNeutrons->GetEntriesFast(); 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;
+}

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:

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<TParticle>& 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<TParticle>& 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)

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<TParticle>& 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;
     }