4444#include < cmath>
4545#include < cstdint>
4646#include < iterator>
47+ #include < limits>
4748#include < memory>
4849#include < ranges>
4950#include < string>
@@ -561,11 +562,12 @@ struct HolderTrack {
561562 static constexpr double truncateNSigmaPid (const double value) { return (!(std::abs (value) < TruncationAbsNSigmaPid) ? -TruncationAbsNSigmaPid : value); }
562563
563564 std::int32_t sign = 0 ;
565+ double p = 0 .;
564566 double pt = 0 .;
565567 double eta = 0 .;
566568 double phi = 0 .;
567- double ptOverQ = 0 .;
568569 double pOverQ = 0 .;
570+ double ptOverQ = 0 .;
569571 double rapidityPi = 0 .;
570572 double rapidityKa = 0 .;
571573 double rapidityPr = 0 .;
@@ -585,11 +587,12 @@ struct HolderTrack {
585587 void clear ()
586588 {
587589 sign = 0 ;
590+ p = 0 .;
588591 pt = 0 .;
589592 eta = 0 .;
590593 phi = 0 .;
591- ptOverQ = 0 .;
592594 pOverQ = 0 .;
595+ ptOverQ = 0 .;
593596 rapidityPi = 0 .;
594597 rapidityKa = 0 .;
595598 rapidityPr = 0 .;
@@ -974,33 +977,48 @@ struct PartNumFluc {
974977
975978 if (doprocessMc.value ) {
976979 if (cfgFlagCalculationPurityPi.value || cfgFlagCalculationPurityKa.value || cfgFlagCalculationPurityPr.value ) {
977- HistogramConfigSpec hcsCalculationPurity (HistType::kTProfile3D , {{20 , 0 ., 100 ., " Centrality (%)" }, {20 , 0 ., 2 ., " #it{p}_{T} (GeV/#it{c})" }, {24 , -1.2 , 1.2 , " #it{#eta}" }});
980+ AxisSpec asQaCentrality (20 , 0 ., 100 ., " Centrality (%)" );
981+ AxisSpec asEta (24 , -1.2 , 1.2 , " #it{#eta}" );
982+ HistogramConfigSpec hcsCalculationPurityP (HistType::kTProfile3D , {asQaCentrality, {35 , 0 ., 3.5 , " #it{p} (GeV/#it{c})" }, asEta});
983+ HistogramConfigSpec hcsCalculationPurityPt (HistType::kTProfile3D , {asQaCentrality, {20 , 0 ., 2 ., " #it{p}_{T} (GeV/#it{c})" }, asEta});
978984
979985 if (cfgFlagCalculationPurityPi.value ) {
980986 LOG (info) << " Enabling pion purity calculation." ;
981987
982- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPiP" , " " , hcsCalculationPurity);
983- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPiM" , " " , hcsCalculationPurity);
984- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiP" , " " , hcsCalculationPurity);
985- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiM" , " " , hcsCalculationPurity);
988+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcPiP" , " " , hcsCalculationPurityP);
989+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcPiM" , " " , hcsCalculationPurityP);
990+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofPiP" , " " , hcsCalculationPurityP);
991+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofPiM" , " " , hcsCalculationPurityP);
992+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPiP" , " " , hcsCalculationPurityPt);
993+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPiM" , " " , hcsCalculationPurityPt);
994+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiP" , " " , hcsCalculationPurityPt);
995+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiM" , " " , hcsCalculationPurityPt);
986996 }
987997
988998 if (cfgFlagCalculationPurityKa.value ) {
989999 LOG (info) << " Enabling kaon purity calculation." ;
9901000
991- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcKaP" , " " , hcsCalculationPurity);
992- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcKaM" , " " , hcsCalculationPurity);
993- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaP" , " " , hcsCalculationPurity);
994- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaM" , " " , hcsCalculationPurity);
1001+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcKaP" , " " , hcsCalculationPurityP);
1002+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcKaM" , " " , hcsCalculationPurityP);
1003+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofKaP" , " " , hcsCalculationPurityP);
1004+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofKaM" , " " , hcsCalculationPurityP);
1005+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcKaP" , " " , hcsCalculationPurityPt);
1006+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcKaM" , " " , hcsCalculationPurityPt);
1007+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaP" , " " , hcsCalculationPurityPt);
1008+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaM" , " " , hcsCalculationPurityPt);
9951009 }
9961010
9971011 if (cfgFlagCalculationPurityPr.value ) {
9981012 LOG (info) << " Enabling (anti)proton purity calculation." ;
9991013
1000- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPrP" , " " , hcsCalculationPurity);
1001- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPrM" , " " , hcsCalculationPurity);
1002- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrP" , " " , hcsCalculationPurity);
1003- hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrM" , " " , hcsCalculationPurity);
1014+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcPrP" , " " , hcsCalculationPurityP);
1015+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcPrM" , " " , hcsCalculationPurityP);
1016+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofPrP" , " " , hcsCalculationPurityP);
1017+ hrCalculationPurity.add (" CalculationPurity/pCentralityPEtaPurityTpcTofPrM" , " " , hcsCalculationPurityP);
1018+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPrP" , " " , hcsCalculationPurityPt);
1019+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcPrM" , " " , hcsCalculationPurityPt);
1020+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrP" , " " , hcsCalculationPurityPt);
1021+ hrCalculationPurity.add (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrM" , " " , hcsCalculationPurityPt);
10041022 }
10051023 }
10061024
@@ -1204,7 +1222,7 @@ struct PartNumFluc {
12041222 return nullptr ;
12051223 }
12061224 }();
1207- return hCentralityPtEtaShiftNSigmaPid ? hCentralityPtEtaShiftNSigmaPid->Interpolate (std::max (hCentralityPtEtaShiftNSigmaPid->GetXaxis ()->GetBinCenter (1 ), std::min (holderEvent.centrality , hCentralityPtEtaShiftNSigmaPid->GetXaxis ()->GetBinCenter (hCentralityPtEtaShiftNSigmaPid->GetNbinsX ()))), std::max (hCentralityPtEtaShiftNSigmaPid->GetYaxis ()->GetBinCenter (1 ), std::min (holderTrack.pt , hCentralityPtEtaShiftNSigmaPid->GetYaxis ()->GetBinCenter (hCentralityPtEtaShiftNSigmaPid->GetNbinsY ()))), std::max (hCentralityPtEtaShiftNSigmaPid->GetZaxis ()->GetBinCenter (1 ), std::min (holderTrack.eta , hCentralityPtEtaShiftNSigmaPid->GetZaxis ()->GetBinCenter (hCentralityPtEtaShiftNSigmaPid->GetNbinsZ ())))) : 0 .;
1225+ return hCentralityPtEtaShiftNSigmaPid ? hCentralityPtEtaShiftNSigmaPid->Interpolate(std::max(std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetXaxis()->GetBinCenter(1), std::numeric_limits<double>::infinity()), std::min(holderEvent.centrality, std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetXaxis()->GetBinCenter(hCentralityPtEtaShiftNSigmaPid->GetNbinsX()), -std::numeric_limits<double>::infinity()))), std::max(std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetYaxis()->GetBinCenter(1), std::numeric_limits<double>::infinity()), std::min(holderTrack.pt, std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetYaxis()->GetBinCenter(hCentralityPtEtaShiftNSigmaPid->GetNbinsY()), -std::numeric_limits<double>::infinity()))), std::max(std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetZaxis()->GetBinCenter(1), std::numeric_limits<double>::infinity()), std::min(holderTrack.eta, std::nextafter(hCentralityPtEtaShiftNSigmaPid->GetZaxis()->GetBinCenter(hCentralityPtEtaShiftNSigmaPid->GetNbinsZ()), -std::numeric_limits<double>::infinity())))) : 0.;
12081226 }
12091227
12101228 template <bool doProcessingMc>
@@ -1381,11 +1399,12 @@ struct PartNumFluc {
13811399 {
13821400 holderTrack.clear ();
13831401 holderTrack.sign = track.sign ();
1402+ holderTrack.p = track.p ();
13841403 holderTrack.pt = track.pt ();
13851404 holderTrack.eta = track.eta ();
13861405 holderTrack.phi = track.phi ();
1406+ holderTrack.pOverQ = holderTrack.p / holderTrack.sign ;
13871407 holderTrack.ptOverQ = holderTrack.pt / holderTrack.sign ;
1388- holderTrack.pOverQ = track.p () / holderTrack.sign ;
13891408 holderTrack.rapidityPi = track.rapidity (constants::physics::MassPiPlus);
13901409 holderTrack.rapidityKa = track.rapidity (constants::physics::MassKPlus);
13911410 holderTrack.rapidityPr = track.rapidity (constants::physics::MassProton);
@@ -2174,9 +2193,11 @@ struct PartNumFluc {
21742193 if (cfgFlagCalculationPurityPi.value ) {
21752194 switch (isPi<false , true >()) {
21762195 case 1 :
2196+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcPiP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiPlus ? 1 . : 0 .);
21772197 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcPiP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiPlus ? 1 . : 0 .);
21782198 break ;
21792199 case -1 :
2200+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcPiM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiMinus ? 1 . : 0 .);
21802201 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcPiM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiMinus ? 1 . : 0 .);
21812202 break ;
21822203 }
@@ -2185,9 +2206,11 @@ struct PartNumFluc {
21852206 if (cfgFlagCalculationPurityKa.value ) {
21862207 switch (isKa<false , true >()) {
21872208 case 1 :
2209+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcKaP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKPlus ? 1 . : 0 .);
21882210 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcKaP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKPlus ? 1 . : 0 .);
21892211 break ;
21902212 case -1 :
2213+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcKaM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKMinus ? 1 . : 0 .);
21912214 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcKaM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKMinus ? 1 . : 0 .);
21922215 break ;
21932216 }
@@ -2196,9 +2219,11 @@ struct PartNumFluc {
21962219 if (cfgFlagCalculationPurityPr.value ) {
21972220 switch (isPr<false , true >()) {
21982221 case 1 :
2222+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcPrP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProton ? 1 . : 0 .);
21992223 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcPrP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProton ? 1 . : 0 .);
22002224 break ;
22012225 case -1 :
2226+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcPrM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProtonBar ? 1 . : 0 .);
22022227 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcPrM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProtonBar ? 1 . : 0 .);
22032228 break ;
22042229 }
@@ -2208,9 +2233,11 @@ struct PartNumFluc {
22082233 if (cfgFlagCalculationPurityPi.value ) {
22092234 switch (isPi<true , true >()) {
22102235 case 1 :
2236+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofPiP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiPlus ? 1 . : 0 .);
22112237 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiPlus ? 1 . : 0 .);
22122238 break ;
22132239 case -1 :
2240+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofPiM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiMinus ? 1 . : 0 .);
22142241 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofPiM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kPiMinus ? 1 . : 0 .);
22152242 break ;
22162243 }
@@ -2219,9 +2246,11 @@ struct PartNumFluc {
22192246 if (cfgFlagCalculationPurityKa.value ) {
22202247 switch (isKa<true , true >()) {
22212248 case 1 :
2249+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofKaP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKPlus ? 1 . : 0 .);
22222250 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKPlus ? 1 . : 0 .);
22232251 break ;
22242252 case -1 :
2253+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofKaM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKMinus ? 1 . : 0 .);
22252254 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofKaM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kKMinus ? 1 . : 0 .);
22262255 break ;
22272256 }
@@ -2230,9 +2259,11 @@ struct PartNumFluc {
22302259 if (cfgFlagCalculationPurityPr.value ) {
22312260 switch (isPr<true , true >()) {
22322261 case 1 :
2262+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofPrP" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProton ? 1 . : 0 .);
22332263 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrP" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProton ? 1 . : 0 .);
22342264 break ;
22352265 case -1 :
2266+ hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPEtaPurityTpcTofPrM" ), holderEvent.centrality , holderTrack.p , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProtonBar ? 1 . : 0 .);
22362267 hrCalculationPurity.fill (HIST (" CalculationPurity/pCentralityPtEtaPurityTpcTofPrM" ), holderEvent.centrality , holderTrack.pt , holderTrack.eta , holderMcParticle.pdgCode == PDG_t::kProtonBar ? 1 . : 0 .);
22372268 break ;
22382269 }
0 commit comments