Merge branch 'AliceO2Group:master' into master

Author: YazhenLin

PWGEM/Dilepton/DataModel/lmeeMLTables.h

@@ -159,31 +159,38 @@ DECLARE_SOA_COLUMN(Chi2MFTMCHMID, chi2MFTMCHMID, float);  //! chi2/ndf of MFT-MC
 DECLARE_SOA_COLUMN(NClustersMFT, nClustersMFT, uint8_t);  //!
 DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool);           //!
 DECLARE_SOA_COLUMN(IsCorrectMatch, isCorrectMatch, bool); //!
+
+DECLARE_SOA_COLUMN(NMFTs, nMFTs, uint16_t); //! number of MFTsa tracks per collision
 } // namespace emmlfwdtrack
 
 DECLARE_SOA_TABLE(EMFwdTracksForML, "AOD", "EMFWDTRKML", //!
-                  o2::soa::Index<>, collision::PosZ, collision::NumContrib, mult::MultFT0C, evsel::NumTracksInTimeRange, evsel::SumAmpFT0CInTimeRange, emmltrack::HadronicRate,
-                  fwdtrack::TrackType,
+                  o2::soa::Index<>, collision::PosZ, /*collision::NumContrib,*/ mult::MultFT0C, /*evsel::NumTracksInTimeRange,*/ evsel::SumAmpFT0CInTimeRange, emmltrack::HadronicRate, emmlfwdtrack::NMFTs,
+                  // fwdtrack::TrackType,
 
                   emmlfwdtrack::Signed1PtMFTatMP, emmlfwdtrack::TglMFTatMP, emmlfwdtrack::PhiMFTatMP,
-                  emmlfwdtrack::Signed1PtErrMFTatMP, emmlfwdtrack::TglErrMFTatMP, emmlfwdtrack::PhiErrMFTatMP,
                   emmlfwdtrack::XMFTatMP, emmlfwdtrack::YMFTatMP,
-                  emmlfwdtrack::XErrMFTatMP, emmlfwdtrack::YErrMFTatMP,
 
                   emmlfwdtrack::Signed1PtMCHMIDatMP, emmlfwdtrack::TglMCHMIDatMP, emmlfwdtrack::PhiMCHMIDatMP,
-                  emmlfwdtrack::Signed1PtErrMCHMIDatMP, emmlfwdtrack::TglErrMCHMIDatMP, emmlfwdtrack::PhiErrMCHMIDatMP,
                   emmlfwdtrack::XMCHMIDatMP, emmlfwdtrack::YMCHMIDatMP,
-                  emmlfwdtrack::XErrMCHMIDatMP, emmlfwdtrack::YErrMCHMIDatMP,
 
-                  // fwdtrack::FwdDcaX, fwdtrack::FwdDcaY,
-                  fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd,
-                  fwdtrack::Chi2MatchMCHMID, fwdtrack::Chi2MatchMCHMFT,
-                  fwdtrack::MFTClusterSizesAndTrackFlags, emmlfwdtrack::Chi2MFTMCHMID, emmlfwdtrack::Chi2MCHMID, emmlfwdtrack::Chi2MFT, emmlfwdtrack::NClustersMFT,
+                  // fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd,
+                  // fwdtrack::Chi2MatchMCHMID,
+                  fwdtrack::Chi2MatchMCHMFT,
+                  // fwdtrack::MFTClusterSizesAndTrackFlags, emmlfwdtrack::Chi2MFTMCHMID, emmlfwdtrack::Chi2MCHMID, emmlfwdtrack::Chi2MFT, emmlfwdtrack::NClustersMFT,
                   mcparticle::PdgCode, emmlfwdtrack::IsPrimary, emmlfwdtrack::IsCorrectMatch, mcparticle::Pt, mcparticle::Eta, mcparticle::Phi);
 
 // iterators
 using EMFwdTrackForML = EMFwdTracksForML::iterator;
 
