Merge branch 'master' into extrawork

Common/TableProducer/CMakeLists.txt

@@ -76,6 +76,11 @@ o2physics_add_dpl_workflow(tracks-extra-converter
                     PUBLIC_LINK_LIBRARIES
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(mft-tracks-converter
+                    SOURCES mftTracksConverter.cxx
+                    PUBLIC_LINK_LIBRARIES
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(fdd-converter
                     SOURCES fddConverter.cxx
                     PUBLIC_LINK_LIBRARIES

Common/TableProducer/mftTracksConverter.cxx

@@ -0,0 +1,69 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file mftTracksConverter.cxx
+/// \brief Converts MFTTracks table from version 000 to 001
+
+/// This task allows for the conversion of the MFTTracks table from version 000 to 001.
+/// The conversion is needed because the table has been extended with the MFTClusterSizesAndTracksFlags column
+/// while nClusters and isCA columns are evaluated dynamically from it.
+/// In the converter a dummy MFTClusterSizesAndTracksFlags is created and filled according to the number of clusters
+
+/// \author L.Micheletti <[email protected]>
+
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+struct MftTracksConverter {
+  Produces<aod::StoredMFTTracks_001> mftTracks_001;
+  void process(aod::MFTTracks_000 const& mftTracks_000)
+  {
+
+    for (const auto& track0 : mftTracks_000) {
+      uint64_t mftClusterSizesAndTrackFlags = 0;
+      int8_t nClusters = track0.nClusters();
+
+      for (int layer = 0; layer < 10; ++layer) {
+        mftClusterSizesAndTrackFlags &= ~(0x3fULL << (layer * 6));
+        mftClusterSizesAndTrackFlags |= (layer < nClusters) ? (1ULL << (layer * 6)) : 0;
+      }
+
+      mftTracks_001(track0.collisionId(),
+                    track0.x(),
+                    track0.y(),
+                    track0.z(),
+                    track0.phi(),
+                    track0.tgl(),
+                    track0.signed1Pt(),
+                    mftClusterSizesAndTrackFlags,
+                    track0.chi2(),
+                    track0.trackTime(),
+                    track0.trackTimeRes());
+    }
+  }
+};
+
+/// Spawn the extended table for MFTTracks001 to avoid the call to the internal spawner and a consequent circular dependency
+struct MFTTracksSpawner {
+  Spawns<aod::MFTTracks_001> mftTracks_001;
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<MftTracksConverter>(cfgc),
+    adaptAnalysisTask<MFTTracksSpawner>(cfgc),
+  };
+}

Common/TableProducer/mftmchMatchingML.cxx

@@ -203,13 +203,24 @@ struct mftmchMatchingML {
 
   void process(aod::Collisions const& collisions, soa::Filtered<aod::FwdTracks> const& fwdtracks, aod::MFTTracks const& mfttracks)
   {
-    for (auto& [fwdtrack, mfttrack] : combinations(CombinationsFullIndexPolicy(fwdtracks, mfttracks))) {
-
+    for (auto& fwdtrack : fwdtracks) {
       if (fwdtrack.trackType() == aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) {
-        if (fwdtrack.has_collision() && mfttrack.has_collision()) {
-          if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) {
-            double result = matchONNX(fwdtrack, mfttrack);
-            if (result > cfgThrScore) {
+        double bestscore = 0;
+        int bestmfttrackid = -1;
+        for (auto& mfttrack : mfttracks) {
+          if (fwdtrack.has_collision() && mfttrack.has_collision()) {
+            if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) {
+              double result = matchONNX(fwdtrack, mfttrack);
+              if (result > cfgThrScore) {
+                bestscore = result;
+                bestmfttrackid = mfttrack.globalIndex();
+              }
+            }
+          }
+        }
+        if (bestmfttrackid != -1) {
+          for (auto& mfttrack : mfttracks) {
+            if (mfttrack.globalIndex() == bestmfttrackid) {
               double mftchi2 = mfttrack.chi2();
               SMatrix5 mftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt());
               std::vector<double> mftv1;
@@ -222,7 +233,7 @@ struct mftmchMatchingML {
               double px = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi());
               double py = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi());
               double pz = fwdtrack.p() * cos(M_PI / 2 - atan(mfttrack.tgl()));
-              fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), fwdtrack.pDca(), fwdtrack.rAtAbsorberEnd(), 0, 0, 0, result, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY);
+              fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), fwdtrack.pDca(), fwdtrack.rAtAbsorberEnd(), 0, 0, 0, bestscore, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY);
             }
           }
         }

DPG/Tasks/AOTTrack/PID/qaPIDTPCSignal.cxx

@@ -35,9 +35,18 @@ struct tpcPidQaSignal {
   Configurable<int> nBinsP{"nBinsP", 3000, "Number of bins for the momentum"};
   Configurable<bool> runPerRunOutput{"runPerRunOutput", false, "Flag to produce run by run output for e.g. lite calibration purposes"};
   Configurable<bool> enabledEdxPerID{"enabledEdxPerID", false, "Flag to produce dE/dx per particle ID histograms"};
+  Configurable<int> trdSelection{"trdSelection", 0, "Flag for the TRD selection: -1 no TRD, 0 no selection, 1 TRD"};
   Configurable<float> fractionOfEvents{"fractionOfEvents", 0.1f, "Downsampling factor for the events for derived data"};
   Configurable<float> minP{"minP", 0.01, "Minimum momentum in range"};
+  Configurable<float> maxEta{"maxEta", 0.8, "Maximum eta in range"};
+  Configurable<float> maxITSChi2{"maxITSChi2", 36, "Maximum chi2 in range"};
+  Configurable<float> maxTPCChi2{"maxTPCChi2", 4, "Maximum chi2 in range"};
   Configurable<float> maxP{"maxP", 20, "Maximum momentum in range"};
+  Configurable<int> minITSCls{"minITSCls", 6, "Minimum number of ITS clusters"};
+  Configurable<int> pidInTracking{"pidInTracking", -1, "PID in tracking"};
+  Configurable<int> minTPCCls{"minTPCCls", 0, "Minimum number of TPC clusters"};
+  Configurable<int> minTPCClsFindable{"minTPCClsFindable", 0, "Minimum number of TPC findable clusters"};
+  Configurable<float> minCrossedRowsOverFindableCls{"minCrossedRowsOverFindableCls", 0.8f, "Minimum number of TPC found/findable clusters"};
   ConfigurableAxis dEdxBins{"dEdxBins", {5000, 0.f, 5000.f}, "Binning in dE/dx"};
   Configurable<float> minTPCNcls{"minTPCNcls", 0.f, "Minimum number or TPC Clusters for tracks"};
   Configurable<float> minNClsCrossedRows{"minNClsCrossedRows", 70.f, "Minimum number or TPC crossed rows for tracks"};
@@ -63,6 +72,25 @@ struct tpcPidQaSignal {
     h->GetXaxis()->SetBinLabel(2, "Passed ev. sel.");
     h->GetXaxis()->SetBinLabel(3, "Passed vtx Z");
 
+    // Event properties
+    h = histos.add<TH1>("trksel", "", kTH1D, {{10, 0.5, 10.5, "Trk. Sel."}});
+    h->GetXaxis()->SetBinLabel(1, "Tracks read");
+    h->GetXaxis()->SetBinLabel(2, Form("Has > %i ITS clusters", minITSCls.value));
+    h->GetXaxis()->SetBinLabel(3, Form("Has > %i TPC clusters findable", minTPCClsFindable.value));
+    h->GetXaxis()->SetBinLabel(4, Form("Has > %f TPC clusters found", minTPCNcls.value));
+    h->GetXaxis()->SetBinLabel(5, Form("Has > %f Found/Findable Ratio", minCrossedRowsOverFindableCls.value));
+    h->GetXaxis()->SetBinLabel(6, Form("Has > %f Xrows", minNClsCrossedRows.value));
+    h->GetXaxis()->SetBinLabel(7, "All PID in trk");
+    if (pidInTracking == -1) {
+      h->GetXaxis()->SetBinLabel(7, Form("PID in trk %i", pidInTracking.value));
+    }
+    h->GetXaxis()->SetBinLabel(8, "no TRD sel.");
+    if (trdSelection == -1) {
+      h->GetXaxis()->SetBinLabel(8, "has no TRD");
+    } else if (trdSelection == 1) {
+      h->GetXaxis()->SetBinLabel(8, "has TRD");
+    }
+
     histos.add("event/vertexz", "", kTH1D, {vtxZAxis});
     hdedx = histos.add<TH3>("event/tpcsignal", "", kTH3D, {pAxis, dedxAxis, chargeAxis});
     if (enabledEdxPerID) {
@@ -74,9 +102,57 @@ struct tpcPidQaSignal {
     histos.print();
   }
 
-  Filter trackFilterEta = (nabs(aod::track::eta) < 0.8f);
-  Filter trackFilterITS = (aod::track::itsChi2NCl < 36.f);
-  Filter trackFilterTPC = ((aod::track::tpcChi2NCl < 4.f));
+  template <typename T>
+  bool isTrackSelected(const T& track)
+  {
+    histos.fill(HIST("trksel"), 1);
+    if (track.itsNCls() < minITSCls) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 2);
+    if (track.tpcNClsFindable() <= minTPCClsFindable) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 3);
+    if (track.tpcNClsFound() < minTPCNcls) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 4);
+    if (track.tpcCrossedRowsOverFindableCls() < minCrossedRowsOverFindableCls) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 5);
+    if (track.tpcNClsCrossedRows() < minNClsCrossedRows) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 6);
+    if (pidInTracking != -1 && (track.pidForTracking() != std::abs(pidInTracking))) {
+      return false;
+    }
+    histos.fill(HIST("trksel"), 7);
+    switch (trdSelection) {
+      case 0:
+        break;
+      case -1:
+        if (track.hasTRD()) {
+          return false;
+        }
+        break;
+      case 1:
+        if (!track.hasTRD()) {
+          return false;
+        }
+        break;
+      default:
+        LOG(fatal) << "Invalid TRD selection";
+    }
+    histos.fill(HIST("trksel"), 8);
+    return true;
+  }
+
+  Filter trackFilterEta = (nabs(aod::track::eta) < maxEta);
+  Filter trackFilterITS = (aod::track::itsChi2NCl < maxITSChi2);
+  Filter trackFilterTPC = ((aod::track::tpcChi2NCl < maxTPCChi2));
   using TrackCandidates = soa::Join<aod::TracksIU, aod::TracksExtra>;
   void processEvSel(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision,
                     soa::Filtered<TrackCandidates> const& tracks,
@@ -100,22 +176,25 @@ struct tpcPidQaSignal {
     histos.fill(HIST("event/vertexz"), collision.posZ());
 
     for (const auto& t : tracks) {
-      if (t.itsNCls() < 6) {
-        continue;
-      }
-      if (t.tpcNClsFindable() <= 0) {
+      if (!isTrackSelected(t)) {
         continue;
       }
-      if (t.tpcNClsFound() < minTPCNcls) {
-        continue;
+      switch (trdSelection) {
+        case 0:
+          break;
+        case -1:
+          if (t.hasTRD()) {
+            continue;
+          }
+          break;
+        case 1:
+          if (!t.hasTRD()) {
+            continue;
+          }
+          break;
+        default:
+          LOG(fatal) << "Invalid TRD selection";
       }
-      if (t.tpcCrossedRowsOverFindableCls() < 0.8f) {
-        continue;
-      }
-      if (t.tpcNClsCrossedRows() < minNClsCrossedRows) {
-        continue;
-      }
-
       if (runPerRunOutput.value == true && lastRun != collision.bc().runNumber()) {
         lastRun = collision.bc().runNumber();
         AxisSpec pAxis{nBinsP, minP, maxP, "#it{p}/|Z| (GeV/#it{c})"};
@@ -134,38 +213,25 @@ struct tpcPidQaSignal {
   }
   PROCESS_SWITCH(tpcPidQaSignal, processEvSel, "Process with event selection", false);
 
-  void processNoEvSel(aod::Collision const& collision,
-                      soa::Filtered<TrackCandidates> const& tracks)
+  void processNoEvSel(soa::Filtered<TrackCandidates> const& tracks, aod::Collisions const& collisions)
   {
     if (fractionOfEvents < 1.f && (static_cast<float>(rand_r(&randomSeed)) / static_cast<float>(RAND_MAX)) > fractionOfEvents) { // Skip events that are not sampled
       return;
     }
 
-    histos.fill(HIST("event/evsel"), 1);
-    histos.fill(HIST("event/evsel"), 2);
-
-    if (abs(collision.posZ()) > 10.f) {
-      return;
-    }
-    histos.fill(HIST("event/evsel"), 3);
-    histos.fill(HIST("event/vertexz"), collision.posZ());
+    histos.fill(HIST("event/evsel"), 1, collisions.size());
 
     for (const auto& t : tracks) {
-      if (t.itsNCls() < 6) {
+      if (!t.has_collision()) {
         continue;
       }
-      if (t.tpcNClsFindable() <= 0) {
+      if (abs(t.collision().posZ()) > 10.f) {
         continue;
       }
-      if (t.tpcNClsFound() < minTPCNcls) {
-        continue;
-      }
-      if (t.tpcCrossedRowsOverFindableCls() < 0.8f) {
-        continue;
-      }
-      if (t.tpcNClsCrossedRows() < minNClsCrossedRows) {
+      if (!isTrackSelected(t)) {
         continue;
       }
+
       hdedx->Fill(t.tpcInnerParam(), t.tpcSignal(), t.sign());
       if (enabledEdxPerID) {
         hdedx_perID[t.pidForTracking()]->Fill(t.tpcInnerParam(), t.tpcSignal(), t.sign());

DPG/Tasks/AOTTrack/perfK0sResolution.cxx

@@ -60,33 +60,63 @@ struct perfK0sResolution {
   Configurable<float> v0lifetime{"v0lifetime", 3., "n ctau"};
   Configurable<float> rapidity{"rapidity", 0.5, "rapidity"};
   Configurable<float> nSigTPC{"nSigTPC", 10., "nSigTPC"};
-  Configurable<int> requireTRDneg{"requireTRDneg", 0, "requireTRDneg"}; // 0: no requirement, >0: TRD only tracks, <0: reject TRD tracks
-  Configurable<int> requireTRDpos{"requireTRDpos", 0, "requireTRDpos"}; // 0: no requirement, >0: TRD only tracks, <0: reject TRD tracks
+  Configurable<int> trdSelectionPos{"trdSelectionPos", 0, "Flag for the TRD selection on positive daughters: -1 no TRD, 0 no selection, 1 TRD"};
+  Configurable<int> trdSelectionNeg{"trdSelectionNeg", 0, "Flag for the TRD selection on negative daughters: -1 no TRD, 0 no selection, 1 TRD"};
   Configurable<bool> eventSelection{"eventSelection", true, "event selection"};
 
   template <typename T1, typename T2, typename C>
   bool acceptV0(const T1& v0, const T2& ntrack, const T2& ptrack, const C& collision)
   {
     // Apply selections on V0
     if (TMath::Abs(v0.yK0Short()) > rapidity)
-      return kFALSE;
+      return false;
     if (v0.v0cosPA() < v0setting_cospa)
-      return kFALSE;
+      return false;
     if (v0.v0radius() < v0setting_radius)
-      return kFALSE;
+      return false;
     if (v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * pid_constants::sMasses[PID::K0] > 2.684 * v0lifetime)
-      return kFALSE;
+      return false;
 
     // Apply selections on V0 daughters
     if (!ntrack.hasTPC() || !ptrack.hasTPC())
-      return kFALSE;
+      return false;
     if (ntrack.tpcNSigmaPi() > nSigTPC || ptrack.tpcNSigmaPi() > nSigTPC)
-      return kFALSE;
-    if ((requireTRDneg > 0 && !ntrack.hasTRD()) || (requireTRDneg < 0 && ntrack.hasTRD()))
-      return kFALSE;
-    if ((requireTRDpos > 0 && !ptrack.hasTRD()) || (requireTRDpos < 0 && ptrack.hasTRD()))
-      return kFALSE;
-    return kTRUE;
+      return false;
+    switch (trdSelectionPos) {
+      case -1:
+        if (ptrack.hasTRD()) {
+          return false;
+        }
+        break;
+      case 0:
+        break;
+      case 1:
+        if (!ptrack.hasTRD()) {
+          return false;
+        }
+        break;
+      default:
+        LOG(fatal) << "Invalid TRD selection for positive daughter";
+        break;
+    }
+    switch (trdSelectionNeg) {
+      case -1:
+        if (ntrack.hasTRD()) {
+          return false;
+        }
+        break;
+      case 0:
+        break;
+      case 1:
+        if (!ntrack.hasTRD()) {
+          return false;
+        }
+        break;
+      default:
+        LOG(fatal) << "Invalid TRD selection for negative daughter";
+        break;
+    }
+    return true;
   }
 
   Filter v0Filter = nabs(aod::v0data::dcapostopv) > v0setting_dcapostopv&& nabs(aod::v0data::dcanegtopv) > v0setting_dcanegtopv&& aod::v0data::dcaV0daughters < v0setting_dcav0dau;