DPMJET and other new processes in STARlight (#1765)

* DPMJET and other new processes * Changing suffix since this is not a normal ini file for o2sim
AliceO2Group · Nov 5, 2024 · 039a4e0 · 039a4e0
1 parent 5869eb1
commit 039a4e0
Show file tree

Hide file tree

Showing 7 changed files with 130 additions and 23 deletions.
diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/TAUTAU.EL3PI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/TAUTAU.EL3PI.DEC
@@ -0,0 +1,7 @@
+Decay tau-
+1.0 e-      anti-nu_e nu_tau                        PHOTOS TAULNUNU; #[Reconstructed PDG2011]
+Enddecay
+Decay tau+
+1.0 pi+     pi+     pi-     anti-nu_tau                  TAUHADNU -0.108 0.775 0.149 1.364 0.400 1.23 0.4; #[Reconstructed PDG2011]
+Enddecay
+End
diff --git a/MC/config/PWGUD/external/generator/DecayTablesEvtGen/TAUTAU.PO3PI.DEC b/MC/config/PWGUD/external/generator/DecayTablesEvtGen/TAUTAU.PO3PI.DEC
@@ -0,0 +1,7 @@
+Decay tau-
+1.0 pi-     pi-     pi+     nu_tau                  TAUHADNU -0.108 0.775 0.149 1.364 0.400 1.23 0.4; #[Reconstructed PDG2011]
+Enddecay
+Decay tau+
+1.0 e+      anti-nu_tau nu_e                        PHOTOS TAULNUNU; #[Reconstructed PDG2011]
+Enddecay
+End
diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlight.C b/MC/config/PWGUD/external/generator/GeneratorStarlight.C
@@ -1,8 +1,11 @@
+R__LOAD_LIBRARY(libDPMJET.so)
+R__LOAD_LIBRARY(libDpmJetLib.so)
 R__LOAD_LIBRARY(libStarlib.so)
 R__ADD_INCLUDE_PATH($STARlight_ROOT/include)
 
 #include "randomgenerator.h"
 #include "upcXevent.h"
+#include "upcevent.h"
 #include "starlight.h"
 #include "inputParameters.h"
 
@@ -96,7 +99,14 @@ class GeneratorStarlight_class : public Generator
     {"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 }, //	
+    {"kIncohUpsilonToEl",    4,  553011,   20, -1.0, -1.0, 0.01, 553, 0 }, //
+    {"kDpmjetSingle",        5,  113,   20, -1.0, -1.0, 0.01, -1, 0 }, //
+    {"kTauLowToEl3Pi",       1,      15,  292,  0.4, 15.0, 0.01, -1, 1 }, // from 0.4 to 15 GeV
+    {"kTauLowToPo3Pi",       1,      15,  292,  0.4, 15.0, 0.01, -1, 1 }, // from 0.4 to 15 GeV
+    {"kTauMediumToEl3Pi",    1,      15,  264,  1.8, 15.0, 0.01, -1, 1 }, // from 1.8 to 15 GeV
+    {"kTauMediumToPo3Pi",    1,      15,  264,  1.8, 15.0, 0.01, -1, 1 }, // from 1.8 to 15 GeV
+    {"kTauHighToEl3Pi",      1,      15,  220,  4.0, 15.0, 0.01, -1, 1 }, // from 4.0 to 15 GeV
+    {"kTauHighToPo3Pi",      1,      15,  220,  4.0, 15.0, 0.01, -1, 1 }, // from 4.0 to 15 GeV
 	};
 
   const int nProcess = sizeof(slConfig)/sizeof(SLConfig);
@@ -146,10 +156,16 @@ class GeneratorStarlight_class : public Generator
   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("XSEC_METHOD   = 0      # Set to 0 to use old method for calculating gamma-gamma luminosity"); //CM
+  setParameter("BSLOPE_DEFINITION = 1");   // 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
+
+  // Photonuclear specific options, energies in Lab frame. These values should be within the range of the values specified in the DPMJet input file (when DPMJet is used)
+  if(slConfig[idx].prod_mode == 5 || slConfig[idx].prod_mode == 6 || slConfig[idx].prod_mode == 7){
+    setParameter("MIN_GAMMA_ENERGY = 1000.0");
+    setParameter("MAX_GAMMA_ENERGY = 600000.0");
+  }
 
   if (not mInputParameters.init()) {
       std::cout << "InitStarLight parameter initialization has failed" << std::endl;
@@ -171,6 +187,11 @@ class GeneratorStarlight_class : public Generator
     return false;
   }
 
+  if (mInputParameters.interactionType() >= 5) {
+    mUpcEvent = mStarLight->produceUpcEvent();
+    mUpcEvent.boost(0.5*(TMath::ACosH(mInputParameters.beam1LorentzGamma()) - TMath::ACosH(mInputParameters.beam2LorentzGamma())));
+  }
+
   mEvent = mStarLight->produceEvent();
   // boost event to the experiment CM frame
   mEvent.boost(0.5*(TMath::ACosH(mInputParameters.beam1LorentzGamma()) - TMath::ACosH(mInputParameters.beam2LorentzGamma())));
@@ -185,17 +206,28 @@ class GeneratorStarlight_class : public Generator
   {    
   int nVtx(0);
   float vtx(0), vty(0), vtz(0), vtt(0);
-  const std::vector<vector3>* slVtx(mEvent.getVertices());
+  const std::vector<vector3>* slVtx;
+  const std::vector<starlightParticle>* slPartArr;
+  int npart = 0;
+
+  if (mInputParameters.interactionType() >= 5) {
+    slVtx = mUpcEvent.getVertices();
+    slPartArr = mUpcEvent.getParticles();
+    npart = mUpcEvent.getParticles()->size();
+  }
+  else{
+    slVtx = mEvent.getVertices();
+    slPartArr = mEvent.getParticles();
+    npart = mEvent.getParticles()->size();
+  }
+
   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<starlightParticle>* slPartArr(mEvent.getParticles());
-  const int npart(mEvent.getParticles()->size());
-
+    }
+  else { // a vertex exits
+      nVtx = slVtx->size();
+    } // end if
+
   if(mPdgMother != -1){ //Reconstruct mother particle for VM processes
   TLorentzVector lmoth;
   TLorentzVector ldaug;
@@ -219,7 +251,7 @@ class GeneratorStarlight_class : public Generator
 	  mParticles.push_back(particle);
 	  o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), false);
   }
