PWGEM/Dilepton: fix dimuonQCMC

PWGEM/Dilepton/Core/DielectronCut.h

@@ -40,6 +40,7 @@ class DielectronCut : public TNamed
  public:
   DielectronCut() = default;
   DielectronCut(const char* name, const char* title) : TNamed(name, title) {}
+  ~DielectronCut() {}
 
   enum class DielectronCuts : int {
     // pair cut

PWGEM/Dilepton/Core/DimuonCut.h

@@ -27,13 +27,18 @@
 #include "Framework/Logger.h"
 #include "Framework/DataTypes.h"
 #include "CommonConstants/PhysicsConstants.h"
+#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h"
+
+using namespace o2::aod::pwgem::dilepton::utils::emtrackutil;
 
 class DimuonCut : public TNamed
 {
  public:
   DimuonCut() = default;
   DimuonCut(const char* name, const char* title) : TNamed(name, title) {}
 
+  ~DimuonCut() {}
+
   enum class DimuonCuts : int {
     // pair cut
     kMass = 0,
@@ -79,36 +84,9 @@ class DimuonCut : public TNamed
     ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon);
     ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
 
-    float dcaX1 = t1.fwdDcaX();
-    float dcaY1 = t1.fwdDcaY();
-    float cXX1 = t1.cXX();
-    float cYY1 = t1.cYY();
-    float cXY1 = t1.cXY();
-    float dcaX2 = t2.fwdDcaX();
-    float dcaY2 = t2.fwdDcaY();
-    float cXX2 = t2.cXX();
-    float cYY2 = t2.cYY();
-    float cXY2 = t2.cXY();
-
-    float dca_xy1 = 999.f;
-    float det1 = cXX1 * cYY1 - cXY1 * cXY1;
-    if (det1 < 0) {
-      dca_xy1 = 999.f;
-    } else {
-      float chi2 = (dcaX1 * dcaX1 * cYY1 + dcaY1 * dcaY1 * cXX1 - 2. * dcaX1 * dcaY1 * cXY1) / det1;
-      dca_xy1 = std::sqrt(std::abs(chi2) / 2.); // in sigma
-    }
-
-    float dca_xy2 = 999.f;
-    float det2 = cXX2 * cYY2 - cXY2 * cXY2;
-    if (det2 < 0) {
-      dca_xy2 = 999.f;
-    } else {
-      float chi2 = (dcaX2 * dcaX2 * cYY2 + dcaY2 * dcaY2 * cXX2 - 2. * dcaX2 * dcaY2 * cXY2) / det2;
-      dca_xy2 = std::sqrt(std::abs(chi2) / 2.); // in sigma
-    }
-
-    float pair_dca_xy = std::sqrt((std::pow(dca_xy1, 2) + std::pow(dca_xy2, 2)) / 2.);
+    float dca_xy_t1 = fwdDcaXYinSigma(t1);
+    float dca_xy_t2 = fwdDcaXYinSigma(t2);
+    float pair_dca_xy = std::sqrt((dca_xy_t1 * dca_xy_t1 + dca_xy_t2 * dca_xy_t2) / 2.);
 
     if (v12.M() < mMinMass || mMaxMass < v12.M()) {
       return false;

PWGEM/Dilepton/Core/EMEventCut.h

@@ -27,6 +27,7 @@ class EMEventCut : public TNamed
  public:
   EMEventCut() = default;
   EMEventCut(const char* name, const char* title) : TNamed(name, title) {}
+  ~EMEventCut() {}
 
   enum class EMEventCuts : int {
     kSel8 = 0,

PWGEM/Dilepton/Core/PhotonHBT.h

@@ -207,6 +207,10 @@ struct PhotonHBT {
     used_photonIds.shrink_to_fit();
     used_dileptonIds.clear();
     used_dileptonIds.shrink_to_fit();
+
+    if (eid_bdt) {
+      delete eid_bdt;
+    }
   }
 
   HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
@@ -298,12 +302,13 @@ struct PhotonHBT {
     const AxisSpec axis_kt{ConfKtBins, "k_{T} (GeV/c)"};
 
     const AxisSpec axis_qinv{30, 0.0, +0.3, "q_{inv} (GeV/c)"};
-    const AxisSpec axis_qout_cms{60, -0.3, +0.3, "q_{out}^{CMS} (GeV/c)"};
-    const AxisSpec axis_qside_cms{60, -0.3, +0.3, "q_{side}^{CMS} (GeV/c)"};
-    const AxisSpec axis_qlong_cms{60, -0.3, +0.3, "q_{long}^{CMS} (GeV/c)"};
-    const AxisSpec axis_qlong_lcms{60, -0.3, +0.3, "q_{long}^{LCMS} (GeV/c)"};
+    const AxisSpec axis_qabs_lcms{30, 0.0, +0.3, "|q|^{LCMS} (GeV/c)"};
+    const AxisSpec axis_qout{60, -0.3, +0.3, "q_{out} (GeV/c)"};   // qout does not change between LAB and LCMS frame
+    const AxisSpec axis_qside{60, -0.3, +0.3, "q_{side} (GeV/c)"}; // qside does not change between LAB and LCMS frame
+    const AxisSpec axis_qlong{60, -0.3, +0.3, "q_{long} (GeV/c)"};
 
-    fRegistry.add("Pair/same/hs", "diphoton correlation", kTHnSparseD, {axis_kt, axis_qinv, axis_qout_cms, axis_qside_cms, axis_qlong_cms, axis_qlong_lcms}, true);
+    fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axis_kt, axis_qinv, axis_qabs_lcms}, true);
+    fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axis_kt, axis_qout, axis_qside, axis_qlong}, true);
     fRegistry.addClone("Pair/same/", "Pair/mix/");
   }
 
