PWGEM/PhotonMeson: implement searching window for PCM

PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx

@@ -60,7 +60,7 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
     list->Add(new TH1F("hNclsITS", "number of ITS clusters", 8, -0.5, 7.5));
     list->Add(new TH1F("hChi2ITS", "chi2/number of ITS clusters", 360, 0, 36));
     list->Add(new TH2F("hXY", "X vs. Y;X;Y", 100, 0, 100, 100, -50, 50));
-    list->Add(new TH2F("hZX", "Z vs. X;Z;X", 1000, -100, 100, 100, 0, 100));
+    list->Add(new TH2F("hZX", "Z vs. X;Z;X", 400, -100, 100, 100, 0, 100));
     list->Add(new TH2F("hZY", "Z vs. Y;Z;Y", 200, -100, 100, 100, -50, 50));
     list->Add(new TH2F("hDCAxyEta", "DCAxy vs. #eta;#eta;DCA_{xy} (cm)", 400, -2, +2, 100, -50, 50));
     list->Add(new TH2F("hDCAxyZ", "DCAxy vs. Z;Z (cm);DCA_{xy} (cm)", 200, -100, +100, 100, -50, 50));
@@ -85,6 +85,7 @@ void o2::aod::emphotonhistograms::DefineHistograms(THashList* list, const char*
     list->Add(new TH2F("hKFChi2vsY", "KF chi2 vs. recalc. conversion point in Y;Y (cm);KF chi2/NDF", 400, -200.0f, 200.0f, 100, 0.f, 100.0f));
     list->Add(new TH2F("hKFChi2vsZ", "KF chi2 vs. recalc. conversion point in Z;Z (cm);KF chi2/NDF", 500, -250.0f, 250.0f, 100, 0.f, 100.0f));
     list->Add(new TH1F("hNgamma", "Number of #gamma candidates per collision", 101, -0.5f, 100.5f));
+    list->Add(new TH2F("hV0R_minTrackX", "V0R vs. trackiu x;R_{xy} of V0 (cm);min TrackIU X (cm)", 100, 0.0f, 100.0f, 200, 0.f, 100.0f));
 
     const int nrxy = 102;
     double rxy[nrxy] = {0.f};

PWGEM/PhotonMeson/Core/V0PhotonCut.h

@@ -116,9 +116,13 @@ class V0PhotonCut : public TNamed
     auto pos = v0.template posTrack_as<TLeg>();
     auto ele = v0.template negTrack_as<TLeg>();
 
-    // float pos_rxy = sqrt(pow(pos.x(),2) + pow(pos.y(),2));
-    // float ele_rxy = sqrt(pow(ele.x(),2) + pow(ele.y(),2));
-    // if (v0.recalculatedVtxR() > std::min(pos_rxy, ele_rxy)) {
+    // if(pos.hasITS() && ele.hasITS()){
+    //   bool ret = sqrt(pow(v0.alpha() / 0.95, 2) + pow(v0.qtarm() / 0.02, 2) ) < 1.0;
+    //   if (!ret){
+    //     return false;
+    //   }
+    // }
+    // if (v0.recalculatedVtxR() > std::min(pos.x(), ele.x()) + 7.f && (pos.x() > 1.f && ele.x() > 1.f)) {
     //   return false;
     // }
 

PWGEM/PhotonMeson/TableProducer/CMakeLists.txt

@@ -53,3 +53,8 @@ o2physics_add_dpl_workflow(produce-meson-calo
                     SOURCES produceMesonCalo.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(check-v0-mc
+                    SOURCES checkV0MC.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)

PWGEM/PhotonMeson/TableProducer/checkV0MC.cxx

@@ -0,0 +1,155 @@
+// 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.
+
+/// \brief check MC true V0s
+/// \author daiki.sekihata@cern.ch
+
+#include <TMath.h>
+#include <TVector2.h>
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Common/Core/trackUtilities.h"
+#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
+#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h"
+#include "PWGEM/PhotonMeson/Utils/MCUtilities.h"
+
+using namespace o2;
+using namespace o2::soa;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+using MyTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksDCA, aod::TracksCovIU>;
+using MyTracksMC = soa::Join<MyTracks, aod::McTrackLabels>;
+
+struct checkV0MC {
+
+  Configurable<float> maxX{"maxX", 83.1, "maximum X (starting point of track X)"};
+  Configurable<float> maxY{"maxY", 20.0, "maximum X (starting point of track X)"};
+  Configurable<float> minpt{"minpt", 0.01, "min pt for track"};
+  Configurable<float> maxeta{"maxeta", 0.9, "eta acceptance"};
+  Configurable<float> dcamin{"dcamin", 0.1, "dcamin"};
+  Configurable<float> dcamax{"dcamax", 1e+10, "dcamax"};
+
+  HistogramRegistry fRegistry{
+    "fRegistry",
+    {
+      {"hDiffCollId", "difference in collision id between ele/pos", {HistType::kTH1F, {{101, -50.5f, +50.5f}}}},
+      {"hOrphanTrack", "Number of orphan track from true photon conversions", {HistType::kTH1F, {{1, 0.5f, +1.5f}}}},
+      {"hV0PhotonCorrX", "V0 photon track iu X correlation;track iu X for ele;track iu X for pos", {HistType::kTH2F, {{100, 0, +100}, {100, 0, +100}}}},
+      {"hV0PhotonCorrY", "V0 photon track iu Y correlation;track iu Y for ele;track iu Y diff", {HistType::kTH2F, {{200, -100, +100}, {400, -20, +20}}}},
+      {"hV0PhotonCorrZ", "V0 photon track iu Z correlation;track iu Z for ele;track iu Z diff", {HistType::kTH2F, {{200, -100, +100}, {400, -20, +20}}}},
+      {"hV0PhotonCorrYvsX", "V0 photon track iu Y correlation;track iu Y for ele;track iu Y diff", {HistType::kTH2F, {{100, 0, +100}, {400, -20, +20}}}},
+      {"hV0PhotonCorrZvsX", "V0 photon track iu Z correlation;track iu Z for ele;track iu Z diff", {HistType::kTH2F, {{100, 0, +100}, {400, -20, +20}}}},
+    },
+  };
+
+  template <typename TTrack>
+  bool isSelected(TTrack const& track)
+  {
+    if (track.pt() < minpt || abs(track.eta()) > maxeta) {
+      return false;
+    }
+    if (abs(track.dcaXY()) < dcamin || dcamax < abs(track.dcaXY())) {
+      return false;
+    }
+    if (!track.hasITS() && !track.hasTPC()) {
+      return false;
+    }
+
+    if (track.hasTPC()) {
+      if (track.tpcSignal() < 40.f || 110.f < track.tpcSignal()) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  Filter trackFilter = o2::aod::track::x < maxX && nabs(o2::aod::track::y) < maxY && o2::aod::track::pt > minpt&& nabs(o2::aod::track::eta) < maxeta&& dcamin < nabs(o2::aod::track::dcaXY) && nabs(o2::aod::track::dcaXY) < dcamax;
+  using MyFilteredTracksMC = soa::Filtered<MyTracksMC>;
+  Partition<MyFilteredTracksMC> posTracks = o2::aod::track::signed1Pt > 0.f;
+  Partition<MyFilteredTracksMC> negTracks = o2::aod::track::signed1Pt < 0.f;
+  Partition<MyFilteredTracksMC> orphan_posTracks = o2::aod::track::signed1Pt > 0.f && o2::aod::track::collisionId < 0;
+  Partition<MyFilteredTracksMC> orphan_negTracks = o2::aod::track::signed1Pt < 0.f && o2::aod::track::collisionId < 0;
+
+  int ndf = 0;
+
+  using MyCollisions = soa::Join<aod::McCollisionLabels, aod::Collisions>;
+  void processMC(MyCollisions const& collisions, aod::McCollisions const&, aod::BCsWithTimestamps const& bcs, MyFilteredTracksMC const& tracks, aod::McParticles const& mcTracks)
+  {
+    LOGF(info, "n pos tracks = %d , n neg tracks = %d", posTracks.size(), negTracks.size());
+    ndf++;
+    if (ndf > 3) { // 3 DFs are enough.
+      return;
+    }
+
+    for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks, posTracks))) {
+      if (!ele.has_mcParticle() || !pos.has_mcParticle()) {
+        continue;
+      }
+
+      if (!isSelected(ele) || !isSelected(pos)) {
+        continue;
+      }
+
+      auto elemc = ele.template mcParticle_as<aod::McParticles>();
+      auto posmc = pos.template mcParticle_as<aod::McParticles>();
+
+      int photonid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcTracks);
+      if (photonid < 0) { // check swap, true electron is reconstructed as positron and vice versa.
+        photonid = FindCommonMotherFrom2Prongs(posmc, elemc, 11, -11, 22, mcTracks);
+      }
+      if (photonid < 0) {
+        continue;
+      }
+      // At this point, there are only true photon conversions.
+      auto photonmc = mcTracks.iteratorAt(photonid);
+      if (!IsPhysicalPrimary(photonmc.mcCollision(), photonmc, mcTracks)) {
+        continue;
+      }
+
+      bool has_coll_ele = ele.has_collision();
+      bool has_coll_pos = pos.has_collision();
+
+      if (!has_coll_ele) {
+        fRegistry.fill(HIST("hOrphanTrack"), 1);
+      }
+      if (!has_coll_pos) {
+        fRegistry.fill(HIST("hOrphanTrack"), 1);
+      }
+
+      if (!has_coll_ele || !has_coll_pos) {
+        continue;
+      }
+
+      auto collision_ele = ele.collision();
+      auto collision_pos = pos.collision();
+      int diff = collision_ele.globalIndex() - collision_pos.globalIndex();
+      fRegistry.fill(HIST("hDiffCollId"), diff);
+
+      fRegistry.fill(HIST("hV0PhotonCorrX"), ele.x(), pos.x());
+      fRegistry.fill(HIST("hV0PhotonCorrY"), ele.y(), pos.y() - ele.y());
+      fRegistry.fill(HIST("hV0PhotonCorrZ"), ele.z(), pos.z() - ele.z());
+      fRegistry.fill(HIST("hV0PhotonCorrYvsX"), ele.x(), pos.y() - ele.y());
+      fRegistry.fill(HIST("hV0PhotonCorrZvsX"), ele.x(), pos.z() - ele.z());
+
+    } // end of pairing loop
+  }
+  PROCESS_SWITCH(checkV0MC, processMC, "process mc truth info", true);
+
+  void processDummy(soa::Join<aod::McCollisionLabels, aod::Collisions> const& collisions) {}
+  PROCESS_SWITCH(checkV0MC, processDummy, "process dummy", false);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<checkV0MC>(cfgc, TaskName{"check-v0-mc"})};
+}