+DECLARE_SOA_TABLE(EMFwdTrackErrsForML, "AOD", "EMFWDTRKERRML", //! Joinable with EMFwdTracksForML
+                  emmlfwdtrack::Signed1PtErrMFTatMP, emmlfwdtrack::TglErrMFTatMP, emmlfwdtrack::PhiErrMFTatMP,
+                  emmlfwdtrack::XErrMFTatMP, emmlfwdtrack::YErrMFTatMP,
+                  emmlfwdtrack::Signed1PtErrMCHMIDatMP, emmlfwdtrack::TglErrMCHMIDatMP, emmlfwdtrack::PhiErrMCHMIDatMP,
+                  emmlfwdtrack::XErrMCHMIDatMP, emmlfwdtrack::YErrMCHMIDatMP);
+
+// iterators
+using EMFwdTrackErrForML = EMFwdTrackErrsForML::iterator;
+
 } // namespace o2::aod
 
 #endif // PWGEM_DILEPTON_DATAMODEL_LMEEMLTABLES_H_

PWGEM/Dilepton/TableProducer/treeCreatorMuonML.cxx

@@ -76,7 +76,8 @@ struct TreeCreatorMuonML {
   using MyMFTTracksMC = soa::Join<o2::aod::MFTTracks, aod::McMFTTrackLabels>;
   using MyMFTTrackMC = MyMFTTracksMC::iterator;
 
-  Produces<aod::EMFwdTracksForML> mltable;
+  Produces<aod::EMFwdTracksForML> trackTable;
+  Produces<aod::EMFwdTrackErrsForML> trackErrTable;
 
   // Configurables
   Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
@@ -88,6 +89,8 @@ struct TreeCreatorMuonML {
   Configurable<std::string> cfgZShiftPath{"cfgZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"};
   Configurable<float> cfgManualZShift{"cfgManualZShift", 0, "manual z-shift for propagation of global muon to PV"};
   Configurable<float> cfgDownSampling{"cfgDownSampling", 1.1, "down sampling for fake matches"};
+  Configurable<float> matchingZ{"matchingZ", -77.5, "z position where matching is performed"};
+  Configurable<float> maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 100.f, "max. chi2 for MCH-MFT matching"};
 
   struct : ConfigurableGroup {
     std::string prefix = "eventCutGroup";
@@ -107,14 +110,19 @@ struct TreeCreatorMuonML {
 
   struct : ConfigurableGroup {
     std::string prefix = "glMuonCutGroup";
-    Configurable<float> minEta{"minEta", -3.6, "min. eta acceptance for MFT-MCH-MID"};
-    Configurable<float> maxEta{"maxEta", -2.5, "max. eta acceptance for MFT-MCH-MID"};
-    Configurable<float> maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 100.f, "max. chi2 for MCH-MFT matching"};
-    Configurable<float> maxChi2{"maxChi2", 20.f, "max. chi2 for global muon"};
+    // Configurable<float> minEta{"minEta", -3.6, "min. eta acceptance for MFT-MCH-MID"};
+    // Configurable<float> maxEta{"maxEta", -2.5, "max. eta acceptance for MFT-MCH-MID"};
+    // Configurable<float> maxChi2{"maxChi2", 1e+10, "max. chi2 for global muon"};
     Configurable<bool> refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"};
-    Configurable<float> matchingZ{"matchingZ", -77.5, "z position where matching is performed"};
   } glMuonCutGroup;
 
+  struct : ConfigurableGroup {
+    std::string prefix = "MFTCutGroup";
+    Configurable<float> minPt{"minPt", 0.04f, "min. pT for MFTsa to reject crazy tracks"};
+    Configurable<float> minEta{"minEta", -4.1f, "min. eta acceptance for MFTsa to reject crazy tracks"};
+    Configurable<float> maxEta{"maxEta", -2.0f, "max. eta acceptance for MFTsa to reject crazy tracks"};
+  } MFTCutGroup;
+
   o2::ccdb::CcdbApi ccdbApi;
   Service<o2::ccdb::BasicCCDBManager> ccdb;
   std::mt19937 engine;
@@ -198,6 +206,8 @@ struct TreeCreatorMuonML {
 
   void addHistograms()
   {
+    fRegistry.add("Event/hCorrFT0CvsMFT", "mult. corr. between FT0C and MFT;multFT0C;N_{MFTsa}", kTH2F, {{600, 0, 60000}, {100, 0, 10000}}, false);
+
     auto hMuonType = fRegistry.add<TH1>("hMuonType", "muon type", kTH1F, {{5, -0.5f, 4.5f}}, false);
     hMuonType->GetXaxis()->SetBinLabel(1, "MFT-MCH-MID (global muon)");
     hMuonType->GetXaxis()->SetBinLabel(2, "MFT-MCH-MID (global muon other match)");
@@ -270,7 +280,7 @@ struct TreeCreatorMuonML {
   }
 
   template <bool withMFTCov, typename TCollision, typename TFwdTrack, typename TFwdTracks, typename TMFTTracks, typename TMFTTracksCov>
-  bool fillFwdTrackTable(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const&, TMFTTracks const&, TMFTTracksCov const& mftCovs, const float hadronicRate)
+  bool fillFwdTrackTable(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const&, TMFTTracks const&, TMFTTracksCov const& mftCovs, const float hadronicRate, const uint16_t nmft)
   {
     if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) {
       return false;
@@ -284,6 +294,10 @@ struct TreeCreatorMuonML {
     float chi2 = fwdtrack.chi2() / (2.f * (mchtrack.nClusters() + mfttrack.nClusters()) - 5.f);
     float chi2mft = mfttrack.chi2() / (2.f * mfttrack.nClusters() - 5.f);
 
+    if (mfttrack.eta() < MFTCutGroup.minEta || MFTCutGroup.maxEta < mfttrack.eta() || mfttrack.pt() < MFTCutGroup.minPt) {
+      return false;
+    }
+
     // auto mcParticle_MFTMCHMID = fwdtrack.template mcParticle_as<aod::McParticles>(); // this is identical to mcParticle_MCHMID
     auto mcParticle_MCHMID = mchtrack.template mcParticle_as<aod::McParticles>(); // this is identical to mcParticle_MFTMCHMID
     auto mcParticle_MFT = mfttrack.template mcParticle_as<aod::McParticles>();
@@ -299,33 +313,33 @@ struct TreeCreatorMuonML {
     if (fwdtrack.chi2MatchMCHMID() < 0.f) { // this should never happen. only for protection.
       return false;
     }
-    if (fwdtrack.chi2MatchMCHMFT() < 0.f || glMuonCutGroup.maxMatchingChi2MCHMFT < fwdtrack.chi2MatchMCHMFT()) { // this should never happen. only for protection.
-      return false;
-    }
-    if (fwdtrack.chi2() < 0.f || glMuonCutGroup.maxChi2 < chi2) { // this should never happen. only for protection.
+    if (fwdtrack.chi2MatchMCHMFT() < 0.f || maxMatchingChi2MCHMFT < fwdtrack.chi2MatchMCHMFT()) {
       return false;
     }
+    // if (fwdtrack.chi2() < 0.f || glMuonCutGroup.maxChi2 < chi2) {
+    //   return false;
+    // }
     if (mfttrack.chi2() < 0.f) { // this should never happen. only for protection.
       return false;
     }
 
-    o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mZShift);
+    o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift);
     float pt = propmuonAtPV.getPt();
     float eta = propmuonAtPV.getEta();
     float phi = propmuonAtPV.getPhi();
     // o2::math_utils::bringTo02Pi(phi);
     phi = RecoDecay::constrainAngle(phi, 0, 1U);
 
-    if (eta < glMuonCutGroup.minEta || glMuonCutGroup.maxEta < eta) {
-      return false;
-    }
+    // if (eta < glMuonCutGroup.minEta || glMuonCutGroup.maxEta < eta) {
+    //   return false;
+    // }
 
     float dcaX = propmuonAtPV.getX() - collision.posX();
     float dcaY = propmuonAtPV.getY() - collision.posY();
     // float dcaZ = propmuonAtPV.getZ() - collision.posZ();
     float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
 
-    o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mZShift);
+    o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift);
     // float ptMatchedMCHMID = propmuonAtPV_Matched.getPt();
     float etaMatchedMCHMID = propmuonAtPV_Matched.getEta();
     float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi();
@@ -335,7 +349,7 @@ struct TreeCreatorMuonML {
       pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta)));
     }
 
-    o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, glMuonCutGroup.matchingZ, mBz, mZShift);
+    o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift);
     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);
