[PWGCF] adding efficiency correction and too for systematic studies

PWGCF/JCorran/Tasks/jEPFlowAnalysis.cxx

@@ -39,7 +39,7 @@
 
 using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Qvectors>;
 using MyCollisionsWithSC = soa::Join<aod::Collisions, aod::EvSels, aod::QvectorsShifteds>;
-using MyTracks = aod::Tracks;
+using MyTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection>;
 
 struct jEPFlowAnalysis {
 
@@ -50,9 +50,17 @@
   Service<o2::ccdb::BasicCCDBManager> ccdb;
   o2::ccdb::CcdbApi ccdbApi;
 
+  struct : ConfigurableGroup {
+    Configurable<std::string> cfgURL{"cfgURL",
+                                     "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"};
+    Configurable<int64_t> nolaterthan{"ccdb-no-later-than",
+                                      std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(),
+                                      "Latest acceptable timestamp of creation for the object"};
+  } cfgCcdbParam;
+
   // Set Configurables here
   struct : ConfigurableGroup {
-    Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "Minimum pT used for track selection."};
+    Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "Minimum pT used for track seletion."};
     Configurable<float> cfgEtaMax{"cfgEtaMax", 1.f, "Maximum eta used for track selection."};
   } cfgTrackCuts;
 
@@ -61,6 +69,18 @@
   Configurable<int> cfgMaxOccupancy{"cfgMaxOccupancy", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
   Configurable<int> cfgMinOccupancy{"cfgMinOccupancy", 0, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
 
+  Configurable<bool> cfgEffCor{"cfgEffCor", false, "flag for efficiency correction"};
+  Configurable<std::string> cfgEffCorDir{"cfgEffCorDir", "Users/n/nmallick/Run3OO/Eff/LHC25h3b_FT0C", "path for efficiency correction"};
+
+  Configurable<bool> cfgSystStudy{"cfgSystStudy", false, "flag for syst study"};
+  Configurable<int> cfgITSNCls{"cfgITSNCls", 5, "minimum number of its clusters"};
+  Configurable<int> cfgTPCNclsCR{"cfgTPCNclsCR", 70, "minimum number of tpc cluster crossed rows"};
+  Configurable<float> cfgTPCChi2{"cfgTPCChi2", 4.0, "maximum TPC chi2"};
+  Configurable<float> cfgITSChi2{"cfgITSChi2", 36.0, "maximum ITS chi2"};
+  Configurable<float> cfgdcaZ{"cfgdcaZ", 2.0, "maximum dca z"};
+  Configurable<float> cfgdcaXY0{"cfgdcaXY0", 0.0105, "maximum constant dca xy"};
+  Configurable<float> cfgdcaXY1{"cfgdcaXY1", 0.035, "maximum pt deepdent dca xy"};
+
   Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
   Configurable<int> cfgnMode{"cfgnMode", 1, "the number of modulations"};
 
@@ -90,6 +110,8 @@
   std::vector<TProfile3D*> shiftprofile{};
   std::string fullCCDBShiftCorrPath;
 
+  THn* effMap = nullptr;
+
   template <typename T>
   int getdetId(const T& name)
   {
@@ -112,8 +134,53 @@
     }
   }
 
+  template <typename Trk>
+  uint8_t trackSel(const Trk& track)
+  {
+    uint8_t tracksel = 0;
+    if (!track.isGlobalTrack()) {
+      tracksel += 1;
+    }
+    if (track.itsNCls() <= cfgITSNCls && cfgSystStudy) {
+      tracksel += 2;
+    }
+    if (track.tpcNClsCrossedRows() <= cfgTPCNclsCR && cfgSystStudy) {
+      tracksel += 4;
+    }
+    if (track.tpcChi2NCl() >= cfgTPCChi2 && cfgSystStudy) {
+      tracksel += 8;
+    }
+    if (track.itsChi2NCl() >= cfgITSChi2 && cfgSystStudy) {
+      tracksel += 16;
+    }
+    if (std::abs(track.dcaZ()) >= cfgdcaZ && cfgSystStudy) {
+      tracksel += 32;
+    }
+    if (std::abs(track.dcaXY()) >= cfgdcaXY0 + cfgdcaXY1 / std::pow(track.pt(), 1.1) && cfgSystStudy) {
+      tracksel += 64;
+    }
+
+    return tracksel;
+  }
+
+  double getEfficiencyCorrection(THn* eff, float eta, float pt, float multiplicity, float posZ)
+  {
+    int effVars[4];
+    effVars[0] = eff->GetAxis(0)->FindBin(eta);
+    effVars[1] = eff->GetAxis(1)->FindBin(pt);
+    effVars[2] = eff->GetAxis(2)->FindBin(multiplicity);
+    effVars[3] = eff->GetAxis(3)->FindBin(posZ);
+    return eff->GetBinContent(effVars);
+  }
+
   void init(InitContext const&)
   {
+    ccdb->setURL(cfgCcdbParam.cfgURL);
+    ccdbApi.init("http://alice-ccdb.cern.ch");
+    ccdb->setCaching(true);
+    ccdb->setLocalObjectValidityChecking();
+    ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());    
+
     detId = getdetId(cfgDetName);
     refAId = getdetId(cfgRefAName);
     refBId = getdetId(cfgRefBName);
@@ -183,6 +250,11 @@
         return;
     }
 