-  if(!mDecayEvtGen){ // Don't import daughters in case of external decayer
+  if(!mDecayEvtGen || mPdgMother == -1){ // Don't import daughters in case of external decayer
   for(int ipart=0;ipart<npart;ipart++) {
     const starlightParticle* slPart(&(slPartArr->at(ipart)));
       if (nVtx < 1) { // No verticies
@@ -232,7 +264,7 @@ class GeneratorStarlight_class : public Generator
       } // end if
       TParticle particle(slPart->getPdgCode(),
 				   1,
-				   0,
+				   (mPdgMother != -1 ? 0 :-1),
 				   -1,
 				   slPart->getFirstDaughter(),
 				   slPart->getLastDaughter(),
@@ -245,7 +277,7 @@ class GeneratorStarlight_class : public Generator
 	  mParticles.push_back(particle);
       o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(mParticles.back(), true);
     }  
-    }
+  }
   return true; 
   }
 
@@ -254,6 +286,7 @@ class GeneratorStarlight_class : public Generator
    inputParameters  mInputParameters;  //   simulation input information.
    randomGenerator  mRandomGenerator;  //   STARLIGHT's own random generator
    upcXEvent        mEvent;            //  object holding STARlight simulated event.
+   upcEvent         mUpcEvent;
    std::string mSelectedConfiguration = "";
    int mPdgMother = -1;
    bool mDecayEvtGen = 0;

diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C b/MC/config/PWGUD/external/generator/GeneratorStarlightToEvtGen.C
@@ -16,18 +16,22 @@ FairGenerator*
   gen->selectConfiguration(configuration);
   gen->setCollisionSystem(energyCM, beam1Z, beam1A, beam2Z, beam2A);
 
-  gen->SetSizePdg(3);
+  gen->SetSizePdg(5);
   gen->AddPdg(443,0);
   gen->AddPdg(100443,1);
   gen->AddPdg(223,2);
-  gen->SetPolarization(1); //Transversal
+  gen->AddPdg(15,3);
+  gen->AddPdg(-15,4);
+  if (configuration.find("kTau") == std::string::npos) 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()));
+  else if (configuration.find("ToEl3Pi") != std::string::npos) gen->SetDecayTable(Form("%s/TAUTAU.EL3PI.DEC",pathO2.Data()));
+  else if (configuration.find("ToPo3Pi") != std::string::npos) gen->SetDecayTable(Form("%s/TAUTAU.PO3PI.DEC",pathO2.Data()));
 
   return gen;
