@@ -80,73 +80,8 @@ enum AssociatedParticleType {
8080 kAssocPartSize
8181};
8282
83- /* enum ParticleOfInterest {
84- Phi = 0,
85- K0S,
86- Pion,
87- ParticleOfInterestSize
88- };*/
89-
90- /*
91- #define LIST_OF_PARTICLES_OF_INTEREST \
92- X(Phi) \
93- X(K0S) \
94- X(Pion) \
95- //X(PionTPC) \
96- //X(PionTPCTOF)
97-
98- enum ParticleOfInterest {
99- #define X(name) name,
100- LIST_OF_PARTICLES_OF_INTEREST
101- #undef X
102- ParticleOfInterestSize
103- };
104-
105- static constexpr std::array<std::string_view, ParticleOfInterestSize> particleOfInterestLabels{
106- #define X(name) #name,
107- LIST_OF_PARTICLES_OF_INTEREST
108- #undef X
109- };
110-
111- static constexpr auto particleOfInterestLabels = std::to_array<std::string_view>({
112- #define X(name) #name,
113- LIST_OF_PARTICLES_OF_INTEREST
114- #undef X
115- });
116- */
117-
11883using EffMapPtr = std::variant<std::shared_ptr<TH2>, std::shared_ptr<TH3>>;
11984
120- /* struct BoundEfficiencyMap {
121- using CoordsTuple = std::tuple<float, float, float>;
122-
123- const TH3* effMap;
124- CoordsTuple coords;
125-
126- BoundEfficiencyMap(const std::shared_ptr<TH3>& effMap, float x, float y, float z) : effMap(effMap.get()), coords(x, y, z) {}
127- BoundEfficiencyMap(const std::shared_ptr<TH3>& effMap, const CoordsTuple& coords) : effMap(effMap.get()), coords(coords) {}
128-
129- float getBinEfficiency() const
130- {
131- if (!effMap) {
132- return 1.0f;
133- }
134-
135- const auto& [x, y, z] = coords;
136- return effMap->GetBinContent(effMap->FindFixBin(x, y, z));
137- }
138-
139- float interpolateEfficiency() const
140- {
141- if (!effMap) {
142- return 1.0f;
143- }
144-
145- const auto& [x, y, z] = coords;
146- return effMap->Interpolate(x, y, z);
147- }
148- };*/
149-
15085struct BoundEfficiencyMap {
15186 using CoordsTuple = std::tuple<float , float , float >;
15287
@@ -156,12 +91,12 @@ struct BoundEfficiencyMap {
15691 BoundEfficiencyMap (const EffMapPtr& effMap, float x, float y, float z) : effMap(effMap), coords(x, y, z) {}
15792 BoundEfficiencyMap (const EffMapPtr& effMap, const CoordsTuple& coords) : effMap(effMap), coords(coords) {}
15893
159- float getBinEfficiency () const
94+ std::pair< float , float > getBinEfficiencyAndError () const
16095 {
16196 return std::visit (
162- [this ](auto && mapPtr) -> float {
97+ [this ](auto && mapPtr) -> std::pair< float , float > {
16398 if (!mapPtr)
164- return 1 .0f ;
99+ return { 1 .0f , 0 . 0f } ;
165100
166101 const auto & [x, y, z] = coords;
167102
@@ -170,9 +105,11 @@ struct BoundEfficiencyMap {
170105
171106 // Compile-time branching: generates the exact correct function call
172107 if constexpr (std::is_same_v<HistoType, TH2>) {
173- return mapPtr->GetBinContent (mapPtr->FindFixBin (y, z)); // 2D case only
108+ int bin = mapPtr->FindFixBin (y, z); // 2D case only
109+ return {mapPtr->GetBinContent (bin), mapPtr->GetBinError (bin)};
174110 } else {
175- return mapPtr->GetBinContent (mapPtr->FindFixBin (x, y, z)); // Full 3D case
111+ int bin = mapPtr->FindFixBin (x, y, z); // Full 3D case
112+ return {mapPtr->GetBinContent (bin), mapPtr->GetBinError (bin)};
176113 }
177114 },
178115 effMap);
@@ -264,7 +201,8 @@ struct PhiStrangenessCorrelation {
264201 struct : ConfigurableGroup {
265202 Configurable<bool > applyEfficiency{" applyEfficiency" , false , " Use efficiency for filling histograms" };
266203 Configurable<bool > useEffInterpolation{" useEffInterpolation" , false , " If true, interpolates efficiency map, else uses bin center" };
267- Configurable<bool > applyPhiEfficiency{" applyPhiEfficiency" , true , " Apply efficiency for Phi candidates" };
204+ Configurable<bool > applyPhiEfficiency{" applyPhiEfficiency" , false , " Apply efficiency for Phi candidates" };
205+ Configurable<bool > propagateEffError{" propagateEffError" , false , " Propagate efficiency error" };
268206 } efficiencyConfigs;
269207
270208 // Configurable for event mixing
@@ -473,14 +411,39 @@ struct PhiStrangenessCorrelation {
473411
474412 // Compute weight based on efficiencies
475413 template <typename ... BoundEffMaps>
476- float computeWeight (const BoundEffMaps&... boundEffMaps)
414+ std::pair< float , float > computeWeightAndError (const BoundEffMaps&... boundEffMaps)
477415 {
478416 if (!efficiencyConfigs.applyEfficiency )
479- return 1 .0f ;
417+ return {1 .0f , 0 .0f };
418+
419+ float effTot = 1 .0f ;
420+ float relErrSqSum = 0 .0f ;
421+
422+ auto processMap = [&](const auto & boundMap) {
423+ auto [eff, err] = boundMap.getBinEfficiencyAndError (); // Unpack efficiency and error from the map
424+
425+ if (efficiencyConfigs.useEffInterpolation ) {
426+ eff = boundMap.interpolateEfficiency ();
427+ // For simplicity, we keep the error from the bin content even when interpolating, but this can be refined if needed
428+ }
429+
430+ effTot *= eff;
480431
481- float totalEfficiency = ((efficiencyConfigs.useEffInterpolation ? boundEffMaps.interpolateEfficiency () : boundEffMaps.getBinEfficiency ()) * ...);
432+ if (eff > 0 .0f ) {
433+ float mapErr = efficiencyConfigs.propagateEffError ? err : 0 .0f ; // Optionally propagate error, otherwise treat as zero
434+ relErrSqSum += (mapErr / eff) * (mapErr / eff);
435+ }
436+ };
437+
438+ (processMap (boundEffMaps), ...); // Fold expression to process all bound efficiency maps
439+
440+ if (effTot <= 0 .0f )
441+ return {1 .0f , 0 .0f };
482442
483- return totalEfficiency <= 0 .0f ? 1 .0f : 1 .0f / totalEfficiency;
443+ float weight = 1 .0f / effTot;
444+ float weightErr = weight * std::sqrt (relErrSqSum); // Propagate relative error to the weight
445+
446+ return {weight, weightErr};
484447 }
485448
486449 float getDeltaPhi (float phiTrigger, float phiAssociated)
@@ -497,6 +460,72 @@ struct PhiStrangenessCorrelation {
497460 return std::distance (binsMult->begin (), it) - 1 ;
498461 }
499462
463+ template <typename THist, typename ... Args>
464+ void customFillHist (auto histId, const std::pair<float , float >& weightPair, Args... coords)
465+ {
466+ auto hist = histos.get <THist>(histId);
467+ if (!hist) {
468+ return ;
469+ }
470+
471+ // Extract the weight and its propagated uncertainty
472+ auto [w, wErr] = weightPair;
473+
474+ // Find the global bin number for the given physical coordinates
475+ int bin = hist->FindFixBin (coords...);
476+
477+ // Retrieve the previous bin content and its absolute error
478+ double prevContent = hist->GetBinContent (bin);
479+ double prevErr = hist->GetBinError (bin);
480+
481+ // Calculate the new content by adding the current weight
482+ double newContent = prevContent + w;
483+
484+ // Error propagation in quadrature:
485+ // prevErr^2 : previous accumulated variance
486+ // (w * w) : statistical fluctuation added by this specific particle count
487+ // (wErr * wErr) : systematic uncertainty added by the efficiency map inaccuracy
488+ double newErr = std::sqrt (prevErr * prevErr + (w * w) + (wErr * wErr));
489+
490+ // Update the histogram bin
491+ hist->SetBinContent (bin, newContent);
492+ hist->SetBinError (bin, newErr);
493+ }
494+
495+ template <typename ... Args>
496+ void customFillTHn (auto histId, const std::pair<float , float >& weightPair, Args... coords)
497+ {
498+ auto hist = histos.get <THnSparse>(histId);
499+ if (!hist) {
500+ return ;
501+ }
502+
503+ // Extract the weight and its propagated uncertainty
504+ auto [w, wErr] = weightPair;
505+
506+ // Create an array of floats for the coordinates to pass to GetBin and SetBinContent
507+ double coordArray[] = {static_cast <double >(coords)...};
508+
509+ // Find the bin number.
510+ // The 'true' flag is mandatory: it allocates the bin in memory if it doesn't exist yet.
511+ long bin = hist->GetBin (coordArray, true );
512+
513+ // Retrieve the previous content and the squared error (variance)
514+ double prevContent = hist->GetBinContent (bin);
515+ double prevErr2 = hist->GetBinError2 (bin);
516+
517+ // Calculate the new content
518+ double newContent = prevContent + w;
519+
520+ // Add the new variances to the accumulated variance
521+ // (w * w) is the statistical term, (wErr * wErr) is the map uncertainty term
522+ double newErr2 = prevErr2 + (w * w) + (wErr * wErr);
523+
524+ // Update the THnSparse bin
525+ hist->SetBinContent (bin, newContent);
526+ hist->SetBinError2 (bin, newErr2); // Note: SetBinError2 takes the squared error directly
527+ }
528+
500529 template <typename TCollision, typename TPhiCands, typename TK0SCands, typename TPionCands>
501530 void processPhiK0SPionSE (TCollision const & collision, TPhiCands const & phiCandidates, TK0SCands const & k0sReduced, TPionCands const & pionTracks)
502531 {
@@ -527,9 +556,11 @@ struct PhiStrangenessCorrelation {
527556 if (efficiencyConfigs.applyEfficiency && efficiencyConfigs.applyPhiEfficiency && phiCand.pt () >= binspTPhi->back ())
528557 continue ;
529558
530- float weightPhi = computeWeight (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()));
559+ // float weightPhi = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()));
560+ auto weightPhi = computeWeightAndError (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()));
531561
532- histos.fill (HIST (" phi/h3PhiData" ), multiplicity, phiCand.pt (), phiCand.m (), weightPhi);
562+ // histos.fill(HIST("phi/h3PhiData"), multiplicity, phiCand.pt(), phiCand.m(), weightPhi);
563+ customFillHist<TH3>(HIST (" phi/h3PhiData" ), weightPhi, multiplicity, phiCand.pt (), phiCand.m ());
533564
534565 auto processCorrelations = [&](auto fillK0S, auto fillPion) {
535566 if (doAssocCorrelations[kK0S ]) {
@@ -538,8 +569,10 @@ struct PhiStrangenessCorrelation {
538569 if (!isK0sValid (k0s))
539570 continue ;
540571
541- float weightPhiK0S = computeWeight (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
542- BoundEfficiencyMap (effMapsAssoc[kK0S ], multiplicity, k0s.pt (), k0s.y ()));
572+ // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
573+ // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y()));
574+ auto weightPhiK0S = computeWeightAndError (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
575+ BoundEfficiencyMap (effMapsAssoc[kK0S ], multiplicity, k0s.pt (), k0s.y ()));
543576 fillK0S (k0s, weightPhiK0S);
544577 }
545578 }
@@ -550,8 +583,10 @@ struct PhiStrangenessCorrelation {
550583 if (!isPionValid (pionTrack))
551584 continue ;
552585
553- float weightPhiPion = computeWeight (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
554- BoundEfficiencyMap (effMapsAssoc[kPion ], multiplicity, pionTrack.pt (), pionTrack.y ()));
586+ // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
587+ // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, pionTrack.pt(), pionTrack.y()));
588+ auto weightPhiPion = computeWeightAndError (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
589+ BoundEfficiencyMap (effMapsAssoc[kPion ], multiplicity, pionTrack.pt (), pionTrack.y ()));
555590 fillPion (pionTrack, weightPhiPion);
556591 }
557592 }
@@ -563,12 +598,17 @@ struct PhiStrangenessCorrelation {
563598 auto piTOFHistID = HIST (" phiPi/h6PhiPiTOFData" );
564599
565600 processCorrelations (
566- [&](const auto & k0s, float w) {
567- histos.fill (k0sHistID, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), phiCand.m (), k0s.m (), w);
601+ // [&](const auto& k0s, float w) {
602+ // histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w);
603+ [&](const auto & k0s, const std::pair<float , float >& w) {
604+ customFillTHn (k0sHistID, w, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), phiCand.m (), k0s.m ());
568605 },
569- [&](const auto & pion, float w) {
570- histos.fill (piTPCHistID, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTPC (), w);
571- histos.fill (piTOFHistID, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTOF (), w);
606+ // [&](const auto& pion, float w) {
607+ // histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w);
608+ // histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w);
609+ [&](const auto & pion, const std::pair<float , float >& w) {
610+ customFillTHn (piTPCHistID, w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTPC ());
611+ customFillTHn (piTOFHistID, w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTOF ());
572612 });
573613 } else if (analysisMode == kDeltaYvsDeltaPhi ) {
574614 auto k0sHistID = std::make_tuple (HIST (" phiK0S/h5PhiK0SDataSignal" ), HIST (" phiK0S/h5PhiK0SDataSideband" ));
@@ -585,11 +625,15 @@ struct PhiStrangenessCorrelation {
585625 // auto piHistID = HIST("phiPi/h5PhiPiData") + HIST(phiMassRegionLabels[i]);
586626
587627 processCorrelations (
588- [&](const auto & k0s, float w) {
589- histos.fill (std::get<i>(k0sHistID), multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), getDeltaPhi (phiCand.phi (), k0s.phi ()), w);
628+ // [&](const auto& k0s, float w) {
629+ // histos.fill(std::get<i>(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w);
630+ [&](const auto & k0s, const std::pair<float , float >& w) {
631+ customFillTHn (std::get<i>(k0sHistID), w, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), getDeltaPhi (phiCand.phi (), k0s.phi ()));
590632 },
591- [&](const auto & pion, float w) {
592- histos.fill (std::get<i>(piHistID), multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), getDeltaPhi (phiCand.phi (), pion.phi ()), w);
633+ // [&](const auto& pion, float w) {
634+ // histos.fill(std::get<i>(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w);
635+ [&](const auto & pion, const std::pair<float , float >& w) {
636+ customFillTHn (std::get<i>(piHistID), w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), getDeltaPhi (phiCand.phi (), pion.phi ()));
593637 });
594638 });
595639 }
@@ -668,17 +712,21 @@ struct PhiStrangenessCorrelation {
668712 continue ;
669713
670714 auto processCorrelations = [&](auto fillK0S) {
671- float weightPhiK0S = computeWeight (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
672- BoundEfficiencyMap (effMapsAssoc[kK0S ], multiplicity, k0s.pt (), k0s.y ()));
715+ // float weightPhiK0S = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
716+ // BoundEfficiencyMap(effMapsAssoc[kK0S], multiplicity, k0s.pt(), k0s.y()));
717+ auto weightPhiK0S = computeWeightAndError (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
718+ BoundEfficiencyMap (effMapsAssoc[kK0S ], multiplicity, k0s.pt (), k0s.y ()));
673719 fillK0S (k0s, weightPhiK0S);
674720 };
675721
676722 if (analysisMode == kMassvsMass ) {
677723 auto k0sHistID = HIST (" phiK0S/h6PhiK0SDataME" );
678724
679725 processCorrelations (
680- [&](const auto & k0s, float w) {
681- histos.fill (k0sHistID, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), phiCand.m (), k0s.m (), w);
726+ // [&](const auto& k0s, float w) {
727+ // histos.fill(k0sHistID, multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), phiCand.m(), k0s.m(), w);
728+ [&](const auto & k0s, const std::pair<float , float >& w) {
729+ customFillTHn (k0sHistID, w, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), phiCand.m (), k0s.m ());
682730 });
683731 } else if (analysisMode == kDeltaYvsDeltaPhi ) {
684732 auto k0sHistID = std::make_tuple (HIST (" phiK0S/h5PhiK0SDataMESignal" ), HIST (" phiK0S/h5PhiK0SDataMESideband" ));
@@ -691,8 +739,10 @@ struct PhiStrangenessCorrelation {
691739 return ;
692740
693741 processCorrelations (
694- [&](const auto & k0s, float w) {
695- histos.fill (std::get<i>(k0sHistID), multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), getDeltaPhi (phiCand.phi (), k0s.phi ()), w);
742+ // [&](const auto& k0s, float w) {
743+ // histos.fill(std::get<i>(k0sHistID), multiplicity, phiCand.pt(), k0s.pt(), phiCand.y() - k0s.y(), getDeltaPhi(phiCand.phi(), k0s.phi()), w);
744+ [&](const auto & k0s, const std::pair<float , float >& w) {
745+ customFillTHn (std::get<i>(k0sHistID), w, multiplicity, phiCand.pt (), k0s.pt (), phiCand.y () - k0s.y (), getDeltaPhi (phiCand.phi (), k0s.phi ()));
696746 });
697747 });
698748 }
@@ -763,8 +813,10 @@ struct PhiStrangenessCorrelation {
763813 continue ;
764814
765815 auto processCorrelations = [&](auto fillPion) {
766- float weightPhiPion = computeWeight (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
767- BoundEfficiencyMap (effMapsAssoc[kPion ], multiplicity, piTrack.pt (), piTrack.y ()));
816+ // float weightPhiPion = computeWeight(BoundEfficiencyMap(effMapPhi, multiplicity, phiCand.pt(), phiCand.y()),
817+ // BoundEfficiencyMap(effMapsAssoc[kPion], multiplicity, piTrack.pt(), piTrack.y()));
818+ auto weightPhiPion = computeWeightAndError (BoundEfficiencyMap (effMapPhi, multiplicity, phiCand.pt (), phiCand.y ()),
819+ BoundEfficiencyMap (effMapsAssoc[kPion ], multiplicity, piTrack.pt (), piTrack.y ()));
768820 fillPion (piTrack, weightPhiPion);
769821 };
770822
@@ -773,9 +825,12 @@ struct PhiStrangenessCorrelation {
773825 auto piTOFHistID = HIST (" phiPi/h6PhiPiTOFDataME" );
774826
775827 processCorrelations (
776- [&](const auto & pion, float w) {
777- histos.fill (piTPCHistID, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTPC (), w);
778- histos.fill (piTOFHistID, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTOF (), w);
828+ // [&](const auto& pion, float w) {
829+ // histos.fill(piTPCHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTPC(), w);
830+ // histos.fill(piTOFHistID, multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), phiCand.m(), pion.nSigmaTOF(), w);
831+ [&](const auto & pion, const std::pair<float , float >& w) {
832+ customFillTHn (piTPCHistID, w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTPC ());
833+ customFillTHn (piTOFHistID, w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), phiCand.m (), pion.nSigmaTOF ());
779834 });
780835 } else if (analysisMode == kDeltaYvsDeltaPhi ) {
781836 auto piHistID = std::make_tuple (HIST (" phiPi/h5PhiPiDataMESignal" ), HIST (" phiPi/h5PhiPiDataMESideband" ));
@@ -788,8 +843,10 @@ struct PhiStrangenessCorrelation {
788843 return ;
789844
790845 processCorrelations (
791- [&](const auto & pion, float w) {
792- histos.fill (std::get<i>(piHistID), multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), getDeltaPhi (phiCand.phi (), pion.phi ()), w);
846+ // [&](const auto& pion, float w) {
847+ // histos.fill(std::get<i>(piHistID), multiplicity, phiCand.pt(), pion.pt(), phiCand.y() - pion.y(), getDeltaPhi(phiCand.phi(), pion.phi()), w);
848+ [&](const auto & pion, const std::pair<float , float >& w) {
849+ customFillTHn (std::get<i>(piHistID), w, multiplicity, phiCand.pt (), pion.pt (), phiCand.y () - pion.y (), getDeltaPhi (phiCand.phi (), pion.phi ()));
793850 });
794851 });
795852 }
0 commit comments