AliceO2Group · ktf · Jul 27, 2021 · Jul 26, 2021 · Jul 26, 2021 · Jul 27, 2021
@@ -103,7 +103,7 @@ struct Alice3LutMaker {
      if (mcParticle.pdgCode() != pdg) {
        continue;
      }
-      if (selPrim.value && !MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Requiring is physical primary
+      if (selPrim.value && !MC::isPhysicalPrimary(mcParticle)) { // Requiring is physical primary
        continue;
      }

@@ -133,7 +133,7 @@ struct Alice3LutMaker {
      if (mcParticle.pdgCode() != pdg) {
        continue;
      }
-      if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Requiring is physical primary
+      if (!MC::isPhysicalPrimary(mcParticle)) { // Requiring is physical primary
        continue;
      }


@@ -118,13 +118,13 @@ struct Alice3SingleParticle {
          continue;
        }
        auto mother = mcParticle.mother0_as<aod::McParticles>();
-        if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+        if (MC::isPhysicalPrimary(mcParticle)) {
          histos.get<TH1>(HIST("primaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
        } else {
          histos.get<TH1>(HIST("secondaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
        }
        if (doPrint) {
-          LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mother.pdgCode() << " and is " << (MC::isPhysicalPrimary<aod::McParticles>(mcParticle) ? "" : "not") << " a primary";
+          LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mother.pdgCode() << " and is " << (MC::isPhysicalPrimary(mcParticle) ? "" : "not") << " a primary";
        }
      }
    }

@@ -152,7 +152,7 @@ struct ftofPidQaMC {
      if (pdgCode != 0 && abs(mcParticle.pdgCode()) != pdgCode) {
        continue;
      }
-      if (useOnlyPhysicsPrimary == 1 && !MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Selecting primaries
+      if (useOnlyPhysicsPrimary == 1 && !MC::isPhysicalPrimary(mcParticle)) { // Selecting primaries
        histos.fill(HIST("p/Sec"), track.p());
        continue;
      }

@@ -167,7 +167,7 @@ struct richPidQaMc {
    if (abs(particle.pdgCode()) == PDGs[pidIndex]) {
      histos.fill(HIST(hnsigmaMC[pidIndex]), track.pt(), nsigma);

-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
        histos.fill(HIST(hnsigmaMCprm[pidIndex]), track.pt(), nsigma);
      } else {
        histos.fill(HIST(hnsigmaMCsec[pidIndex]), track.pt(), nsigma);
@@ -250,7 +250,7 @@ struct richPidQaMc {
      }
      histos.fill(HIST(hnsigma[pid_type]), track.pt(), nsigma);
      histos.fill(HIST(hdelta[pid_type]), track.p(), delta);
-      if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(mcParticle)) { // Selecting primaries
        histos.fill(HIST(hnsigmaprm[pid_type]), track.pt(), nsigma);
        histos.fill(HIST("p/Prim"), track.p());
      } else {

@@ -20,7 +20,6 @@ o2_target_root_dictionary(AnalysisCore
                      include/AnalysisCore/TrackSelection.h
                      include/AnalysisCore/TrackSelectionDefaults.h
                      include/AnalysisCore/TriggerAliases.h
-                      include/AnalysisCore/MC.h
              LINKDEF src/AnalysisCoreLinkDef.h)

 o2_add_executable(merger

@@ -13,6 +13,7 @@
 #define MC_H

 #include "Framework/Logger.h"
+#include "Framework/AnalysisDataModel.h"

 #include "TPDGCode.h"

@@ -61,7 +62,7 @@ bool isStable(int pdg)

 // Ported from AliRoot AliStack::IsPhysicalPrimary
 template <typename Particles, typename Particle>
-bool isPhysicalPrimary(Particle& particle)
+bool isPhysicalPrimary(Particle const& particle)
 {
  // Test if a particle is a physical primary according to the following definition:
  // Particles produced in the collision including products of strong and
@@ -169,6 +170,12 @@ bool isPhysicalPrimary(Particle& particle)
  }
 }

+// Short hand for the standard type
+bool isPhysicalPrimary(o2::aod::McParticle const& particle)
+{
+  return isPhysicalPrimary<o2::aod::McParticles>(particle);
+}
+
 }; // namespace MC

 #endif
@@ -135,7 +135,7 @@ struct pidTOFTaskQA {
    if (abs(particle.pdgCode()) == PDGs[pidIndex]) {
      histos.fill(HIST(hnsigmaMC[pidIndex]), track.pt(), nsigma);

-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
        histos.fill(HIST(hnsigmaMCprm[pidIndex]), track.pt(), nsigma);
      } else {
        histos.fill(HIST(hnsigmaMCsec[pidIndex]), track.pt(), nsigma);
@@ -190,7 +190,7 @@ struct pidTOFTaskQA {
      histos.fill(HIST(hnsigma[pid_type]), t.pt(), nsigma);
      histos.fill(HIST("event/tofbeta"), t.p(), t.beta());
      const auto particle = t.mcParticle();
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { // Selecting primaries
+      if (MC::isPhysicalPrimary(particle)) { // Selecting primaries
        histos.fill(HIST(hnsigmaprm[pid_type]), t.pt(), nsigma);
        histos.fill(HIST("event/tofbetaPrm"), t.p(), t.beta());
      } else {

@@ -70,7 +70,7 @@ struct NucleiSpectraEfficiencyGen {
    spectra.add("histGenPt", "generated particles", HistType::kTH1F, {ptAxis});
  }

-  void process(aod::McCollision const& mcCollision, aod::McParticles& mcParticles)
+  void process(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles)
  {
    //
    // loop over generated particles and fill generated particles
@@ -79,7 +79,7 @@ struct NucleiSpectraEfficiencyGen {
      if (mcParticleGen.pdgCode() != -1000020030) {
        continue;
      }
-      if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticleGen)) {
+      if (!MC::isPhysicalPrimary(mcParticleGen)) {
        continue;
      }
      if (abs(mcParticleGen.y()) > 0.5) {

@@ -152,7 +152,7 @@ struct raacharged {

  Configurable<bool> isMC{"isMC", 1, "0 - data, 1 - MC"};

-  void process(soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::Cents>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection> const& tracks, aod::McParticles& mcParticles)
+  void process(soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::Cents>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection> const& tracks, aod::McParticles const& mcParticles)
  {
    if (!collision.alias()[kINT7]) {
      return;
@@ -179,7 +179,7 @@ struct raacharged {
        continue;
      }
      const auto particle = track.mcParticle();
-      if (MC::isPhysicalPrimary<aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
        mcInfoVal = 0.0;
      } else {
        mcInfoVal = 1.0;
@@ -198,7 +198,7 @@ struct raacharged {
        if (abs(mcParticle.eta()) > 0.8) {
          continue;
        }
-        if (!MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+        if (!MC::isPhysicalPrimary(mcParticle)) {
          continue;
        }


@@ -90,7 +90,7 @@ struct TrackCheckTaskEvSel {
  //Filters
  Filter collfilter = nabs(aod::collision::posZ) < cfgCutVZ;
  void process(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>>::iterator const& col,
-               soa::Join<aod::Tracks, aod::TracksExtra, aod::McTrackLabels>& tracks, aod::McParticles& mcParticles)
+               soa::Join<aod::Tracks, aod::TracksExtra, aod::McTrackLabels>& tracks, aod::McParticles const& mcParticles)
  {

    //event selection
@@ -110,7 +110,7 @@ struct TrackCheckTaskEvSel {
        const auto particle = track.mcParticle();
        int pdgcode = particle.pdgCode();

-        if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { //is primary?
+        if (MC::isPhysicalPrimary(particle)) { //is primary?
          isPrimary = true;
        }

@@ -208,7 +208,7 @@ struct TrackCheckTaskEvSelTrackSel {
        const auto particle = track.mcParticle();
        int pdgcode = particle.pdgCode();

-        if (MC::isPhysicalPrimary<aod::McParticles>(particle)) { //is primary?
+        if (MC::isPhysicalPrimary(particle)) { //is primary?
          isPrimary = true;
        }
        //Calculate y

@@ -229,7 +229,7 @@ struct QaTrackingEfficiency {
        return true;
      }
      histos.fill(h, 5);
-      if ((selPrim == 1) && (!MC::isPhysicalPrimary<o2::aod::McParticles>(p))) { // Requiring is physical primary
+      if ((selPrim == 1) && (!MC::isPhysicalPrimary(p))) { // Requiring is physical primary
        return true;
      }
      histos.fill(h, 6);

@@ -203,7 +203,7 @@ struct QaTrackingKine {
      if (pdgCodeSel != 0 && particle.pdgCode() != pdgCodeSel) { // Checking PDG code
        continue;
      }
-      if (MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
        histos.fill(HIST("trackingPrm/pt"), t.pt());
        histos.fill(HIST("trackingPrm/eta"), t.eta());
        histos.fill(HIST("trackingPrm/phi"), t.phi());
@@ -223,7 +223,7 @@ struct QaTrackingKine {
      histos.fill(HIST("particle/pt"), particle.pt());
      histos.fill(HIST("particle/eta"), particle.eta());
      histos.fill(HIST("particle/phi"), particle.phi());
-      if (MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (MC::isPhysicalPrimary(particle)) {
        histos.fill(HIST("particlePrm/pt"), particle.pt());
        histos.fill(HIST("particlePrm/eta"), particle.eta());
        histos.fill(HIST("particlePrm/phi"), particle.phi());
@@ -376,7 +376,7 @@ struct QaTrackingResolution {
      if (pdgCodeSel != 0 && particle.pdgCode() != pdgCodeSel) {
        continue;
      }
-      if (useOnlyPhysicsPrimary && !MC::isPhysicalPrimary<o2::aod::McParticles>(particle)) {
+      if (useOnlyPhysicsPrimary && !MC::isPhysicalPrimary(particle)) {
        continue;
      }
      const double deltaPt = track.pt() - particle.pt();

@@ -38,14 +38,14 @@ struct AccessMcData {
  OutputObj<TH1F> etaH{TH1F("eta", "eta", 102, -2.01, 2.01)};

  // group according to McCollisions
-  void process(aod::McCollision const& mcCollision, aod::McParticles& mcParticles)
+  void process(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles)
  {
    // access MC truth information with mcCollision() and mcParticle() methods
    LOGF(info, "MC. vtx-z = %f", mcCollision.posZ());
    LOGF(info, "First: %d | Length: %d", mcParticles.begin().index(), mcParticles.size());
    int count = 0;
    for (auto& mcParticle : mcParticles) {
-      if (MC::isPhysicalPrimary<aod::McParticles>(mcParticle)) {
+      if (MC::isPhysicalPrimary(mcParticle)) {
        phiH->Fill(mcParticle.phi());
        etaH->Fill(mcParticle.eta());
        count++;
@@ -61,7 +61,8 @@ struct AccessMcTruth {
  OutputObj<TH1F> phiDiff{TH1F("phiDiff", ";phi_{MC} - phi_{Rec}", 100, -M_PI, M_PI)};

  // group according to reconstructed Collisions
-  void process(soa::Join<aod::Collisions, aod::McCollisionLabels>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions)
+  void process(soa::Join<aod::Collisions, aod::McCollisionLabels>::iterator const& collision, soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks,
+               aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions)
  {
    // access MC truth information with mcCollision() and mcParticle() methods
    LOGF(info, "vtx-z (data) = %f | vtx-z (MC) = %f", collision.posZ(), collision.mcCollision().posZ());
@@ -84,11 +85,34 @@ struct AccessMcTruth {
  }
 };

+// Loop over MCColisions and get corresponding collisions (there can be more than one)
+// For each of them get the corresponding tracks
+struct LoopOverMcMatched {
+  OutputObj<TH1F> etaDiff{TH1F("etaDiff", ";eta_{MC} - eta_{Rec}", 100, -2, 2)};
+  void process(aod::McCollision const& mcCollision, soa::Join<aod::McCollisionLabels, aod::Collisions> const& collisions,
+               soa::Join<aod::Tracks, aod::McTrackLabels> const& tracks, aod::McParticles const& mcParticles)
+  {
+    // access MC truth information with mcCollision() and mcParticle() methods
+    LOGF(info, "MC collision at vtx-z = %f with %d mc particles and %d reconstructed collisions", mcCollision.posZ(), mcParticles.size(), collisions.size());
+    for (auto& collision : collisions) {
+      LOGF(info, "  Reconstructed collision at vtx-z = %f", collision.posZ());
+
+      // NOTE this will be replaced by a improved grouping in the future
+      auto groupedTracks = tracks.sliceBy(aod::track::collisionId, collision.globalIndex());
+      LOGF(info, "  which has %d tracks", groupedTracks.size());
+      for (auto& track : groupedTracks) {
+        etaDiff->Fill(track.mcParticle().eta() - track.eta());
+      }
+    }
+  }
+};
+
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
  return WorkflowSpec{
    adaptAnalysisTask<VertexDistribution>(cfgc),
    adaptAnalysisTask<AccessMcData>(cfgc),
    adaptAnalysisTask<AccessMcTruth>(cfgc),
+    adaptAnalysisTask<LoopOverMcMatched>(cfgc),
  };
 }