Merge branch 'AliceO2Group:master' into master

Author: ljoergen

Common/TableProducer/fwdtrackPropagation.cxx

@@ -61,6 +61,7 @@ struct FwdTrackPropagation {
   Configurable<float> maxEtaSA{"maxEtaSA", -2.5, "max. eta acceptance for MCH-MID"};
   Configurable<float> minEtaGL{"minEtaGL", -3.6, "min. eta acceptance for MFT-MCH-MID"};
   Configurable<float> maxEtaGL{"maxEtaGL", -2.5, "max. eta acceptance for MFT-MCH-MID"};
+  Configurable<float> minRabsGL{"minRabsGL", 27.6, "min. R at absorber end for global muon (min. eta = -3.6)"}; // std::tan(2.f * std::atan(std::exp(- -3.6)) ) * -505.
   Configurable<float> minRabs{"minRabs", 17.6, "min. R at absorber end"};
   Configurable<float> midRabs{"midRabs", 26.5, "middle R at absorber end for pDCA cut"};
   Configurable<float> maxRabs{"maxRabs", 89.5, "max. R at absorber end"};
@@ -69,11 +70,8 @@ struct FwdTrackPropagation {
   Configurable<float> maxPDCAforSmallR{"maxPDCAforSmallR", 594.f, "max. pDCA for small R at absorber end"};
   Configurable<float> maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 50.f, "max. chi2 for MCH-MFT matching"};
   Configurable<float> maxChi2SA{"maxChi2SA", 1e+6, "max. chi2 for standalone muon"};
-  Configurable<float> maxChi2GL{"maxChi2GL", 1e+6, "max. chi2 for global muon"};
+  Configurable<float> maxChi2GL{"maxChi2GL", 50.f, "max. chi2 for global muon"};
   Configurable<bool> refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"};
-  Configurable<bool> applyDEtaDPhi{"cfgApplyDEtaDPhi", false, "flag to apply deta-dphi elliptic cut"};
-  Configurable<float> minDEta{"minDEta", 0.1, "min deta between MFT-MCH-MID and its attached MID-MCH at PV"};
-  Configurable<float> minDPhi{"minDPhi", 0.1, "min dphi between MFT-MCH-MID and its attached MID-MCH at PV"};
 
   HistogramRegistry fRegistry{"fRegistry"};
   static constexpr std::string_view muon_types[5] = {"MFTMCHMID/", "MFTMCHMIDOtherMatch/", "MFTMCH/", "MCHMID/", "MCH/"};
@@ -175,6 +173,9 @@ struct FwdTrackPropagation {
       if (chi2 < 0.f || maxChi2GL < chi2) {
         return false;
       }
+      if (rAtAbsorberEnd < minRabsGL || maxRabs < rAtAbsorberEnd) {
+        return false;
+      }
     } else if (trackType == static_cast<uint8_t>(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) {
       if (eta < minEtaSA || maxEtaSA < eta) {
         return false;
@@ -189,89 +190,16 @@ struct FwdTrackPropagation {
     return true;
   }
 
-  template <typename TFwdTracks, typename TMFTTracks, typename TCollision, typename TTarget, typename TCandidates>
-  bool isBestMatch(TCollision const& collision, TTarget const& target, TCandidates const& candidates)
+  template <typename TCollision, typename TFwdTrack, typename TFwdTracks, typename TMFTTracks>
+  void fillFwdTrackTable(TCollision const& collision, TFwdTrack fwdtrack, TFwdTracks const&, TMFTTracks const&, const bool isAmbiguous)
   {
-    std::map<int64_t, float> map_chi2MFTMCH;
-    for (const auto& matchedtrack : candidates) { // MFT-MCH-MID or MFT-MCH
-      if (matchedtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) {
-        continue;
-      }
-      o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(matchedtrack, collision, propagationPoint::kToVertex);
-      float eta = propmuonAtPV.getEta();
-      float phi = propmuonAtPV.getPhi();
-      if (refitGlobalMuon) {
-        const auto& mfttrack = matchedtrack.template matchMFTTrack_as<TMFTTracks>();
-        eta = mfttrack.eta();
-        phi = mfttrack.phi();
-      }
-      o2::math_utils::bringTo02Pi(phi);
-      if (eta < minEtaGL || maxEtaGL < eta) {
-        continue;
-      }
-
-      const auto& mchtrack = matchedtrack.template matchMCHTrack_as<TFwdTracks>(); // MCH-MID
-      o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToVertex);
-      float etaMatchedMCHMID = propmuonAtPV_Matched.getEta();
-      float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi();
-      o2::math_utils::bringTo02Pi(phiMatchedMCHMID);
-
-      float deta = etaMatchedMCHMID - eta;
-      float dphi = phiMatchedMCHMID - phi;
-      o2::math_utils::bringToPMPi(dphi);
-      if (applyDEtaDPhi && std::sqrt(std::pow(deta / minDEta, 2) + std::pow(dphi / minDPhi, 2)) > 1.f) {
-        continue;
-      }
-
-      if (matchedtrack.chi2() < 0.f || maxChi2GL < matchedtrack.chi2()) {
-        continue;
-      }
-
-      float rAtAbsorberEnd = matchedtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack
-      if (rAtAbsorberEnd < minRabs || maxRabs < rAtAbsorberEnd) {
-        continue;
-      }
-      o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(matchedtrack, collision, propagationPoint::kToDCA);
-      float dcaX = propmuonAtDCA.getX() - collision.posX();
-      float dcaY = propmuonAtDCA.getY() - collision.posY();
-      float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
-      if (maxDCAxy < dcaXY) {
-        continue;
-      }
-
-      o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToDCA);
-      float dcaX_Matched = propmuonAtDCA_Matched.getX() - collision.posX();
-      float dcaY_Matched = propmuonAtDCA_Matched.getY() - collision.posY();
-      float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched);
-      float pDCA = mchtrack.p() * dcaXY_Matched;
-
-      if (rAtAbsorberEnd < midRabs ? pDCA > maxPDCAforSmallR : pDCA > maxPDCAforLargeR) {
-        continue;
-      }
-
-      map_chi2MFTMCH[matchedtrack.globalIndex()] = matchedtrack.chi2MatchMCHMFT();
-    }
-    if (map_chi2MFTMCH.begin()->first != target.globalIndex()) { // search for minimum matching chi2
-      map_chi2MFTMCH.clear();
-      return false;
-    }
-    map_chi2MFTMCH.clear();
+    if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && (fwdtrack.chi2MatchMCHMFT() > maxMatchingChi2MCHMFT || fwdtrack.chi2() > maxChi2GL)) {
+      return;
+    } // Users have to decide the best match between MFT and MCH-MID at analysis level. The same global muon is repeatedly stored.
 
-    if (target.chi2MatchMCHMFT() > maxMatchingChi2MCHMFT) {
-      return false;
+    if (fwdtrack.chi2MatchMCHMID() < 0.f) { // this should never happen. only for protection.
+      return;
     }
-    return true;
-  }
-
-  template <typename TCollision, typename TFwdTrack, typename TFwdTracks, typename TMFTTracks>
-  void fillFwdTrackTable(TCollision const& collision, TFwdTrack fwdtrack, TFwdTracks const& fwdtracks, TMFTTracks const&, const bool isAmbiguous)
-  {
-    if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) {
-      const auto& matchedGlobalTracks = fwdtracks.sliceBy(perMFTTrack, fwdtrack.matchMFTTrackId()); // MFT-MCH-MID or MFT-MCH
-      if (!isBestMatch<TFwdTracks, TMFTTracks>(collision, fwdtrack, matchedGlobalTracks)) {
-        return;
-      }
-    } // find the best match between MFT and MCH-MID
 
     o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, collision, propagationPoint::kToVertex);
     o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, collision, propagationPoint::kToDCA);

EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx

@@ -208,6 +208,12 @@ struct HeavyNeutralMesonFilter {
     false,
     "Evt sel: check for offline selection"};
 
+  Configurable<bool> ConfDoEMCShift{"ConfDoEMCShift", false, "Apply SM-wise shift in eta and phi to EMCal clusters to align with TPC tracks"};
+  Configurable<std::vector<float>> ConfEMCEtaShift{"ConfEMCEtaShift", {0.f}, "values for SM-wise shift in eta to be added to EMCal clusters to align with TPC tracks"};
+  Configurable<std::vector<float>> ConfEMCPhiShift{"ConfEMCPhiShift", {0.f}, "values for SM-wise shift in phi to be added to EMCal clusters to align with TPC tracks"};
+  std::array<float, 20> EMCEtaShift = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::array<float, 20> EMCPhiShift = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
   template <typename T>
   bool isSelectedTrack(T const& track, CFTrigger::FemtoPartners partSpecies)
   {
@@ -760,6 +766,14 @@ struct HeavyNeutralMesonFilter {
     mHistManager.add("etaprimep/fetaprimePtVskstar_EMC", "Same Event distribution", HistType::kTH1F, {{8000, 0, 8}});
     mHistManager.add("etaprimep/fProtonPtVskstar_EMC", "Same Event distribution", HistType::kTH1F, {{8000, 0, 8}});
     mHistManager.add("etaprimep/fAntiProtonPtVskstar_EMC", "Same Event distribution", HistType::kTH1F, {{8000, 0, 8}});
+
+    if (ConfDoEMCShift.value) {
+      for (unsigned short iSM = 0; iSM < 20; iSM++) {
+        EMCEtaShift[iSM] = ConfEMCEtaShift.value[iSM];
+        EMCPhiShift[iSM] = ConfEMCPhiShift.value[iSM];
+        LOG(info) << "SM-wise shift in eta/phi for SM " << iSM << ": " << EMCEtaShift[iSM] << " / " << EMCPhiShift[iSM];
+      }
+    }
   }
   Preslice<aod::V0PhotonsKF> perCollision_pcm = aod::v0photonkf::collisionId;
   Preslice<aod::SkimEMCClusters> perCollision_emc = aod::skimmedcluster::collisionId;
@@ -817,10 +831,11 @@ struct HeavyNeutralMesonFilter {
     mHistManager.fill(HIST("Event/nClustersVsV0s"), clustersInThisCollision.size(), v0sInThisCollision.size());
     mHistManager.fill(HIST("Event/nTracks"), tracksWithItsPid.size());
 
-    hnmutilities::reconstructGGs(clustersInThisCollision, v0sInThisCollision, vGGs);
+    std::vector<hnmutilities::Photon> vGammas;
+    hnmutilities::storeGammasInVector(clustersInThisCollision, v0sInThisCollision, vGammas, EMCEtaShift, EMCPhiShift);
+    hnmutilities::reconstructGGs(vGammas, vGGs);
+    vGammas.clear();
     processGGs(vGGs);
-    // hnmutilities::reconstructHeavyNeutralMesons(tracks, vGGs, vHNMs);
-    // processHNMs(vHNMs);
 
     bool isProton = false;
     bool isDeuteron = false;

EventFiltering/Zorro.cxx

@@ -113,7 +113,7 @@ void Zorro::populateExternalHists(int runNumber, TH2* ZorroHisto, TH2* ToiHisto)
   }
   if (!ToiHisto) {
     LOGF(info, "TOI histogram not set, creating a new one");
-    ToiHisto = new TH2D("TOI", "TOI", 1, -0.5, 0.5, mTOIs.size(), -0.5, mTOIs.size() - 0.5);
+    ToiHisto = new TH2D("TOI", "TOI", 1, -0.5, 0.5, mTOIs.size() * 2, -0.5, mTOIs.size() * 2 - 0.5);
   }
   // if it is the first run, initialize the histogram
   if (mRunNumberHistos.size() == 0) {
@@ -126,17 +126,22 @@ void Zorro::populateExternalHists(int runNumber, TH2* ZorroHisto, TH2* ToiHisto)
       ZorroHisto->GetYaxis()->SetBinLabel(i + 2 + mTOIs.size(), Form("%s scalers", mTOIs[i].data()));
     }
     // TOI histogram
-    ToiHisto->SetBins(1, -0.5, 0.5, mTOIs.size(), -0.5, mTOIs.size() - 0.5);
+    ToiHisto->SetBins(1, -0.5, 0.5, mTOIs.size() * 2, -0.5, mTOIs.size() * 2 - 0.5);
     ToiHisto->GetXaxis()->SetBinLabel(1, Form("%d", runNumber));
     for (size_t i{0}; i < mTOIs.size(); ++i) {
-      ToiHisto->GetYaxis()->SetBinLabel(i + 1, mTOIs[i].data());
+      ToiHisto->GetYaxis()->SetBinLabel(i * 2 + 1, mTOIs[i].data());
+      ToiHisto->GetYaxis()->SetBinLabel(i * 2 + 2, Form("%s AnalysedTriggers", mTOIs[i].data()));
     }
   }
   if (mInspectedTVX) {
     ZorroHisto->Fill(Form("%d", runNumber), "inspected TVX", mInspectedTVX->GetBinContent(1));
     ZorroHisto->SetBinError(mRunNumberHistos.size() + 1, 1, mInspectedTVX->GetBinError(1));
   }
   if (mSelections) {
+    mAnalysedTriggers = new TH1D("AnalysedTriggers", "", mSelections->GetNbinsX() - 2, -0.5, mSelections->GetNbinsX() - 2.5);
+    for (int iBin{2}; iBin < mSelections->GetNbinsX(); ++iBin) { // Exclude first and last bins as they are total number of analysed and selected events, respectively
+      mAnalysedTriggers->GetXaxis()->SetBinLabel(iBin - 1, mSelections->GetXaxis()->GetBinLabel(iBin));
+    }
     for (size_t i{0}; i < mTOIs.size(); ++i) {
       int bin = findBin(mSelections, mTOIs[i]);
       ZorroHisto->Fill(Form("%d", runNumber), Form("%s selections", mTOIs[i].data()), mSelections->GetBinContent(bin));
@@ -248,6 +253,11 @@ bool Zorro::isSelected(uint64_t bcGlobalId, uint64_t tolerance, TH2* ToiHisto)
     if (mTOIidx[i] < 0) {
       continue;
     } else if (mLastResult.test(mTOIidx[i])) {
+      if (ToiHisto && mAnalysedTriggers) {
+        int binX = ToiHisto->GetXaxis()->FindBin(Form("%d", mRunNumber));
+        int binY = ToiHisto->GetYaxis()->FindBin(Form("%s AnalysedTriggers", mTOIs[i].data()));
+        ToiHisto->SetBinContent(binX, binY, mAnalysedTriggers->GetBinContent(mAnalysedTriggers->GetXaxis()->FindBin(mTOIs[i].data())));
+      }
       mTOIcounts[i] += (lastSelectedIdx != mLastSelectedIdx); /// Avoid double counting
       if (mAnalysedTriggersOfInterest && lastSelectedIdx != mLastSelectedIdx) {
         mAnalysedTriggersOfInterest->Fill(i);

PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

@@ -130,6 +130,7 @@ struct FemtoUniverseProducerTask {
   Configurable<bool> confIsForceGRP{"confIsForceGRP", false, "Set true if the magnetic field configuration is not available in the usual CCDB directory (e.g. for Run 2 converted data or unanchorad Monte Carlo)"};
 
   Configurable<bool> confDoSpher{"confDoSpher", false, "Calculate sphericity. If false sphericity will take value of 2."};
+  Configurable<bool> confStoreMCmothers{"confStoreMCmothers", false, "MC truth: Fill with not only primary particles and store mothers' PDG in tempFitVar."};
   Configurable<bool> confFillCollExt{"confFillCollExt", false, "Option to fill collision extended table"};
 
   /// Event cuts
@@ -675,20 +676,23 @@ struct FemtoUniverseProducerTask {
   }
 
   template <typename ParticleType>
-  void fillDebugParticleMC(ParticleType const& particle)
+  int32_t getMotherPDG(ParticleType particle)
   {
     auto motherparticlesMC = particle.template mothers_as<aod::McParticles>();
     if (!motherparticlesMC.empty()) {
       auto motherparticleMC = motherparticlesMC.front();
-      if (particle.isPhysicalPrimary())
-        outputDebugPartsMC(0);
-      else
-        outputDebugPartsMC(motherparticleMC.pdgCode());
+      return particle.isPhysicalPrimary() ? 0 : motherparticleMC.pdgCode();
     } else {
-      outputDebugPartsMC(9999);
+      return 9999;
     }
   }
 
+  template <typename ParticleType>
+  void fillDebugParticleMC(ParticleType const& particle)
+  {
+    outputDebugPartsMC(getMotherPDG(particle));
+  }
+
   template <typename ParticleType>
   void fillMCParticle(ParticleType const& particle, o2::aod::femtouniverseparticle::ParticleType fdparttype)
   {
@@ -1691,7 +1695,7 @@ struct FemtoUniverseProducerTask {
             if (pdgCode == 333) { // && (recoMcIds && recoMcIds->get().contains(particle.globalIndex()))) { // ATTENTION: all Phi mesons are NOT primary particles
               pass = true;
             } else {
-              if (particle.isPhysicalPrimary() || (confActivateSecondaries && recoMcIds && recoMcIds->get().contains(particle.globalIndex())))
+              if (confStoreMCmothers || particle.isPhysicalPrimary() || (confActivateSecondaries && recoMcIds && recoMcIds->get().contains(particle.globalIndex())))
                 pass = true;
             }
           }
@@ -1711,6 +1715,8 @@ struct FemtoUniverseProducerTask {
       // auto cutContainer = trackCuts.getCutContainer<aod::femtouniverseparticle::CutContainerType>(track);
       // instead of the bitmask, the PDG of the particle is stored as uint32_t
 
+      int32_t variablePDG = confStoreMCmothers ? getMotherPDG(particle) : particle.pdgCode();
+
       // now the table is filled
       if constexpr (resolveDaughs) {
         tmpIDtrack.push_back(particle.globalIndex());
@@ -1724,7 +1730,7 @@ struct FemtoUniverseProducerTask {
                     aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
                     0,
                     pdgCode,
-                    pdgCode,
+                    variablePDG,
                     childIDs,
                     0,
                     0);
@@ -1770,6 +1776,9 @@ struct FemtoUniverseProducerTask {
             }
           }
         }
+
+        int32_t variablePDG = confStoreMCmothers ? getMotherPDG(particle) : particle.pdgCode();
+
         if (!confIsActivateCascade) {
           outputParts(outputCollision.lastIndex(),
                       particle.pt(),
@@ -1778,7 +1787,7 @@ struct FemtoUniverseProducerTask {
                       aod::femtouniverseparticle::ParticleType::kMCTruthTrack,
                       0,
                       static_cast<uint32_t>(particle.pdgCode()),
-                      particle.pdgCode(),
+                      variablePDG,
                       childIDs,
                       0,
                       0);

PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx

@@ -248,6 +248,8 @@ struct FemtoUniversePairTaskTrackV0Extended {
     registryMCtruth.add("minus/MCtruthPiPt", "MC truth pions;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
     registryMCtruth.add("minus/MCtruthPrPt", "MC truth protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}});
 
+    registryMCtruth.add("motherParticle", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}});
+
     // MC reco
     registryMCreco.add("plus/MCrecoLambda", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
     registryMCreco.add("plus/MCrecoLambdaChildPr", "MC reco Lambdas;#it{p}_{T} (GeV/c); #eta", {HistType::kTH2F, {{500, 0, 5}, {400, -1.0, 1.0}}});
@@ -1013,6 +1015,39 @@ struct FemtoUniversePairTaskTrackV0Extended {
   }
   PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Extended, processPairFractionsV0, "Process MC data to obtain pair fractions for V0V0 pairs", false);
 
+  void processPairFractionsMCTruthV0(FilteredFDCollisions const& cols, FemtoFullParticles const& /*parts*/)
+  {
+    ColumnBinningPolicy<aod::collision::PosZ, aod::femtouniversecollision::MultNtr> colBinningMult{{confVtxBins, confMultBins}, true};
+    ColumnBinningPolicy<aod::collision::PosZ, aod::femtouniversecollision::MultV0M> colBinningCent{{confVtxBins, confMultBins}, true};
+
+    auto mixedCollProcessFunc = [&](auto& collision1, auto& collision2) -> void {
+      auto groupPartsOne = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache);
+      auto groupPartsTwo = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache);
+
+      for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) {
+        int pdgCode1 = static_cast<int>(p1.pidCut());
+        if ((confV0Type1 == 0 && pdgCode1 != kLambda0) || (confV0Type1 == 1 && pdgCode1 != kLambda0Bar))
+          continue;
+        int pdgCode2 = static_cast<int>(p2.pidCut());
+        if ((confV0Type2 == 0 && pdgCode2 != kLambda0) || (confV0Type2 == 1 && pdgCode2 != kLambda0Bar))
+          continue;
+
+        registryMCtruth.fill(HIST("motherParticle"), p1.tempFitVar(), p2.tempFitVar());
+      }
+    };
+
+    if (confUseCent) {
+      for (const auto& [collision1, collision2] : soa::selfCombinations(colBinningCent, confNEventsMix, -1, cols, cols)) {
+        mixedCollProcessFunc(collision1, collision2);
+      }
+    } else {
+      for (const auto& [collision1, collision2] : soa::selfCombinations(colBinningMult, confNEventsMix, -1, cols, cols)) {
+        mixedCollProcessFunc(collision1, collision2);
+      }
+    }
+  }
+  PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Extended, processPairFractionsMCTruthV0, "Process MC data to obtain pair fractions for V0V0 MC truth pairs", false);
+
   void processMCReco(FemtoRecoParticles const& parts, aod::FdMCParticles const& mcparts)
   {
     for (const auto& part : parts) {