[Common] ITS PID: Use truncated mean to suppress dedx tails (AliceO2Group#9334)

Co-authored-by: Francesco Mazzaschi <[email protected]>

Author: fmazzasc

Common/DataModel/PIDResponseITS.h

@@ -34,19 +34,24 @@ namespace o2::aod
 struct ITSResponse {
   static float averageClusterSize(uint32_t itsClusterSizes)
   {
-    float average = 0;
+    float sum = 0;
     int nclusters = 0;
-
+    int max = 0;
     for (int layer = 0; layer < 7; layer++) {
-      if ((itsClusterSizes >> (layer * 4)) & 0xf) {
+      int clsize = (itsClusterSizes >> (layer * 4)) & 0xf;
+      if (clsize > 0) {
         nclusters++;
-        average += (itsClusterSizes >> (layer * 4)) & 0xf;
+        sum += clsize;
+        if (clsize > max) {
+          max = clsize;
+        }
       }
     }
     if (nclusters == 0) {
       return 0;
     }
-    return average / nclusters;
+    // truncated mean
+    return (sum - max) / (nclusters - 1);
   };
 
   template <o2::track::PID::ID id>
@@ -67,8 +72,10 @@ struct ITSResponse {
     static constexpr float inverseMass = 1. / o2::track::pid_constants::sMasses[id];
     // static constexpr float charge = static_cast<float>(o2::track::pid_constants::sCharges[id]);
     const float bg = momentum * inverseMass;
-    float relRes = mResolutionParams[0] * std::erf((bg - mResolutionParams[1]) / mResolutionParams[2]);
-    return relRes;
+    if (id == o2::track::PID::Helium3 || id == o2::track::PID::Alpha) {
+      return mResolutionParamsZ2[0] * std::erf((bg - mResolutionParamsZ2[1]) / mResolutionParamsZ2[2]);
+    }
+    return mResolutionParams[0] * std::erf((bg - mResolutionParams[1]) / mResolutionParams[2]);
   }
 
   template <o2::track::PID::ID id>
@@ -87,7 +94,10 @@ struct ITSResponse {
     return nSigmaITS<id>(track.itsClusterSizes(), track.p(), track.eta());
   }
 
-  static void setParameters(float p0, float p1, float p2, float p0_Z2, float p1_Z2, float p2_Z2, float p0_res, float p1_res, float p2_res)
+  static void setParameters(float p0, float p1, float p2,
+                            float p0_Z2, float p1_Z2, float p2_Z2,
+                            float p0_res, float p1_res, float p2_res,
+                            float p0_res_Z2, float p1_res_Z2, float p2_res_Z2)
   {
     if (mIsInitialized) {
       LOG(fatal) << "ITSResponse parameters already initialized";
@@ -102,19 +112,24 @@ struct ITSResponse {
     mResolutionParams[0] = p0_res;
     mResolutionParams[1] = p1_res;
     mResolutionParams[2] = p2_res;
+    mResolutionParamsZ2[0] = p0_res_Z2;
+    mResolutionParamsZ2[1] = p1_res_Z2;
+    mResolutionParamsZ2[2] = p2_res_Z2;
   }
 
  private:
   static std::array<float, 3> mITSRespParams;
   static std::array<float, 3> mITSRespParamsZ2;
   static std::array<float, 3> mResolutionParams;
+  static std::array<float, 3> mResolutionParamsZ2;
   static bool mIsInitialized;
 };
 
-std::array<float, 3> ITSResponse::mITSRespParams = {1.1576, 1.684, 1.9453};
-std::array<float, 3> ITSResponse::mITSRespParamsZ2 = {2.8752, 1.1246, 5.0259};
+std::array<float, 3> ITSResponse::mITSRespParams = {1.18941, 1.53792, 1.69961};
+std::array<float, 3> ITSResponse::mITSRespParamsZ2 = {2.35117, 1.80347, 5.14355};
 // relative resolution is modelled with an erf function: [0]*TMath::Erf((x-[1])/[2])
-std::array<float, 3> ITSResponse::mResolutionParams = {0.2431, -0.3293, 1.533};
+std::array<float, 3> ITSResponse::mResolutionParams = {1.94669e-01, -2.08616e-01, 1.30753};
+std::array<float, 3> ITSResponse::mResolutionParamsZ2 = {8.74371e-02, -1.82804, 5.06449e-01};
 bool ITSResponse::mIsInitialized = false;
 
 namespace pidits

Common/TableProducer/PID/pidITS.cxx

@@ -43,17 +43,16 @@ using namespace o2::track;
 MetadataHelper metadataInfo;
 
 static constexpr int nCases = 2;
-static constexpr int nParameters = 9;
+static constexpr int nParameters = 12;
 static const std::vector<std::string> casesNames{"Data", "MC"};
 static const std::vector<std::string> parameterNames{"RespITSPar1", "RespITSPar2", "RespITSPar3",
                                                      "RespITSPar1_Z2", "RespITSPar2_Z2", "RespITSPar3_Z2",
-                                                     "ResolutionPar1", "ResolutionPar2", "ResolutionPar3"};
-static constexpr float defaultParameters[nCases][nParameters]{{0.903, 2.014, 2.440,
-                                                               2.8752, 1.1246, 5.0259,
-                                                               0.2431, -0.3293, 1.533},
-                                                              {0.903, 2.014, 2.440,
-                                                               2.8752, 1.1246, 5.0259,
-                                                               0.2431, -0.3293, 1.533}};
+                                                     "ResolutionPar1", "ResolutionPar2", "ResolutionPar3",
+                                                     "ResolutionPar1_Z2", "ResolutionPar2_Z2", "ResolutionPar3_Z2"};
+
+static constexpr float defaultParameters[nCases][nParameters] = {
+  {1.18941, 1.53792, 1.69961, 2.35117, 1.80347, 5.14355, 1.94669e-01, -2.08616e-01, 1.30753, 8.74371e-02, -1.82804, 5.06449e-01},
+  {1.18941, 1.53792, 1.69961, 2.35117, 1.80347, 5.14355, 1.94669e-01, -2.08616e-01, 1.30753, 8.74371e-02, -1.82804, 5.06449e-01}};
 
 /// Task to produce the ITS PID information for each particle species
 /// The parametrization is: [p0/(bg)**p1 + p2] being bg = p/m. Different parametrizations are used for He3 and Alpha particles.
@@ -90,7 +89,10 @@ struct itsPid {
                                           itsParams->get(key, "RespITSPar3_Z2"),
                                           itsParams->get(key, "ResolutionPar1"),
                                           itsParams->get(key, "ResolutionPar2"),
-                                          itsParams->get(key, "ResolutionPar3"));
+                                          itsParams->get(key, "ResolutionPar3"),
+                                          itsParams->get(key, "ResolutionPar1_Z2"),
+                                          itsParams->get(key, "ResolutionPar2_Z2"),
+                                          itsParams->get(key, "ResolutionPar3_Z2"));
     }
   }