Skip to content

Commit 656233c

Browse files
njacazionjacazio
authored
Update onTheFlyRichPid.cxx
1 parent bddbaea commit 656233c

File tree

1 file changed

+58
-63
lines changed

1 file changed

+58
-63
lines changed

‎ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx‎

Lines changed: 58 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -170,15 +170,9 @@ struct OnTheFlyRichPid {
170170
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
171171

172172
/// Flag unphysical and unavailable values (must be negative)
173-
float error_value = -1000;
173+
static constexpr float kErrorValue = -1000;
174174

175175
// Variables projective/hybrid layout
176-
int mNumberSectors = bRichNumberOfSectors; // 21;
177-
float mTileLength = bRichPhotodetectorOtherModuleLength; // 18.4; // [cm]
178-
float mTileLengthCentral = bRichPhotodetectorCentralModuleHalfLength; // 18.4 / 2.0; // [cm]
179-
float mProjectiveLengthInner = mTileLengthCentral;
180-
float mRadiusProjIn = bRichRadiatorInnerRadius; // 85.; // [cm]
181-
float mRadiusProjOut = bRichPhotodetectorOuterRadius; // 112.; // [cm]
182176
std::vector<TVector3> det_centers;
183177
std::vector<TVector3> rad_centers;
184178
std::vector<float> angle_centers;
@@ -192,13 +186,7 @@ struct OnTheFlyRichPid {
192186
// Update projective geometry
193187
void updateProjectiveParameters()
194188
{
195-
mNumberSectors = bRichNumberOfSectors;
196-
mTileLength = bRichPhotodetectorOtherModuleLength;
197-
mTileLengthCentral = bRichPhotodetectorCentralModuleHalfLength;
198-
mProjectiveLengthInner = mTileLengthCentral;
199-
mRadiusProjIn = bRichRadiatorInnerRadius;
200-
mRadiusProjOut = bRichPhotodetectorOuterRadius;
201-
const int number_of_sectors_in_z = mNumberSectors;
189+
const int number_of_sectors_in_z = bRichNumberOfSectors;
202190
det_centers.resize(number_of_sectors_in_z);
203191
rad_centers.resize(number_of_sectors_in_z);
204192
angle_centers.resize(number_of_sectors_in_z);
@@ -207,10 +195,10 @@ struct OnTheFlyRichPid {
207195
aerogel_rindex.resize(number_of_sectors_in_z);
208196
photodetrctor_length.resize(number_of_sectors_in_z);
209197
gap_thickness.resize(number_of_sectors_in_z);
210-
float square_size_barrel_cylinder = 2.0 * mTileLengthCentral;
211-
float square_size_z = mTileLength;
212-
float R_min = mRadiusProjIn;
213-
float R_max = mRadiusProjOut;
198+
float square_size_barrel_cylinder = 2.0 * bRichPhotodetectorCentralModuleHalfLength;
199+
float square_size_z = bRichPhotodetectorOtherModuleLength;
200+
float R_min = bRichRadiatorInnerRadius;
201+
float R_max = bRichPhotodetectorOuterRadius;
214202
std::vector<float> theta_bi;
215203
std::vector<float> R0_tilt;
216204
std::vector<float> z0_tilt;
@@ -225,7 +213,8 @@ struct OnTheFlyRichPid {
225213
l_detector_z.resize(number_of_sectors_in_z);
226214

227215
// Odd number of sectors
228-
if (number_of_sectors_in_z % 2 != 0) {
216+
static constexpr int kTwo = 2;
217+
if (number_of_sectors_in_z % kTwo != 0) {
229218
int i_central_mirror = static_cast<int>((number_of_sectors_in_z) / 2.0);
230219
float m_val = std::tan(0.0);
231220
theta_bi[i_central_mirror] = std::atan(m_val);
@@ -237,7 +226,7 @@ struct OnTheFlyRichPid {
237226
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)));
238227
theta_max[i_central_mirror] = o2::constants::math::PIHalf - std::atan(t);
239228
theta_min[i_central_mirror] = o2::constants::math::PIHalf + std::atan(t);
240-
mProjectiveLengthInner = R_min * t;
229+
bRichPhotodetectorCentralModuleHalfLength = R_min * t;
241230
aerogel_rindex[i_central_mirror] = bRichRefractiveIndexSector[0];
242231
for (int i = i_central_mirror + 1; i < number_of_sectors_in_z; i++) {
243232
float par_a = t;
@@ -264,14 +253,14 @@ struct OnTheFlyRichPid {
264253
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);
265254
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]);
266255
aerogel_rindex[2 * i_central_mirror - i] = bRichRefractiveIndexSector[i - i_central_mirror];
267-
mProjectiveLengthInner = R_min * t; // <-- At the end of the loop this will be the maximum Z
256+
bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z
268257
}
269258
} else { // Even number of sectors
270259
float two_half_gap = 1.0;
271260
int i_central_mirror = static_cast<int>((number_of_sectors_in_z) / 2.0);
272261
float m_val = std::tan(0.0);
273262
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)));
274-
mProjectiveLengthInner = R_min * t;
263+
bRichPhotodetectorCentralModuleHalfLength = R_min * t;
275264
for (int i = i_central_mirror; i < number_of_sectors_in_z; i++) {
276265
float par_a = t;
277266
float par_b = 2.0 * R_max / square_size_z;
@@ -297,7 +286,7 @@ struct OnTheFlyRichPid {
297286
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);
298287
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]);
299288
aerogel_rindex[2 * i_central_mirror - i - 1] = bRichRefractiveIndexSector[i - i_central_mirror];
300-
mProjectiveLengthInner = R_min * t; // <-- At the end of the loop this will be the maximum Z
289+
bRichPhotodetectorCentralModuleHalfLength = R_min * t; // <-- At the end of the loop this will be the maximum Z
301290
}
302291
}
303292
// Coordinate radiali layer considerati
@@ -326,7 +315,7 @@ struct OnTheFlyRichPid {
326315
{
327316
pRandomNumberGenerator.SetSeed(0); // fully randomize
328317

329-
if (magneticField.value < 0.0001f) {
318+
if (magneticField.value < o2::constants::math::Epsilon) {
330319
LOG(info) << "Getting the magnetic field from the on-the-fly tracker task";
331320
if (!getTaskOptionValue(initContext, "on-the-fly-tracker", magneticField, false)) {
332321
LOG(fatal) << "Could not get Bz from on-the-fly-tracker task";
@@ -487,7 +476,7 @@ struct OnTheFlyRichPid {
487476
int findSector(const float eta)
488477
{
489478
float polar = 2.0 * std::atan(std::exp(-eta));
490-
for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) {
479+
for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) {
491480
if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) {
492481
return j_sec;
493482
}
@@ -508,7 +497,7 @@ struct OnTheFlyRichPid {
508497
const float z_sec_rich_squared = z_sec_rich * z_sec_rich;
509498
return (R_sec_rich_squared + z_sec_rich_squared) / (R_sec_rich + z_sec_rich / std::tan(polar));
510499
} else {
511-
return error_value;
500+
return kErrorValue;
512501
}
513502
}
514503

@@ -528,7 +517,8 @@ struct OnTheFlyRichPid {
528517
const float meanNumberofDetectedPhotons = 230. * std::sin(angle) * std::sin(angle) * bRichRadiatorThickness;
529518

530519
// Require at least 3 photons on average for real angle reconstruction
531-
if (meanNumberofDetectedPhotons > 3.) {
520+
static constexpr float kMinPhotons = 3.f;
521+
if (meanNumberofDetectedPhotons > kMinPhotons) {
532522
return true; // Particle is above the threshold and enough photons
533523
}
534524
return false; // Particle is above the threshold, but not enough photons
@@ -573,7 +563,7 @@ struct OnTheFlyRichPid {
573563
}
574564
} else {
575565
// std::cout << "Unable to interpolate. Target x value is outside the range of available data." << std::endl;
576-
return error_value;
566+
return kErrorValue;
577567
}
578568
}
579569

@@ -587,14 +577,14 @@ struct OnTheFlyRichPid {
587577
const float polar = 2.0 * std::atan(std::exp(-eta));
588578
int i_sector = 0;
589579
bool flag_sector = false;
590-
for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) {
580+
for (int j_sec = 0; j_sec < bRichNumberOfSectors; j_sec++) {
591581
if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) {
592582
flag_sector = true;
593583
i_sector = j_sec;
594584
}
595585
}
596586
if (flag_sector == false) {
597-
return error_value; // <-- Returning negative value
587+
return kErrorValue; // <-- Returning negative value
598588
}
599589
float R_sec_rich = rad_centers[i_sector].X();
600590
float z_sec_rich = rad_centers[i_sector].Z();
@@ -627,8 +617,8 @@ struct OnTheFlyRichPid {
627617
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)
628618
{
629619
// Check if input angle is error value
630-
if (theta_c <= error_value + 1)
631-
return error_value;
620+
if (theta_c <= kErrorValue + 1)
621+
return kErrorValue;
632622
// Parametrization variables (notation from https://doi.org/10.1016/0168-9002(94)90532-0)
633623
const float phi_c = 0.;
634624
const float theta_p = 0.;
@@ -700,8 +690,8 @@ struct OnTheFlyRichPid {
700690
// 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))));
701691
// Considering "exact" resolution (eta by eta)
702692
float n_photons = N0 * multiplicity_spectrum_factor * fractionPhotonsProjectiveRICH(eta, tile_z_length, radius);
703-
if (n_photons <= error_value + 1)
704-
return error_value;
693+
if (n_photons <= kErrorValue + 1)
694+
return kErrorValue;
705695
// Ring angular resolution
706696
float ring_angular_resolution = single_photon_angular_resolution / std::sqrt(n_photons);
707697
return ring_angular_resolution;
@@ -774,7 +764,7 @@ struct OnTheFlyRichPid {
774764
for (const auto& track : tracks) {
775765

776766
auto fillDummyValues = [&]() {
777-
upgradeRich(error_value, error_value, error_value, error_value, error_value);
767+
upgradeRich(kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue);
778768
upgradeRichSignal(false, false, false, false, false, false);
779769
};
780770

@@ -788,7 +778,7 @@ struct OnTheFlyRichPid {
788778
const auto& mcParticle = track.mcParticle();
789779
o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg);
790780

791-
// float xPv = error_value;
781+
// float xPv = kErrorValue;
792782
if (o2track.propagateToDCA(mcPvVtx, magneticField)) {
793783
// xPv = o2track.getX();
794784
}
@@ -807,7 +797,7 @@ struct OnTheFlyRichPid {
807797
continue;
808798
}
809799

810-
float expectedAngleBarrelRich = error_value;
800+
float expectedAngleBarrelRich = kErrorValue;
811801
const bool expectedAngleBarrelRichOk = CherenkovAngle(o2track.getP(), pdgInfo->Mass(), aerogel_rindex[i_sector], expectedAngleBarrelRich);
812802
if (!expectedAngleBarrelRichOk) {
813803
fillDummyValues();
@@ -817,7 +807,7 @@ struct OnTheFlyRichPid {
817807
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]);
818808
const float projectiveRadiatorRadius = radiusRipple(o2track.getEta(), i_sector);
819809
bool flagReachesRadiator = false;
820-
if (projectiveRadiatorRadius > error_value + 1.) {
810+
if (projectiveRadiatorRadius > kErrorValue + 1.) {
821811
flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, magneticField);
822812
}
823813
/// DISCLAIMER: Exact extrapolation of angular resolution would require track propagation
@@ -837,7 +827,12 @@ struct OnTheFlyRichPid {
837827

838828
// Straight to Nsigma
839829
static constexpr int kNspecies = 5;
840-
float nSigmaBarrelRich[kNspecies] = {error_value, error_value, error_value, error_value, error_value};
830+
static constexpr int kEl = 0;
831+
static constexpr int kMu = 1;
832+
static constexpr int kPi = 2;
833+
static constexpr int kKa = 3;
834+
static constexpr int kPr = 4;
835+
float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue};
841836
bool signalBarrelRich[kNspecies] = {false, false, false, false, false};
842837
float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies];
843838
static constexpr int lpdg_array[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
@@ -849,7 +844,7 @@ struct OnTheFlyRichPid {
849844

850845
for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses
851846

852-
float hypothesisAngleBarrelRich = error_value;
847+
float hypothesisAngleBarrelRich = kErrorValue;
853848
const bool hypothesisAngleBarrelRichOk = CherenkovAngle(recoTrack.getP(), masses[ii], aerogel_rindex[i_sector], hypothesisAngleBarrelRich);
854849
signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons
855850

@@ -867,42 +862,42 @@ struct OnTheFlyRichPid {
867862
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]);
868863
barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso);
869864
if (doQAplots &&
870-
hypothesisAngleBarrelRich > error_value + 1. &&
871-
measuredAngleBarrelRich > error_value + 1. &&
872-
barrelRICHAngularResolution > error_value + 1. &&
865+
hypothesisAngleBarrelRich > kErrorValue + 1. &&
866+
measuredAngleBarrelRich > kErrorValue + 1. &&
867+
barrelRICHAngularResolution > kErrorValue + 1. &&
873868
flagReachesRadiator) {
874869
switch (mcParticle.pdgCode()) {
875-
case lpdg_array[0]: // Electron
876-
case -lpdg_array[0]: // Positron
877-
if (ii == 0) {
870+
case lpdg_array[kEl]: // Electron
871+
case -lpdg_array[kEl]: // Positron
872+
if (ii == kEl) {
878873
histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
879874
histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
880875
}
881876
break;
882-
case lpdg_array[1]: // Muon
883-
case -lpdg_array[1]: // AntiMuon
884-
if (ii == 1) {
877+
case lpdg_array[kMu]: // Muon
878+
case -lpdg_array[kMu]: // AntiMuon
879+
if (ii == kMu) {
885880
histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
886881
histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
887882
}
888883
break;
889-
case lpdg_array[2]: // Pion
890-
case -lpdg_array[2]: // AntiPion
891-
if (ii == 2) {
884+
case lpdg_array[kPi]: // Pion
885+
case -lpdg_array[kPi]: // AntiPion
886+
if (ii == kPi) {
892887
histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
893888
histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
894889
}
895890
break;
896-
case lpdg_array[3]: // Kaon
897-
case -lpdg_array[3]: // AntiKaon
898-
if (ii == 3) {
891+
case lpdg_array[kKa]: // Kaon
892+
case -lpdg_array[kKa]: // AntiKaon
893+
if (ii == kKa) {
899894
histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
900895
histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
901896
}
902897
break;
903-
case lpdg_array[4]: // Proton
904-
case -lpdg_array[4]: // AntiProton
905-
if (ii == 4) {
898+
case lpdg_array[kPr]: // Proton
899+
case -lpdg_array[kPr]: // AntiProton
900+
if (ii == kPr) {
906901
histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso);
907902
histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso);
908903
}
@@ -916,9 +911,9 @@ struct OnTheFlyRichPid {
916911
/// DISCLAIMER: here tracking is accounted only for momentum value, but not for track parameters at impact point on the
917912
/// RICH radiator, since exact resolution would require photon generation and transport to photodetector.
918913
/// Effects are expected to be negligible (a few tenths of a milliradian) but further studies are required !
919-
if (hypothesisAngleBarrelRich > error_value + 1. &&
920-
measuredAngleBarrelRich > error_value + 1. &&
921-
barrelRICHAngularResolution > error_value + 1. &&
914+
if (hypothesisAngleBarrelRich > kErrorValue + 1. &&
915+
measuredAngleBarrelRich > kErrorValue + 1. &&
916+
barrelRICHAngularResolution > kErrorValue + 1. &&
922917
flagReachesRadiator) {
923918
deltaThetaBarrelRich[ii] = hypothesisAngleBarrelRich - measuredAngleBarrelRich;
924919
nSigmaBarrelRich[ii] = deltaThetaBarrelRich[ii] / barrelTotalAngularReso;
@@ -927,8 +922,8 @@ struct OnTheFlyRichPid {
927922

928923
// Fill histograms
929924
if (doQAplots) {
930-
if (measuredAngleBarrelRich > error_value + 1. &&
931-
barrelRICHAngularResolution > error_value + 1. &&
925+
if (measuredAngleBarrelRich > kErrorValue + 1. &&
926+
barrelRICHAngularResolution > kErrorValue + 1. &&
932927
flagReachesRadiator) {
933928
histos.fill(HIST("h2dAngleVsMomentumBarrelRICH"), recoTrack.getP(), measuredAngleBarrelRich);
934929
histos.fill(HIST("h2dAngleVsEtaBarrelRICH"), recoTrack.getEta(), measuredAngleBarrelRich);

0 commit comments

Comments
 (0)