1414// /
1515// / \author M. Hemmer, marvin.hemmer@cern.ch
1616
17- #include < numbers>
18- #include < iterator>
19- #include < array>
20- #include < string>
21- #include < map>
22- #include < algorithm>
23- #include < vector>
24- #include < tuple>
25- #include < utility>
26- #include < cmath>
27-
28- #include " Math/Vector4D.h"
29- #include " Math/Vector3D.h"
30- #include " Math/LorentzRotation.h"
31- #include " Math/Rotation3D.h"
32- #include " Math/AxisAngle.h"
33-
34- #include " CCDB/BasicCCDBManager.h"
35- #include " Framework/AnalysisTask.h"
36- #include " Framework/ASoAHelpers.h"
37- #include " Framework/HistogramRegistry.h"
38- #include " Framework/runDataProcessing.h"
39-
40- #include " Common/Core/EventPlaneHelper.h"
41- #include " Common/Core/RecoDecay.h"
42- #include " Common/DataModel/Qvectors.h"
43- #include " CommonConstants/MathConstants.h"
44-
45- #include " DetectorsBase/GeometryManager.h"
46- #include " DataFormatsEMCAL/Constants.h"
47- #include " EMCALBase/Geometry.h"
48- #include " EMCALCalib/BadChannelMap.h"
49-
50- #include " PWGEM/Dilepton/Utils/EMTrackUtilities.h"
5117#include " PWGEM/PhotonMeson/Core/EMCPhotonCut.h"
5218#include " PWGEM/PhotonMeson/Core/EMPhotonEventCut.h"
5319#include " PWGEM/PhotonMeson/DataModel/gammaTables.h"
54- #include " PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h"
55- #include " PWGEM/PhotonMeson/Utils/PairUtilities.h"
5620#include " PWGEM/PhotonMeson/Utils/EventHistograms.h"
57- #include " PWGEM/PhotonMeson/Utils/NMHistograms.h"
21+ #include " PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h"
22+
23+ #include " Common/CCDB/TriggerAliases.h"
24+ #include " Common/Core/EventPlaneHelper.h"
25+ #include " Common/Core/RecoDecay.h"
26+ #include " Common/DataModel/Centrality.h"
27+ #include " Common/DataModel/EventSelection.h"
28+
29+ #include < CCDB/BasicCCDBManager.h>
30+ #include < CommonConstants/MathConstants.h>
31+ #include < EMCALBase/Geometry.h>
32+ #include < EMCALBase/GeometryBase.h>
33+ #include < EMCALCalib/BadChannelMap.h>
34+ #include < Framework/ASoA.h>
35+ #include < Framework/ASoAHelpers.h>
36+ #include < Framework/AnalysisDataModel.h>
37+ #include < Framework/AnalysisHelpers.h>
38+ #include < Framework/AnalysisTask.h>
39+ #include < Framework/BinningPolicy.h>
40+ #include < Framework/Configurable.h>
41+ #include < Framework/Expressions.h>
42+ #include < Framework/GroupedCombinations.h>
43+ #include < Framework/HistogramRegistry.h>
44+ #include < Framework/HistogramSpec.h>
45+ #include < Framework/InitContext.h>
46+ #include < Framework/OutputObjHeader.h>
47+ #include < Framework/SliceCache.h>
48+ #include < Framework/runDataProcessing.h>
49+
50+ #include < Math/GenVector/AxisAngle.h>
51+ #include < Math/GenVector/Rotation3D.h>
52+ #include < Math/Vector4D.h> // IWYU pragma: keep
53+ #include < TH1.h>
54+ #include < TString.h>
55+
56+ #include < array>
57+ #include < cmath>
58+ #include < cstdint>
59+ #include < numbers>
60+ #include < string>
61+ #include < string_view>
62+ #include < tuple>
63+ #include < utility>
64+ #include < vector>
5865
5966using namespace o2 ;
6067using namespace o2 ::aod;
6168using namespace o2 ::framework;
6269using namespace o2 ::framework::expressions;
6370using namespace o2 ::soa;
64- using namespace o2 ::aod::pwgem::photonmeson::photonpair;
6571using namespace o2 ::aod::pwgem::photon;
66- using namespace o2 ::aod::pwgem::dilepton::utils;
67-
68- enum QvecEstimator { FT0M = 0 ,
69- FT0A = 1 ,
70- FT0C,
71- TPCPos,
72- TPCNeg,
73- TPCTot };
74-
75- enum CentralityEstimator { None = 0 ,
76- CFT0A = 1 ,
77- CFT0C,
78- CFT0M,
79- NCentralityEstimators
72+
73+ enum QvecEstimator {
74+ FT0M = 0 ,
75+ FT0A = 1 ,
76+ FT0C,
77+ TPCPos,
78+ TPCNeg,
79+ TPCTot
80+ };
81+
82+ enum CentralityEstimator {
83+ None = 0 ,
84+ CFT0A = 1 ,
85+ CFT0C,
86+ CFT0M,
87+ NCentralityEstimators
88+ };
89+
90+ enum Harmonics {
91+ kNone = 0 ,
92+ kDirect = 1 ,
93+ kElliptic = 2 ,
94+ kTriangluar = 3 ,
95+ kQuadrangular = 4 ,
96+ kPentagonal = 5 ,
97+ kHexagonal = 6 ,
98+ kHeptagonal = 7 ,
99+ kOctagonal = 8
80100};
81101
82102struct TaskPi0FlowEMC {
@@ -98,6 +118,8 @@ struct TaskPi0FlowEMC {
98118 Configurable<bool > cfgDoM02{" cfgDoM02" , false , " Flag to enable flow vs M02 for single photons" };
99119 Configurable<bool > cfgDoReverseScaling{" cfgDoReverseScaling" , false , " Flag to reverse the scaling that is possibly applied during NonLin" };
100120 Configurable<bool > cfgDoPlaneQA{" cfgDoPlaneQA" , false , " Flag to enable QA plots comparing in and out of plane" };
121+ Configurable<float > cfgMaxQVector{" cfgMaxQVector" , 20 .f , " Maximum allowed absolute QVector value." };
122+ Configurable<float > cfgMaxAsymmetry{" cfgMaxAsymmetry" , 0 .1f , " Maximum allowed asymmetry for photon pairs used in calibration." };
101123
102124 // configurable axis
103125 ConfigurableAxis thnConfigAxisInvMass{" thnConfigAxisInvMass" , {400 , 0.0 , 0.8 }, " " };
@@ -205,6 +227,9 @@ struct TaskPi0FlowEMC {
205227 std::vector<int8_t > lookupTable1D;
206228 float epsilon = 1 .e-8 ;
207229
230+ // static constexpr
231+ static constexpr int64_t NMinPhotonRotBkg = 3 ;
232+
208233 // To access the 1D array
209234 inline int getIndex (int iEta, int iPhi)
210235 {
@@ -265,7 +290,7 @@ struct TaskPi0FlowEMC {
265290
266291 void init (InitContext&)
267292 {
268- if (harmonic != 2 && harmonic != 3 ) {
293+ if (harmonic != kElliptic && harmonic != kTriangluar ) {
269294 LOG (info) << " Harmonic was set to " << harmonic << " but can only be 2 or 3!" ;
270295 }
271296
@@ -492,65 +517,65 @@ struct TaskPi0FlowEMC {
492517
493518 switch (detector) {
494519 case QvecEstimator::FT0M:
495- if (harmonic == 2 ) {
520+ if (harmonic == kElliptic ) {
496521 xQVec = collision.q2xft0m ();
497522 yQVec = collision.q2yft0m ();
498- } else if (harmonic == 3 ) {
523+ } else if (harmonic == kTriangluar ) {
499524 xQVec = collision.q3xft0m ();
500525 yQVec = collision.q3yft0m ();
501526 }
502527 break ;
503528 case QvecEstimator::FT0A:
504- if (harmonic == 2 ) {
529+ if (harmonic == kElliptic ) {
505530 xQVec = collision.q2xft0a ();
506531 yQVec = collision.q2yft0a ();
507- } else if (harmonic == 3 ) {
532+ } else if (harmonic == kTriangluar ) {
508533 xQVec = collision.q3xft0a ();
509534 yQVec = collision.q3yft0a ();
510535 }
511536 break ;
512537 case QvecEstimator::FT0C:
513- if (harmonic == 2 ) {
538+ if (harmonic == kElliptic ) {
514539 xQVec = collision.q2xft0c ();
515540 yQVec = collision.q2yft0c ();
516- } else if (harmonic == 3 ) {
541+ } else if (harmonic == kTriangluar ) {
517542 xQVec = collision.q3xft0c ();
518543 yQVec = collision.q3yft0c ();
519544 }
520545 break ;
521546 case QvecEstimator::TPCPos:
522- if (harmonic == 2 ) {
547+ if (harmonic == kElliptic ) {
523548 xQVec = collision.q2xbpos ();
524549 yQVec = collision.q2ybpos ();
525- } else if (harmonic == 3 ) {
550+ } else if (harmonic == kTriangluar ) {
526551 xQVec = collision.q3xbpos ();
527552 yQVec = collision.q3ybpos ();
528553 }
529554 break ;
530555 case QvecEstimator::TPCNeg:
531- if (harmonic == 2 ) {
556+ if (harmonic == kElliptic ) {
532557 xQVec = collision.q2xbneg ();
533558 yQVec = collision.q2ybneg ();
534- } else if (harmonic == 3 ) {
559+ } else if (harmonic == kTriangluar ) {
535560 xQVec = collision.q3xbneg ();
536561 yQVec = collision.q3ybneg ();
537562 }
538563 break ;
539564 case QvecEstimator::TPCTot:
540- if (harmonic == 2 ) {
565+ if (harmonic == kElliptic ) {
541566 xQVec = collision.q2xbtot ();
542567 yQVec = collision.q2ybtot ();
543- } else if (harmonic == 3 ) {
568+ } else if (harmonic == kTriangluar ) {
544569 xQVec = collision.q3xbtot ();
545570 yQVec = collision.q3ybtot ();
546571 }
547572 break ;
548573 default :
549574 LOG (warning) << " Q vector estimator not valid. Falling back to FT0M" ;
550- if (harmonic == 2 ) {
575+ if (harmonic == kElliptic ) {
551576 xQVec = collision.q2xft0m ();
552577 yQVec = collision.q2yft0m ();
553- } else if (harmonic == 3 ) {
578+ } else if (harmonic == kTriangluar ) {
554579 xQVec = collision.q3xft0m ();
555580 yQVec = collision.q3yft0m ();
556581 }
@@ -565,7 +590,7 @@ struct TaskPi0FlowEMC {
565590 {
566591 bool isgood = true ;
567592 for (const auto & QVec : QVecs) {
568- if (std::fabs (QVec) > 20 . f ) {
593+ if (std::fabs (QVec) > cfgMaxQVector ) {
569594 isgood = false ;
570595 break ;
571596 }
@@ -660,7 +685,7 @@ struct TaskPi0FlowEMC {
660685 void rotationBackground (const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const & photons_coll, unsigned int ig1, unsigned int ig2, CollsWithQvecs::iterator const & collision)
661686 {
662687 // if less than 3 clusters are present skip event since we need at least 3 clusters
663- if (photons_coll.size () < 3 ) {
688+ if (photons_coll.size () < NMinPhotonRotBkg ) {
664689 return ;
665690 }
666691
@@ -759,7 +784,7 @@ struct TaskPi0FlowEMC {
759784 void rotationBackgroundCalib (const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const & photons_coll, unsigned int ig1, unsigned int ig2, CollsWithQvecs::iterator const & collision)
760785 {
761786 // if less than 3 clusters are present skip event since we need at least 3 clusters
762- if (photons_coll.size () < 3 ) {
787+ if (photons_coll.size () < NMinPhotonRotBkg ) {
763788 return ;
764789 }
765790 float cent = getCentrality (collision);
@@ -794,7 +819,7 @@ struct TaskPi0FlowEMC {
794819 }
795820 ROOT::Math::PtEtaPhiMVector photon3 (photon.pt (), photon.eta (), photon.phi (), 0 .);
796821 if (iCellIDPhoton1 >= 0 ) {
797- if (std::fabs ((photon1.E () - photon3.E ()) / (photon1.E () + photon3.E ()) < 0.1 )) { // only use symmetric decays
822+ if (std::fabs ((photon1.E () - photon3.E ()) / (photon1.E () + photon3.E ()) < cfgMaxAsymmetry )) { // only use symmetric decays
798823 ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3;
799824 float openingAngle1 = std::acos (photon1.Vect ().Dot (photon3.Vect ()) / (photon1.P () * photon3.P ()));
800825
@@ -812,7 +837,7 @@ struct TaskPi0FlowEMC {
812837 }
813838 }
814839 if (iCellIDPhoton2 >= 0 ) {
815- if (std::fabs ((photon2.E () - photon3.E ()) / (photon2.E () + photon3.E ()) < 0.1 )) { // only use symmetric decays
840+ if (std::fabs ((photon2.E () - photon3.E ()) / (photon2.E () + photon3.E ()) < cfgMaxAsymmetry )) { // only use symmetric decays
816841 ROOT::Math::PtEtaPhiMVector mother2 = photon2 + photon3;
817842 float openingAngle2 = std::acos (photon2.Vect ().Dot (photon3.Vect ()) / (photon2.P () * photon3.P ()));
818843
@@ -1138,7 +1163,7 @@ struct TaskPi0FlowEMC {
11381163 float yQVecBNeg = -999 .f ;
11391164 float xQVecBTot = -999 .f ;
11401165 float yQVecBTot = -999 .f ;
1141- if (harmonic == 2 ) {
1166+ if (harmonic == kElliptic ) {
11421167 xQVecFT0a = collision.q2xft0a ();
11431168 yQVecFT0a = collision.q2yft0a ();
11441169 xQVecFT0c = collision.q2xft0c ();
@@ -1151,7 +1176,7 @@ struct TaskPi0FlowEMC {
11511176 yQVecBNeg = collision.q2ybneg ();
11521177 xQVecBTot = collision.q2xbtot ();
11531178 yQVecBTot = collision.q2ybtot ();
1154- } else if (harmonic == 3 ) {
1179+ } else if (harmonic == kTriangluar ) {
11551180 xQVecFT0a = collision.q3xft0a ();
11561181 yQVecFT0a = collision.q3yft0a ();
11571182 xQVecFT0c = collision.q3xft0c ();
@@ -1203,7 +1228,7 @@ struct TaskPi0FlowEMC {
12031228 registry.fill (HIST (" epReso/hEpResoFT0mTPCneg" ), centrality, std::cos (harmonic * getDeltaPsiInRange (epFT0m, epBNegs)));
12041229 registry.fill (HIST (" epReso/hEpResoFT0mTPCtot" ), centrality, std::cos (harmonic * getDeltaPsiInRange (epFT0m, epBTots)));
12051230 registry.fill (HIST (" epReso/hEpResoTPCposTPCneg" ), centrality, std::cos (harmonic * getDeltaPsiInRange (epBPoss, epBNegs)));
1206- for (int n = 1 ; n <= 8 ; n++) {
1231+ for (int n = 1 ; n <= kOctagonal ; n++) {
12071232 registry.fill (HIST (" epReso/hEpCosCoefficientsFT0c" ), centrality, n, std::cos (n * epFT0c));
12081233 registry.fill (HIST (" epReso/hEpSinCoefficientsFT0c" ), centrality, n, std::sin (n * epFT0c));
12091234 registry.fill (HIST (" epReso/hEpCosCoefficientsFT0a" ), centrality, n, std::cos (n * epFT0a));
@@ -1323,7 +1348,7 @@ struct TaskPi0FlowEMC {
13231348 registry.fill (HIST (" hClusterCuts" ), 5 );
13241349 continue ;
13251350 }
1326- if (std::fabs ((v1.E () - v2.E ()) / (v1.E () + v2.E ()) < 0.1 )) { // only use symmetric decays
1351+ if (std::fabs ((v1.E () - v2.E ()) / (v1.E () + v2.E ()) < cfgMaxAsymmetry )) { // only use symmetric decays
13271352 registry.fill (HIST (" hClusterCuts" ), 6 );
13281353 registry.fill (HIST (" hSparseCalibSE" ), vMeson.M (), vMeson.E () / 2 ., getCentrality (collision));
13291354 }
0 commit comments