Merge pull request #3578 from stahlleiton/STARlight_11_12_2023

Update Starlight generator to v326

bin/Starlight/gen_fragment.py

@@ -1,25 +1,26 @@
 import FWCore.ParameterSet.Config as cms
 
 externalLHEProducer = cms.EDProducer("ExternalLHEProducer",
-    args = cms.vstring(ADD_PATH_TO_TARBALL),
+    args = cms.vstring(),
     nEvents = cms.untracked.uint32(1),
     numberOfParameters = cms.uint32(1),
     outputFile = cms.string('cmsgrid_final.lhe'),
     scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh')
 )
 
-
 generator = cms.EDFilter("Pythia8HadronizerFilter",
-                         maxEventsToPrint = cms.untracked.int32(1),
-                         pythiaPylistVerbosity = cms.untracked.int32(1),
-                         filterEfficiency = cms.untracked.double(1.0),
-                         pythiaHepMCVerbosity = cms.untracked.bool(False),
-                         comEnergy = cms.double(ADD_BEAM_ENERGY),
-                         PythiaParameters = cms.PSet(
-                                parameterSets = cms.vstring('skip_hadronization'),
-                                skip_hadronization = cms.vstring('ProcessLevel:all = off',
-                                        'Check:event = off')
-                        )
+    PythiaParameters = cms.PSet(
+        skip_hadronization = cms.vstring(
+            'ProcessLevel:all = off',
+            'Check:event = off'
+        ),
+        parameterSets = cms.vstring('skip_hadronization')
+    ),
+    comEnergy = cms.double(5360.0),
+    filterEfficiency = cms.untracked.double(1.0),
+    maxEventsToPrint = cms.untracked.int32(1),
+    pythiaHepMCVerbosity = cms.untracked.bool(False),
+    pythiaPylistVerbosity = cms.untracked.int32(1)
 )
 
 ProductionFilterSequence = cms.Sequence(generator)

bin/Starlight/gridpack_generation.sh

