From bddbaeac79887666ea5568ff20401d2026e74e08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 18:22:00 +0200 Subject: [PATCH 1/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 529 +++++++++++-------- 1 file changed, 301 insertions(+), 228 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index e3c740966d1..c61db969688 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -30,41 +30,45 @@ /// \since May 22, 2024 /// -#include -#include -#include -#include -#include +#include "TableHelper.h" -#include +#include "ALICE3/Core/DelphesO2TrackSmearer.h" +#include "ALICE3/Core/TrackUtilities.h" +#include "ALICE3/DataModel/OTFRICH.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CcdbApi.h" +#include "CommonConstants/GeomConstants.h" +#include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "CommonUtils/NameConf.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsVertexing/HelixHelper.h" +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "Framework/RunningWorkflowInfo.h" #include "Framework/HistogramRegistry.h" #include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/ASoAHelpers.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/trackUtilities.h" -#include "ALICE3/Core/TrackUtilities.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/DCA.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "CommonUtils/NameConf.h" -#include "CCDB/CcdbApi.h" -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "CommonConstants/GeomConstants.h" -#include "CommonConstants/PhysicsConstants.h" +#include "ReconstructionDataFormats/PID.h" + #include "TRandom3.h" -#include "TVector3.h" #include "TString.h" -#include "ALICE3/DataModel/OTFRICH.h" -#include "DetectorsVertexing/HelixHelper.h" -#include "TableHelper.h" -#include "ALICE3/Core/DelphesO2TrackSmearer.h" +#include "TVector3.h" +#include + +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -72,12 +76,13 @@ using namespace o2::constants::math; struct OnTheFlyRichPid { Produces upgradeRich; + Produces upgradeRichSignal; // necessary for particle charges Service pdg; // master setting: magnetic field - Configurable dBz{"dBz", 20, "magnetic field (kilogauss)"}; + Configurable magneticField{"magneticField", 0, "magnetic field (kilogauss) if 0, taken from the tracker task"}; // add rich-specific configurables here Configurable bRichNumberOfSectors{"bRichNumberOfSectors", 21, "barrel RICH number of sectors"}; @@ -230,19 +235,19 @@ struct OnTheFlyRichPid { l_aerogel_z[i_central_mirror] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_barrel_cylinder / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_barrel_cylinder); T_r_plus_g[i_central_mirror] = R_max - R_min; float t = std::tan(std::atan(m_val) + std::atan(square_size_barrel_cylinder / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_barrel_cylinder * m_val))); - theta_max[i_central_mirror] = M_PI / 2.0 - std::atan(t); - theta_min[i_central_mirror] = M_PI / 2.0 + std::atan(t); + theta_max[i_central_mirror] = o2::constants::math::PIHalf - std::atan(t); + theta_min[i_central_mirror] = o2::constants::math::PIHalf + std::atan(t); mProjectiveLengthInner = R_min * t; aerogel_rindex[i_central_mirror] = bRichRefractiveIndexSector[0]; for (int i = i_central_mirror + 1; i < number_of_sectors_in_z; i++) { float par_a = t; float par_b = 2.0 * R_max / square_size_z; m_val = (std::sqrt(par_a * par_a * par_b * par_b + par_b * par_b - 1.0) + par_a * par_b * par_b) / (par_b * par_b - 1.0); - theta_min[i] = M_PI / 2.0 - std::atan(t); - theta_max[2 * i_central_mirror - i] = M_PI / 2.0 + std::atan(t); + theta_min[i] = o2::constants::math::PIHalf - std::atan(t); + theta_max[2 * i_central_mirror - i] = o2::constants::math::PIHalf + std::atan(t); t = std::tan(std::atan(m_val) + std::atan(square_size_z / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_z * m_val))); - theta_max[i] = M_PI / 2.0 - std::atan(t); - theta_min[2 * i_central_mirror - i] = M_PI / 2.0 + std::atan(t); + theta_max[i] = o2::constants::math::PIHalf - std::atan(t); + theta_min[2 * i_central_mirror - i] = o2::constants::math::PIHalf + std::atan(t); // Forward sectors theta_bi[i] = std::atan(m_val); R0_tilt[i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); @@ -271,11 +276,11 @@ struct OnTheFlyRichPid { float par_a = t; float par_b = 2.0 * R_max / square_size_z; m_val = (std::sqrt(par_a * par_a * par_b * par_b + par_b * par_b - 1.0) + par_a * par_b * par_b) / (par_b * par_b - 1.0); - theta_min[i] = M_PI / 2.0 - std::atan(t); - theta_max[2 * i_central_mirror - i - 1] = M_PI / 2.0 + std::atan(t); + theta_min[i] = o2::constants::math::PIHalf - std::atan(t); + theta_max[2 * i_central_mirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); t = std::tan(std::atan(m_val) + std::atan(square_size_z / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_z * m_val))); - theta_max[i] = M_PI / 2.0 - std::atan(t); - theta_min[2 * i_central_mirror - i - 1] = M_PI / 2.0 + std::atan(t); + theta_max[i] = o2::constants::math::PIHalf - std::atan(t); + theta_min[2 * i_central_mirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); // Forward sectors theta_bi[i] = std::atan(m_val); R0_tilt[i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); @@ -321,6 +326,14 @@ struct OnTheFlyRichPid { { pRandomNumberGenerator.SetSeed(0); // fully randomize + if (magneticField.value < 0.0001f) { + LOG(info) << "Getting the magnetic field from the on-the-fly tracker task"; + if (!getTaskOptionValue(initContext, "on-the-fly-tracker", magneticField, false)) { + LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; + } + LOG(info) << "Bz = " << magneticField.value << " T"; + } + // Check if inheriting the LUT configuration auto configLutPath = [&](Configurable& lut) { if (lut.value != "inherit") { @@ -376,9 +389,10 @@ struct OnTheFlyRichPid { histos.add("h2dAngularResolutionVsEtaBarrelRICH", "h2dAngularResolutionVsEtaBarrelRICH", kTH2F, {axisEta, axisRingAngularResolution}); histos.add("hSectorID", "hSectorID", kTH1F, {axisSector}); - std::string particle_names1[5] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; - std::string particle_names2[5] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; - for (int i_true = 0; i_true < 5; i_true++) { + const int kNspec = 5; // electron, muon, pion, kaon, proton + std::string particle_names1[kNspec] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; + std::string particle_names2[kNspec] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; + for (int i_true = 0; i_true < kNspec; i_true++) { std::string name_title_barrel_track_res = "h2dBarrelAngularResTrack" + particle_names2[i_true] + "VsP"; std::string name_title_barrel_total_res = "h2dBarrelAngularResTotal" + particle_names2[i_true] + "VsP"; const AxisSpec axisTrackAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Track angular resolution - " + particle_names1[i_true] + " (mrad)"}; @@ -386,8 +400,8 @@ struct OnTheFlyRichPid { histos.add(name_title_barrel_track_res.c_str(), name_title_barrel_track_res.c_str(), kTH2F, {axisMomentum, axisTrackAngularRes}); histos.add(name_title_barrel_total_res.c_str(), name_title_barrel_total_res.c_str(), kTH2F, {axisMomentum, axisTotalAngularRes}); } - for (int i_true = 0; i_true < 5; i_true++) { - for (int i_hyp = 0; i_hyp < 5; i_hyp++) { + for (int i_true = 0; i_true < kNspec; i_true++) { + for (int i_hyp = 0; i_hyp < kNspec; i_hyp++) { std::string name_title = "h2dBarrelNsigmaTrue" + particle_names2[i_true] + "Vs" + particle_names2[i_hyp] + "Hypothesis"; if (i_true == i_hyp) { const AxisSpec axisNsigmaCorrect{static_cast(nBinsNsigmaCorrectSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"}; @@ -451,7 +465,7 @@ struct OnTheFlyRichPid { /// \param track the input track /// \param radius the radius of the layer you're calculating the length to /// \param magneticField the magnetic field to use when propagating - bool checkMagfieldLimit(o2::track::TrackParCov track, float radius, float magneticField) + bool checkMagfieldLimit(o2::track::TrackParCov track, const float radius, const float magneticField) { o2::math_utils::CircleXYf_t trcCircle; float sna, csa; @@ -470,7 +484,7 @@ struct OnTheFlyRichPid { /// returns sector hit by the track (if any), -1 otherwise /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) - int findSector(float eta) + int findSector(const float eta) { float polar = 2.0 * std::atan(std::exp(-eta)); for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) { @@ -484,13 +498,15 @@ struct OnTheFlyRichPid { /// returns radiator radius in cm at considered eta (accounting for sector inclination) /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) /// \param i_sector the index of the track RICH sector - float radiusRipple(float eta, int i_sector) + float radiusRipple(const float eta, const int i_sector) { - float polar = 2.0 * std::atan(std::exp(-eta)); if (i_sector > 0) { - float R_sec_rich = rad_centers[i_sector].X(); - float z_sec_rich = rad_centers[i_sector].Z(); - return (std::pow(R_sec_rich, 2) + std::pow(z_sec_rich, 2)) / (R_sec_rich + z_sec_rich / std::tan(polar)); + const float polar = 2.0 * std::atan(std::exp(-eta)); + const float R_sec_rich = rad_centers[i_sector].X(); + const float z_sec_rich = rad_centers[i_sector].Z(); + const float R_sec_rich_squared = R_sec_rich * R_sec_rich; + const float z_sec_rich_squared = z_sec_rich * z_sec_rich; + return (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar)); } else { return error_value; } @@ -500,25 +516,24 @@ struct OnTheFlyRichPid { /// \param momentum the momentum of the tarck /// \param mass the mass of the particle /// \param n the refractive index of the considered sector - float CherenkovAngle(float momentum, float mass, float n) + /// \param angle the output angle (passed by reference) + /// \return true if particle is above the threshold and enough photons, false otherwise + bool CherenkovAngle(const float momentum, const float mass, const float n, float& angle) { - // Check if particle is above the threshold - if (momentum > mass / std::sqrt(n * n - 1.0)) { + if (momentum > mass / std::sqrt(n * n - 1.0)) { // Check if particle is above the threshold // Calculate angle - float angle = std::acos(std::sqrt(momentum * momentum + mass * mass) / (momentum * n)); + angle = std::acos(std::sqrt(momentum * momentum + mass * mass) / (momentum * n)); // Mean number of detected photons (SiPM PDE is accountd in multiplicative factor!!!) - float meanNumberofDetectedPhotons = 230. * std::sin(angle) * std::sin(angle) * bRichRadiatorThickness; + const float meanNumberofDetectedPhotons = 230. * std::sin(angle) * std::sin(angle) * bRichRadiatorThickness; // Require at least 3 photons on average for real angle reconstruction if (meanNumberofDetectedPhotons > 3.) { - return angle; - } else { - return error_value; + return true; // Particle is above the threshold and enough photons } - } else { - return error_value; + return false; // Particle is above the threshold, but not enough photons } + return false; // Particle is below the threshold } /// returns linear interpolation @@ -538,13 +553,13 @@ struct OnTheFlyRichPid { float AngularResolution(float eta) { // Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS) - float eta_sampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; - float res_ring_sampling_with_abs_walls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; - float res_ring_sampling_without_abs_walls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; - int size_res_vector = sizeof(eta_sampling) / sizeof(eta_sampling[0]); + static constexpr float eta_sampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; + static constexpr float res_ring_sampling_with_abs_walls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; + static constexpr float res_ring_sampling_without_abs_walls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; + static constexpr int size_res_vector = sizeof(eta_sampling) / sizeof(eta_sampling[0]); // Use binary search to find the lower and upper indices - int lowerIndex = std::lower_bound(eta_sampling, eta_sampling + size_res_vector, eta) - eta_sampling - 1; - int upperIndex = lowerIndex + 1; + const int lowerIndex = std::lower_bound(eta_sampling, eta_sampling + size_res_vector, eta) - eta_sampling - 1; + const int upperIndex = lowerIndex + 1; if (lowerIndex >= 0 && upperIndex < size_res_vector) { // Resolution interpolation if (bRichFlagAbsorbingWalls) { @@ -569,7 +584,7 @@ struct OnTheFlyRichPid { /// \param radius the nominal radius of the Cherenkov ring float fractionPhotonsProjectiveRICH(float eta, float tile_z_length, float radius) { - float polar = 2.0 * std::atan(std::exp(-eta)); + const float polar = 2.0 * std::atan(std::exp(-eta)); int i_sector = 0; bool flag_sector = false; for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) { @@ -585,14 +600,16 @@ struct OnTheFlyRichPid { float z_sec_rich = rad_centers[i_sector].Z(); // float R_sec_tof = det_centers[i_sector].X(); // float z_sec_tof = det_centers[i_sector].Z(); - float radius_ripple = (std::pow(R_sec_rich, 2) + std::pow(z_sec_rich, 2)) / (R_sec_rich + z_sec_rich / std::tan(polar)); - float z_ripple = radius_ripple / std::tan(polar); - float absZ = std::hypot(radius_ripple - R_sec_rich, z_ripple - z_sec_rich); + const float R_sec_rich_squared = R_sec_rich * R_sec_rich; + const float z_sec_rich_squared = z_sec_rich * z_sec_rich; + const float radius_ripple = (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar)); + const float z_ripple = radius_ripple / std::tan(polar); + const float absZ = std::hypot(radius_ripple - R_sec_rich, z_ripple - z_sec_rich); float fraction = 1.; if (tile_z_length / 2. - absZ < radius) { - fraction = fraction - (1. / M_PI) * std::acos((tile_z_length / 2. - absZ) / radius); + fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tile_z_length / 2. - absZ) / radius); if (tile_z_length / 2. + absZ < radius) { - fraction = fraction - (1. / M_PI) * std::acos((tile_z_length / 2. + absZ) / radius); + fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tile_z_length / 2. + absZ) / radius); } } return fraction; @@ -607,73 +624,80 @@ struct OnTheFlyRichPid { /// \param pixel_size the SiPM pixel size in cm /// \param theta_c the Cherenkov angle for the considered track /// \param tile_z_length the Z-length of the photodetector plane of the considered sector - float extract_ring_angular_resolution(float eta, float n, float n_g, float T_r, float T_g, float pixel_size, float theta_c, float tile_z_length) + float extract_ring_angular_resolution(const float eta, const float n, const float n_g, const float T_r, const float T_g, const float pixel_size, const float theta_c, const float tile_z_length) { // Check if input angle is error value if (theta_c <= error_value + 1) return error_value; // Parametrization variables (notation from https://doi.org/10.1016/0168-9002(94)90532-0) - float phi_c = 0.; - float theta_p = 0.; - float sin_theta_c = std::sin(theta_c); - float cos_theta_c = std::cos(theta_c); - float x_p = std::sin(theta_p); - float z_p = std::cos(theta_p); - float sin_phi = std::sin(phi_c); - float cos_phi = std::cos(phi_c); - // float ze = T_r / (2. * z_p); - float az = z_p * cos_theta_c - x_p * sin_theta_c * cos_phi; - float e3z = std::sqrt(az * az + (n_g / n) * (n_g / n) - 1.); - float Z = T_g; - float alpha = e3z / az; - float etac = e3z * n * n; - float k = T_r / (2. * Z); - float m = 1. / (n * n); - float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); + const float phi_c = 0.; + const float theta_p = 0.; + const float sin_theta_c = std::sin(theta_c); + const float cos_theta_c = std::cos(theta_c); + const float x_p = std::sin(theta_p); + const float z_p = std::cos(theta_p); + const float sin_phi = std::sin(phi_c); + const float cos_phi = std::cos(phi_c); + // const float ze = T_r / (2. * z_p); + const float az = z_p * cos_theta_c - x_p * sin_theta_c * cos_phi; + const float e3z = std::sqrt(az * az + (n_g / n) * (n_g / n) - 1.); + const float Z = T_g; + const float alpha = e3z / az; + const float etac = e3z * n * n; + const float k = T_r / (2. * Z); + const float m = 1. / (n * n); + const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); // Derivative d(theta_c)/dx - float temp1 = etac / Z; - float temp2 = alpha * e3z * cos_phi; - float temp3 = x_p * sin_theta_c * sin_phi * sin_phi; - float d_theta_x = temp1 * (temp2 - temp3); + const float temp1 = etac / Z; + const float temp2 = alpha * e3z * cos_phi; + const float temp3 = x_p * sin_theta_c * sin_phi * sin_phi; + const float d_theta_x = temp1 * (temp2 - temp3); // Derivative d(theta_c)/dy - float temp4 = etac * sin_phi / Z; - float temp5 = cos_theta_c - z_p * (1 - m) / az; - float d_theta_y = temp4 * temp5; + const float temp4 = etac * sin_phi / Z; + const float temp5 = cos_theta_c - z_p * (1 - m) / az; + const float d_theta_y = temp4 * temp5; // Derivative d(theta_c)/dze - float temp8 = etac * sin_theta_c; - float temp9 = Z * (1.0 + k * alpha * alpha * alpha * n * n); - float temp10 = alpha * alpha * (1.0 - x_p * x_p * sin_phi * sin_phi) + lambda * x_p * x_p * sin_phi * sin_phi; - float d_theta_ze = temp8 * temp10 / temp9; + const float temp8 = etac * sin_theta_c; + const float temp9 = Z * (1.0 + k * alpha * alpha * alpha * n * n); + const float temp10 = alpha * alpha * (1.0 - x_p * x_p * sin_phi * sin_phi) + lambda * x_p * x_p * sin_phi * sin_phi; + const float d_theta_ze = temp8 * temp10 / temp9; // Derivative d(theta_c)/dn - float d_theta_n = z_p / (n * az * sin_theta_c); + const float d_theta_n = z_p / (n * az * sin_theta_c); // RMS wavelength (using Sellmeier formula with measured aerogel parameters) - float a0 = 0.0616; - float w0 = 56.5; - float n0 = 1.0307; - float w_mean = 436.0; // 450.0;//484.9;//426.0; //450;//485;//450; // nm //450 426 493.5 // 426.0 - float scaling = n / n0; - float temp11 = a0 * w0 * w0 * w_mean; - float temp12 = std::pow(w_mean * w_mean - w0 * w0, 1.5); - float temp13 = std::sqrt(w_mean * w_mean * (1.0 + a0) - w0 * w0); - float d_n_w = temp11 / (temp12 * temp13) * scaling; - float sigma_w = 115.0; + const float a0 = 0.0616; + const float w0 = 56.5; + const float n0 = 1.0307; + const float w_mean = 436.0; // 450.0;//484.9;//426.0; //450;//485;//450; // nm //450 426 493.5 // 426.0 + const float scaling = n / n0; + const float temp11 = a0 * w0 * w0 * w_mean; + const float temp12 = std::pow(w_mean * w_mean - w0 * w0, 1.5); + const float temp13 = std::sqrt(w_mean * w_mean * (1.0 + a0) - w0 * w0); + const float d_n_w = temp11 / (temp12 * temp13) * scaling; + const float sigma_w = 115.0; // RMS x y - float sigma_theta_x = pixel_size / std::sqrt(12.0); - float sigma_theta_y = pixel_size / std::sqrt(12.0); + const float sigma_theta_x = pixel_size / std::sqrt(12.0); + const float sigma_theta_y = pixel_size / std::sqrt(12.0); // RMS emission point - float sigma_theta_ze = T_r / (z_p * std::sqrt(12.0)); + const float sigma_theta_ze = T_r / (z_p * std::sqrt(12.0)); // Single photon angular resolution - float res_x = d_theta_x * sigma_theta_x; - float res_y = d_theta_y * sigma_theta_y; - float res_ze = d_theta_ze * sigma_theta_ze; - float res_n = d_theta_n * d_n_w * sigma_w; - float single_photon_angular_resolution = std::sqrt(std::pow(res_x, 2) + std::pow(res_y, 2) + std::pow(res_ze, 2) + std::pow(res_n, 2)); + const float res_x = d_theta_x * sigma_theta_x; + const float res_y = d_theta_y * sigma_theta_y; + const float res_ze = d_theta_ze * sigma_theta_ze; + const float res_n = d_theta_n * d_n_w * sigma_w; + const float res_x_squared = res_x * res_x; + const float res_y_squared = res_y * res_y; + const float res_ze_squared = res_ze * res_ze; + const float res_n_squared = res_n * res_n; + const float single_photon_angular_resolution = std::sqrt(res_x_squared + res_y_squared + res_ze_squared + res_n_squared); // Fraction of photons in rings (to account loss close to sector boundary in projective geometry) - float radius = (T_r / 2.0) * std::tan(theta_c) + T_g * std::tan(std::asin((n / n_g) * std::sin(theta_c))); - float N0 = 24. * T_r / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) - float multiplicity_spectrum_factor = std::pow(std::sin(theta_c), 2.) / std::pow(std::sin(std::acos(1. / 1.03)), 2.); // scale multiplicity w.r.t. N = 1.03 at saturation + const float sin_theta_c_squared = sin_theta_c * sin_theta_c; + const float radius = (T_r / 2.0) * std::tan(theta_c) + T_g * std::tan(std::asin((n / n_g) * sin_theta_c)); + const float N0 = 24. * T_r / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) + const float sin_angle = std::sin(std::acos(1. / 1.03)); + const float sin_angle_squared = sin_angle * sin_angle; + const float multiplicity_spectrum_factor = sin_theta_c_squared / sin_angle_squared; // scale multiplicity w.r.t. N = 1.03 at saturation // Considering average resolution (integrated over the sector) - // float n_photons = (tile_z_length / 2.0 > radius) ? N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(M_PI*tile_z_length)) : N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(M_PI*tile_z_length) - (2.0/(tile_z_length*M_PI))*(-(tile_z_length/(2.0))*std::acos(tile_z_length/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tile_z_length/(2.0*radius),2.0)))); + // float n_photons = (tile_z_length / 2.0 > radius) ? N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length)) : N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length) - (2.0/(tile_z_length*o2::constants::math::PI))*(-(tile_z_length/(2.0))*std::acos(tile_z_length/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tile_z_length/(2.0*radius),2.0)))); // Considering "exact" resolution (eta by eta) float n_photons = N0 * multiplicity_spectrum_factor * fractionPhotonsProjectiveRICH(eta, tile_z_length, radius); if (n_photons <= error_value + 1) @@ -690,25 +714,28 @@ struct OnTheFlyRichPid { /// \param track_pt_resolution the absolute resolution on eta /// \param mass the mass of the particle /// \param refractive_index the refractive index of the radiator - double calculate_track_angular_resolution_advanced(float pt, float eta, float track_pt_resolution, float track_eta_resolution, float mass, float refractive_index) + double calculate_track_angular_resolution_advanced(const float pt, const float eta, const float track_pt_resolution, const float track_eta_resolution, const float mass, const float refractive_index) { // Compute tracking contribution to timing using the error propagation formula // Uses light speed in m/ps, magnetic field in T (*0.1 for conversion kGauss -> T) - double a0 = mass * mass; - double a1 = refractive_index; - double dtheta_on_dpt = a0 / (pt * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)) * std::sqrt(std::pow(pt * std::cosh(eta), 2) * (std::pow(a1, 2) - 1.0) - a0)); - double dtheta_on_deta = (a0 * std::tanh(eta)) / (std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)) * std::sqrt(std::pow(pt * std::cosh(eta), 2) * (std::pow(a1, 2) - 1.0) - a0)); - double track_angular_resolution = std::hypot(std::fabs(dtheta_on_dpt) * track_pt_resolution, std::fabs(dtheta_on_deta) * track_eta_resolution); + const double a0 = mass * mass; + const double a1 = refractive_index; + const double a1_squared = a1 * a1; + const float pt_cosh_eta = pt * std::cosh(eta); + const float pt_cosh_eta_squared = pt_cosh_eta * pt_cosh_eta; + const double dtheta_on_dpt = a0 / (pt * std::sqrt(a0 + pt_cosh_eta_squared) * std::sqrt(pt_cosh_eta_squared * (a1_squared - 1.0) - a0)); + const double dtheta_on_deta = (a0 * std::tanh(eta)) / (std::sqrt(a0 + pt_cosh_eta_squared) * std::sqrt(pt_cosh_eta_squared * (a1_squared - 1.0) - a0)); + const double track_angular_resolution = std::hypot(std::fabs(dtheta_on_dpt) * track_pt_resolution, std::fabs(dtheta_on_deta) * track_eta_resolution); return track_angular_resolution; } - void process(soa::Join::iterator const& collision, soa::Join const& tracks, aod::McParticles const&, aod::McCollisions const&) + void process(soa::Join::iterator const& collision, + soa::Join const& tracks, + aod::McParticles const&, + aod::McCollisions const&) { - // for (int i = 0; i < 1000; i++) - // std::cout << "Prova" << std::endl; - - o2::dataformats::VertexBase pvVtx({collision.posX(), collision.posY(), collision.posZ()}, - {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}); + const o2::dataformats::VertexBase pvVtx({collision.posX(), collision.posY(), collision.posZ()}, + {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}); std::array mcPvCov = {0.}; o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, mcPvCov); @@ -746,44 +773,52 @@ struct OnTheFlyRichPid { for (const auto& track : tracks) { - float nSigmaBarrelRich[5] = {error_value, error_value, error_value, error_value, error_value}; + auto fillDummyValues = [&]() { + upgradeRich(error_value, error_value, error_value, error_value, error_value); + upgradeRichSignal(false, false, false, false, false, false); + }; // first step: find precise arrival time (if any) // --- convert track into perfect track if (!track.has_mcParticle()) { // should always be OK but check please - upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]); + fillDummyValues(); continue; } - auto mcParticle = track.mcParticle(); + const auto& mcParticle = track.mcParticle(); o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg); // float xPv = error_value; - if (o2track.propagateToDCA(mcPvVtx, dBz)) { + if (o2track.propagateToDCA(mcPvVtx, magneticField)) { // xPv = o2track.getX(); } // get particle to calculate Cherenkov angle and resolution auto pdgInfo = pdg->GetParticle(mcParticle.pdgCode()); if (pdgInfo == nullptr) { - upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]); + fillDummyValues(); continue; } // find track bRICH sector int i_sector = findSector(o2track.getEta()); if (i_sector < 0) { - upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]); + fillDummyValues(); continue; } - float expectedAngleBarrelRich = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogel_rindex[i_sector]); + float expectedAngleBarrelRich = error_value; + const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogel_rindex[i_sector], expectedAngleBarrelRich); + if (!expectedAngleBarrelRichOk) { + fillDummyValues(); + continue; // Particle is below the threshold or not enough photons + } // float barrelRICHAngularResolution = AngularResolution(o2track.getEta()); - float barrelRICHAngularResolution = extract_ring_angular_resolution(o2track.getEta(), aerogel_rindex[i_sector], bRichGapRefractiveIndex, bRichRadiatorThickness, gap_thickness[i_sector], bRICHPixelSize, expectedAngleBarrelRich, photodetrctor_length[i_sector]); - float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), i_sector); + const float barrelRICHAngularResolution = extract_ring_angular_resolution(o2track.getEta(), aerogel_rindex[i_sector], bRichGapRefractiveIndex, bRichRadiatorThickness, gap_thickness[i_sector], bRICHPixelSize, expectedAngleBarrelRich, photodetrctor_length[i_sector]); + const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), i_sector); bool flagReachesRadiator = false; if (projectiveRadiatorRadius > error_value + 1.) { - flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, dBz); + flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, magneticField); } /// DISCLAIMER: Exact extrapolation of angular resolution would require track propagation /// to the RICH radiator (accounting sector inclination) in terms of (R,z). @@ -791,62 +826,89 @@ struct OnTheFlyRichPid { /// Discrepancies may be negligible, but would be more rigorous if propagation tool is available // Smear with expected resolutions - float measuredAngleBarrelRich = pRandomNumberGenerator.Gaus(expectedAngleBarrelRich, barrelRICHAngularResolution); + const float measuredAngleBarrelRich = pRandomNumberGenerator.Gaus(expectedAngleBarrelRich, barrelRICHAngularResolution); // Now we calculate the expected Cherenkov angle following certain mass hypotheses // and the (imperfect!) reconstructed track parametrizations auto recoTrack = getTrackParCov(track); - if (recoTrack.propagateToDCA(pvVtx, dBz)) { + if (recoTrack.propagateToDCA(pvVtx, magneticField)) { // xPv = recoTrack.getX(); } // Straight to Nsigma - float deltaThetaBarrelRich[5]; //, nSigmaBarrelRich[5]; - int lpdg_array[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; - float masses[5]; - - for (int ii = 0; ii < 5; ii++) { - // nSigmaBarrelRich[ii] = error_value; - - auto pdgInfoThis = pdg->GetParticle(lpdg_array[ii]); - masses[ii] = pdgInfoThis->Mass(); - float hypothesisAngleBarrelRich = CherenkovAngle(recoTrack.getP(), masses[ii], aerogel_rindex[i_sector]); + static constexpr int kNspecies = 5; + float nSigmaBarrelRich[kNspecies] = {error_value, error_value, error_value, error_value, error_value}; + bool signalBarrelRich[kNspecies] = {false, false, false, false, false}; + float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies]; + static constexpr int lpdg_array[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static constexpr float masses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron], + o2::track::pid_constants::sMasses[o2::track::PID::Muon], + o2::track::pid_constants::sMasses[o2::track::PID::Pion], + o2::track::pid_constants::sMasses[o2::track::PID::Kaon], + o2::track::pid_constants::sMasses[o2::track::PID::Proton]}; + + for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses + + float hypothesisAngleBarrelRich = error_value; + const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogel_rindex[i_sector], hypothesisAngleBarrelRich); + signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons // Evaluate total sigma (layer + tracking resolution) float barrelTotalAngularReso = barrelRICHAngularResolution; if (flagIncludeTrackAngularRes) { - double pt_resolution = std::pow(recoTrack.getP() / std::cosh(recoTrack.getEta()), 2) * std::sqrt(recoTrack.getSigma1Pt2()); + const float transverseMomentum = recoTrack.getP() / std::cosh(recoTrack.getEta()); + double pt_resolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2()); double eta_resolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2()); if (flagRICHLoadDelphesLUTs) { - pt_resolution = mSmearer.getAbsPtRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta())); - eta_resolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta())); + pt_resolution = mSmearer.getAbsPtRes(lpdg_array[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + eta_resolution = mSmearer.getAbsEtaRes(lpdg_array[ii], dNdEta, recoTrack.getEta(), transverseMomentum); } // cout << endl << "Pt resolution: " << pt_resolution << ", Eta resolution: " << eta_resolution << endl << endl; - float barrelTrackAngularReso = calculate_track_angular_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], aerogel_rindex[i_sector]); + const float barrelTrackAngularReso = calculate_track_angular_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], aerogel_rindex[i_sector]); barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso); - if (doQAplots && hypothesisAngleBarrelRich > error_value + 1. && measuredAngleBarrelRich > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) { - float momentum = recoTrack.getP(); - // float pseudorapidity = recoTrack.getEta(); - // float transverse_momentum = momentum / std::cosh(pseudorapidity); - if (ii == 0 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) { - histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), momentum, 1000.0 * barrelTrackAngularReso); - histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), momentum, 1000.0 * barrelTotalAngularReso); - } - if (ii == 1 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) { - histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), momentum, 1000.0 * barrelTrackAngularReso); - histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), momentum, 1000.0 * barrelTotalAngularReso); - } - if (ii == 2 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) { - histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), momentum, 1000.0 * barrelTrackAngularReso); - histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), momentum, 1000.0 * barrelTotalAngularReso); - } - if (ii == 3 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) { - histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), momentum, 1000.0 * barrelTrackAngularReso); - histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), momentum, 1000.0 * barrelTotalAngularReso); - } - if (ii == 4 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) { - histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), momentum, 1000.0 * barrelTrackAngularReso); - histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), momentum, 1000.0 * barrelTotalAngularReso); + if (doQAplots && + hypothesisAngleBarrelRich > error_value + 1. && + measuredAngleBarrelRich > error_value + 1. && + barrelRICHAngularResolution > error_value + 1. && + flagReachesRadiator) { + switch (mcParticle.pdgCode()) { + case lpdg_array[0]: // Electron + case -lpdg_array[0]: // Positron + if (ii == 0) { + histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case lpdg_array[1]: // Muon + case -lpdg_array[1]: // AntiMuon + if (ii == 1) { + histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case lpdg_array[2]: // Pion + case -lpdg_array[2]: // AntiPion + if (ii == 2) { + histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case lpdg_array[3]: // Kaon + case -lpdg_array[3]: // AntiKaon + if (ii == 3) { + histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case lpdg_array[4]: // Proton + case -lpdg_array[4]: // AntiProton + if (ii == 4) { + histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + default: + break; } } } @@ -854,7 +916,10 @@ struct OnTheFlyRichPid { /// DISCLAIMER: here tracking is accounted only for momentum value, but not for track parameters at impact point on the /// RICH radiator, since exact resolution would require photon generation and transport to photodetector. /// Effects are expected to be negligible (a few tenths of a milliradian) but further studies are required ! - if (hypothesisAngleBarrelRich > error_value + 1. && measuredAngleBarrelRich > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) { + if (hypothesisAngleBarrelRich > error_value + 1. && + measuredAngleBarrelRich > error_value + 1. && + barrelRICHAngularResolution > error_value + 1. && + flagReachesRadiator) { deltaThetaBarrelRich[ii] = hypothesisAngleBarrelRich - measuredAngleBarrelRich; nSigmaBarrelRich[ii] = deltaThetaBarrelRich[ii] / barrelTotalAngularReso; } @@ -862,57 +927,65 @@ struct OnTheFlyRichPid { // Fill histograms if (doQAplots) { - float momentum = recoTrack.getP(); - float pseudorapidity = recoTrack.getEta(); - float barrelRichTheta = measuredAngleBarrelRich; - - if (barrelRichTheta > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) { - histos.fill(HIST("h2dAngleVsMomentumBarrelRICH"), momentum, barrelRichTheta); - histos.fill(HIST("h2dAngleVsEtaBarrelRICH"), pseudorapidity, barrelRichTheta); - histos.fill(HIST("h2dAngularResolutionVsMomentumBarrelRICH"), momentum, 1000 * barrelRICHAngularResolution); - histos.fill(HIST("h2dAngularResolutionVsEtaBarrelRICH"), pseudorapidity, 1000 * barrelRICHAngularResolution); + if (measuredAngleBarrelRich > error_value + 1. && + barrelRICHAngularResolution > error_value + 1. && + flagReachesRadiator) { + histos.fill(HIST("h2dAngleVsMomentumBarrelRICH"), recoTrack.getP(), measuredAngleBarrelRich); + histos.fill(HIST("h2dAngleVsEtaBarrelRICH"), recoTrack.getEta(), measuredAngleBarrelRich); + histos.fill(HIST("h2dAngularResolutionVsMomentumBarrelRICH"), recoTrack.getP(), 1000 * barrelRICHAngularResolution); + histos.fill(HIST("h2dAngularResolutionVsEtaBarrelRICH"), recoTrack.getEta(), 1000 * barrelRICHAngularResolution); histos.fill(HIST("hSectorID"), i_sector); - if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) { - histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), momentum, nSigmaBarrelRich[0]); - histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]); - histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), momentum, nSigmaBarrelRich[2]); - histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]); - histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), momentum, nSigmaBarrelRich[4]); - } - if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) { - histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), momentum, nSigmaBarrelRich[0]); - histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]); - histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), momentum, nSigmaBarrelRich[2]); - histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]); - histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), momentum, nSigmaBarrelRich[4]); - } - if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) { - histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), momentum, nSigmaBarrelRich[0]); - histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]); - histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), momentum, nSigmaBarrelRich[2]); - histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]); - histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), momentum, nSigmaBarrelRich[4]); - } - if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) { - histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), momentum, nSigmaBarrelRich[0]); - histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]); - histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), momentum, nSigmaBarrelRich[2]); - histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]); - histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), momentum, nSigmaBarrelRich[4]); - } - if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) { - histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), momentum, nSigmaBarrelRich[0]); - histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]); - histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), momentum, nSigmaBarrelRich[2]); - histos.fill(HIST("h2dBarrelNsigmaTrueProtVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]); - histos.fill(HIST("h2dBarrelNsigmaTrueProtVsProtHypothesis"), momentum, nSigmaBarrelRich[4]); + switch (mcParticle.pdgCode()) { + case lpdg_array[0]: // Electron + case -lpdg_array[0]: // Positron + histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); + histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); + histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); + histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); + histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + break; + case lpdg_array[1]: // Muon + case -lpdg_array[1]: // AntiMuon + histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); + histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); + histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); + histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); + histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + break; + case lpdg_array[2]: // Pion + case -lpdg_array[2]: // AntiPion + histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); + histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); + histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); + histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); + histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + break; + case lpdg_array[3]: // Kaon + case -lpdg_array[3]: // AntiKaon + histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); + histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); + histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); + histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); + histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + break; + case lpdg_array[4]: // Proton + case -lpdg_array[4]: // AntiProton + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + break; + default: + break; } } } // Sigmas have been fully calculated. Please populate the NSigma helper table (once per track) upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]); + upgradeRichSignal(expectedAngleBarrelRichOk, signalBarrelRich[0], signalBarrelRich[1], signalBarrelRich[2], signalBarrelRich[3], signalBarrelRich[4]); } } }; From 656233c2460a7e0d5e481a18a76ff6c344013f63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 18:50:22 +0200 Subject: [PATCH 2/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 121 +++++++++---------- 1 file changed, 58 insertions(+), 63 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index c61db969688..abd9e75201f 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -170,15 +170,9 @@ struct OnTheFlyRichPid { HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; /// Flag unphysical and unavailable values (must be negative) - float error_value = -1000; + static constexpr float kErrorValue = -1000; // Variables projective/hybrid layout - int mNumberSectors = bRichNumberOfSectors; // 21; - float mTileLength = bRichPhotodetectorOtherModuleLength; // 18.4; // [cm] - float mTileLengthCentral = bRichPhotodetectorCentralModuleHalfLength; // 18.4 / 2.0; // [cm] - float mProjectiveLengthInner = mTileLengthCentral; - float mRadiusProjIn = bRichRadiatorInnerRadius; // 85.; // [cm] - float mRadiusProjOut = bRichPhotodetectorOuterRadius; // 112.; // [cm] std::vector det_centers; std::vector rad_centers; std::vector angle_centers; @@ -192,13 +186,7 @@ struct OnTheFlyRichPid { // Update projective geometry void updateProjectiveParameters() { - mNumberSectors = bRichNumberOfSectors; - mTileLength = bRichPhotodetectorOtherModuleLength; - mTileLengthCentral = bRichPhotodetectorCentralModuleHalfLength; - mProjectiveLengthInner = mTileLengthCentral; - mRadiusProjIn = bRichRadiatorInnerRadius; - mRadiusProjOut = bRichPhotodetectorOuterRadius; - const int number_of_sectors_in_z = mNumberSectors; + const int number_of_sectors_in_z = bRichNumberOfSectors; det_centers.resize(number_of_sectors_in_z); rad_centers.resize(number_of_sectors_in_z); angle_centers.resize(number_of_sectors_in_z); @@ -207,10 +195,10 @@ struct OnTheFlyRichPid { aerogel_rindex.resize(number_of_sectors_in_z); photodetrctor_length.resize(number_of_sectors_in_z); gap_thickness.resize(number_of_sectors_in_z); - float square_size_barrel_cylinder = 2.0 * mTileLengthCentral; - float square_size_z = mTileLength; - float R_min = mRadiusProjIn; - float R_max = mRadiusProjOut; + float square_size_barrel_cylinder = 2.0 * bRichPhotodetectorCentralModuleHalfLength; + float square_size_z = bRichPhotodetectorOtherModuleLength; + float R_min = bRichRadiatorInnerRadius; + float R_max = bRichPhotodetectorOuterRadius; std::vector theta_bi; std::vector R0_tilt; std::vector z0_tilt; @@ -225,7 +213,8 @@ struct OnTheFlyRichPid { l_detector_z.resize(number_of_sectors_in_z); // Odd number of sectors - if (number_of_sectors_in_z % 2 != 0) { + static constexpr int kTwo = 2; + if (number_of_sectors_in_z % kTwo != 0) { int i_central_mirror = static_cast((number_of_sectors_in_z) / 2.0); float m_val = std::tan(0.0); theta_bi[i_central_mirror] = std::atan(m_val); @@ -237,7 +226,7 @@ struct OnTheFlyRichPid { float t = std::tan(std::atan(m_val) + std::atan(square_size_barrel_cylinder / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_barrel_cylinder * m_val))); theta_max[i_central_mirror] = o2::constants::math::PIHalf - std::atan(t); theta_min[i_central_mirror] = o2::constants::math::PIHalf + std::atan(t); - mProjectiveLengthInner = R_min * t; + bRichPhotodetectorCentralModuleHalfLength = R_min * t; aerogel_rindex[i_central_mirror] = bRichRefractiveIndexSector[0]; for (int i = i_central_mirror + 1; i < number_of_sectors_in_z; i++) { float par_a = t; @@ -264,14 +253,14 @@ struct OnTheFlyRichPid { l_aerogel_z[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); T_r_plus_g[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); aerogel_rindex[2 * i_central_mirror - i] = bRichRefractiveIndexSector[i - i_central_mirror]; - mProjectiveLengthInner = R_min * t; // <-- At the end of the loop this will be the maximum Z + bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z } } else { // Even number of sectors float two_half_gap = 1.0; int i_central_mirror = static_cast((number_of_sectors_in_z) / 2.0); float m_val = std::tan(0.0); float t = std::tan(std::atan(m_val) + std::atan(two_half_gap / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - two_half_gap * m_val))); - mProjectiveLengthInner = R_min * t; + bRichPhotodetectorCentralModuleHalfLength = R_min * t; for (int i = i_central_mirror; i < number_of_sectors_in_z; i++) { float par_a = t; float par_b = 2.0 * R_max / square_size_z; @@ -297,7 +286,7 @@ struct OnTheFlyRichPid { l_aerogel_z[2 * i_central_mirror - i - 1] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); T_r_plus_g[2 * i_central_mirror - i - 1] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); aerogel_rindex[2 * i_central_mirror - i - 1] = bRichRefractiveIndexSector[i - i_central_mirror]; - mProjectiveLengthInner = R_min * t; // <-- At the end of the loop this will be the maximum Z + bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z } } // Coordinate radiali layer considerati @@ -326,7 +315,7 @@ struct OnTheFlyRichPid { { pRandomNumberGenerator.SetSeed(0); // fully randomize - if (magneticField.value < 0.0001f) { + if (magneticField.value < o2::constants::math::Epsilon) { LOG(info) << "Getting the magnetic field from the on-the-fly tracker task"; if (!getTaskOptionValue(initContext, "on-the-fly-tracker", magneticField, false)) { LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; @@ -487,7 +476,7 @@ struct OnTheFlyRichPid { int findSector(const float eta) { float polar = 2.0 * std::atan(std::exp(-eta)); - for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) { + for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) { if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) { return j_sec; } @@ -508,7 +497,7 @@ struct OnTheFlyRichPid { const float z_sec_rich_squared = z_sec_rich * z_sec_rich; return (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar)); } else { - return error_value; + return kErrorValue; } } @@ -528,7 +517,8 @@ struct OnTheFlyRichPid { const float meanNumberofDetectedPhotons = 230. * std::sin(angle) * std::sin(angle) * bRichRadiatorThickness; // Require at least 3 photons on average for real angle reconstruction - if (meanNumberofDetectedPhotons > 3.) { + static constexpr float kMinPhotons = 3.f; + if (meanNumberofDetectedPhotons > kMinPhotons) { return true; // Particle is above the threshold and enough photons } return false; // Particle is above the threshold, but not enough photons @@ -573,7 +563,7 @@ struct OnTheFlyRichPid { } } else { // std::cout << "Unable to interpolate. Target x value is outside the range of available data." << std::endl; - return error_value; + return kErrorValue; } } @@ -587,14 +577,14 @@ struct OnTheFlyRichPid { const float polar = 2.0 * std::atan(std::exp(-eta)); int i_sector = 0; bool flag_sector = false; - for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) { + for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) { if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) { flag_sector = true; i_sector = j_sec; } } if (flag_sector == false) { - return error_value; // <-- Returning negative value + return kErrorValue; // <-- Returning negative value } float R_sec_rich = rad_centers[i_sector].X(); float z_sec_rich = rad_centers[i_sector].Z(); @@ -627,8 +617,8 @@ struct OnTheFlyRichPid { float extract_ring_angular_resolution(const float eta, const float n, const float n_g, const float T_r, const float T_g, const float pixel_size, const float theta_c, const float tile_z_length) { // Check if input angle is error value - if (theta_c <= error_value + 1) - return error_value; + if (theta_c <= kErrorValue + 1) + return kErrorValue; // Parametrization variables (notation from https://doi.org/10.1016/0168-9002(94)90532-0) const float phi_c = 0.; const float theta_p = 0.; @@ -700,8 +690,8 @@ struct OnTheFlyRichPid { // float n_photons = (tile_z_length / 2.0 > radius) ? N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length)) : N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length) - (2.0/(tile_z_length*o2::constants::math::PI))*(-(tile_z_length/(2.0))*std::acos(tile_z_length/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tile_z_length/(2.0*radius),2.0)))); // Considering "exact" resolution (eta by eta) float n_photons = N0 * multiplicity_spectrum_factor * fractionPhotonsProjectiveRICH(eta, tile_z_length, radius); - if (n_photons <= error_value + 1) - return error_value; + if (n_photons <= kErrorValue + 1) + return kErrorValue; // Ring angular resolution float ring_angular_resolution = single_photon_angular_resolution / std::sqrt(n_photons); return ring_angular_resolution; @@ -774,7 +764,7 @@ struct OnTheFlyRichPid { for (const auto& track : tracks) { auto fillDummyValues = [&]() { - upgradeRich(error_value, error_value, error_value, error_value, error_value); + upgradeRich(kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue); upgradeRichSignal(false, false, false, false, false, false); }; @@ -788,7 +778,7 @@ struct OnTheFlyRichPid { const auto& mcParticle = track.mcParticle(); o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg); - // float xPv = error_value; + // float xPv = kErrorValue; if (o2track.propagateToDCA(mcPvVtx, magneticField)) { // xPv = o2track.getX(); } @@ -807,7 +797,7 @@ struct OnTheFlyRichPid { continue; } - float expectedAngleBarrelRich = error_value; + float expectedAngleBarrelRich = kErrorValue; const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogel_rindex[i_sector], expectedAngleBarrelRich); if (!expectedAngleBarrelRichOk) { fillDummyValues(); @@ -817,7 +807,7 @@ struct OnTheFlyRichPid { const float barrelRICHAngularResolution = extract_ring_angular_resolution(o2track.getEta(), aerogel_rindex[i_sector], bRichGapRefractiveIndex, bRichRadiatorThickness, gap_thickness[i_sector], bRICHPixelSize, expectedAngleBarrelRich, photodetrctor_length[i_sector]); const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), i_sector); bool flagReachesRadiator = false; - if (projectiveRadiatorRadius > error_value + 1.) { + if (projectiveRadiatorRadius > kErrorValue + 1.) { flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, magneticField); } /// DISCLAIMER: Exact extrapolation of angular resolution would require track propagation @@ -837,7 +827,12 @@ struct OnTheFlyRichPid { // Straight to Nsigma static constexpr int kNspecies = 5; - float nSigmaBarrelRich[kNspecies] = {error_value, error_value, error_value, error_value, error_value}; + static constexpr int kEl = 0; + static constexpr int kMu = 1; + static constexpr int kPi = 2; + static constexpr int kKa = 3; + static constexpr int kPr = 4; + float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue}; bool signalBarrelRich[kNspecies] = {false, false, false, false, false}; float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies]; static constexpr int lpdg_array[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; @@ -849,7 +844,7 @@ struct OnTheFlyRichPid { for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses - float hypothesisAngleBarrelRich = error_value; + float hypothesisAngleBarrelRich = kErrorValue; const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogel_rindex[i_sector], hypothesisAngleBarrelRich); signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons @@ -867,42 +862,42 @@ struct OnTheFlyRichPid { const float barrelTrackAngularReso = calculate_track_angular_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], aerogel_rindex[i_sector]); barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso); if (doQAplots && - hypothesisAngleBarrelRich > error_value + 1. && - measuredAngleBarrelRich > error_value + 1. && - barrelRICHAngularResolution > error_value + 1. && + hypothesisAngleBarrelRich > kErrorValue + 1. && + measuredAngleBarrelRich > kErrorValue + 1. && + barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { switch (mcParticle.pdgCode()) { - case lpdg_array[0]: // Electron - case -lpdg_array[0]: // Positron - if (ii == 0) { + case lpdg_array[kEl]: // Electron + case -lpdg_array[kEl]: // Positron + if (ii == kEl) { histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[1]: // Muon - case -lpdg_array[1]: // AntiMuon - if (ii == 1) { + case lpdg_array[kMu]: // Muon + case -lpdg_array[kMu]: // AntiMuon + if (ii == kMu) { histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[2]: // Pion - case -lpdg_array[2]: // AntiPion - if (ii == 2) { + case lpdg_array[kPi]: // Pion + case -lpdg_array[kPi]: // AntiPion + if (ii == kPi) { histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[3]: // Kaon - case -lpdg_array[3]: // AntiKaon - if (ii == 3) { + case lpdg_array[kKa]: // Kaon + case -lpdg_array[kKa]: // AntiKaon + if (ii == kKa) { histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[4]: // Proton - case -lpdg_array[4]: // AntiProton - if (ii == 4) { + case lpdg_array[kPr]: // Proton + case -lpdg_array[kPr]: // AntiProton + if (ii == kPr) { histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } @@ -916,9 +911,9 @@ struct OnTheFlyRichPid { /// DISCLAIMER: here tracking is accounted only for momentum value, but not for track parameters at impact point on the /// RICH radiator, since exact resolution would require photon generation and transport to photodetector. /// Effects are expected to be negligible (a few tenths of a milliradian) but further studies are required ! - if (hypothesisAngleBarrelRich > error_value + 1. && - measuredAngleBarrelRich > error_value + 1. && - barrelRICHAngularResolution > error_value + 1. && + if (hypothesisAngleBarrelRich > kErrorValue + 1. && + measuredAngleBarrelRich > kErrorValue + 1. && + barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { deltaThetaBarrelRich[ii] = hypothesisAngleBarrelRich - measuredAngleBarrelRich; nSigmaBarrelRich[ii] = deltaThetaBarrelRich[ii] / barrelTotalAngularReso; @@ -927,8 +922,8 @@ struct OnTheFlyRichPid { // Fill histograms if (doQAplots) { - if (measuredAngleBarrelRich > error_value + 1. && - barrelRICHAngularResolution > error_value + 1. && + if (measuredAngleBarrelRich > kErrorValue + 1. && + barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { histos.fill(HIST("h2dAngleVsMomentumBarrelRICH"), recoTrack.getP(), measuredAngleBarrelRich); histos.fill(HIST("h2dAngleVsEtaBarrelRICH"), recoTrack.getEta(), measuredAngleBarrelRich); From 933b09f17ce4785168d3e7c85fe760d072d49351 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 22:59:17 +0200 Subject: [PATCH 3/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 564 +++++++++---------- 1 file changed, 282 insertions(+), 282 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index abd9e75201f..ee7fd93e881 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -173,140 +173,140 @@ struct OnTheFlyRichPid { static constexpr float kErrorValue = -1000; // Variables projective/hybrid layout - std::vector det_centers; - std::vector rad_centers; - std::vector angle_centers; - std::vector theta_min; - std::vector theta_max; + std::vector detCenters; + std::vector radCenters; + std::vector angleCenters; + std::vector thetaMin; + std::vector thetaMax; std::vector bRichRefractiveIndexSector; - std::vector aerogel_rindex; - std::vector photodetrctor_length; - std::vector gap_thickness; + std::vector aerogelRindex; + std::vector photodetrctorLength; + std::vector gapThickness; // Update projective geometry void updateProjectiveParameters() { - const int number_of_sectors_in_z = bRichNumberOfSectors; - det_centers.resize(number_of_sectors_in_z); - rad_centers.resize(number_of_sectors_in_z); - angle_centers.resize(number_of_sectors_in_z); - theta_min.resize(number_of_sectors_in_z); - theta_max.resize(number_of_sectors_in_z); - aerogel_rindex.resize(number_of_sectors_in_z); - photodetrctor_length.resize(number_of_sectors_in_z); - gap_thickness.resize(number_of_sectors_in_z); - float square_size_barrel_cylinder = 2.0 * bRichPhotodetectorCentralModuleHalfLength; - float square_size_z = bRichPhotodetectorOtherModuleLength; - float R_min = bRichRadiatorInnerRadius; - float R_max = bRichPhotodetectorOuterRadius; - std::vector theta_bi; - std::vector R0_tilt; - std::vector z0_tilt; - std::vector T_r_plus_g; - std::vector l_aerogel_z; - std::vector l_detector_z; - theta_bi.resize(number_of_sectors_in_z); - R0_tilt.resize(number_of_sectors_in_z); - z0_tilt.resize(number_of_sectors_in_z); - T_r_plus_g.resize(number_of_sectors_in_z); - l_aerogel_z.resize(number_of_sectors_in_z); - l_detector_z.resize(number_of_sectors_in_z); + const int numberOfSectorsInZ = bRichNumberOfSectors; + detCenters.resize(numberOfSectorsInZ); + radCenters.resize(numberOfSectorsInZ); + angleCenters.resize(numberOfSectorsInZ); + thetaMin.resize(numberOfSectorsInZ); + thetaMax.resize(numberOfSectorsInZ); + aerogelRindex.resize(numberOfSectorsInZ); + photodetrctorLength.resize(numberOfSectorsInZ); + gapThickness.resize(numberOfSectorsInZ); + float squareSizeBarrelCylinder = 2.0 * bRichPhotodetectorCentralModuleHalfLength; + float squareSizeZ = bRichPhotodetectorOtherModuleLength; + float rMin = bRichRadiatorInnerRadius; + float rMax = bRichPhotodetectorOuterRadius; + std::vector thetaBi; + std::vector r0Tilt; + std::vector z0Tilt; + std::vector tRPlusG; + std::vector lAerogelZ; + std::vector lDetectorZ; + thetaBi.resize(numberOfSectorsInZ); + r0Tilt.resize(numberOfSectorsInZ); + z0Tilt.resize(numberOfSectorsInZ); + tRPlusG.resize(numberOfSectorsInZ); + lAerogelZ.resize(numberOfSectorsInZ); + lDetectorZ.resize(numberOfSectorsInZ); // Odd number of sectors static constexpr int kTwo = 2; - if (number_of_sectors_in_z % kTwo != 0) { - int i_central_mirror = static_cast((number_of_sectors_in_z) / 2.0); - float m_val = std::tan(0.0); - theta_bi[i_central_mirror] = std::atan(m_val); - R0_tilt[i_central_mirror] = R_max; - z0_tilt[i_central_mirror] = R0_tilt[i_central_mirror] * std::tan(theta_bi[i_central_mirror]); - l_detector_z[i_central_mirror] = square_size_barrel_cylinder; - l_aerogel_z[i_central_mirror] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_barrel_cylinder / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_barrel_cylinder); - T_r_plus_g[i_central_mirror] = R_max - R_min; - float t = std::tan(std::atan(m_val) + std::atan(square_size_barrel_cylinder / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_barrel_cylinder * m_val))); - theta_max[i_central_mirror] = o2::constants::math::PIHalf - std::atan(t); - theta_min[i_central_mirror] = o2::constants::math::PIHalf + std::atan(t); - bRichPhotodetectorCentralModuleHalfLength = R_min * t; - aerogel_rindex[i_central_mirror] = bRichRefractiveIndexSector[0]; - for (int i = i_central_mirror + 1; i < number_of_sectors_in_z; i++) { - float par_a = t; - float par_b = 2.0 * R_max / square_size_z; - m_val = (std::sqrt(par_a * par_a * par_b * par_b + par_b * par_b - 1.0) + par_a * par_b * par_b) / (par_b * par_b - 1.0); - theta_min[i] = o2::constants::math::PIHalf - std::atan(t); - theta_max[2 * i_central_mirror - i] = o2::constants::math::PIHalf + std::atan(t); - t = std::tan(std::atan(m_val) + std::atan(square_size_z / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_z * m_val))); - theta_max[i] = o2::constants::math::PIHalf - std::atan(t); - theta_min[2 * i_central_mirror - i] = o2::constants::math::PIHalf + std::atan(t); + if (numberOfSectorsInZ % kTwo != 0) { + int iCentralMirror = static_cast((numberOfSectorsInZ) / 2.0); + float mVal = std::tan(0.0); + thetaBi[iCentralMirror] = std::atan(mVal); + r0Tilt[iCentralMirror] = rMax; + z0Tilt[iCentralMirror] = r0Tilt[iCentralMirror] * std::tan(thetaBi[iCentralMirror]); + lDetectorZ[iCentralMirror] = squareSizeBarrelCylinder; + lAerogelZ[iCentralMirror] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeBarrelCylinder / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeBarrelCylinder); + tRPlusG[iCentralMirror] = rMax - rMin; + float t = std::tan(std::atan(mVal) + std::atan(squareSizeBarrelCylinder / (2.0 * rMax * std::sqrt(1.0 + mVal * mVal) - squareSizeBarrelCylinder * mVal))); + thetaMax[iCentralMirror] = o2::constants::math::PIHalf - std::atan(t); + thetaMin[iCentralMirror] = o2::constants::math::PIHalf + std::atan(t); + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; + aerogelRindex[iCentralMirror] = bRichRefractiveIndexSector[0]; + for (int i = iCentralMirror + 1; i < numberOfSectorsInZ; i++) { + float parA = t; + float parB = 2.0 * rMax / squareSizeZ; + mVal = (std::sqrt(parA * parA * parB * parB + parB * parB - 1.0) + parA * parB * parB) / (parB * parB - 1.0); + thetaMin[i] = o2::constants::math::PIHalf - std::atan(t); + thetaMax[2 * iCentralMirror - i] = o2::constants::math::PIHalf + std::atan(t); + t = std::tan(std::atan(mVal) + std::atan(squareSizeZ / (2.0 * rMax * std::sqrt(1.0 + mVal * mVal) - squareSizeZ * mVal))); + thetaMax[i] = o2::constants::math::PIHalf - std::atan(t); + thetaMin[2 * iCentralMirror - i] = o2::constants::math::PIHalf + std::atan(t); // Forward sectors - theta_bi[i] = std::atan(m_val); - R0_tilt[i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); - z0_tilt[i] = R0_tilt[i] * std::tan(theta_bi[i]); - l_detector_z[i] = square_size_z; - l_aerogel_z[i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); - T_r_plus_g[i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); - aerogel_rindex[i] = bRichRefractiveIndexSector[i - i_central_mirror]; + thetaBi[i] = std::atan(mVal); + r0Tilt[i] = rMax - squareSizeZ / 2.0 * std::sin(std::atan(mVal)); + z0Tilt[i] = r0Tilt[i] * std::tan(thetaBi[i]); + lDetectorZ[i] = squareSizeZ; + lAerogelZ[i] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); + tRPlusG[i] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); + aerogelRindex[i] = bRichRefractiveIndexSector[i - iCentralMirror]; // Backword sectors - theta_bi[2 * i_central_mirror - i] = -std::atan(m_val); - R0_tilt[2 * i_central_mirror - i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); - z0_tilt[2 * i_central_mirror - i] = -R0_tilt[i] * std::tan(theta_bi[i]); - l_detector_z[2 * i_central_mirror - i] = square_size_z; - l_aerogel_z[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); - T_r_plus_g[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); - aerogel_rindex[2 * i_central_mirror - i] = bRichRefractiveIndexSector[i - i_central_mirror]; - bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z + thetaBi[2 * iCentralMirror - i] = -std::atan(mVal); + r0Tilt[2 * iCentralMirror - i] = rMax - squareSizeZ / 2.0 * std::sin(std::atan(mVal)); + z0Tilt[2 * iCentralMirror - i] = -r0Tilt[i] * std::tan(thetaBi[i]); + lDetectorZ[2 * iCentralMirror - i] = squareSizeZ; + lAerogelZ[2 * iCentralMirror - i] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); + tRPlusG[2 * iCentralMirror - i] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); + aerogelRindex[2 * iCentralMirror - i] = bRichRefractiveIndexSector[i - iCentralMirror]; + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum Z } } else { // Even number of sectors - float two_half_gap = 1.0; - int i_central_mirror = static_cast((number_of_sectors_in_z) / 2.0); - float m_val = std::tan(0.0); - float t = std::tan(std::atan(m_val) + std::atan(two_half_gap / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - two_half_gap * m_val))); - bRichPhotodetectorCentralModuleHalfLength = R_min * t; - for (int i = i_central_mirror; i < number_of_sectors_in_z; i++) { - float par_a = t; - float par_b = 2.0 * R_max / square_size_z; - m_val = (std::sqrt(par_a * par_a * par_b * par_b + par_b * par_b - 1.0) + par_a * par_b * par_b) / (par_b * par_b - 1.0); - theta_min[i] = o2::constants::math::PIHalf - std::atan(t); - theta_max[2 * i_central_mirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); - t = std::tan(std::atan(m_val) + std::atan(square_size_z / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_z * m_val))); - theta_max[i] = o2::constants::math::PIHalf - std::atan(t); - theta_min[2 * i_central_mirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); + float twoHalfGap = 1.0; + int iCentralMirror = static_cast((numberOfSectorsInZ) / 2.0); + float mVal = std::tan(0.0); + float t = std::tan(std::atan(mVal) + std::atan(twoHalfGap / (2.0 * rMax * std::sqrt(1.0 + mVal * mVal) - twoHalfGap * mVal))); + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; + for (int i = iCentralMirror; i < numberOfSectorsInZ; i++) { + float parA = t; + float parB = 2.0 * rMax / squareSizeZ; + mVal = (std::sqrt(parA * parA * parB * parB + parB * parB - 1.0) + parA * parB * parB) / (parB * parB - 1.0); + thetaMin[i] = o2::constants::math::PIHalf - std::atan(t); + thetaMax[2 * iCentralMirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); + t = std::tan(std::atan(mVal) + std::atan(squareSizeZ / (2.0 * rMax * std::sqrt(1.0 + mVal * mVal) - squareSizeZ * mVal))); + thetaMax[i] = o2::constants::math::PIHalf - std::atan(t); + thetaMin[2 * iCentralMirror - i - 1] = o2::constants::math::PIHalf + std::atan(t); // Forward sectors - theta_bi[i] = std::atan(m_val); - R0_tilt[i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); - z0_tilt[i] = R0_tilt[i] * std::tan(theta_bi[i]); - l_detector_z[i] = square_size_z; - l_aerogel_z[i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); - T_r_plus_g[i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); - aerogel_rindex[i] = bRichRefractiveIndexSector[i - i_central_mirror]; + thetaBi[i] = std::atan(mVal); + r0Tilt[i] = rMax - squareSizeZ / 2.0 * std::sin(std::atan(mVal)); + z0Tilt[i] = r0Tilt[i] * std::tan(thetaBi[i]); + lDetectorZ[i] = squareSizeZ; + lAerogelZ[i] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); + tRPlusG[i] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); + aerogelRindex[i] = bRichRefractiveIndexSector[i - iCentralMirror]; // Backword sectors - theta_bi[2 * i_central_mirror - i - 1] = -std::atan(m_val); - R0_tilt[2 * i_central_mirror - i - 1] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val)); - z0_tilt[2 * i_central_mirror - i - 1] = -R0_tilt[i] * std::tan(theta_bi[i]); - l_detector_z[2 * i_central_mirror - i - 1] = square_size_z; - l_aerogel_z[2 * i_central_mirror - i - 1] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z); - T_r_plus_g[2 * i_central_mirror - i - 1] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]); - aerogel_rindex[2 * i_central_mirror - i - 1] = bRichRefractiveIndexSector[i - i_central_mirror]; - bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z + thetaBi[2 * iCentralMirror - i - 1] = -std::atan(mVal); + r0Tilt[2 * iCentralMirror - i - 1] = rMax - squareSizeZ / 2.0 * std::sin(std::atan(mVal)); + z0Tilt[2 * iCentralMirror - i - 1] = -r0Tilt[i] * std::tan(thetaBi[i]); + lDetectorZ[2 * iCentralMirror - i - 1] = squareSizeZ; + lAerogelZ[2 * iCentralMirror - i - 1] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); + tRPlusG[2 * iCentralMirror - i - 1] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); + aerogelRindex[2 * iCentralMirror - i - 1] = bRichRefractiveIndexSector[i - iCentralMirror]; + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum Z } } // Coordinate radiali layer considerati - std::vector R0_detector; - std::vector R0_aerogel; - R0_detector.resize(number_of_sectors_in_z); - R0_aerogel.resize(number_of_sectors_in_z); - for (int i = 0; i < number_of_sectors_in_z; i++) { - R0_detector[i] = R0_tilt[i]; // + (detector_thickness / 2.0) * std::cos(theta_bi[i]) NEGLIGIBLE - det_centers[i].SetXYZ(R0_detector[i], 0, R0_detector[i] * std::tan(theta_bi[i])); - R0_aerogel[i] = R0_tilt[i] - (T_r_plus_g[i] - bRichRadiatorThickness / 2.0) * std::cos(theta_bi[i]); - rad_centers[i].SetXYZ(R0_aerogel[i], 0, R0_aerogel[i] * std::tan(theta_bi[i])); - angle_centers[i] = theta_bi[i]; - photodetrctor_length[i] = l_detector_z[i]; - gap_thickness[i] = (det_centers[i] - rad_centers[i]).Mag() - bRichRadiatorThickness / 2.0; + std::vector r0Detector; + std::vector r0Aerogel; + r0Detector.resize(numberOfSectorsInZ); + r0Aerogel.resize(numberOfSectorsInZ); + for (int i = 0; i < numberOfSectorsInZ; i++) { + r0Detector[i] = r0Tilt[i]; // + (detector_thickness / 2.0) * std::cos(thetaBi[i]) NEGLIGIBLE + detCenters[i].SetXYZ(r0Detector[i], 0, r0Detector[i] * std::tan(thetaBi[i])); + r0Aerogel[i] = r0Tilt[i] - (tRPlusG[i] - bRichRadiatorThickness / 2.0) * std::cos(thetaBi[i]); + radCenters[i].SetXYZ(r0Aerogel[i], 0, r0Aerogel[i] * std::tan(thetaBi[i])); + angleCenters[i] = thetaBi[i]; + photodetrctorLength[i] = lDetectorZ[i]; + gapThickness[i] = (detCenters[i] - radCenters[i]).Mag() - bRichRadiatorThickness / 2.0; } // DEBUG // std::cout << std::endl << std::endl; - // for (int i = 0; i < number_of_sectors_in_z; i++) { - // std::cout << "Sector " << i << ": Gap = " << gap_thickness[i] << " cm, Index aerogel = " << aerogel_rindex[i] << ", Gap thickness = " << gap_thickness[i] << " cm" << std::endl; + // for (int i = 0; i < numberOfSectorsInZ; i++) { + // std::cout << "Sector " << i << ": Gap = " << gapThickness[i] << " cm, Index aerogel = " << aerogelRindex[i] << ", Gap thickness = " << gapThickness[i] << " cm" << std::endl; // } // std::cout << std::endl << std::endl; } @@ -379,25 +379,25 @@ struct OnTheFlyRichPid { histos.add("hSectorID", "hSectorID", kTH1F, {axisSector}); const int kNspec = 5; // electron, muon, pion, kaon, proton - std::string particle_names1[kNspec] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; - std::string particle_names2[kNspec] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; - for (int i_true = 0; i_true < kNspec; i_true++) { - std::string name_title_barrel_track_res = "h2dBarrelAngularResTrack" + particle_names2[i_true] + "VsP"; - std::string name_title_barrel_total_res = "h2dBarrelAngularResTotal" + particle_names2[i_true] + "VsP"; - const AxisSpec axisTrackAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Track angular resolution - " + particle_names1[i_true] + " (mrad)"}; - const AxisSpec axisTotalAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Total angular resolution - " + particle_names1[i_true] + " (mrad)"}; - histos.add(name_title_barrel_track_res.c_str(), name_title_barrel_track_res.c_str(), kTH2F, {axisMomentum, axisTrackAngularRes}); - histos.add(name_title_barrel_total_res.c_str(), name_title_barrel_total_res.c_str(), kTH2F, {axisMomentum, axisTotalAngularRes}); + std::string particleNames1[kNspec] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; + std::string particleNames2[kNspec] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; + for (int iTrue = 0; iTrue < kNspec; iTrue++) { + std::string nameTitleBarrelTrackRes = "h2dBarrelAngularResTrack" + particleNames2[iTrue] + "VsP"; + std::string nameTitleBarrelTotalRes = "h2dBarrelAngularResTotal" + particleNames2[iTrue] + "VsP"; + const AxisSpec axisTrackAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Track angular resolution - " + particleNames1[iTrue] + " (mrad)"}; + const AxisSpec axisTotalAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Total angular resolution - " + particleNames1[iTrue] + " (mrad)"}; + histos.add(nameTitleBarrelTrackRes.c_str(), nameTitleBarrelTrackRes.c_str(), kTH2F, {axisMomentum, axisTrackAngularRes}); + histos.add(nameTitleBarrelTotalRes.c_str(), nameTitleBarrelTotalRes.c_str(), kTH2F, {axisMomentum, axisTotalAngularRes}); } - for (int i_true = 0; i_true < kNspec; i_true++) { - for (int i_hyp = 0; i_hyp < kNspec; i_hyp++) { - std::string name_title = "h2dBarrelNsigmaTrue" + particle_names2[i_true] + "Vs" + particle_names2[i_hyp] + "Hypothesis"; - if (i_true == i_hyp) { - const AxisSpec axisNsigmaCorrect{static_cast(nBinsNsigmaCorrectSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"}; - histos.add(name_title.c_str(), name_title.c_str(), kTH2F, {axisMomentum, axisNsigmaCorrect}); + for (int iTrue = 0; iTrue < kNspec; iTrue++) { + for (int iHyp = 0; iHyp < kNspec; iHyp++) { + std::string nameTitle = "h2dBarrelNsigmaTrue" + particleNames2[iTrue] + "Vs" + particleNames2[iHyp] + "Hypothesis"; + if (iTrue == iHyp) { + const AxisSpec axisNsigmaCorrect{static_cast(nBinsNsigmaCorrectSpecies), -10.0f, +10.0f, "N#sigma - True " + particleNames1[iTrue] + " vs " + particleNames1[iHyp] + " hypothesis"}; + histos.add(nameTitle.c_str(), nameTitle.c_str(), kTH2F, {axisMomentum, axisNsigmaCorrect}); } else { - const AxisSpec axisNsigmaWrong{static_cast(nBinsNsigmaWrongSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"}; - histos.add(name_title.c_str(), name_title.c_str(), kTH2F, {axisMomentum, axisNsigmaWrong}); + const AxisSpec axisNsigmaWrong{static_cast(nBinsNsigmaWrongSpecies), -10.0f, +10.0f, "N#sigma - True " + particleNames1[iTrue] + " vs " + particleNames1[iHyp] + " hypothesis"}; + histos.add(nameTitle.c_str(), nameTitle.c_str(), kTH2F, {axisMomentum, axisNsigmaWrong}); } } } @@ -476,9 +476,9 @@ struct OnTheFlyRichPid { int findSector(const float eta) { float polar = 2.0 * std::atan(std::exp(-eta)); - for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) { - if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) { - return j_sec; + for (int jSector = 0; jSector < bRichNumberOfSectors; jSector++) { + if (polar > thetaMax[jSector] && polar < thetaMin[jSector]) { + return jSector; } } return -1; @@ -486,16 +486,16 @@ struct OnTheFlyRichPid { /// returns radiator radius in cm at considered eta (accounting for sector inclination) /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) - /// \param i_sector the index of the track RICH sector - float radiusRipple(const float eta, const int i_sector) + /// \param iSecor the index of the track RICH sector + float radiusRipple(const float eta, const int iSecor) { - if (i_sector > 0) { + if (iSecor > 0) { const float polar = 2.0 * std::atan(std::exp(-eta)); - const float R_sec_rich = rad_centers[i_sector].X(); - const float z_sec_rich = rad_centers[i_sector].Z(); - const float R_sec_rich_squared = R_sec_rich * R_sec_rich; - const float z_sec_rich_squared = z_sec_rich * z_sec_rich; - return (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar)); + const float rSecRich = radCenters[iSecor].X(); + const float zSecRich = radCenters[iSecor].Z(); + const float rSecRichSquared = rSecRich * rSecRich; + const float zSecRichSquared = zSecRich * zSecRich; + return (rSecRichSquared + zSecRichSquared) / (rSecRich + zSecRich / std::tan(polar)); } else { return kErrorValue; } @@ -543,21 +543,21 @@ struct OnTheFlyRichPid { float AngularResolution(float eta) { // Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS) - static constexpr float eta_sampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; - static constexpr float res_ring_sampling_with_abs_walls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; - static constexpr float res_ring_sampling_without_abs_walls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; - static constexpr int size_res_vector = sizeof(eta_sampling) / sizeof(eta_sampling[0]); + static constexpr float etaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; + static constexpr float resRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; + static constexpr float resRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; + static constexpr int sizeResVector = sizeof(etaSampling) / sizeof(etaSampling[0]); // Use binary search to find the lower and upper indices - const int lowerIndex = std::lower_bound(eta_sampling, eta_sampling + size_res_vector, eta) - eta_sampling - 1; + const int lowerIndex = std::lower_bound(etaSampling, etaSampling + sizeResVector, eta) - etaSampling - 1; const int upperIndex = lowerIndex + 1; - if (lowerIndex >= 0 && upperIndex < size_res_vector) { + if (lowerIndex >= 0 && upperIndex < sizeResVector) { // Resolution interpolation if (bRichFlagAbsorbingWalls) { - float interpolatedResRing = interpolate(eta, eta_sampling[lowerIndex], eta_sampling[upperIndex], res_ring_sampling_with_abs_walls[lowerIndex], res_ring_sampling_with_abs_walls[upperIndex]); + float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithAbsWalls[lowerIndex], resRingSamplingWithAbsWalls[upperIndex]); // std::cout << "Interpolated y value: " << interpolatedY << std::endl; return interpolatedResRing; } else { - float interpolatedResRing = interpolate(eta, eta_sampling[lowerIndex], eta_sampling[upperIndex], res_ring_sampling_without_abs_walls[lowerIndex], res_ring_sampling_without_abs_walls[upperIndex]); + float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithoutAbsWalls[lowerIndex], resRingSamplingWithoutAbsWalls[upperIndex]); // std::cout << "Interpolated y value: " << interpolatedY << std::endl; return interpolatedResRing; } @@ -570,36 +570,36 @@ struct OnTheFlyRichPid { /// To account border effects in bRICH /// Approximation for analytic calculation: track along projective sector normal (reasonable) /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) - /// \param tile_z_length the Z-length of the photodetector plane of the considered sector + /// \param tileZlength the Z-length of the photodetector plane of the considered sector /// \param radius the nominal radius of the Cherenkov ring - float fractionPhotonsProjectiveRICH(float eta, float tile_z_length, float radius) + float fractionPhotonsProjectiveRICH(float eta, float tileZlength, float radius) { const float polar = 2.0 * std::atan(std::exp(-eta)); - int i_sector = 0; - bool flag_sector = false; - for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) { - if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) { - flag_sector = true; - i_sector = j_sec; + int iSecor = 0; + bool flagSector = false; + for (int jSector = 0; jSector < bRichNumberOfSectors; jSector++) { + if (polar > thetaMax[jSector] && polar < thetaMin[jSector]) { + flagSector = true; + iSecor = jSector; } } - if (flag_sector == false) { + if (flagSector == false) { return kErrorValue; // <-- Returning negative value } - float R_sec_rich = rad_centers[i_sector].X(); - float z_sec_rich = rad_centers[i_sector].Z(); - // float R_sec_tof = det_centers[i_sector].X(); - // float z_sec_tof = det_centers[i_sector].Z(); - const float R_sec_rich_squared = R_sec_rich * R_sec_rich; - const float z_sec_rich_squared = z_sec_rich * z_sec_rich; - const float radius_ripple = (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar)); - const float z_ripple = radius_ripple / std::tan(polar); - const float absZ = std::hypot(radius_ripple - R_sec_rich, z_ripple - z_sec_rich); + float rSecRich = radCenters[iSecor].X(); + float zSecRich = radCenters[iSecor].Z(); + // float rSecTof = detCenters[iSecor].X(); + // float zSecTof = detCenters[iSecor].Z(); + const float rSecRichSquared = rSecRich * rSecRich; + const float zSecRichSquared = zSecRich * zSecRich; + const float radiusRipple = (rSecRichSquared + zSecRichSquared) / (rSecRich + zSecRich / std::tan(polar)); + const float zRipple = radiusRipple / std::tan(polar); + const float absZ = std::hypot(radiusRipple - rSecRich, zRipple - zSecRich); float fraction = 1.; - if (tile_z_length / 2. - absZ < radius) { - fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tile_z_length / 2. - absZ) / radius); - if (tile_z_length / 2. + absZ < radius) { - fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tile_z_length / 2. + absZ) / radius); + if (tileZlength / 2. - absZ < radius) { + fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tileZlength / 2. - absZ) / radius); + if (tileZlength / 2. + absZ < radius) { + fraction = fraction - (1. / o2::constants::math::PI) * std::acos((tileZlength / 2. + absZ) / radius); } } return fraction; @@ -608,115 +608,115 @@ struct OnTheFlyRichPid { /// returns refined ring angular resolution /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) /// \param n the refractive index of the considered sector - /// \param n_g the refractive index of the gas filling the RICH proximity gap - /// \param T_r the radiator thickness in cm - /// \param T_g the proximity gap thickness of the considered sector in cm - /// \param pixel_size the SiPM pixel size in cm - /// \param theta_c the Cherenkov angle for the considered track - /// \param tile_z_length the Z-length of the photodetector plane of the considered sector - float extract_ring_angular_resolution(const float eta, const float n, const float n_g, const float T_r, const float T_g, const float pixel_size, const float theta_c, const float tile_z_length) + /// \param nGas the refractive index of the gas filling the RICH proximity gap + /// \param thicknessRad the radiator thickness in cm + /// \param thicknessGas the proximity gap thickness of the considered sector in cm + /// \param pixelSize the SiPM pixel size in cm + /// \param thetaCherenkov the Cherenkov angle for the considered track + /// \param tileZlength the Z-length of the photodetector plane of the considered sector + float extractRingAngularResolution(const float eta, const float n, const float nGas, const float thicknessRad, const float thicknessGas, const float pixelSize, const float thetaCherenkov, const float tileZlength) { // Check if input angle is error value - if (theta_c <= kErrorValue + 1) + if (thetaCherenkov <= kErrorValue + 1) return kErrorValue; // Parametrization variables (notation from https://doi.org/10.1016/0168-9002(94)90532-0) - const float phi_c = 0.; - const float theta_p = 0.; - const float sin_theta_c = std::sin(theta_c); - const float cos_theta_c = std::cos(theta_c); - const float x_p = std::sin(theta_p); - const float z_p = std::cos(theta_p); - const float sin_phi = std::sin(phi_c); - const float cos_phi = std::cos(phi_c); - // const float ze = T_r / (2. * z_p); - const float az = z_p * cos_theta_c - x_p * sin_theta_c * cos_phi; - const float e3z = std::sqrt(az * az + (n_g / n) * (n_g / n) - 1.); - const float Z = T_g; + const float phiC = 0.; + const float thetaP = 0.; + const float sinThetaCherenkov = std::sin(thetaCherenkov); + const float cosThetaCherenkov = std::cos(thetaCherenkov); + const float xP = std::sin(thetaP); + const float zP = std::cos(thetaP); + const float sinPhi = std::sin(phiC); + const float cosPhi = std::cos(phiC); + // const float ze = thicknessRad / (2. * zP); + const float az = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi; + const float e3z = std::sqrt(az * az + (nGas / n) * (nGas / n) - 1.); + const float Z = thicknessGas; const float alpha = e3z / az; const float etac = e3z * n * n; - const float k = T_r / (2. * Z); + const float k = thicknessRad / (2. * Z); const float m = 1. / (n * n); const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); - // Derivative d(theta_c)/dx + // Derivative d(thetaCherenkov)/dx const float temp1 = etac / Z; - const float temp2 = alpha * e3z * cos_phi; - const float temp3 = x_p * sin_theta_c * sin_phi * sin_phi; - const float d_theta_x = temp1 * (temp2 - temp3); - // Derivative d(theta_c)/dy - const float temp4 = etac * sin_phi / Z; - const float temp5 = cos_theta_c - z_p * (1 - m) / az; - const float d_theta_y = temp4 * temp5; - // Derivative d(theta_c)/dze - const float temp8 = etac * sin_theta_c; + const float temp2 = alpha * e3z * cosPhi; + const float temp3 = xP * sinThetaCherenkov * sinPhi * sinPhi; + const float dThetaX = temp1 * (temp2 - temp3); + // Derivative d(thetaCherenkov)/dy + const float temp4 = etac * sinPhi / Z; + const float temp5 = cosThetaCherenkov - zP * (1 - m) / az; + const float dThetaY = temp4 * temp5; + // Derivative d(thetaCherenkov)/dze + const float temp8 = etac * sinThetaCherenkov; const float temp9 = Z * (1.0 + k * alpha * alpha * alpha * n * n); - const float temp10 = alpha * alpha * (1.0 - x_p * x_p * sin_phi * sin_phi) + lambda * x_p * x_p * sin_phi * sin_phi; - const float d_theta_ze = temp8 * temp10 / temp9; - // Derivative d(theta_c)/dn - const float d_theta_n = z_p / (n * az * sin_theta_c); + const float temp10 = alpha * alpha * (1.0 - xP * xP * sinPhi * sinPhi) + lambda * xP * xP * sinPhi * sinPhi; + const float dThetaZe = temp8 * temp10 / temp9; + // Derivative d(thetaCherenkov)/dn + const float dThetaN = zP / (n * az * sinThetaCherenkov); // RMS wavelength (using Sellmeier formula with measured aerogel parameters) const float a0 = 0.0616; const float w0 = 56.5; const float n0 = 1.0307; - const float w_mean = 436.0; // 450.0;//484.9;//426.0; //450;//485;//450; // nm //450 426 493.5 // 426.0 + const float wMean = 436.0; // 450.0;//484.9;//426.0; //450;//485;//450; // nm //450 426 493.5 // 426.0 const float scaling = n / n0; - const float temp11 = a0 * w0 * w0 * w_mean; - const float temp12 = std::pow(w_mean * w_mean - w0 * w0, 1.5); - const float temp13 = std::sqrt(w_mean * w_mean * (1.0 + a0) - w0 * w0); - const float d_n_w = temp11 / (temp12 * temp13) * scaling; - const float sigma_w = 115.0; + const float temp11 = a0 * w0 * w0 * wMean; + const float temp12 = std::pow(wMean * wMean - w0 * w0, 1.5); + const float temp13 = std::sqrt(wMean * wMean * (1.0 + a0) - w0 * w0); + const float dNw = temp11 / (temp12 * temp13) * scaling; + const float sigmaW = 115.0; // RMS x y - const float sigma_theta_x = pixel_size / std::sqrt(12.0); - const float sigma_theta_y = pixel_size / std::sqrt(12.0); + const float sigmaThetaX = pixelSize / std::sqrt(12.0); + const float sigmaThetaY = pixelSize / std::sqrt(12.0); // RMS emission point - const float sigma_theta_ze = T_r / (z_p * std::sqrt(12.0)); + const float sigmaThetaZe = thicknessRad / (zP * std::sqrt(12.0)); // Single photon angular resolution - const float res_x = d_theta_x * sigma_theta_x; - const float res_y = d_theta_y * sigma_theta_y; - const float res_ze = d_theta_ze * sigma_theta_ze; - const float res_n = d_theta_n * d_n_w * sigma_w; - const float res_x_squared = res_x * res_x; - const float res_y_squared = res_y * res_y; - const float res_ze_squared = res_ze * res_ze; - const float res_n_squared = res_n * res_n; - const float single_photon_angular_resolution = std::sqrt(res_x_squared + res_y_squared + res_ze_squared + res_n_squared); + const float resX = dThetaX * sigmaThetaX; + const float resY = dThetaY * sigmaThetaY; + const float resZe = dThetaZe * sigmaThetaZe; + const float resN = dThetaN * dNw * sigmaW; + const float resXsquared = resX * resX; + const float resYsquared = resY * resY; + const float resZesquared = resZe * resZe; + const float resNsquared = resN * resN; + const float singlePhotonAngularResolution = std::sqrt(resXsquared + resYsquared + resZesquared + resNsquared); // Fraction of photons in rings (to account loss close to sector boundary in projective geometry) - const float sin_theta_c_squared = sin_theta_c * sin_theta_c; - const float radius = (T_r / 2.0) * std::tan(theta_c) + T_g * std::tan(std::asin((n / n_g) * sin_theta_c)); - const float N0 = 24. * T_r / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) - const float sin_angle = std::sin(std::acos(1. / 1.03)); - const float sin_angle_squared = sin_angle * sin_angle; - const float multiplicity_spectrum_factor = sin_theta_c_squared / sin_angle_squared; // scale multiplicity w.r.t. N = 1.03 at saturation + const float sinThetaCherenkovSquared = sinThetaCherenkov * sinThetaCherenkov; + const float radius = (thicknessRad / 2.0) * std::tan(thetaCherenkov) + thicknessGas * std::tan(std::asin((n / nGas) * sinThetaCherenkov)); + const float N0 = 24. * thicknessRad / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) + const float sinAngle = std::sin(std::acos(1. / 1.03)); + const float sinAngleSquared = sinAngle * sinAngle; + const float multiplicitySpectrumFactor = sinThetaCherenkovSquared / sinAngleSquared; // scale multiplicity w.r.t. N = 1.03 at saturation // Considering average resolution (integrated over the sector) - // float n_photons = (tile_z_length / 2.0 > radius) ? N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length)) : N0 * multiplicity_spectrum_factor * (1.-(2.0*radius)/(o2::constants::math::PI*tile_z_length) - (2.0/(tile_z_length*o2::constants::math::PI))*(-(tile_z_length/(2.0))*std::acos(tile_z_length/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tile_z_length/(2.0*radius),2.0)))); + // float nPhotons = (tileZlength / 2.0 > radius) ? N0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength)) : N0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength) - (2.0/(tileZlength*o2::constants::math::PI))*(-(tileZlength/(2.0))*std::acos(tileZlength/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tileZlength/(2.0*radius),2.0)))); // Considering "exact" resolution (eta by eta) - float n_photons = N0 * multiplicity_spectrum_factor * fractionPhotonsProjectiveRICH(eta, tile_z_length, radius); - if (n_photons <= kErrorValue + 1) + const float nPhotons = N0 * multiplicitySpectrumFactor * fractionPhotonsProjectiveRICH(eta, tileZlength, radius); + if (nPhotons <= kErrorValue + 1) return kErrorValue; // Ring angular resolution - float ring_angular_resolution = single_photon_angular_resolution / std::sqrt(n_photons); - return ring_angular_resolution; + const float ringAngularResolution = singlePhotonAngularResolution / std::sqrt(nPhotons); + return ringAngularResolution; } /// returns track angular resolution /// \param pt the transverse momentum of the tarck /// \param eta the pseudorapidity of the tarck - /// \param track_pt_resolution the absolute resolution on pt - /// \param track_pt_resolution the absolute resolution on eta + /// \param trackPtResolution the absolute resolution on pt + /// \param trackPtResolution the absolute resolution on eta /// \param mass the mass of the particle - /// \param refractive_index the refractive index of the radiator - double calculate_track_angular_resolution_advanced(const float pt, const float eta, const float track_pt_resolution, const float track_eta_resolution, const float mass, const float refractive_index) + /// \param refractiveIndex the refractive index of the radiator + double calculateTrackAngularResolutionAdvanced(const float pt, const float eta, const float trackPtResolution, const float trackEtaResolution, const float mass, const float refractiveIndex) { // Compute tracking contribution to timing using the error propagation formula // Uses light speed in m/ps, magnetic field in T (*0.1 for conversion kGauss -> T) const double a0 = mass * mass; - const double a1 = refractive_index; - const double a1_squared = a1 * a1; - const float pt_cosh_eta = pt * std::cosh(eta); - const float pt_cosh_eta_squared = pt_cosh_eta * pt_cosh_eta; - const double dtheta_on_dpt = a0 / (pt * std::sqrt(a0 + pt_cosh_eta_squared) * std::sqrt(pt_cosh_eta_squared * (a1_squared - 1.0) - a0)); - const double dtheta_on_deta = (a0 * std::tanh(eta)) / (std::sqrt(a0 + pt_cosh_eta_squared) * std::sqrt(pt_cosh_eta_squared * (a1_squared - 1.0) - a0)); - const double track_angular_resolution = std::hypot(std::fabs(dtheta_on_dpt) * track_pt_resolution, std::fabs(dtheta_on_deta) * track_eta_resolution); - return track_angular_resolution; + const double a1 = refractiveIndex; + const double a1Squared = a1 * a1; + const float ptCoshEta = pt * std::cosh(eta); + const float ptCoshEtaSquared = ptCoshEta * ptCoshEta; + const double dThetaOndPt = a0 / (pt * std::sqrt(a0 + ptCoshEtaSquared) * std::sqrt(ptCoshEtaSquared * (a1Squared - 1.0) - a0)); + const double dThetaOndEta = (a0 * std::tanh(eta)) / (std::sqrt(a0 + ptCoshEtaSquared) * std::sqrt(ptCoshEtaSquared * (a1Squared - 1.0) - a0)); + const double trackAngularResolution = std::hypot(std::fabs(dThetaOndPt) * trackPtResolution, std::fabs(dThetaOndEta) * trackEtaResolution); + return trackAngularResolution; } void process(soa::Join::iterator const& collision, @@ -791,21 +791,21 @@ struct OnTheFlyRichPid { } // find track bRICH sector - int i_sector = findSector(o2track.getEta()); - if (i_sector < 0) { + int iSecor = findSector(o2track.getEta()); + if (iSecor < 0) { fillDummyValues(); continue; } float expectedAngleBarrelRich = kErrorValue; - const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogel_rindex[i_sector], expectedAngleBarrelRich); + const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogelRindex[iSecor], expectedAngleBarrelRich); if (!expectedAngleBarrelRichOk) { fillDummyValues(); continue; // Particle is below the threshold or not enough photons } // float barrelRICHAngularResolution = AngularResolution(o2track.getEta()); - const float barrelRICHAngularResolution = extract_ring_angular_resolution(o2track.getEta(), aerogel_rindex[i_sector], bRichGapRefractiveIndex, bRichRadiatorThickness, gap_thickness[i_sector], bRICHPixelSize, expectedAngleBarrelRich, photodetrctor_length[i_sector]); - const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), i_sector); + const float barrelRICHAngularResolution = extractRingAngularResolution(o2track.getEta(), aerogelRindex[iSecor], bRichGapRefractiveIndex, bRichRadiatorThickness, gapThickness[iSecor], bRICHPixelSize, expectedAngleBarrelRich, photodetrctorLength[iSecor]); + const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), iSecor); bool flagReachesRadiator = false; if (projectiveRadiatorRadius > kErrorValue + 1.) { flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, magneticField); @@ -835,7 +835,7 @@ struct OnTheFlyRichPid { float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue}; bool signalBarrelRich[kNspecies] = {false, false, false, false, false}; float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies]; - static constexpr int lpdg_array[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static constexpr int lpdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; static constexpr float masses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron], o2::track::pid_constants::sMasses[o2::track::PID::Muon], o2::track::pid_constants::sMasses[o2::track::PID::Pion], @@ -845,21 +845,21 @@ struct OnTheFlyRichPid { for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses float hypothesisAngleBarrelRich = kErrorValue; - const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogel_rindex[i_sector], hypothesisAngleBarrelRich); + const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons // Evaluate total sigma (layer + tracking resolution) float barrelTotalAngularReso = barrelRICHAngularResolution; if (flagIncludeTrackAngularRes) { const float transverseMomentum = recoTrack.getP() / std::cosh(recoTrack.getEta()); - double pt_resolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2()); - double eta_resolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2()); + double ptResolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2()); + double etaResolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2()); if (flagRICHLoadDelphesLUTs) { - pt_resolution = mSmearer.getAbsPtRes(lpdg_array[ii], dNdEta, recoTrack.getEta(), transverseMomentum); - eta_resolution = mSmearer.getAbsEtaRes(lpdg_array[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + ptResolution = mSmearer.getAbsPtRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + etaResolution = mSmearer.getAbsEtaRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); } - // cout << endl << "Pt resolution: " << pt_resolution << ", Eta resolution: " << eta_resolution << endl << endl; - const float barrelTrackAngularReso = calculate_track_angular_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], aerogel_rindex[i_sector]); + // cout << endl << "Pt resolution: " << ptResolution << ", Eta resolution: " << etaResolution << endl << endl; + const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, masses[ii], aerogelRindex[iSecor]); barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso); if (doQAplots && hypothesisAngleBarrelRich > kErrorValue + 1. && @@ -867,36 +867,36 @@ struct OnTheFlyRichPid { barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { switch (mcParticle.pdgCode()) { - case lpdg_array[kEl]: // Electron - case -lpdg_array[kEl]: // Positron + case lpdgArray[kEl]: // Electron + case -lpdgArray[kEl]: // Positron if (ii == kEl) { histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[kMu]: // Muon - case -lpdg_array[kMu]: // AntiMuon + case lpdgArray[kMu]: // Muon + case -lpdgArray[kMu]: // AntiMuon if (ii == kMu) { histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[kPi]: // Pion - case -lpdg_array[kPi]: // AntiPion + case lpdgArray[kPi]: // Pion + case -lpdgArray[kPi]: // AntiPion if (ii == kPi) { histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[kKa]: // Kaon - case -lpdg_array[kKa]: // AntiKaon + case lpdgArray[kKa]: // Kaon + case -lpdgArray[kKa]: // AntiKaon if (ii == kKa) { histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdg_array[kPr]: // Proton - case -lpdg_array[kPr]: // AntiProton + case lpdgArray[kPr]: // Proton + case -lpdgArray[kPr]: // AntiProton if (ii == kPr) { histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); @@ -929,43 +929,43 @@ struct OnTheFlyRichPid { histos.fill(HIST("h2dAngleVsEtaBarrelRICH"), recoTrack.getEta(), measuredAngleBarrelRich); histos.fill(HIST("h2dAngularResolutionVsMomentumBarrelRICH"), recoTrack.getP(), 1000 * barrelRICHAngularResolution); histos.fill(HIST("h2dAngularResolutionVsEtaBarrelRICH"), recoTrack.getEta(), 1000 * barrelRICHAngularResolution); - histos.fill(HIST("hSectorID"), i_sector); + histos.fill(HIST("hSectorID"), iSecor); switch (mcParticle.pdgCode()) { - case lpdg_array[0]: // Electron - case -lpdg_array[0]: // Positron + case lpdgArray[0]: // Electron + case -lpdgArray[0]: // Positron histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdg_array[1]: // Muon - case -lpdg_array[1]: // AntiMuon + case lpdgArray[1]: // Muon + case -lpdgArray[1]: // AntiMuon histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdg_array[2]: // Pion - case -lpdg_array[2]: // AntiPion + case lpdgArray[2]: // Pion + case -lpdgArray[2]: // AntiPion histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdg_array[3]: // Kaon - case -lpdg_array[3]: // AntiKaon + case lpdgArray[3]: // Kaon + case -lpdgArray[3]: // AntiKaon histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdg_array[4]: // Proton - case -lpdg_array[4]: // AntiProton + case lpdgArray[4]: // Proton + case -lpdgArray[4]: // AntiProton histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); From aed822640871e9709cc95825b204cc2523b785e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 23:05:30 +0200 Subject: [PATCH 4/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 23 +++++++++++--------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index ee7fd93e881..86146b91f98 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -41,7 +41,10 @@ #include "CCDB/BasicCCDBManager.h" #include "CCDB/CcdbApi.h" #include "CommonConstants/GeomConstants.h" +<<<<<<< Updated upstream +======= #include "CommonConstants/MathConstants.h" +>>>>>>> Stashed changes #include "CommonConstants/PhysicsConstants.h" #include "CommonUtils/NameConf.h" #include "DataFormatsCalibration/MeanVertexObject.h" @@ -932,40 +935,40 @@ struct OnTheFlyRichPid { histos.fill(HIST("hSectorID"), iSecor); switch (mcParticle.pdgCode()) { - case lpdgArray[0]: // Electron - case -lpdgArray[0]: // Positron + case lpdgArray[kEl]: // Electron + case -lpdgArray[kEl]: // Positron histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[1]: // Muon - case -lpdgArray[1]: // AntiMuon + case lpdgArray[kMu]: // Muon + case -lpdgArray[kMu]: // AntiMuon histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[2]: // Pion - case -lpdgArray[2]: // AntiPion + case lpdgArray[kPi]: // Pion + case -lpdgArray[kPi]: // AntiPion histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[3]: // Kaon - case -lpdgArray[3]: // AntiKaon + case lpdgArray[kka]: // Kaon + case -lpdgArray[kka]: // AntiKaon histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[4]: // Proton - case -lpdgArray[4]: // AntiProton + case lpdgArray[kPr]: // Proton + case -lpdgArray[kPr]: // AntiProton histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); From ce21fdbfe432cce37fb07e79c2accdd6c8158d73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 23:10:44 +0200 Subject: [PATCH 5/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 31 ++++++++++---------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index 5ebf627b4fc..73b5ef3d643 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -41,6 +41,7 @@ #include "CCDB/BasicCCDBManager.h" #include "CCDB/CcdbApi.h" #include "CommonConstants/GeomConstants.h" +#include "CommonConstants/MathConstants.h" #include "CommonConstants/PhysicsConstants.h" #include "CommonUtils/NameConf.h" #include "DataFormatsCalibration/MeanVertexObject.h" @@ -506,7 +507,7 @@ struct OnTheFlyRichPid { /// \param n the refractive index of the considered sector /// \param angle the output angle (passed by reference) /// \return true if particle is above the threshold and enough photons, false otherwise - bool CherenkovAngle(const float momentum, const float mass, const float n, float& angle) + bool cherenkovAngle(const float momentum, const float mass, const float n, float& angle) { if (momentum > mass / std::sqrt(n * n - 1.0)) { // Check if particle is above the threshold // Calculate angle @@ -539,7 +540,7 @@ struct OnTheFlyRichPid { /// returns angular resolution for considered track eta /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin) - float AngularResolution(float eta) + float angularResolution(float eta) { // Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS) static constexpr float etaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; @@ -628,22 +629,22 @@ struct OnTheFlyRichPid { const float sinPhi = std::sin(phiC); const float cosPhi = std::cos(phiC); // const float ze = thicknessRad / (2. * zP); - const float az = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi; - const float e3z = std::sqrt(az * az + (nGas / n) * (nGas / n) - 1.); + const float aZ = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi; + const float e3Z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.); const float Z = thicknessGas; - const float alpha = e3z / az; - const float etac = e3z * n * n; + const float alpha = e3Z / aZ; + const float etac = e3Z * n * n; const float k = thicknessRad / (2. * Z); const float m = 1. / (n * n); const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); // Derivative d(thetaCherenkov)/dx const float temp1 = etac / Z; - const float temp2 = alpha * e3z * cosPhi; + const float temp2 = alpha * e3Z * cosPhi; const float temp3 = xP * sinThetaCherenkov * sinPhi * sinPhi; const float dThetaX = temp1 * (temp2 - temp3); // Derivative d(thetaCherenkov)/dy const float temp4 = etac * sinPhi / Z; - const float temp5 = cosThetaCherenkov - zP * (1 - m) / az; + const float temp5 = cosThetaCherenkov - zP * (1 - m) / aZ; const float dThetaY = temp4 * temp5; // Derivative d(thetaCherenkov)/dze const float temp8 = etac * sinThetaCherenkov; @@ -651,7 +652,7 @@ struct OnTheFlyRichPid { const float temp10 = alpha * alpha * (1.0 - xP * xP * sinPhi * sinPhi) + lambda * xP * xP * sinPhi * sinPhi; const float dThetaZe = temp8 * temp10 / temp9; // Derivative d(thetaCherenkov)/dn - const float dThetaN = zP / (n * az * sinThetaCherenkov); + const float dThetaN = zP / (n * aZ * sinThetaCherenkov); // RMS wavelength (using Sellmeier formula with measured aerogel parameters) const float a0 = 0.0616; const float w0 = 56.5; @@ -681,14 +682,14 @@ struct OnTheFlyRichPid { // Fraction of photons in rings (to account loss close to sector boundary in projective geometry) const float sinThetaCherenkovSquared = sinThetaCherenkov * sinThetaCherenkov; const float radius = (thicknessRad / 2.0) * std::tan(thetaCherenkov) + thicknessGas * std::tan(std::asin((n / nGas) * sinThetaCherenkov)); - const float N0 = 24. * thicknessRad / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) + const float n0 = 24. * thicknessRad / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) const float sinAngle = std::sin(std::acos(1. / 1.03)); const float sinAngleSquared = sinAngle * sinAngle; const float multiplicitySpectrumFactor = sinThetaCherenkovSquared / sinAngleSquared; // scale multiplicity w.r.t. N = 1.03 at saturation // Considering average resolution (integrated over the sector) - // float nPhotons = (tileZlength / 2.0 > radius) ? N0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength)) : N0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength) - (2.0/(tileZlength*o2::constants::math::PI))*(-(tileZlength/(2.0))*std::acos(tileZlength/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tileZlength/(2.0*radius),2.0)))); + // float nPhotons = (tileZlength / 2.0 > radius) ? n0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength)) : n0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength) - (2.0/(tileZlength*o2::constants::math::PI))*(-(tileZlength/(2.0))*std::acos(tileZlength/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tileZlength/(2.0*radius),2.0)))); // Considering "exact" resolution (eta by eta) - const float nPhotons = N0 * multiplicitySpectrumFactor * fractionPhotonsProjectiveRICH(eta, tileZlength, radius); + const float nPhotons = n0 * multiplicitySpectrumFactor * fractionPhotonsProjectiveRICH(eta, tileZlength, radius); if (nPhotons <= kErrorValue + 1) return kErrorValue; // Ring angular resolution @@ -797,12 +798,12 @@ struct OnTheFlyRichPid { } float expectedAngleBarrelRich = kErrorValue; - const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogelRindex[iSecor], expectedAngleBarrelRich); + const bool expectedAngleBarrelRichOk = cherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogelRindex[iSecor], expectedAngleBarrelRich); if (!expectedAngleBarrelRichOk) { fillDummyValues(); continue; // Particle is below the threshold or not enough photons } - // float barrelRICHAngularResolution = AngularResolution(o2track.getEta()); + // float barrelRICHAngularResolution = angularResolution(o2track.getEta()); const float barrelRICHAngularResolution = extractRingAngularResolution(o2track.getEta(), aerogelRindex[iSecor], bRichGapRefractiveIndex, bRichRadiatorThickness, gapThickness[iSecor], bRICHPixelSize, expectedAngleBarrelRich, photodetrctorLength[iSecor]); const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), iSecor); bool flagReachesRadiator = false; @@ -844,7 +845,7 @@ struct OnTheFlyRichPid { for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses float hypothesisAngleBarrelRich = kErrorValue; - const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); + const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), masses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons // Evaluate total sigma (layer + tracking resolution) From f5ddf80187aa72c5f7ba2f80ffd1e33eefd9d30e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 23:14:40 +0200 Subject: [PATCH 6/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 103 ++++++++----------- 1 file changed, 41 insertions(+), 62 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index 73b5ef3d643..67d3721791f 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -105,27 +105,6 @@ struct OnTheFlyRichPid { Configurable flagIncludeTrackAngularRes{"flagIncludeTrackAngularRes", true, "flag to include or exclude track time resolution"}; Configurable multiplicityEtaRange{"multiplicityEtaRange", 0.800000012, "eta range to compute the multiplicity"}; Configurable flagRICHLoadDelphesLUTs{"flagRICHLoadDelphesLUTs", false, "flag to load Delphes LUTs for tracking correction (use recoTrack parameters if false)"}; - /*Configurable bRichRefractiveIndexSector0AndNMinus1{"bRichRefractiveIndexSector0AndMinus1", 1.03, "barrel RICH refractive index sector 0 and N-1"}; - Configurable bRichRefractiveIndexSector1AndNMinus2{"bRichRefractiveIndexSector1AndMinus2", 1.03, "barrel RICH refractive index sector 1 and N-2"}; - Configurable bRichRefractiveIndexSector2AndNMinus3{"bRichRefractiveIndexSector2AndMinus3", 1.03, "barrel RICH refractive index sector 2 and N-3"}; - Configurable bRichRefractiveIndexSector3AndNMinus4{"bRichRefractiveIndexSector3AndMinus4", 1.03, "barrel RICH refractive index sector 3 and N-4"}; - Configurable bRichRefractiveIndexSector4AndNMinus5{"bRichRefractiveIndexSector4AndMinus5", 1.03, "barrel RICH refractive index sector 4 and N-5"}; - Configurable bRichRefractiveIndexSector5AndNMinus6{"bRichRefractiveIndexSector5AndMinus6", 1.03, "barrel RICH refractive index sector 5 and N-6"}; - Configurable bRichRefractiveIndexSector6AndNMinus7{"bRichRefractiveIndexSector6AndMinus7", 1.03, "barrel RICH refractive index sector 6 and N-7"}; - Configurable bRichRefractiveIndexSector7AndNMinus8{"bRichRefractiveIndexSector7AndMinus8", 1.03, "barrel RICH refractive index sector 7 and N-8"}; - Configurable bRichRefractiveIndexSector8AndNMinus9{"bRichRefractiveIndexSector8AndMinus9", 1.03, "barrel RICH refractive index sector 8 and N-9"}; - Configurable bRichRefractiveIndexSector9AndNMinus10{"bRichRefractiveIndexSector9AndMinus10", 1.03, "barrel RICH refractive index sector 9 and N-10"}; - Configurable bRichRefractiveIndexSector10AndNMinus11{"bRichRefractiveIndexSector10AndMinus11", 1.03, "barrel RICH refractive index sector 10 and N-11"}; - Configurable bRichRefractiveIndexSector11AndNMinus12{"bRichRefractiveIndexSector11AndMinus12", 1.03, "barrel RICH refractive index sector 11 and N-12"}; - Configurable bRichRefractiveIndexSector12AndNMinus13{"bRichRefractiveIndexSector12AndMinus13", 1.03, "barrel RICH refractive index sector 12 and N-13"}; - Configurable bRichRefractiveIndexSector13AndNMinus14{"bRichRefractiveIndexSector13AndMinus14", 1.03, "barrel RICH refractive index sector 13 and N-14"}; - Configurable bRichRefractiveIndexSector14AndNMinus15{"bRichRefractiveIndexSector14AndMinus15", 1.03, "barrel RICH refractive index sector 14 and N-15"}; - Configurable bRichRefractiveIndexSector15AndNMinus16{"bRichRefractiveIndexSector15AndMinus16", 1.03, "barrel RICH refractive index sector 15 and N-16"}; - Configurable bRichRefractiveIndexSector16AndNMinus17{"bRichRefractiveIndexSector16AndMinus17", 1.03, "barrel RICH refractive index sector 16 and N-17"}; - Configurable bRichRefractiveIndexSector17AndNMinus18{"bRichRefractiveIndexSector17AndMinus18", 1.03, "barrel RICH refractive index sector 17 and N-18"}; - Configurable bRichRefractiveIndexSector18AndNMinus19{"bRichRefractiveIndexSector18AndMinus19", 1.03, "barrel RICH refractive index sector 18 and N-19"}; - Configurable bRichRefractiveIndexSector19AndNMinus20{"bRichRefractiveIndexSector19AndMinus20", 1.03, "barrel RICH refractive index sector 19 and N-20"}; - Configurable bRichRefractiveIndexSector20AndNMinus21{"bRichRefractiveIndexSector20AndMinus21", 1.03, "barrel RICH refractive index sector 20 and N-21"};*/ Configurable bRichRefractiveIndexSector0{"bRichRefractiveIndexSector0", 1.03, "barrel RICH refractive index central(s)"}; // central(s) Configurable bRichRefractiveIndexSector1{"bRichRefractiveIndexSector1", 1.03, "barrel RICH refractive index central(s)-1 and central(s)+1"}; // central(s)-1 and central(s)+1 Configurable bRichRefractiveIndexSector2{"bRichRefractiveIndexSector2", 1.03, "barrel RICH refractive index central(s)-2 and central(s)+2"}; // central(s)-2 and central(s)+2 @@ -543,21 +522,21 @@ struct OnTheFlyRichPid { float angularResolution(float eta) { // Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS) - static constexpr float etaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; - static constexpr float resRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; - static constexpr float resRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; - static constexpr int sizeResVector = sizeof(etaSampling) / sizeof(etaSampling[0]); + static constexpr float kEtaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; + static constexpr float kResRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; + static constexpr float kResRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; + static constexpr int sizeResVector = sizeof(kEtaSampling) / sizeof(kEtaSampling[0]); // Use binary search to find the lower and upper indices - const int lowerIndex = std::lower_bound(etaSampling, etaSampling + sizeResVector, eta) - etaSampling - 1; + const int lowerIndex = std::lower_bound(kEtaSampling, kEtaSampling + sizeResVector, eta) - kEtaSampling - 1; const int upperIndex = lowerIndex + 1; if (lowerIndex >= 0 && upperIndex < sizeResVector) { // Resolution interpolation if (bRichFlagAbsorbingWalls) { - float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithAbsWalls[lowerIndex], resRingSamplingWithAbsWalls[upperIndex]); + float interpolatedResRing = interpolate(eta, kEtaSampling[lowerIndex], kEtaSampling[upperIndex], kResRingSamplingWithAbsWalls[lowerIndex], kResRingSamplingWithAbsWalls[upperIndex]); // std::cout << "Interpolated y value: " << interpolatedY << std::endl; return interpolatedResRing; } else { - float interpolatedResRing = interpolate(eta, etaSampling[lowerIndex], etaSampling[upperIndex], resRingSamplingWithoutAbsWalls[lowerIndex], resRingSamplingWithoutAbsWalls[upperIndex]); + float interpolatedResRing = interpolate(eta, kEtaSampling[lowerIndex], kEtaSampling[upperIndex], kResRingSamplingWithoutAbsWalls[lowerIndex], kResRingSamplingWithoutAbsWalls[upperIndex]); // std::cout << "Interpolated y value: " << interpolatedY << std::endl; return interpolatedResRing; } @@ -630,16 +609,16 @@ struct OnTheFlyRichPid { const float cosPhi = std::cos(phiC); // const float ze = thicknessRad / (2. * zP); const float aZ = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi; - const float e3Z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.); + const float e3z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.); const float Z = thicknessGas; - const float alpha = e3Z / aZ; - const float etac = e3Z * n * n; + const float alpha = e3z / aZ; + const float etac = e3z * n * n; const float k = thicknessRad / (2. * Z); const float m = 1. / (n * n); const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); // Derivative d(thetaCherenkov)/dx const float temp1 = etac / Z; - const float temp2 = alpha * e3Z * cosPhi; + const float temp2 = alpha * e3z * cosPhi; const float temp3 = xP * sinThetaCherenkov * sinPhi * sinPhi; const float dThetaX = temp1 * (temp2 - temp3); // Derivative d(thetaCherenkov)/dy @@ -835,17 +814,17 @@ struct OnTheFlyRichPid { float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue}; bool signalBarrelRich[kNspecies] = {false, false, false, false, false}; float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies]; - static constexpr int lpdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; - static constexpr float masses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron], - o2::track::pid_constants::sMasses[o2::track::PID::Muon], - o2::track::pid_constants::sMasses[o2::track::PID::Pion], - o2::track::pid_constants::sMasses[o2::track::PID::Kaon], - o2::track::pid_constants::sMasses[o2::track::PID::Proton]}; + static constexpr int kPdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static constexpr float kMasses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron], + o2::track::pid_constants::sMasses[o2::track::PID::Muon], + o2::track::pid_constants::sMasses[o2::track::PID::Pion], + o2::track::pid_constants::sMasses[o2::track::PID::Kaon], + o2::track::pid_constants::sMasses[o2::track::PID::Proton]}; for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses float hypothesisAngleBarrelRich = kErrorValue; - const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), masses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); + const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), kMasses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons // Evaluate total sigma (layer + tracking resolution) @@ -855,11 +834,11 @@ struct OnTheFlyRichPid { double ptResolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2()); double etaResolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2()); if (flagRICHLoadDelphesLUTs) { - ptResolution = mSmearer.getAbsPtRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); - etaResolution = mSmearer.getAbsEtaRes(lpdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + ptResolution = mSmearer.getAbsPtRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + etaResolution = mSmearer.getAbsEtaRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); } // cout << endl << "Pt resolution: " << ptResolution << ", Eta resolution: " << etaResolution << endl << endl; - const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, masses[ii], aerogelRindex[iSecor]); + const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, kMasses[ii], aerogelRindex[iSecor]); barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso); if (doQAplots && hypothesisAngleBarrelRich > kErrorValue + 1. && @@ -867,36 +846,36 @@ struct OnTheFlyRichPid { barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { switch (mcParticle.pdgCode()) { - case lpdgArray[kEl]: // Electron - case -lpdgArray[kEl]: // Positron + case kPdgArray[kEl]: // Electron + case -kPdgArray[kEl]: // Positron if (ii == kEl) { histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdgArray[kMu]: // Muon - case -lpdgArray[kMu]: // AntiMuon + case kPdgArray[kMu]: // Muon + case -kPdgArray[kMu]: // AntiMuon if (ii == kMu) { histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdgArray[kPi]: // Pion - case -lpdgArray[kPi]: // AntiPion + case kPdgArray[kPi]: // Pion + case -kPdgArray[kPi]: // AntiPion if (ii == kPi) { histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdgArray[kKa]: // Kaon - case -lpdgArray[kKa]: // AntiKaon + case kPdgArray[kKa]: // Kaon + case -kPdgArray[kKa]: // AntiKaon if (ii == kKa) { histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case lpdgArray[kPr]: // Proton - case -lpdgArray[kPr]: // AntiProton + case kPdgArray[kPr]: // Proton + case -kPdgArray[kPr]: // AntiProton if (ii == kPr) { histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); @@ -932,40 +911,40 @@ struct OnTheFlyRichPid { histos.fill(HIST("hSectorID"), iSecor); switch (mcParticle.pdgCode()) { - case lpdgArray[kEl]: // Electron - case -lpdgArray[kEl]: // Positron + case kPdgArray[kEl]: // Electron + case -kPdgArray[kEl]: // Positron histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[kMu]: // Muon - case -lpdgArray[kMu]: // AntiMuon + case kPdgArray[kMu]: // Muon + case -kPdgArray[kMu]: // AntiMuon histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[kPi]: // Pion - case -lpdgArray[kPi]: // AntiPion + case kPdgArray[kPi]: // Pion + case -kPdgArray[kPi]: // AntiPion histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[kka]: // Kaon - case -lpdgArray[kka]: // AntiKaon + case kPdgArray[kka]: // Kaon + case -kPdgArray[kka]: // AntiKaon histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case lpdgArray[kPr]: // Proton - case -lpdgArray[kPr]: // AntiProton + case kPdgArray[kPr]: // Proton + case -kPdgArray[kPr]: // AntiProton histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); From 72e92a97b6f4ffb8612d26ad37429dd401559802 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sun, 13 Jul 2025 23:17:38 +0200 Subject: [PATCH 7/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index 67d3721791f..afbc5b4fab1 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -232,7 +232,7 @@ struct OnTheFlyRichPid { lAerogelZ[2 * iCentralMirror - i] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); tRPlusG[2 * iCentralMirror - i] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); aerogelRindex[2 * iCentralMirror - i] = bRichRefractiveIndexSector[i - iCentralMirror]; - bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum Z + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum z } } else { // Even number of sectors float twoHalfGap = 1.0; @@ -265,7 +265,7 @@ struct OnTheFlyRichPid { lAerogelZ[2 * iCentralMirror - i - 1] = std::sqrt(1.0 + mVal * mVal) * rMin * squareSizeZ / (std::sqrt(1.0 + mVal * mVal) * rMax - mVal * squareSizeZ); tRPlusG[2 * iCentralMirror - i - 1] = std::sqrt(1.0 + mVal * mVal) * (rMax - rMin) - mVal / 2.0 * (squareSizeZ + lAerogelZ[i]); aerogelRindex[2 * iCentralMirror - i - 1] = bRichRefractiveIndexSector[i - iCentralMirror]; - bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum Z + bRichPhotodetectorCentralModuleHalfLength.value = rMin * t; // <-- At the end of the loop this will be the maximum z } } // Coordinate radiali layer considerati @@ -525,11 +525,11 @@ struct OnTheFlyRichPid { static constexpr float kEtaSampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000}; static constexpr float kResRingSamplingWithAbsWalls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165}; static constexpr float kResRingSamplingWithoutAbsWalls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165}; - static constexpr int sizeResVector = sizeof(kEtaSampling) / sizeof(kEtaSampling[0]); + static constexpr int kSizeResVector = sizeof(kEtaSampling) / sizeof(kEtaSampling[0]); // Use binary search to find the lower and upper indices - const int lowerIndex = std::lower_bound(kEtaSampling, kEtaSampling + sizeResVector, eta) - kEtaSampling - 1; + const int lowerIndex = std::lower_bound(kEtaSampling, kEtaSampling + kSizeResVector, eta) - kEtaSampling - 1; const int upperIndex = lowerIndex + 1; - if (lowerIndex >= 0 && upperIndex < sizeResVector) { + if (lowerIndex >= 0 && upperIndex < kSizeResVector) { // Resolution interpolation if (bRichFlagAbsorbingWalls) { float interpolatedResRing = interpolate(eta, kEtaSampling[lowerIndex], kEtaSampling[upperIndex], kResRingSamplingWithAbsWalls[lowerIndex], kResRingSamplingWithAbsWalls[upperIndex]); @@ -610,24 +610,24 @@ struct OnTheFlyRichPid { // const float ze = thicknessRad / (2. * zP); const float aZ = zP * cosThetaCherenkov - xP * sinThetaCherenkov * cosPhi; const float e3z = std::sqrt(aZ * aZ + (nGas / n) * (nGas / n) - 1.); - const float Z = thicknessGas; + const float z = thicknessGas; const float alpha = e3z / aZ; const float etac = e3z * n * n; - const float k = thicknessRad / (2. * Z); + const float k = thicknessRad / (2. * z); const float m = 1. / (n * n); const float lambda = (1. + k * alpha * alpha * alpha) / (1. + k * alpha); // Derivative d(thetaCherenkov)/dx - const float temp1 = etac / Z; + const float temp1 = etac / z; const float temp2 = alpha * e3z * cosPhi; const float temp3 = xP * sinThetaCherenkov * sinPhi * sinPhi; const float dThetaX = temp1 * (temp2 - temp3); // Derivative d(thetaCherenkov)/dy - const float temp4 = etac * sinPhi / Z; + const float temp4 = etac * sinPhi / z; const float temp5 = cosThetaCherenkov - zP * (1 - m) / aZ; const float dThetaY = temp4 * temp5; // Derivative d(thetaCherenkov)/dze const float temp8 = etac * sinThetaCherenkov; - const float temp9 = Z * (1.0 + k * alpha * alpha * alpha * n * n); + const float temp9 = z * (1.0 + k * alpha * alpha * alpha * n * n); const float temp10 = alpha * alpha * (1.0 - xP * xP * sinPhi * sinPhi) + lambda * xP * xP * sinPhi * sinPhi; const float dThetaZe = temp8 * temp10 / temp9; // Derivative d(thetaCherenkov)/dn From 28fdd551113babc460df75b133d16fe9f538f3aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Mon, 14 Jul 2025 06:13:00 +0200 Subject: [PATCH 8/8] Update onTheFlyRichPid.cxx --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index afbc5b4fab1..8900d87b331 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -661,14 +661,14 @@ struct OnTheFlyRichPid { // Fraction of photons in rings (to account loss close to sector boundary in projective geometry) const float sinThetaCherenkovSquared = sinThetaCherenkov * sinThetaCherenkov; const float radius = (thicknessRad / 2.0) * std::tan(thetaCherenkov) + thicknessGas * std::tan(std::asin((n / nGas) * sinThetaCherenkov)); - const float n0 = 24. * thicknessRad / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) + const float n0Photons = 24. * thicknessRad / 2.; // photons for N = 1.03 at saturation ( 24/2 factor per radiator cm ) const float sinAngle = std::sin(std::acos(1. / 1.03)); const float sinAngleSquared = sinAngle * sinAngle; const float multiplicitySpectrumFactor = sinThetaCherenkovSquared / sinAngleSquared; // scale multiplicity w.r.t. N = 1.03 at saturation // Considering average resolution (integrated over the sector) - // float nPhotons = (tileZlength / 2.0 > radius) ? n0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength)) : n0 * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength) - (2.0/(tileZlength*o2::constants::math::PI))*(-(tileZlength/(2.0))*std::acos(tileZlength/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tileZlength/(2.0*radius),2.0)))); + // float nPhotons = (tileZlength / 2.0 > radius) ? n0Photons * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength)) : n0Photons * multiplicitySpectrumFactor * (1.-(2.0*radius)/(o2::constants::math::PI*tileZlength) - (2.0/(tileZlength*o2::constants::math::PI))*(-(tileZlength/(2.0))*std::acos(tileZlength/(2.0*radius)) + radius*std::sqrt(1.-std::pow(tileZlength/(2.0*radius),2.0)))); // Considering "exact" resolution (eta by eta) - const float nPhotons = n0 * multiplicitySpectrumFactor * fractionPhotonsProjectiveRICH(eta, tileZlength, radius); + const float nPhotons = n0Photons * multiplicitySpectrumFactor * fractionPhotonsProjectiveRICH(eta, tileZlength, radius); if (nPhotons <= kErrorValue + 1) return kErrorValue; // Ring angular resolution @@ -935,8 +935,8 @@ struct OnTheFlyRichPid { histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case kPdgArray[kka]: // Kaon - case -kPdgArray[kka]: // AntiKaon + case kPdgArray[kKa]: // Kaon + case -kPdgArray[kKa]: // AntiKaon histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]);