@@ -356,7 +370,7 @@ struct TreeCreatorMuonML {
     if constexpr (withMFTCov) {
       auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]);
       o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update
-      mftsaAtMP.propagateToZhelix(glMuonCutGroup.matchingZ, mBz);                                   // propagated to matching plane
+      mftsaAtMP.propagateToZhelix(matchingZ, mBz);                                                  // propagated to matching plane
       xMFTatMP = mftsaAtMP.getX();
       yMFTatMP = mftsaAtMP.getY();
       xErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2X());
@@ -369,7 +383,7 @@ struct TreeCreatorMuonML {
       tglErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Tanl());
       phiErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Phi());
 
-      auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, glMuonCutGroup.matchingZ, mBz, mZShift); // propagated to matching plane
+      auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, mZShift); // propagated to matching plane
       xMCHMIDatMP = muonAtMP.getX();
       yMCHMIDatMP = muonAtMP.getY();
       xErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2X());
@@ -387,24 +401,24 @@ struct TreeCreatorMuonML {
     float dphi = phiMatchedMCHMID - phi;
     o2::math_utils::bringToPMPi(dphi);
 
-    mltable(collision.posZ(), collision.numContrib(), collision.multFT0C(), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), hadronicRate,
-            fwdtrack.trackType(),
+    trackTable(collision.posZ(), /*collision.numContrib(),*/ collision.multFT0C(), /*collision.trackOccupancyInTimeRange(),*/ collision.ft0cOccupancyInTimeRange(), hadronicRate, nmft,
+               // fwdtrack.trackType(),
 
-            signed1PtMFTatMP, tglMFTatMP, phiMFTatMP,
-            signed1PtErrMFTatMP, tglErrMFTatMP, phiErrMFTatMP,
-            xMFTatMP, yMFTatMP,
-            xErrMFTatMP, yErrMFTatMP,
+               signed1PtMFTatMP, tglMFTatMP, phiMFTatMP,
+               xMFTatMP, yMFTatMP,
 
-            signed1PtMCHMIDatMP, tglMCHMIDatMP, phiMCHMIDatMP,
-            signed1PtErrMCHMIDatMP, tglErrMCHMIDatMP, phiErrMCHMIDatMP,
-            xMCHMIDatMP, yMCHMIDatMP,
-            xErrMCHMIDatMP, yErrMCHMIDatMP,
+               signed1PtMCHMIDatMP, tglMCHMIDatMP, phiMCHMIDatMP,
+               xMCHMIDatMP, yMCHMIDatMP,
 
-            // dcaX, dcaY,
-            fwdtrack.nClusters(), pDCA, rAtAbsorberEnd, fwdtrack.chi2MatchMCHMID(), fwdtrack.chi2MatchMCHMFT(),
-            mfttrack.mftClusterSizesAndTrackFlags(), chi2, mchtrack.chi2(), chi2mft, mfttrack.nClusters(),
-            pdgCode, isPrimary, isMatched,
-            mcParticle_MCHMID.pt(), mcParticle_MCHMID.eta(), mcParticle_MCHMID.phi());
+               // fwdtrack.chi2MatchMCHMID(),
+               fwdtrack.chi2MatchMCHMFT(),
+               pdgCode, isPrimary, isMatched,
+               mcParticle_MCHMID.pt(), mcParticle_MCHMID.eta(), mcParticle_MCHMID.phi());
+
+    trackErrTable(signed1PtErrMFTatMP, tglErrMFTatMP, phiErrMFTatMP,
+                  xErrMFTatMP, yErrMFTatMP,
+                  signed1PtErrMCHMIDatMP, tglErrMCHMIDatMP, phiErrMCHMIDatMP,
+                  xErrMCHMIDatMP, yErrMCHMIDatMP);
 
     fRegistry.fill(HIST("hMuonType"), fwdtrack.trackType());
     fRegistry.fill(HIST("MFTMCHMID/hPt"), pt);
