Merge branch 'AliceO2Group:master' into master

Author: shreyasiacharya

DATA/common/setenv_calib.sh

@@ -263,7 +263,10 @@ if [[ -z ${CALIBDATASPEC_BARREL_TF:-} ]]; then
   if [[ $CALIB_PRIMVTX_MEANVTX == 1 ]]; then add_semicolon_separated CALIBDATASPEC_BARREL_TF "pvtx:GLO/PVTX/0"; fi
 
   # ITS
-  if [[ $CALIB_ITS_DEADMAP_TIME == 1 ]]; then add_semicolon_separated CALIBDATASPEC_BARREL_TF "itsChipStatus:ITS/CHIPSSTATUS/0"; fi
+  if [[ $CALIB_ITS_DEADMAP_TIME == 1 ]]; then
+		add_semicolon_separated CALIBDATASPEC_BARREL_TF "itsChipStatus:ITS/CHIPSSTATUS/0"
+		add_semicolon_separated CALIBDATASPEC_BARREL_TF "itsErrorInfo:ITS/ErrorInfo/0"
+  fi
 
   # MFT
   if [[ $CALIB_MFT_DEADMAP_TIME == 1 ]]; then add_semicolon_separated CALIBDATASPEC_BARREL_TF "mftChipStatus:MFT/CHIPSSTATUS/0"; fi

GRID/utils/fetch_output_onfailure.sh

@@ -25,11 +25,11 @@ RecycleBase=""
 for ((i = 0; i < job_number; i++)); do
   jobid=${FAILEDSUBJOBIDS[i]}
   if [ ! "${RecycleBase}" ]; then
-    RecycleOutputDir=$(alien.py ps --trace ${jobid} | awk '/Going to uploadOutputFiles/' | sed 's/.*outputDir=//' | sed 's/)//')
+    RecycleOutputDir=$(alien.py ps --trace ${jobid} | awk '/Going to uploadOutputFiles/' | sed 's/.*outputDir=//' | sed 's/)//' | sort -u | head -n1)
     # /alice/cern.ch/user/a/aliprod/recycle/alien-job-2974093751
     RecycleBase=${RecycleOutputDir%-${jobid}} # Removes the ${jobid} and yields the recycle base path
   fi
-  $(alien.py registerOutput ${jobid}) 2> /dev/null
+  alien.py registerOutput ${jobid} > /dev/null 2>&1
 done
 
 # wait a bit to allow propagation of "registerOutput"

GRID/utils/querymasterjob.sh

@@ -29,5 +29,5 @@ erroredjobs=${tokens[8]}
 echo "MasterID $masterid"
 echo "TotalJobs $totaljobs"
 echo "DoneJobs $donejobs"
-echo "RunningJobs $donejobs"
+echo "RunningJobs $runningjobs"
 echo "ErroredJobs $erroredjobs"

MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C

@@ -170,7 +170,8 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
         mPythiaGun.moreDecays();
         std::array<int, 2> dausToCoal = {-1, -1};
         std::vector<int> pdgShortLivedResos = {313, 2224, 102134};
-        bool isResoFound{false};
+        std::map<int, int> statusResoDecay = {{313, 95}, {2224, 96}, {102134, 97}}; // do not use 94, it is used by default for no resonances
+        int whichReso{0};
         int idxCharmNucleus{-1};
         for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
           auto part = mPythiaGun.event[iPart];
@@ -180,17 +181,18 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
           }
           auto mother = part.mother1();
           // if we find a resonance, we remove it, otherwise we prevent the coalescence of daughters from resonances and daughters from charmed nucleus directly
-          if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end() && mother >= idxCharmNucleus) {
+          auto resoIt = std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg);
+          if (resoIt != pdgShortLivedResos.end() && mother >= idxCharmNucleus) {
             // we need to change the indices of the daughter particles to point to the charmed nucleus
             auto dauList = part.daughterList();
             for (auto const& dau : dauList) {
               mPythiaGun.event[dau].mother1(idxCharmNucleus);
             }
             mPythiaGun.event.remove(iPart, iPart, true);
-            isResoFound=true;
+            whichReso=*resoIt;
           }
         }
