Merge pull request #775 in B2/basf2 from bugfix/BIIDP-5236-pxd-hot-de…

…ad-channel-calibration-bugged-in-proc13_chunk1_rel06-00 to release/06-00

* commit 'b7a866d815314a770bf58401abfaba713142bb1d':
  Check PXDHits pointer before processing.
  Check PXDClusterCounter pointer before processing.
  Added empty tag for beam energy in prompt calibration.
  Bugfix for PXD efficiency estimation and improving plotting
  Bugfix for PXD calibration plotting and adding support for different track particles in gain calibration.
  Style fix

calibration/scripts/prompt/__init__.py

@@ -21,6 +21,7 @@
                       "Beam Energy": {"No Beam": "No Beam",
                                       "4S": "4S",
                                       "Continuum": "Continuum",
+                                      "": "",
                                       "Scan": "Scan"},
                       "Run Type": {"beam": "beam",
                                    "cosmic": "cosmic",

calibration/scripts/prompt/calibrations/caf_pxd.py

@@ -37,6 +37,7 @@
                                        INPUT_DATA_FILTERS["Beam Energy"]["4S"],
                                        INPUT_DATA_FILTERS["Beam Energy"]["Continuum"],
                                        INPUT_DATA_FILTERS["Beam Energy"]["Scan"],
+                                       INPUT_DATA_FILTERS["Beam Energy"][""],
                                        INPUT_DATA_FILTERS["Run Type"]["physics"],
                                        INPUT_DATA_FILTERS["Data Quality Tag"]["Good Or Recoverable"]],
                                    "cosmic": [INPUT_DATA_FILTERS["Run Type"]["cosmic"]]},

calibration/scripts/prompt/calibrations/caf_pxd_gain.py

@@ -34,6 +34,7 @@
                                        INPUT_DATA_FILTERS["Beam Energy"]["4S"],
                                        INPUT_DATA_FILTERS["Beam Energy"]["Continuum"],
                                        INPUT_DATA_FILTERS["Beam Energy"]["Scan"],