@@ -4,7 +4,7 @@ create_setup(){
     echo "Using scram_arch: ${SCRAM_ARCH}"
     echo "Using cmssw_version: ${CMSSW_VERSION}"
 
-    GENDIR=${PRODDIR}/tarball
+    GENDIR=${PRODDIR}/tarball_$RANDOM
     [[ -d "${GENDIR}" ]] && rm -rf ${GENDIR}
     mkdir -p ${GENDIR} && cd ${GENDIR}
 
@@ -23,35 +23,33 @@ create_setup(){
 
 install_starlight(){
     DPMJET=dpmjet3.0-5
-    STARLIGHT=starlight_r317
+    STARLIGHT=STARlight-REV_326
     cd ${WORKDIR}
 
     echo "Downloading "${DPMJET}
     export DPMJETDIR=${WORKDIR}/dpmjet_v${DPMJET//[!0-9]/}
-    wget --no-verbose --no-check-certificate http://cms-project-generators.web.cern.ch/cms-project-generators/starlight/${DPMJET}.tar
+    wget --no-verbose --no-check-certificate https://cms-project-generators.web.cern.ch/cms-project-generators/starlight/${DPMJET}.tar
     tar -xf ${DPMJET}.tar && mv ${DPMJET} ${DPMJETDIR}
     rm -f ${DPMJET}.tar
 
     echo "Downloading "${STARLIGHT}
     STARLIGHT_VER=${STARLIGHT//[!0-9]/}
     STARLIGHTDIR=${WORKDIR}/starlight_v${STARLIGHT_VER}
-    wget --no-verbose --no-check-certificate http://cms-project-generators.web.cern.ch/cms-project-generators/starlight/${STARLIGHT}.tar
-    mkdir -p ${STARLIGHTDIR} && tar -xf ${STARLIGHT}.tar -C ${STARLIGHTDIR}
-    [[ -d "${STARLIGHTDIR}/trunk" ]] && mv ${STARLIGHTDIR}/trunk/* ${STARLIGHTDIR} && rm -rf ${STARLIGHTDIR}/trunk
-    rm -f ${STARLIGHT}.tar
+    wget --no-verbose --no-check-certificate https://cms-project-generators.web.cern.ch/cms-project-generators/starlight/${STARLIGHT}.tar.gz
+    tar -xzf ${STARLIGHT}.tar.gz && mv ${STARLIGHT} ${STARLIGHTDIR}
+    rm -f ${STARLIGHT}.tar.gz
 
     echo "Patching "${DPMJET}" and "${STARLIGHT}
     patch -ufZs -p1 -i ${PRODDIR}/patches/dpmjet.patch -d ${DPMJETDIR}
-    patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_default.patch -d ${STARLIGHTDIR}
     patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_pythia.patch -d ${STARLIGHTDIR}
-    if [ "$STARLIGHT_VER" -lt 317 ]; then
-        patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_randomgenerator.patch -d ${STARLIGHTDIR}
-    fi
+    patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_varnotused.patch -d ${STARLIGHTDIR}
+    patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_xsec.patch -d ${STARLIGHTDIR}
+    patch -ufZs -p1 -i ${PRODDIR}/patches/starlight_nuclearpar.patch -d ${STARLIGHTDIR}
 
     echo "Compiling ${DPMJET}"
     cd ${DPMJETDIR}
     rm -f fpe.o
-    make -j8
+    make -j $(nproc)
 
     echo "Compiling ${STARLIGHT}"
     cd ${STARLIGHTDIR}
@@ -60,12 +58,12 @@ install_starlight(){
     cp ${STARLIGHTDIR}/config/my.input ./
     export PYTHIADIR=$(scram tool tag pythia8 PYTHIA8_BASE)
     cmake ${STARLIGHTDIR} -DENABLE_DPMJET=ON -DENABLE_PYTHIA=ON
-    make -j8
+    make -j $(nproc)
 
     echo "Compiling macros"
     cp -r ${PRODDIR}/macros ${WORKDIR}/
     cd ${WORKDIR}/macros/
-    make -j8
+    make -j $(nproc)
 
     #Set installation parameters
     sed -i 's/SCRAM_ARCH_VERSION_REPLACE/'${SCRAM_ARCH}'/g' ${WORKDIR}/runcmsgrid.sh
@@ -76,13 +74,12 @@ install_starlight(){
 
 make_tarball(){
     #Set tarball name
-    CONFIG=${1}
-    prefix=${CONFIG##*/} ; prefix=${prefix%%.*} ; prefix=${prefix#*starlight_}
+    prefix=${INPUTFILE##*/} ; prefix=${prefix%%.*} ; prefix=${prefix#*starlight_}
     TARBALL=starlight_${prefix}_${SCRAM_ARCH}_${CMSSW_VERSION}_tarball.tgz
 
     echo "Creating tarball"
     cd ${WORKDIR}
-    cp ${CONFIG} ${STARLIGHTDIR}/build/slight.in
+    cp ${INPUTFILE} ${STARLIGHTDIR}/build/slight.in
     tar -czf ${TARBALL} ${STARLIGHTDIR##*/} ${DPMJETDIR##*/} macros runcmsgrid.sh
     mv ${WORKDIR}/${TARBALL} ${PRODDIR}/
     echo "Tarball created successfully at ${PRODDIR}/${TARBALL}"
@@ -118,6 +115,8 @@ else
         SCRAM_ARCH=slc7_amd64_gcc11
     elif [[ $SYSTEM_RELEASE == *"release 8"* ]]; then
         SCRAM_ARCH=el8_amd64_gcc11
+    elif [[ $SYSTEM_RELEASE == *"release 9"* ]]; then
+        SCRAM_ARCH=el9_amd64_gcc11
     else
         echo "No default scram_arch for current OS"
         exit 1
@@ -132,9 +131,11 @@ else
     if [[ $SYSTEM_RELEASE == *"release 6"* ]]; then
         CMSSW_VERSION=CMSSW_10_3_5
     elif [[ $SYSTEM_RELEASE == *"release 7"* ]]; then
-        CMSSW_VERSION=CMSSW_13_0_7
+        CMSSW_VERSION=CMSSW_13_0_17
     elif [[ $SYSTEM_RELEASE == *"release 8"* ]]; then
-        CMSSW_VERSION=CMSSW_13_0_7
+        CMSSW_VERSION=CMSSW_13_0_17
+    elif [[ $SYSTEM_RELEASE == *"release 9"* ]]; then
+        CMSSW_VERSION=CMSSW_13_0_17
     else
         echo "No default CMSSW for current OS"
         exit 1
@@ -149,7 +150,7 @@ create_setup
 install_starlight
 
 #Create tarball
-make_tarball ${INPUTFILE}
+make_tarball
 
 #Clean up
 echo "Removing "${GENDIR}

bin/Starlight/macros/Makefile

@@ -1,6 +1,6 @@
 # compiler flags:
 CFLAGS  = $(shell root-config --cflags)
-LIBS    = $(shell root-config --evelibs) -lGenVector
+LIBS    = -L$(shell root-config --libdir) -lEG -lCore
 
 # The build target
 TARGET = convert_SL2LHE

bin/Starlight/macros/convert_SL2LHE.cpp

@@ -36,6 +36,12 @@ void convert_SL2LHE(const std::string& inFileName, const double& beamE1, const d
   if (not inFile.is_open())
     throw std::logic_error("[convert_SL2LHE] Failed to open input file: "+inFileName);
 
+  // Extract cross section
+  double fidxsec(1), totxsec(3);
+  std::ifstream xsecFile(inFileName.substr(0, inFileName.rfind("/")+1)+"xsec.out");
+  if (xsecFile.is_open())
+    xsecFile >> fidxsec >> totxsec;
+
   // Create output file
   const auto outFileName = inFileName.substr(0, inFileName.rfind(".")).substr(inFileName.rfind("/")+1) + ".lhe";
   std::ofstream outFile(outFileName);
@@ -51,16 +57,18 @@ void convert_SL2LHE(const std::string& inFileName, const double& beamE1, const d
 
   // Put generic initialization-level info since STARLIGHT doesn't save this
   outFile << "<init>" << std::endl;
-  //beam particle 1, 2, beam energy 1, 2, author group, beam 1, 2, PDFSET beam 1, 2,
-  outFile << "22 " << "22 " << beamE1 << " " << beamE2 << " 0 " << "0 " << "0 " << "0 " << "3 " << "1" << std::endl;
-  outFile << "1.0 " << "0.0 " << "3.0 " << "81" << std::endl;
+  outFile << std::fixed << std::setprecision(8) << std::scientific;
+  //LHE format: https://arxiv.org/pdf/hep-ph/0109068.pdf
+  //beam pdg id (1, 2), beam energy [GeV] (1, 2), PDF author group (1, 2), PDF set id (1, 2), weight strategy, # subprocesses
+  outFile << "2212 2212 " << beamE1 << " " << beamE2 << " 0 0 0 0 3 1" << std::endl;
+  //cross section [pb], cross section stat. unc. [pb], maximum event weight, subprocess id
+  outFile << fidxsec << " " << 0.0 << " " << totxsec << " 81" << std::endl;
   outFile << "</init>" << std::endl;
 
   int nEvt(0);
   std::string line, label;
   const auto momPdgId = getMomPdgID(chnId);
   std::unique_ptr<TDatabasePDG> dataPDG(TDatabasePDG::Instance());
-  outFile.precision(10);
 
   // Read input file
   while (getline(inFile, line)) {
@@ -71,15 +79,19 @@ void convert_SL2LHE(const std::string& inFileName, const double& beamE1, const d
         label != "EVENT:")
       throw std::logic_error("[convert_SL2LHE] Failed to parse event line: "+line);
 
-    std::vector<std::string> oLines;
-    oLines.reserve(nTrk + nVtx);
+    // Particle tuple: pdg id, status, mother index, 4-momentum
+    std::vector<std::tuple<int, int, int, ROOT::Math::PxPyPzMVector>> parV;
+    parV.reserve(nVtx + nTrk);
+    double scale(-1);
 
     while (iVtx < nVtx && getline(inFile, line)) {
       // Read vertex line
       std::istringstream stream(line);
       if (not (stream >> label))
         throw std::logic_error("[convert_SL2LHE] Failed to parse line: "+line);
-      else if (label != "VERTEX:")
+      if (label == "GAMMAENERGIES:")
+        stream >> scale;
+      if (label != "VERTEX:")
         continue;
 
       double a;
@@ -88,10 +100,10 @@ void convert_SL2LHE(const std::string& inFileName, const double& beamE1, const d
       if (not (stream >> a >> a >> a >> a >> vtxN >> b >> iMom >> nDau) || vtxN != iVtx+1)
         throw std::logic_error("[convert_SL2LHE] Failed to parse vertex line: "+line);
 
-      ROOT::Math::PxPyPzMVector momP(0,0,0,0);
       const bool addMom = (momPdgId > 0 && iMom == 0);
       const auto& momI = addMom ? jMom : iMom;
-      std::vector<std::string> pLines(nDau + addMom);
+      if (addMom)
+        parV.emplace_back(momPdgId, 2, 0, ROOT::Math::PxPyPzMVector(0,0,0,0));
 
       // Read track lines
       for (int iDau=0; iDau<nDau; iDau++, iTrk++) {
@@ -109,33 +121,28 @@ void convert_SL2LHE(const std::string& inFileName, const double& beamE1, const d
         const auto& par = dataPDG->GetParticle(pdgId);
         if (not par)
           throw std::logic_error("[convert_SL2LHE] Invalid PDG ID for track line: "+line);
-        ROOT::Math::PxPyPzMVector dauP(px, py, pz, par->Mass());
+        parV.emplace_back(pdgId, 1, momI, ROOT::Math::PxPyPzMVector(px, py, pz, par->Mass()));
         if (addMom)
-          momP += dauP;
-
-        std::ostringstream oLine;
-        oLine << pdgId << " 1 " << momI << " " << momI << " 0 0 " << px << " " << py << " " << pz << " " << dauP.E() << " " << dauP.M() << " 0.0 9.0" << std::endl;
-        pLines[iDau + addMom] = oLine.str();
+          std::get<3>(parV[momI-1]) += std::get<3>(parV.back());
       }
 
-      if (addMom) {
-        std::ostringstream oLine;
-        oLine << momPdgId << " 2 0 0 0 0 " << momP.Px() << " " << momP.Py() << " " << momP.Pz() << " " << momP.E() << " " << momP.M() << " 0.0 9.0" << std::endl;
-        pLines[0] = oLine.str();
+      if (addMom)
         jMom += nDau+1;
-      }
-
-      oLines.insert(oLines.end(), pLines.begin(), pLines.end());
       iVtx++;
     }
-
     if (iTrk != nTrk)
       throw std::logic_error("[convert_SL2LHE] Failed to find all tracks in event");
 
     outFile << "<event>" << std::endl;
-    outFile << oLines.size() << " 81" << " 1.0 -1.0 -1.0 -1.0" << std::endl;
-    for (const auto& oLine : oLines)
-        outFile << oLine;
+    outFile << std::fixed << std::setprecision(8) << std::scientific;
+    //# particles, subprocess id, event weight, event scale, alpha_em, alpha_s
+    outFile << parV.size() << " 81 " << 1.0 << " " << scale << " " << -1.0 << " " << -1.0 << std::endl;
+    outFile << std::fixed << std::setprecision(10) << std::scientific;
+    for (const auto& pV : parV) {
+      const auto& [pdgId, status, momI, p] = pV;
+      //particle: pdg id, status, mother index (1, 2), color flow tag (1, 2), (px, py, pz, energy, mass [GeV]), proper lifetime [mm], spin
+      outFile << pdgId << " " << status << " " << momI << " 0 0 0 " << p.Px() << " " << p.Py() << " " << p.Pz() << " " << p.E() << " " << p.M() << " 0.0000e+00 9.0000e+00" << std::endl;
+    }
     outFile << "</event>" << std::endl;
     nEvt++;
   }