@@ -376,6 +381,7 @@ struct PhotonHBT {
     }
   }
 
+  o2::ml::OnnxModel* eid_bdt = nullptr;
   void DefineDileptonCut()
   {
     fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut");
@@ -414,7 +420,7 @@ struct PhotonHBT {
     fDielectronCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl);
 
     if (dielectroncuts.cfg_pid_scheme == static_cast<int>(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut
-      o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel();
+      eid_bdt = new o2::ml::OnnxModel();
       if (dielectroncuts.loadModelsFromCCDB) {
         ccdbApi.init(ccdburl);
         std::map<std::string, std::string> metadata;
@@ -453,14 +459,12 @@ struct PhotonHBT {
     ROOT::Math::PtEtaPhiMVector k12 = 0.5 * (v1 + v2);
     float qinv = -q12.M();
     float kt = k12.Pt();
-    float qlong_cms = q12.Pz();
 
-    ROOT::Math::XYZVector q_3d = q12.Vect();                                   // 3D q vector
+    // ROOT::Math::XYZVector q_3d = q12.Vect();                                   // 3D q vector
     ROOT::Math::XYZVector uv_out(k12.Px() / k12.Pt(), k12.Py() / k12.Pt(), 0); // unit vector for out. i.e. parallel to kt
     ROOT::Math::XYZVector uv_long(0, 0, 1);                                    // unit vector for long, beam axis
     ROOT::Math::XYZVector uv_side = uv_out.Cross(uv_long);                     // unit vector for side
-    float qout_cms = q_3d.Dot(uv_out);
-    float qside_cms = q_3d.Dot(uv_side);
+    // float qlong_lab = q_3d.Dot(uv_long);
 
     // longitudinally co-moving system (LCMS)
     ROOT::Math::PxPyPzEVector v1_cartesian(v1.Px(), v1.Py(), v1.Pz(), v1.E());
@@ -469,17 +473,28 @@ struct PhotonHBT {
     float beta_z = (v1 + v2).Pz() / (v1 + v2).E();
     ROOT::Math::Boost bst_z(0, 0, -beta_z); // Boost supports only PxPyPzEVector
     ROOT::Math::PxPyPzEVector q12_lcms = bst_z(q12_cartesian);
-    float qlong_lcms = q12_lcms.Pz();
+    ROOT::Math::XYZVector q_3d_lcms = q12_lcms.Vect(); // 3D q vector in LCMS
+    float qout_lcms = q_3d_lcms.Dot(uv_out);
+    float qside_lcms = q_3d_lcms.Dot(uv_side);
+    float qlong_lcms = q_3d_lcms.Dot(uv_long);
+    float qabs_lcms = q_3d_lcms.R();
 
     // ROOT::Math::PxPyPzEVector v1_lcms_cartesian = bst_z(v1_cartesian);
     // ROOT::Math::PxPyPzEVector v2_lcms_cartesian = bst_z(v2_cartesian);
     // ROOT::Math::PxPyPzEVector q12_lcms_cartesian = bst_z(q12_cartesian);
-    // LOGF(info, "q12.Pz() = %f, q12_cartesian.Pz() = %f",q12.Pz(), q12_cartesian.Pz());
-    // LOGF(info, "v1.Pz() = %f, v2.Pz() = %f",v1.Pz(), v2.Pz());
-    // LOGF(info, "v1_lcms_cartesian.Pz() = %f, v2_lcms_cartesian.Pz() = %f",v1_lcms_cartesian.Pz(), v2_lcms_cartesian.Pz());
+    // LOGF(info, "q12.Pz() = %f, q12_cartesian.Pz() = %f", q12.Pz(), q12_cartesian.Pz());
+    // LOGF(info, "v1.Pz() = %f, v2.Pz() = %f", v1.Pz(), v2.Pz());
+    // LOGF(info, "v1_lcms_cartesian.Pz() = %f, v2_lcms_cartesian.Pz() = %f", v1_lcms_cartesian.Pz(), v2_lcms_cartesian.Pz());
     // LOGF(info, "q12_lcms_cartesian.Pz() = %f", q12_lcms_cartesian.Pz());
-
-    fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs"), kt, qinv, qout_cms, qside_cms, qlong_cms, qlong_lcms);
+    // LOGF(info, "q_3d_lcms.Dot(uv_out) = %f, q_3d_lcms.Dot(uv_side) = %f, q_3d.Dot(uv_out) = %f, q_3d.Dot(uv_side) = %f", q_3d_lcms.Dot(uv_out), q_3d_lcms.Dot(uv_side), q_3d.Dot(uv_out), q_3d.Dot(uv_side));
+    // LOGF(info, "q12_lcms.Pz() = %f, q_3d_lcms.Dot(uv_long) = %f", q12_lcms.Pz(), q_3d_lcms.Dot(uv_long));
+    // ROOT::Math::PxPyPzEVector q12_lcms_tmp = bst_z(v1_cartesian) - bst_z(v2_cartesian);
+    // LOGF(info, "q12_lcms.Px() = %f, q12_lcms.Py() = %f, q12_lcms.Pz() = %f, q12_lcms_tmp.Px() = %f, q12_lcms_tmp.Py() = %f, q12_lcms_tmp.Pz() = %f", q12_lcms.Px(), q12_lcms.Py(), q12_lcms.Pz(), q12_lcms_tmp.Px(), q12_lcms_tmp.Py(), q12_lcms_tmp.Pz());
+    // float qabs_lcms_tmp = q12_lcms.P();
+    // LOGF(info, "qabs_lcms = %f, qabs_lcms_tmp = %f", qabs_lcms, qabs_lcms_tmp);
+
+    fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), kt, qinv, qabs_lcms);
+    fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), kt, qout_lcms, qside_lcms, qlong_lcms);
   }
 
   template <typename TCollisions, typename TPhotons1, typename TPhotons2, typename TSubInfos1, typename TSubInfos2, typename TPreslice1, typename TPreslice2, typename TCut1, typename TCut2>

PWGEM/Dilepton/DataModel/dileptonTables.h

@@ -9,6 +9,8 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+#include <string>
+#include <unordered_map>
 #include "Common/Core/RecoDecay.h"
 #include "Framework/AnalysisDataModel.h"
 #include "Common/DataModel/PIDResponse.h"
@@ -25,9 +27,39 @@
 namespace o2::aod
 {
 
+namespace pwgem::dilepton::swt
+{
+enum class swtAliases : int { // software trigger aliases for EM
+  kUnDef = 0,
+  kHighTrackMult = 1,
+  kHighFt0Mult,
+  kSingleE,
+  kLMeeIMR,
+  kLMeeHMR,
+  kDiElectron,
+  kSingleMuLow,
+  kSingleMuHigh,
+  kDiMuon,
+  kNaliases
+};
+
+const std::unordered_map<std::string, int> aliasLabels = {
+  {"fHighTrackMult", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kHighTrackMult)},
+  {"fHighFt0Mult", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kHighFt0Mult)},
+  {"fSingleE", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleE)},
+  {"fLMeeIMR", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kLMeeIMR)},
+  {"fLMeeHMR", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kLMeeHMR)},
+  {"fDiElectron", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kDiElectron)},
+  {"fSingleMuLow", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleMuLow)},
+  {"fSingleMuHigh", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleMuHigh)},
+  {"fDiMuon", static_cast<int>(o2::aod::pwgem::dilepton::swt::swtAliases::kDiMuon)},
+};
+} // namespace pwgem::dilepton::swt
+
 namespace emevent
 {
 DECLARE_SOA_COLUMN(CollisionId, collisionId, int);
+DECLARE_SOA_BITMAP_COLUMN(SWTAlias, swtalias, 16); //! Bitmask of fired trigger aliases (see above for definitions)
 DECLARE_SOA_COLUMN(NeeULS, neeuls, int);
 DECLARE_SOA_COLUMN(NeeLSpp, neelspp, int);
 DECLARE_SOA_COLUMN(NeeLSmm, neelsmm, int);
@@ -134,6 +166,10 @@ DECLARE_SOA_TABLE(EMEventsQvec, "AOD", "EMEVENTQVEC", //!   event q vector table
                   emevent::EP3BTot<emevent::Q3xBTot, emevent::Q3yBTot>);
 using EMEventQvec = EMEventsQvec::iterator;
 
+DECLARE_SOA_TABLE(EMSWTriggerBits, "AOD", "EMSWTRIGGERBIT", //!
+                  emevent::SWTAlias);
+using EMSWTriggerBit = EMSWTriggerBits::iterator;
+
 DECLARE_SOA_TABLE(EMEventsNee, "AOD", "EMEVENTNEE", emevent::NeeULS, emevent::NeeLSpp, emevent::NeeLSmm); // joinable to EMEvents
 using EMEventNee = EMEventsNee::iterator;
 

PWGEM/Dilepton/TableProducer/CMakeLists.txt

@@ -47,7 +47,7 @@ o2physics_add_dpl_workflow(skimmer-secondary-electron
 
 o2physics_add_dpl_workflow(create-emevent-dilepton
                     SOURCES createEMEventDilepton.cxx
-                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(associate-mc-info-dilepton

PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx

@@ -24,6 +24,7 @@
 #include "DataFormatsParameters/GRPObject.h"
 #include "DataFormatsParameters/GRPMagField.h"
 #include "CCDB/BasicCCDBManager.h"
+#include "EventFiltering/Zorro.h"
 
 #include "PWGEM/Dilepton/DataModel/dileptonTables.h"
 #include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
@@ -50,6 +51,7 @@ struct CreateEMEventDilepton {
   Produces<o2::aod::EMEventsMult> event_mult;
   Produces<o2::aod::EMEventsCent> event_cent;
   Produces<o2::aod::EMEventsQvec> event_qvec;
+  Produces<o2::aod::EMSWTriggerBits> emswtbit;
 
   enum class EMEventType : int {
     kEvent = 0,
@@ -64,6 +66,8 @@ struct CreateEMEventDilepton {
   Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
   Configurable<double> d_bz_input{"d_bz", -999, "bz field, -999 is automatic"};
   Configurable<bool> applyEveSel_at_skimming{"applyEveSel_at_skimming", false, "flag to apply minimal event selection at the skimming level"};
+  Configurable<bool> enable_swt{"enable_swt", false, "flag to process skimmed data (swt triggered)"};
+  Configurable<std::string> cfg_swt_names{"cfg_swt_names", "fHighTrackMult,fHighFt0Mult", "comma-separated software trigger names"}; // !trigger names have to be pre-registered in dileptonTable.h for bit operation!
 
   HistogramRegistry registry{"registry"};
   void init(o2::framework::InitContext&)
@@ -73,6 +77,14 @@ struct CreateEMEventDilepton {
     hEventCounter->GetXaxis()->SetBinLabel(2, "sel8");
   }
 
+  ~CreateEMEventDilepton()
+  {
+    swt_names.clear();
+    swt_names.shrink_to_fit();
+  }
+
+  Zorro zorro;
+  std::vector<std::string> swt_names;
   int mRunNumber;
   float d_bz;
   Service<o2::ccdb::BasicCCDBManager> ccdb;
@@ -84,6 +96,16 @@ struct CreateEMEventDilepton {
       return;
     }
 
+    if (enable_swt) {
+      LOGF(info, "enable software triggers : %s", cfg_swt_names.value.data());
+      zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), cfg_swt_names.value);
+      std::stringstream tokenizer(cfg_swt_names.value);
+      std::string token;
+      while (std::getline(tokenizer, token, ',')) {
+        swt_names.emplace_back(token);
+      }
+    }
+
     // In case override, don't proceed, please - no CCDB access required
     if (d_bz_input > -990) {
       d_bz = d_bz_input;
@@ -137,13 +159,27 @@ struct CreateEMEventDilepton {
         continue;
       }
 
+      if (enable_swt) {
+        if (zorro.isSelected(bc.globalBC())) { // triggered event
+          uint16_t trigger_bitmap = 0;
+          for (auto& swtname : swt_names) {
+            LOGF(info, "swtname = %s , swt index = %d", swtname.data(), o2::aod::pwgem::dilepton::swt::aliasLabels.at(swtname));
+            trigger_bitmap |= BIT(o2::aod::pwgem::dilepton::swt::aliasLabels.at(swtname));
+          } // end of desired trigger names loop
+          LOGF(info, "trigger_bitmap = %d", trigger_bitmap);
+          emswtbit(trigger_bitmap);
+        } else { // rejected
+          continue;
+        }
+      }
+      // LOGF(info, "collision.multNTracksPV() = %d, collision.multFT0A() = %f, collision.multFT0C() = %f", collision.multNTracksPV(), collision.multFT0A(), collision.multFT0C());
+
       registry.fill(HIST("hEventCounter"), 1);
 
       if (collision.sel8()) {
         registry.fill(HIST("hEventCounter"), 2);
       }
 
-      // uint64_t tag = collision.selection_raw();
       event(collision.globalIndex(), bc.runNumber(), bc.globalBC(), collision.alias_raw(), collision.selection_raw(), bc.timestamp(),
             collision.posX(), collision.posY(), collision.posZ(),
             collision.numContrib(), collision.trackOccupancyInTimeRange());

PWGEM/Dilepton/Tasks/dielectronQC.cxx

@@ -277,6 +277,10 @@ struct dielectronQC {
 
     used_trackIds.clear();
     used_trackIds.shrink_to_fit();
+
+    if (eid_bdt) {
+      delete eid_bdt;
+    }
   }
 
   void addhistograms()
@@ -412,6 +416,7 @@ struct dielectronQC {
     fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax);
   }
 
+  o2::ml::OnnxModel* eid_bdt = nullptr;
   void DefineDileptonCut()
   {
     fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut");
@@ -448,7 +453,7 @@ struct dielectronQC {
     fDielectonCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl);
 
     if (dielectroncuts.cfg_pid_scheme == static_cast<int>(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut
-      o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel();
+      eid_bdt = new o2::ml::OnnxModel();
       if (dielectroncuts.loadModelsFromCCDB) {
         ccdbApi.init(ccdburl);
         std::map<std::string, std::string> metadata;