+                                       INPUT_DATA_FILTERS["Beam Energy"][""],
                                        INPUT_DATA_FILTERS["Run Type"]["physics"],
                                        INPUT_DATA_FILTERS["Data Quality Tag"]["Good"]]},
                                expert_config={

pxd/calibration/src/PXDAnalyticGainCalibrationAlgorithm.cc

@@ -54,6 +54,17 @@ CalibrationAlgorithm::EResult PXDAnalyticGainCalibrationAlgorithm::calibrate()
   // Get counter histograms
   auto cluster_counter = getObjectPtr<TH1I>("PXDClusterCounter");
 
+  // Check if there is any PXD cluster
+  if (cluster_counter == nullptr) {
+    B2WARNING("No PXD cluster reconstructed!");
+    if (not forceContinue)
+      return c_NotEnoughData;
+    else {
+      B2WARNING("Skip processing.");
+      return c_OK;
+    }
+  }
+
   // Extract number of sensors from counter histograms
   auto nSensors = getNumberOfSensors(cluster_counter);
 

pxd/calibration/src/PXDHotPixelMaskCalibrationAlgorithm.cc

@@ -96,6 +96,14 @@ CalibrationAlgorithm::EResult PXDHotPixelMaskCalibrationAlgorithm::calibrate()
   auto collector_pxdhits = getObjectPtr<TH1I>("PXDHits");
   auto collector_pxdhitcounts = getObjectPtr<TH1I>("PXDHitCounts");
 
+  // Check if there is any PXD hit
+  if (!collector_pxdhits) {
+    if (forceContinueMasking)
+      return c_OK;
+    else
+      return c_NotEnoughData;
+  }
+
   // We should have some minimum number of events
   auto nevents = collector_pxdhits->GetEntries();
   if (nevents < minEvents) {

pxd/modules/pxdPerformanceVariablesCollector/include/PXDPerformanceVariablesCollectorModule.h

@@ -166,7 +166,7 @@ namespace Belle2 {
     PXDGainMapPar m_gainMap;
 
     /** Track struct for holding required variables */
-    PXD::Track_t track_struct;
+    PXD::Track_t m_track_struct;
     /** Name of the particle list for gain calibration */
     std::string m_PList4GainName = "";
     /** Name of the particle list for efficiency study */

pxd/modules/pxdPerformanceVariablesCollector/src/PXDPerformanceVariablesCollectorModule.cc

@@ -216,6 +216,7 @@ void PXDPerformanceVariablesCollectorModule::collect() // Do your event() stuff
   StoreObjPtr<ParticleList> particles4Eff(m_PList4EffName);
   if (particles4Eff.isValid() && particles4Eff->getListSize() == 1)
     m_selected4Eff = true;
+  const Particle* vpho4eff = particles4Eff->getParticle(0);
 
   collectDeltaIP(); // using ParticleList(m_PList4ResName)
 
@@ -228,14 +229,15 @@ void PXDPerformanceVariablesCollectorModule::collect() // Do your event() stuff
 
   for (auto const& particle : *particles) {
     const Track* trackPtr = particle.getTrack();
+    auto mass = particle.getMass();
     if (!trackPtr) return;
     auto recoTrackPtr = trackPtr->getRelated<RecoTrack>("");
     if (!recoTrackPtr) return;
     auto pxdIntercepts = recoTrackPtr->getRelationsTo<PXDIntercept>("");
     for (auto const& pxdIntercept : pxdIntercepts) {
       TrackCluster_t trackCluster;
       // Function setValues() also returns a recoTrack pointer
-      if (!trackCluster.setValues(pxdIntercept, "", "PXDClustersFromTracks"))
+      if (!trackCluster.setValues(pxdIntercept, "", "PXDClustersFromTracks", mass))
         continue;
 
       auto const& cluster = trackCluster.cluster;
@@ -244,8 +246,13 @@ void PXDPerformanceVariablesCollectorModule::collect() // Do your event() stuff
 
 
       // Collect info for efficiency study
-      if (m_selected4Eff && intersection.inside)
-        collectEfficiencyVariables(trackCluster);
+      if (m_selected4Eff && intersection.inside) {
+        // check if the particle is vpho4eff's daughter
+        for (auto const& daughter : vpho4eff->getDaughters()) {
+          if (particle.isCopyOf(daughter))
+            collectEfficiencyVariables(trackCluster);
+        }  // end of vpho4eff daughter loop
+      }
 
       // Collect info for gain calibration
       // Check for valid cluster and intersection
@@ -284,12 +291,12 @@ void PXDPerformanceVariablesCollectorModule::collectDeltaIP()
   const TrackFitResult* tr0 = part0->getTrack()->getTrackFitResultWithClosestMass(Const::pion);
   const TrackFitResult* tr1 = part1->getTrack()->getTrackFitResultWithClosestMass(Const::pion);
 
-  track_struct.setTrackVariables(tr0, vertex);
-  auto d0p_0 = track_struct.d0p;
-  auto z0p_0 = track_struct.z0p;
-  track_struct.setTrackVariables(tr1, vertex);
-  auto d0p_1 = track_struct.d0p;
-  auto z0p_1 = track_struct.z0p;
+  m_track_struct.setTrackVariables(tr0, vertex);
+  auto d0p_0 = m_track_struct.d0p;
+  auto z0p_0 = m_track_struct.z0p;
+  m_track_struct.setTrackVariables(tr1, vertex);
+  auto d0p_1 = m_track_struct.d0p;
+  auto z0p_1 = m_track_struct.z0p;
 
   m_deltaD0oSqrt2 = (d0p_0 + d0p_1) / sqrt(2.);
   m_deltaZ0oSqrt2 = (z0p_0 - z0p_1) / sqrt(2.);

pxd/scripts/pxd/calibration/__init__.py

@@ -176,11 +176,15 @@ def gain_calibration(input_files, cal_name="PXDGainCalibration",
       **kwargs: Additional configuration to support extentions without changing scripts in calibration folder.
         Supported options are listed below:
 
-        "collector_prefix" is a string indicating which collector to be used for gain calibration. The supported
+        "collector_prefix": a string indicating which collector to be used for gain calibration. The supported
           collectors are:
             PXDPerformanceVariablesCollector (default),
             PXDPerformanceCollector(using RAVE package for vertexing, obsolete)
         "useClusterPosition": Flag to use cluster postion rather than track point to group pixels for calibration.
+        "particle_type": Particle type assigned to tracks. "e" by default.
+        "track_cuts_4gain": Track cuts used for gain calibration.
+        "track_cuts_4eff": Track cuts used for efficiency study.
+        "track_cuts_4res": Track cuts used for resolution study.
 
     Return:
       A caf.framework.Calibration obj.
@@ -201,6 +205,10 @@ def gain_calibration(input_files, cal_name="PXDGainCalibration",
     useClusterPosition = kwargs.get("useClusterPosition", False)
     if not isinstance(useClusterPosition, bool):
         raise ValueError("useClusterPosition has to be a boolean!")
+    particle_type = kwargs.get("particle_type", "e")  # rely on modular analysis for value check
+    track_cuts_4gain = kwargs.get("track_cuts_4gain", "p > 1.0")  # see above
+    track_cuts_4eff = kwargs.get("track_cuts_4eff", "pt > 2.0")  # see above
+    track_cuts_4res = kwargs.get("track_cuts_4res", "Note2019")  # NOTE-TE-2019-018
 
     # Create basf2 path
 
@@ -226,35 +234,37 @@ def gain_calibration(input_files, cal_name="PXDGainCalibration",
     else:  # the default PXDPerformanceVariablesCollector
         import modularAnalysis as ana
         import vertex
-        from variables import variables as vm
         # Particle list for gain calibration
-        ana.fillParticleList('e+:gain', "p > 1.0", path=main)
+        p = particle_type
+        ana.fillParticleList(f'{p}+:gain', track_cuts_4gain, path=main)
 
         # Particle list for event selection in efficiency monitoring
         # nSVDHits > 5 doesn't help, firstSVDLayer == 3 for < 0.1% improvement?
-        ana.fillParticleList('e+:eff', "pt > 2.0", path=main)
+        ana.fillParticleList(f'{p}+:eff', track_cuts_4eff, path=main)
         # Mass cut (9.5 < M < 11.5) below can improve efficiency by ~ 1%
-        ana.reconstructDecay('vpho:eff -> e+:eff e-:eff', '9.5<M<11.5', path=main)
+        ana.reconstructDecay(f'vpho:eff -> {p}+:eff {p}-:eff', '9.5<M<11.5', path=main)
         # < 0.1% improvement by using kfit pvalue >= 0.01 after mass cut
         # vertex.kFit('vpho:eff', conf_level=0.01, fit_type="fourC", daughtersUpdate=False, path=main)
 
         # Particle list for studying impact parameter resolution
         # Alias dosn't work with airflow implementation
         # vm.addAlias("pBetaSinTheta3o2", "formula(pt * (1./(1. + tanLambda**2)**0.5)**0.5)")
         # vm.addAlias("absLambda", "abs(atan(tanLambda))")
-        mySelection = 'pt>1.0 and abs(dz)<1.0 and dr<0.3'
-        mySelection += ' and nCDCHits>20 and nSVDHits>=8 and nPXDHits>=1'
-        mySelection += ' and [abs(atan(tanLambda)) < 0.5]'
-        mySelection += ' and [formula(pt * (1./(1. + tanLambda**2)**0.5)**0.5) > 2.0]'
-        # mySelection += ' and [absLambda<0.5]'
-        # mySelection += ' and [pBetaSinTheta3o2>2.0]'
-        ana.fillParticleList('e+:res', mySelection, path=main)
-        ana.reconstructDecay('vpho:res -> e+:res e-:res', '9.5<M<11.5', path=main)
+        track_cuts_4res_note2019 = 'pt>1.0 and abs(dz)<1.0 and dr<0.3'
+        track_cuts_4res_note2019 += ' and nCDCHits>20 and nSVDHits>=8 and nPXDHits>=1'
+        track_cuts_4res_note2019 += ' and [abs(atan(tanLambda)) < 0.5]'
+        track_cuts_4res_note2019 += ' and [formula(pt * (1./(1. + tanLambda**2)**0.5)**0.5) > 2.0]'
+        # track_cuts_4res_note2019 += ' and [absLambda<0.5]'
+        # track_cuts_4res_note2019 += ' and [pBetaSinTheta3o2>2.0]'
+        if track_cuts_4res == "Note2019":
+            track_cuts_4res = track_cuts_4res_note2019
+        ana.fillParticleList(f'{p}+:res', track_cuts_4res, path=main)
+        ana.reconstructDecay(f'vpho:res -> {p}+:res {p}-:res', '9.5<M<11.5', path=main)
         # Remove multiple candidate events
         ana.applyCuts('vpho:res', 'nParticlesInList(vpho:res)==1', path=main)
         vertex.kFit('vpho:res', conf_level=0.0, path=main)
 
-        collector.param("PList4GainName", "e+:gain")
+        collector.param("PList4GainName", f"{p}+:gain")
         collector.param("PList4EffName", "vpho:eff")
         collector.param("PList4ResName", "vpho:res")
         collector.param("maskedDistance", 3)

pxd/scripts/pxd/calibration/calibration_checker.py

@@ -221,7 +221,7 @@ def fill_graphs(self):
             self.sum_graphs["hot_dead"].append(current_run, sum_hotdeadfraction / nSensors)
 
         nSensors = nSensors - len(dead_checker.dbObj.getDeadSensorMap())
-        if current_run == saved_run_occ:
+        if current_run == saved_run_occ and nSensors > 0:
             self.sum_graphs["occ_masked"].append(current_run, sum_occupancymasked / nSensors)
             self.sum_graphs["occ_no_mask"].append(current_run, sum_occupancynomask / nSensors)
 

pxd/scripts/pxd/utils/plots.py

@@ -114,6 +114,7 @@ def df_plot_errorbar(df, x, y, yerr_low, yerr_up, ax=None, *args, **kwargs):
     ax.legend()
     ax.set_xlabel(x)
     ax.set_ylabel(y)
+    return ax
 
 
 # Extend pandas.DataFrame
@@ -181,7 +182,7 @@ def plot_efficiency_map(tree, exp_runs=[], num_name="hTrackClustersLayer1", den_
         exp_runs = sorted(exp_runs)
         use_exp_run_tuple = False
     count = 0
-    h_num, h_den = None, None
+    h_eff, h_num, h_den = None, None, None
     for i_evt in range(tree.GetEntries()):
         tree.GetEntry(i_evt)
         exp = getattr(tree, "exp")
@@ -199,14 +200,36 @@ def plot_efficiency_map(tree, exp_runs=[], num_name="hTrackClustersLayer1", den_
             h_den.Add(getattr(tree, den_name))
         count += 1
     if count:
-        h_num.Divide(h_den)
-        h_num.SetTitle(title)
-        h_num.SetStats(False)
-        h_num.Draw("colz")
+        h_temp = h_den.Clone()
+        h_temp.Divide(h_den)  # h_temp bins filled to 1 if there is any counts in h_den
+        h_eff = h_num.Clone()
+        h_eff.Add(h_temp, 1e-9)  # Added 1e-9 which shouldn't bias eff
+        h_eff.Divide(h_den)
+        h_eff.SetTitle(title)
+        h_eff.SetStats(True)
+        # default_opt_stat = ROOT.gStyle.GetOptStat()
+        ROOT.gStyle.SetOptStat(10)
+        h_eff.Draw("colz")
         canvas.Draw()
+        s = h_eff.GetListOfFunctions().FindObject("stats")
+        # print(s.GetX1NDC(), s.GetX2NDC())
+        s.SetX1NDC(0.7)
+        s.SetX2NDC(0.9)
+        s.SetY1NDC(0.9)
+        s.SetY2NDC(0.95)
+        s.SetFillColorAlpha(0, 0.1)
+        s.SetLineColorAlpha(0, 0)
+        canvas.Update()
         if save_to:
+            h_eff.GetZaxis().SetRangeUser(0.9, 1)
+            canvas.Update()
+            canvas.Print(save_to+".above90.png")
+            h_eff.GetZaxis().SetRangeUser(0, 1)
+            canvas.Update()
             canvas.Print(save_to)
-    return h_num
+
+        # ROOT.gStyle.SetOptStat(default_opt_stat)
+    return h_eff, h_num, h_den
 
 
 def plot_in_module_efficiency(df, pxdid=1052, figsize=(12, 16), alpha=0.7, save_to="",
@@ -303,14 +326,16 @@ def plot_module_efficiencies_in_DHHs(df, eff_var="eff_sel", phase="early_phase3"
         }
     for dhh, pxdid_list in dhh_modules_dic.items():
         plt.figure(figsize=figsize)
+        ymin = 1.0
         for pxdid in pxdid_list:
             df.query(f"{eff_var}>0&pxdid=={pxdid}&{eff_var}_err_low<0.01").errorbar(
                 y=eff_var, x="run", yerr_low=f"{eff_var}_err_low", yerr_up=f"{eff_var}_err_up", label=f"{pxdid}", alpha=0.7)
-            plt.title(dhh + " efficiency")
-            ymin, ymax = plt.ylim()
-            plt.ylim(ymin, 1.0)
-            if save_to:
-                plt.savefig(dhh + "_" + save_to)
+            ymin = min(plt.ylim()[0], ymin)
+        plt.ylabel("Efficiency (selected)")
+        plt.title(dhh + " efficiency")
+        plt.ylim(ymin, 1.0)
+        if save_to:
+            plt.savefig(dhh + "_" + save_to)
 
 
 if __name__ == '__main__':

pxd/utilities/include/PXDPerformanceStructs.h

@@ -61,9 +61,11 @@ namespace Belle2 {
       /** Update values from a PXDCluster.
        * @param pxdIntercept a PXDIntercept object.
        * @param recoTracksName Name of RecoTrack collection
+       * @param mass Mass of the impinging particle
        * @return the pointer of the related RecoTrack object.
        */
-      RecoTrack* setValues(const PXDIntercept& pxdIntercept, const std::string& recoTracksName = "");
+      RecoTrack* setValues(const PXDIntercept& pxdIntercept, const std::string& recoTracksName = "",
+                           const double& mass = Const::electronMass);
 
       float x;         /**< Global position in x. */
       float y;         /**< Global position in y. */
@@ -86,11 +88,13 @@ namespace Belle2 {
        * @param pxdIntercept a PXDIntercept object.
        * @param recoTracksName Name of RecoTrack collection
        * @param pxdTrackClustersName Name of track matched PXDClusters
+       * @param mass Mass of the impinging particle
        * @return the pointer of the related RecoTrack object.
        */
       RecoTrack* setValues(const PXDIntercept& pxdIntercept,
                            const std::string& recoTracksName = "",
-                           const std::string& pxdTrackClustersName = "PXDClustersFromTracks");
+                           const std::string& pxdTrackClustersName = "PXDClustersFromTracks",
+                           const double& mass = Const::electronMass);
 
       bool usedInTrack;        /**< True if the cluster is used in tracking */
       float dU;                /**< Residual (meas - prediction) in U. */

pxd/utilities/src/PXDPerformanceStructs.cc

@@ -25,7 +25,7 @@ namespace Belle2 {
       posV = pxdCluster.getV();
     }
 
-    RecoTrack* TrackPoint_t::setValues(const PXDIntercept& pxdIntercept, const std::string& recoTracksName)
+    RecoTrack* TrackPoint_t::setValues(const PXDIntercept& pxdIntercept, const std::string& recoTracksName, const double& mass)
     {
       // Construct VxdID from its baseType (unsigned short)
       VxdID sensorID(pxdIntercept.getSensorID());
@@ -60,17 +60,18 @@ namespace Belle2 {
       auto ADUToEnergy = PXD::PXDGainCalibrator::getInstance().getADUToEnergy(sensorID,
                          sensorInfo.getUCellID(pxdIntercept.getCoorU()),
                          sensorInfo.getVCellID(pxdIntercept.getCoorV()));
-      chargeMPV = getDeltaP(intersec_p.Mag(), length) / ADUToEnergy;
+      chargeMPV = getDeltaP(intersec_p.Mag(), length, mass) / ADUToEnergy;
 
       // Return pointer of the relatd RecoTrack for accessing additional info.
       return recoTracks[0];
     }
 
     RecoTrack* TrackCluster_t::setValues(const PXDIntercept& pxdIntercept,
                                          const std::string& recoTracksName,
-                                         const std::string& pxdTrackClustersName)
+                                         const std::string& pxdTrackClustersName,
+                                         const double& mass)
     {
-      auto recoTrackPtr = intersection.setValues(pxdIntercept, recoTracksName);
+      auto recoTrackPtr = intersection.setValues(pxdIntercept, recoTracksName, mass);
       // sensor ID from intersectioon
       //VxdID sensorID(pxdIntercept.getSensorID());
       if (!recoTrackPtr) return nullptr; // return nullptr