-        if (isResoFound) { // we have to reset all the particles as daughters of the charm nucleus
+        if (whichReso > 0) { // we have to reset all the particles as daughters of the charm nucleus
           std::vector<int> idxDausCharmNucleus{};
           for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
             auto mother = mPythiaGun.event[iPart].mother1();
@@ -242,6 +244,9 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
 
         // we try the coalescence here, if successful we copy particles in the pythia event and we move to the next charm nucleus
         isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<unsigned int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.);
+        if (whichReso > 0) {
+          mPythiaGun.event[idxCharmNucleus].status(statusResoDecay[whichReso]);
+        }
         if (isCoalSuccess) {
           restoreEnergyConservation(mPythiaGun.event, idxCharmNucleus);
           int offset = mPythia.event.size(); // we need to rescale the indices of mothers and daughters, accounting for the particles that are already appended to the event

MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C

@@ -1,7 +1,12 @@
+R__LOAD_LIBRARY(EvtGen)
+R__ADD_INCLUDE_PATH($EVTGEN_ROOT/include)
+
 #include "FairGenerator.h"
 #include "Generators/GeneratorPythia8.h"
 #include "Generators/GeneratorPythia8Param.h"
+#include "EvtGen/EvtGen.hh"
 #include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/EvtGen.h"
 #include "TRandom.h"
 #include "TDatabasePDG.h"
 #include <fairlogger/Logger.h>
@@ -33,6 +38,9 @@ public:
     mHadronPdgList = hadronPdgList;
     mPartPdgToReplaceList = partPdgToReplaceList;
     mFreqReplaceList = freqReplaceList;
+    mEvtGen = nullptr;
+    mUseEvtGen = false;
+    mEvtGenDecTable = "";
     // Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0, LambdaC(2625), LambdaC(2595), LambdaC(2860), LambdaC(2880), LambdaC(2940), ThetaC(3100)
     mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326, 4124, 14122, 24124, 24126, 4125, 9422111};
     mCustomPartMasses[30433] = 2.714f;
@@ -114,6 +122,10 @@ public:
       pdgToReplace.push_back(mPartPdgToReplaceList[iRepl].at(0));
     }
 
+    if (mUseEvtGen) {
+      mEvtGen = new Pythia8::EvtGenDecays(&mPythia, mEvtGenDecTable.data(), gSystem->ExpandPathName("$EVTGEN_ROOT/share/EvtGen/evt.pdl"));
+    }
+
     return o2::eventgen::GeneratorPythia8::Init();
   }
 
@@ -135,6 +147,14 @@ public:
   {
     return mUsedSeed;
   };
+  void setUseEvtGenDecayer(std::string evtGenDecTable = "") {
+    mUseEvtGen = true;
+    if (evtGenDecTable.empty()) {
+      mEvtGenDecTable = gSystem->ExpandPathName("$EVTGEN_ROOT/share/EvtGen/DECAY.DEC");
+    } else {
+      mEvtGenDecTable = gSystem->ExpandPathName(evtGenDecTable.data());
+    }
+  }
 
 protected:
   //__________________________________________________________________
@@ -167,6 +187,10 @@ protected:
       {
         if (GeneratorPythia8::generateEvent())
         {
+          if (mUseEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
           genOk = selectEvent();
         }
       }
@@ -179,6 +203,12 @@ protected:
       while (!genOk)
       {
         genOk = GeneratorPythia8::generateEvent();
+        if (genOk) {
+          if (mUseEvtGen) {
+            mEvtGen->decay();
+            checkConsistency();
+          }
+        }
       }
       notifySubGenerator(0);
     }
@@ -197,6 +227,8 @@ protected:
 
     for (auto iPart{0}; iPart < mPythia.event.size(); ++iPart)
     {
+
+
       // search for Q-Qbar mother with at least one Q in rapidity window
       if (!isGoodAtPartonLevel)
       {
@@ -352,11 +384,22 @@ protected:
     return true;
   }
 
+  void checkConsistency() {
+    for (int iPart{1}; iPart<mPythia.event.size(); ++iPart) {
+      // verify if all particles of decay chain are in the TDatabasePDG
+      // taken from https://github.com/AliceO2Group/O2DPG/blob/master/MC/config/PWGDQ/EvtGen/GeneratorEvtGen.C
+      if (!TDatabasePDG::Instance()->GetParticle(abs(mPythia.event[iPart].id()))) {
+        // std::cout << "Particle code non known in TDatabasePDG " << mPythia.event[iPart].id() << " - set pdg = 89" << std::endl;
+        mPythia.event[iPart].id(89);
+      }
+    }
+  }
+
   private:
   // Interface to override import particles
   Pythia8::Event mOutputEvent;
 
-  // Properties of selection
+// Properties of selection
   int mQuarkPdg;
   float mQuarkRapidityMin;
   float mQuarkRapidityMax;
@@ -379,6 +422,12 @@ protected:
 
   // Control alternate trigger on different hadrons
   std::vector<int> mHadronPdgList = {};
+
+  // EVTGEN decayer from PYTHIA8 plugins
+  Pythia8::EvtGenDecays *mEvtGen;
+  bool mUseEvtGen;
+  std::string mEvtGenDecTable;
+
 };
 
 // Predefined generators:
@@ -446,3 +495,20 @@ FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1
 
   return myGen;
 }