@@ -432,6 +446,7 @@ struct TreeCreatorMuonML {
 
   SliceCache cache;
   Preslice<aod::FwdTracks> perCollision = o2::aod::fwdtrack::collisionId;
+  Preslice<aod::MFTTracks> perCollision_MFT = o2::aod::fwdtrack::collisionId;
 
   std::unordered_map<int, int> map_mfttrackcovs;
   void processWithMFTCov(MyCollisionsMC const& collisions, aod::BCsWithTimestamps const&, MyFwdTracksMC const& fwdtracks, MyMFTTracksMC const& mfttracks, aod::MFTTracksCov const& mftCovs, aod::McParticles const&, aod::McCollisions const&)
@@ -458,6 +473,10 @@ struct TreeCreatorMuonML {
 
       float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSourceForCptFetcher) * 1.e-3; // kHz
 
+      auto mfttracks_per_collision = mfttracks.sliceBy(perCollision_MFT, collision.globalIndex());
+      uint16_t nmft = mfttracks_per_collision.size();
+      fRegistry.fill(HIST("Event/hCorrFT0CvsMFT"), collision.multFT0C(), mfttracks_per_collision.size());
+
       auto fwdtracks_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex());
       for (const auto& fwdtrack : fwdtracks_coll) {
         if (!fwdtrack.has_mcParticle()) {
@@ -467,7 +486,7 @@ struct TreeCreatorMuonML {
           continue;
         }
 
-        fillFwdTrackTable<true>(collision, fwdtrack, fwdtracks, mfttracks, mftCovs, hadronicRate);
+        fillFwdTrackTable<true>(collision, fwdtrack, fwdtracks, mfttracks, mftCovs, hadronicRate, nmft);
 
       } // end of fwdtrack loop
     } // end of collision loop
@@ -496,6 +515,10 @@ struct TreeCreatorMuonML {
 
       float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSourceForCptFetcher) * 1.e-3; // kHz
 
+      auto mfttracks_per_collision = mfttracks.sliceBy(perCollision_MFT, collision.globalIndex());
+      uint16_t nmft = mfttracks_per_collision.size();
+      fRegistry.fill(HIST("Event/hCorrFT0CvsMFT"), collision.multFT0C(), mfttracks_per_collision.size());
+
       auto fwdtracks_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex());
       for (const auto& fwdtrack : fwdtracks_coll) {
         if (!fwdtrack.has_mcParticle()) {
@@ -505,7 +528,7 @@ struct TreeCreatorMuonML {
           continue;
         }
 
-        fillFwdTrackTable<false>(collision, fwdtrack, fwdtracks, mfttracks, nullptr, hadronicRate);
+        fillFwdTrackTable<false>(collision, fwdtrack, fwdtracks, mfttracks, nullptr, hadronicRate, nmft);
 
       } // end of fwdtrack loop
     } // end of collision loop