+    if (cfgEffCor) {
+      auto bc = coll.bc_as<aod::BCsWithTimestamps>();
+      effMap = ccdb->getForTimeStamp<THnT<float>>(cfgEffCorDir, bc.timestamp());
+    }
+
     float cent = coll.cent();
     epFlowHistograms.fill(HIST("hCentrality"), cent);
     epFlowHistograms.fill(HIST("hVertex"), coll.posZ());
@@ -219,6 +291,11 @@
       float weight = 1.0;
 
       for (const auto& track : tracks) {
+        if (trackSel(track)) continue;
+        if (cfgEffCor) {
+          weight /= getEfficiencyCorrection(effMap, track.eta(), track.pt(), cent, coll.posZ()); 
+        }
+
         float vn = std::cos((i + 2) * (track.phi() - eps[0]));
         float vnSin = std::sin((i + 2) * (track.phi() - eps[0]));
 

PWGLF/Tasks/Resonances/heptaquark.cxx

@@ -68,6 +68,7 @@ struct heptaquark {
 
   Configurable<float> cfgSoftFraction{"cfgSoftFraction", 0.01, "Minimum allowed softest fraction"};
   Configurable<float> cfgCollinear{"cfgCollinear", 0.98, "Maximum allowed collinear selection"};
+  Configurable<float> cfgCosPoint{"cfgCosPoint", 0.95, "Minimum pointing angle selection"};
 
   ConfigurableAxis massAxis{"massAxis", {600, 2.8, 3.4}, "Invariant mass axis"};
   ConfigurableAxis ptAxis{"ptAxis", {VARIABLE_WIDTH, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0, 100.0}, "Transverse momentum bins"};
@@ -97,6 +98,18 @@ struct heptaquark {
     histos.add("hDalitzRot", "hDalitzRot", {HistType::kTHnSparseF, {massPPAxis, massPLAxis, massAxis, ptAxis, {2, -0.5f, 1.5f}, centAxis}});
   }
 
+  template <typename HQRow>
+  static inline TLorentzVector makeP4FromHQRow(HQRow const& hq)
+  {
+    const double px = hq.hqPx();
+    const double py = hq.hqPy();
+    const double pz = hq.hqPz();
+    const double m  = hq.hqMass();
+    TLorentzVector v;
+    v.SetXYZM(px, py, pz, m);
+    return v;
+  }
+
   double massLambda = o2::constants::physics::MassLambda;
   double massPr = o2::constants::physics::MassProton;
   double massPi = o2::constants::physics::MassPionCharged;
@@ -185,37 +198,52 @@ struct heptaquark {
     return false;
   }
 
-  template <typename HQ1, typename HQ2, typename HQ3>
-  int selectHQ(HQ1 const& hq1, HQ2 const& hq2, HQ3 const& hq3)
+  template <typename HQRow1, typename HQRow2, typename HQRow3, typename ColRow>
+  int selectHQ(HQRow1 const& hq1r, HQRow2 const& hq2r, HQRow3 const& hq3r, ColRow const& col)
   {
     int selection = 0;
-    if (hq1.Pt() < cfgMinPhiPt || hq2.Pt() < cfgMinPhiPt || hq3.Pt() < cfgMinLambdaPt)
+
+    auto hq1 = makeP4FromHQRow(hq1r);
+    auto hq2 = makeP4FromHQRow(hq2r);
+    auto hq3 = makeP4FromHQRow(hq3r);
+
+    if (hq1.Pt() < cfgMinPhiPt || hq2.Pt() < cfgMinPhiPt || hq3.Pt() < cfgMinLambdaPt) {
       selection += 1;
+    }
 
-    double sumE = hq1.E() + hq2.E() + hq3.E();
-    double emin = std::min({hq1.E(), hq2.E(), hq3.E()});
-    double fmin = emin / std::max(1e-9, sumE);
-    if (fmin < cfgSoftFraction)
+    const double sumE = hq1.E() + hq2.E() + hq3.E();
+    const double emin = std::min({hq1.E(), hq2.E(), hq3.E()});
+    const double fmin = emin / std::max(1e-9, sumE);
+    if (fmin < cfgSoftFraction) {
       selection += 2;
+    }
 
     auto ex = hq1 + hq2 + hq3;
     TVector3 boost = -ex.BoostVector();
+
     auto hqphipair_boost = hq1 + hq2;
-    auto hqlambda_boost = hq3;
+    auto hqlambda_boost  = hq3;
     hqphipair_boost.Boost(boost);
     hqlambda_boost.Boost(boost);
-    double cosHel = hqlambda_boost.Vect().Dot(hqphipair_boost.Vect()) / (hqlambda_boost.Vect().Mag() * hqphipair_boost.Vect().Mag());
-    if (std::abs(cosHel) > cfgCollinear)
+
+    const double denom = (hqlambda_boost.Vect().Mag() * hqphipair_boost.Vect().Mag());
+    const double cosHel = (denom > 0.) ? (hqlambda_boost.Vect().Dot(hqphipair_boost.Vect()) / denom) : 1.0;
+    if (std::abs(cosHel) > cfgCollinear) {
       selection += 4;
-    /*
-        ROOT::Math::XYZVector rPV(col.posX(), col.posY(), col.posZ());
-        ROOT::Math::XYZVector rSV(hq3.hqx(), hq3.hqy(), hq3.hqz());
-        ROOT::Math::XYZVector L = rSV - rPV;
-        ROOT::Math::XYZVector exMom(ex.Px(), ex.Py(), ex.Pz());
-        double cosPoint = L.Dot(exMom) / (L.R() * pEx.R() + 1e-9);
-        if (cosPoint < cfgCosPoint)
-          return 8;
-    */
+    }
+
+    ROOT::Math::XYZVector rPV(col.posX(), col.posY(), col.posZ());
+    ROOT::Math::XYZVector rSV(hq3r.hqx(), hq3r.hqy(), hq3r.hqz());
+
+    ROOT::Math::XYZVector L = rSV - rPV;
+    ROOT::Math::XYZVector exMom(ex.Px(), ex.Py(), ex.Pz());
+
+    const double denom2 = (L.R() * exMom.R() + 1e-9);
+    const double cosPoint = L.Dot(exMom) / denom2;
+    if (cosPoint < cfgCosPoint) {
+      selection += 8;
+    }
+
     return selection;
   }
 
@@ -344,7 +372,7 @@ struct heptaquark {
           HQ12 = HQ1 + HQ2;
           HQ13 = HQ1 + HQ3;
 
-          if (cfgSelectHQ && selectHQ(HQ1, HQ2, HQ3))
+          if (cfgSelectHQ && selectHQ(hqtrackd1, hqtrackd2, hqtrackd3, collision))
             continue;
 
           histos.fill(HIST("h_InvMass_same"), exotic.M(), exotic.Pt(), collision.centrality());