+
+// Beauty-enriched with EvtGen decayer
+FairGenerator *GeneratorPythia8GapTriggeredBeautyWithEvtGen(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {}, std::string decayTable = "")
+{
+  auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{5}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
+  auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+  myGen->setUsedSeed(seed);
+  myGen->readString("Random:setSeed on");
+  myGen->readString("Random:seed " + std::to_string(seed));
+  myGen->setQuarkRapidity(yQuarkMin, yQuarkMax);
+  if (hadronPdgList.size() != 0)
+  {
+    myGen->setHadronRapidity(yHadronMin, yHadronMax);
+  }
+  myGen->setUseEvtGenDecayer(decayTable);
+  return myGen;
+}

MC/config/PWGHF/ini/GeneratorHFTrigger_Bforced.ini

@@ -1,4 +1,4 @@
-#NEV_TEST> 20
+#NEV_TEST> 50
 [GeneratorExternal]
 fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
 funcName=GeneratorPythia8GapTriggeredBeauty(3, -5, 5)

MC/config/PWGHF/ini/GeneratorHF_D2H_bbbar_B2JPsi_gap5_Mode2.ini

@@ -0,0 +1,9 @@
+#NEV_TEST> 20
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredBeautyWithEvtGen(5, -1.5, 1.5, -1.5, 1.5, {}, {}, {}, "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/decayer/EVTGEN_DECAY_B2JPSI.DEC")
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_Mode2_noforcedecays.cfg
+includePartonEvent=false

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.ini

@@ -1,4 +1,4 @@
-#NEV_TEST> 40
+#NEV_TEST> 100
 ### The external generator derives from GeneratorPythia8.
 [GeneratorExternal]
 fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C

MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_B2JPsi_gap5_Mode2.C

@@ -0,0 +1,113 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuark{5};
+    float ratioTrigger = 1./5; // one event triggered out of 5
+
+    std::vector<int> checkPdgHadron{443, 511, 521, 531, 5122};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {443, {{-11, 11}, {-13, 13}}}, // J/psi
+        {511, {{313, 443}}}, // B0
+        {521, {{321, 443}}}, // B+
+        {531, {{333, 443}}}, // Bs0
+        {5122, {{321, 443, 2212}}} // Lb0
+    };
+
+    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");
+    std::vector<o2::MCTrack> *tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+    o2::dataformats::MCEventHeader *eventHeader = nullptr;
+    tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+
+    int nEventsMB{}, nEventsInj{};
+    int nSignals{}, nSignalGoodDecay{};
+    int nJPsiToEE{}, nJPsiToMuMu{};
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+
+        // check subgenerator information
+        if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
+            bool isValid = false;
+            int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
+            if (subGeneratorId == 0) {
+                nEventsMB++;
+            } else if (subGeneratorId == checkPdgQuark) {
+                nEventsInj++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
+                nSignals++; // count signal PDG
+
+                std::vector<int> pdgsDecay{};
+                std::vector<int> pdgsDecayAntiPart{};
+                for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                    auto pdgDau = tracks->at(j).GetPdgCode();
+                    if (pdg == 443 && pdgDau == 22) {
+                        continue;
+                    }
+                    pdgsDecay.push_back(pdgDau);
+                    if (pdgDau != 333) { // phi is antiparticle of itself
+                        pdgsDecayAntiPart.push_back(-pdgDau);
+                    } else {
+                        pdgsDecayAntiPart.push_back(pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+                std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
+
+                for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
+                    if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
+                        nSignalGoodDecay++;
+                        if (pdg == 443) {
+                            if (decay == std::vector{-11, 11}) {
+                                nJPsiToEE++;
+                            } else if (decay == std::vector{-13, 13}) {
+                                nJPsiToMuMu++;
+                            }
+                        }
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# MB events: " << nEventsMB << "\n";
+    std::cout << Form("# events injected with %d quark pair: ", checkPdgQuark) << nEventsInj << "\n";
+    std::cout <<"# signal hadrons: " << nSignals << "\n";
+    std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
+    std::cout <<"# J/Psi decaying to e+e-: " << nJPsiToEE << "\n";
+    std::cout <<"# J/Psi decaying to mu+mu-: " << nJPsiToMuMu << "\n";
+
+    if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small
+        std::cerr << "Number of generated MB events different than expected\n";
+        return 1;
+    }
+    if (nEventsInj < nEvents * ratioTrigger * 0.95 || nEventsInj > nEvents * ratioTrigger * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuark << " different than expected\n";
+        return 1;
+    }
+
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 1.0f;
+    if (1 - fracForcedDecays > 0.5 + uncFracForcedDecays) { // at least 50% in the main decay channels (we also have correlated backgrounds, mostly from other B decays)
+        std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
+        return 1;
+    }
+
+    return 0;
+}