-}
+}
diff --git a/MC/config/PWGUD/ini/dpmjet_PbPb5360.txt b/MC/config/PWGUD/ini/dpmjet_PbPb5360.txt
@@ -0,0 +1,45 @@
+**********************************************************************
+* Example for a DTUNUC input file.
+* Uncomment the input-cards according to your requirements.
+*
+* Format: A10,6E10.0,A8
+*        (except for the section enclosed by "PHOINPUT" and "ENDINPUT"
+*         which is format-free)
+*         lines starting with "*" are comment lines
+**********************************************************************
+*
+* projectile / target / Energy
+* -------------------
+*        1         2         3         4         5         6         7
+*23456789012345678901234567890123456789012345678901234567890123456789012345678
+PROJPAR          0.0                                                  PHOTON
+TARPAR         208.0      82.0
+ENERGY        1000.0  600000.0
+*ENERGY         100.0
+* Initialize the random number generator
+RNDMINIT        55.0     101.0      15.0      73.0              
+*
+*
+* PHOJET-specific input
+* ---------------------
+* The following lines control the event-generation with PHOJET for
+* individual photon/nucleon-nucleon collisions.
+* For details see the PHOJET-manual available at
+*        http://lepton.bartol.udel.edu/~eng/phojet.html
+* Any options explained in the PHOJET-manual can be used in between
+* the "PHOINPUT" and "ENDINPUT" cards.
+PHOINPUT
+PROCESS           1 0 1 1 1 1 1 1
+ENDINPUT
+*
+* Output
+* ------
+*   some default output (particle multiplicities etc.)
+HISTOGRAM      101.0     102.0
+*
+* Start of event generation
+* -------------------------
+*START         5000.0       0.0
+START         100.0       0.0
+STOP
+*...+....1....+....2....+....3....+....4....+....5....+....6....+....7...
diff --git a/MC/config/PWGUD/ini/makeStarlightConfig.py b/MC/config/PWGUD/ini/makeStarlightConfig.py
@@ -17,7 +17,7 @@
 parser.add_argument('--rapidity', default='cent', choices=['cent_rap', 'muon_rap', 'cent_eta', 'muon_eta'],
                     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'],
+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', 'kDpmjetSingle',  'kTauLowToEl3Pi', 'kTauLowToPo3Pi', 'kTauMediumToEl3Pi', 'kTauMediumToPo3Pi' ,'kTauHighToEl3Pi' ,'kTauHighToPo3Pi'],
                     help='Process to switch on')
 
 
@@ -68,7 +68,7 @@
 
 ### 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 :
+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 or 'kTau' 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 ,pZ,pA,tZ,tA))
 else:
@@ -78,7 +78,7 @@
 ###Trigger
 fout.write('[TriggerExternal] \n')
 fout.write('fileName = ${O2DPG_ROOT}/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C \n')
-if 'kTwoGamma' in args.process:
+if 'kTwoGamma' in args.process or 'kTau' in args.process:
     if args.rapidity == 'cent_eta':
         fout.write('funcName = selectDirectPartInAcc(-0.9,0.9) \n')
     if args.rapidity == 'muon_eta':

diff --git a/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C b/MC/config/PWGUD/trigger/selectParticlesInAcceptance.C
@@ -24,13 +24,24 @@ o2::eventgen::Trigger selectDaughterPartInAcc(double etaMin = -1., double etaMax
     for (const auto& particle : particles) {
       if (particle.GetFirstMother() == -1)
         if ((particle.Y() < etaMin) || (particle.Y() > etaMax)) return kFALSE;
-	  if (particle.GetFirstMother() != -1 && particle.GetFirstDaughter() == -1 && particle.GetPdgCode() != 22)
+	  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)
 		  if ((particle.Eta() < etaMin) || (particle.Eta() > etaMax)) return kFALSE; 
     }
     return kTRUE;  
   };
 }
 
+o2::eventgen::Trigger selectDileptonsInAcc(double etaMin = -1., double etaMax = -1.)
+{
+  return [etaMin, etaMax](const std::vector<TParticle>& particles) -> bool {
+    for (const auto& particle : particles) {
+	  if (particle.GetFirstMother() != -1 && particle.GetFirstDaughter() == -1 && (TMath::Abs(particle.GetPdgCode()) == 11 || TMath::Abs(particle.GetPdgCode() == 13)))
+		  if ((particle.Eta() < etaMin) || (particle.Eta() > etaMax)) return kFALSE;
+    }
+    return kTRUE;
+  };
+}
+
 o2::eventgen::Trigger selectDirectPartInAcc(double etaMin = -1., double etaMax = -1.)
 {
   return [etaMin, etaMax](const std::vector<TParticle>& particles) -> bool {
@@ -41,4 +52,4 @@ o2::eventgen::Trigger selectDirectPartInAcc(double etaMin = -1., double etaMax =
     }
     return kTRUE;  
   };
-}
+}