PWGEM/PhotonMeson/TableProducer/createEMReducedEvent.cxx

@@ -46,7 +46,7 @@ struct createEMReducedEvent {
 
   HistogramRegistry registry{"registry"};
 
-  Preslice<aod::V0Photons> perCollision_pcm = aod::v0photon::collisionId;
+  Preslice<aod::V0PhotonsKF> perCollision_pcm = aod::v0photon::collisionId;
   Preslice<aod::PHOSClusters> perCollision_phos = aod::skimmedcluster::collisionId;
   Preslice<aod::SkimEMCClusters> perCollision_emc = aod::skimmedcluster::collisionId;
 
@@ -130,7 +130,7 @@ struct createEMReducedEvent {
 
   } // end of process
 
-  void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons)
+  void process_PCM(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons)
   {
     process<kPCM>(collisions, bcs, v0photons, nullptr, nullptr);
   }
@@ -142,12 +142,12 @@ struct createEMReducedEvent {
   {
     process<kEMC>(collisions, bcs, nullptr, nullptr, emcclusters);
   }
-  void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters)
+  void process_PCM_PHOS(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters)
   {
     const uint8_t sysflag = kPCM | kPHOS;
     process<sysflag>(collisions, bcs, v0photons, phosclusters, nullptr);
   }
-  void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::SkimEMCClusters const& emcclusters)
+  void process_PCM_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::SkimEMCClusters const& emcclusters)
   {
     const uint8_t sysflag = kPCM | kEMC;
     process<sysflag>(collisions, bcs, v0photons, nullptr, emcclusters);
@@ -157,7 +157,7 @@ struct createEMReducedEvent {
     const uint8_t sysflag = kPHOS | kEMC;
     process<sysflag>(collisions, bcs, nullptr, phosclusters, emcclusters);
   }
-  void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0Photons const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
+  void process_PCM_PHOS_EMC(MyCollisions const& collisions, aod::BCs const& bcs, aod::V0PhotonsKF const& v0photons, aod::PHOSClusters const& phosclusters, aod::SkimEMCClusters const& emcclusters)
   {
     const uint8_t sysflag = kPCM | kPHOS | kEMC;
     process<sysflag>(collisions, bcs, v0photons, phosclusters, emcclusters);