Skip to content

Commit 696da6f

Browse files
mhemmer-cernmhemmer-cern
committed
[PWGEM,PWGEM-36] taskPi0FlowEMC: Add additional EMCal calib
- Add option to use additional calibration for EMCal clusters. This calibration aims to shift the reconstructed Pi0 mass to the PDG value.
1 parent fbdd1f1 commit 696da6f

File tree

1 file changed

+60
-43
lines changed

1 file changed

+60
-43
lines changed

‎PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx‎

Lines changed: 60 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -50,12 +50,14 @@
5050
#include <Math/GenVector/AxisAngle.h>
5151
#include <Math/GenVector/Rotation3D.h>
5252
#include <Math/Vector4D.h> // IWYU pragma: keep
53+
#include <TF1.h>
5354
#include <TH1.h>
5455
#include <TString.h>
5556

5657
#include <array>
5758
#include <cmath>
5859
#include <cstdint>
60+
#include <memory>
5961
#include <numbers>
6062
#include <string>
6163
#include <string_view>
@@ -100,6 +102,8 @@ enum Harmonics {
100102
};
101103

102104
struct TaskPi0FlowEMC {
105+
static constexpr float MinEnergy = 0.7f;
106+
103107
// configurable for flow
104108
Configurable<int> harmonic{"harmonic", 2, "harmonic number"};
105109
Configurable<int> qvecDetector{"qvecDetector", 0, "Detector for Q vector estimation (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5)"};
@@ -129,6 +133,7 @@ struct TaskPi0FlowEMC {
129133
ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {8, -1., 1.}, ""};
130134
ConfigurableAxis thnConfigAxisScalarProd{"thnConfigAxisScalarProd", {100, -5., 5.}, ""};
131135
ConfigurableAxis thnConfigAxisM02{"thnConfigAxisM02", {200, 0., 5.}, ""};
136+
ConfigurableAxis thnConfigAxisEnergyCalib{"thnConfigAxisEnergyCalib", {200, 0., 20.}, ""};
132137

133138
EMPhotonEventCut fEMEventCut;
134139
struct : ConfigurableGroup {
@@ -195,6 +200,7 @@ struct TaskPi0FlowEMC {
195200
Configurable<std::string> cfgSpresoPath{"cfgSpresoPath", "Users/m/mhemmer/EM/Flow/Resolution", "Path to SP resolution file"};
196201
Configurable<int> cfgApplySPresolution{"cfgApplySPresolution", 0, "Apply resolution correction"};
197202
Configurable<bool> doEMCalCalib{"doEMCalCalib", 0, "Produce output for EMCal calibration"};
203+
Configurable<bool> cfgEnableNonLin{"cfgEnableNonLin", false, "flag to turn extra non linear energy calibration on/off"};
198204
} correctionConfig;
199205

200206
SliceCache cache;
@@ -209,6 +215,7 @@ struct TaskPi0FlowEMC {
209215
using EMCalPhotons = soa::Join<aod::EMCEMEventIds, aod::SkimEMCClusters>;
210216
using FilteredCollsWithQvecs = soa::Filtered<soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent, aod::EMEventsQvec>>;
211217
using CollsWithQvecs = soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent, aod::EMEventsQvec>;
218+
using Colls = soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent>;
212219

213220
Preslice<EMCalPhotons> perCollisionEMC = aod::emccluster::emeventId;
214221

@@ -230,6 +237,9 @@ struct TaskPi0FlowEMC {
230237
// static constexpr
231238
static constexpr int64_t NMinPhotonRotBkg = 3;
232239

240+
// Usage when cfgEnableNonLin is enabled
241+
std::unique_ptr<TF1> fEMCalCorrectionFactor; // ("fEMCalCorrectionFactor","(1 + [0]/x + [1]/x^2) / (1 + [2]/x)", 0.3, 100.);
242+
233243
// To access the 1D array
234244
inline int getIndex(int iEta, int iPhi)
235245
{
@@ -308,10 +318,10 @@ struct TaskPi0FlowEMC {
308318
const AxisSpec thnAxisM02{thnConfigAxisM02, "M_{02}"};
309319
const AxisSpec thAxisTanThetaPhi{mesonConfig.thConfigAxisTanThetaPhi, "atan(#Delta#theta/#Delta#varphi)"};
310320
const AxisSpec thAxisClusterEnergy{thnConfigAxisPt, "#it{E} (GeV)"};
321+
const AxisSpec thAxisEnergyCalib{thnConfigAxisEnergyCalib, "#it{E}_{clus} (GeV)"};
311322
const AxisSpec thAxisAlpha{100, -1., +1, "#alpha"};
312323
const AxisSpec thAxisMult{1000, 0., +1000, "#it{N}_{ch}"};
313324
const AxisSpec thAxisEnergy{1000, 0., 100., "#it{E}_{clus} (GeV)"};
314-
const AxisSpec thAxisEnergyCalib{400, 0., 20., "#it{E}_{clus} (GeV)"};
315325
const AxisSpec thAxisTime{1500, -600, 900, "#it{t}_{cl} (ns)"};
316326
const AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"};
317327
const AxisSpec thAxisPhi{500, 0, 2 * 3.14159, "phi"};
@@ -439,6 +449,9 @@ struct TaskPi0FlowEMC {
439449

440450
LOG(info) << "thnConfigAxisInvMass.value[1] = " << thnConfigAxisInvMass.value[1] << " thnConfigAxisInvMass.value.back() = " << thnConfigAxisInvMass.value.back();
441451
LOG(info) << "thnConfigAxisPt.value[1] = " << thnConfigAxisPt.value[1] << " thnConfigAxisPt.value.back() = " << thnConfigAxisPt.value.back();
452+
453+
fEMCalCorrectionFactor = std::make_unique<TF1>("fEMCalCorrectionFactor", "(1 + [0]/x + [1]/x^2) / (1 + [2]/x)", 0.3, 100.);
454+
fEMCalCorrectionFactor->SetParameters(-5.33426e-01, 1.40144e-02, -5.24434e-01);
442455
}; // end init
443456

444457
/// Change radians to degree
@@ -475,7 +488,8 @@ struct TaskPi0FlowEMC {
475488

476489
/// Get the centrality
477490
/// \param collision is the collision with the centrality information
478-
float getCentrality(CollsWithQvecs::iterator const& collision)
491+
template <typename TCollision>
492+
float getCentrality(TCollision const& collision)
479493
{
480494
float cent = -999.;
481495
switch (centEstimator) {
@@ -728,7 +742,11 @@ struct TaskPi0FlowEMC {
728742
if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelBackground.value) {
729743
continue;
730744
}
731-
ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.);
745+
float energyCorrectionFactor = 1.f;
746+
if (correctionConfig.cfgEnableNonLin.value) {
747+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(photon.e() > MinEnergy ? photon.e() : MinEnergy);
748+
}
749+
ROOT::Math::PtEtaPhiMVector photon3(energyCorrectionFactor * photon.pt(), photon.eta(), photon.phi(), 0.);
732750
if (iCellIDPhoton1 >= 0) {
733751
ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3;
734752
float openingAngle1 = std::acos(photon1.Vect().Dot(photon3.Vect()) / (photon1.P() * photon3.P()));
@@ -780,8 +798,8 @@ struct TaskPi0FlowEMC {
780798
}
781799

782800
/// \brief Calculate background using rotation background method for calib
783-
template <typename TPhotons>
784-
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)
801+
template <typename TPhotons, typename TCollision>
802+
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, TCollision const& collision)
785803
{
786804
// if less than 3 clusters are present skip event since we need at least 3 clusters
787805
if (photons_coll.size() < NMinPhotonRotBkg) {
@@ -817,7 +835,11 @@ struct TaskPi0FlowEMC {
817835
if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelBackground.value) {
818836
continue;
819837
}
820-
ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.);
838+
float energyCorrectionFactor = 1.f;
839+
if (correctionConfig.cfgEnableNonLin.value) {
840+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(photon.e() > MinEnergy ? photon.e() : MinEnergy);
841+
}
842+
ROOT::Math::PtEtaPhiMVector photon3(energyCorrectionFactor * photon.pt(), photon.eta(), photon.phi(), 0.);
821843
if (iCellIDPhoton1 >= 0) {
822844
if (std::fabs((photon1.E() - photon3.E()) / (photon1.E() + photon3.E()) < cfgMaxAsymmetry)) { // only use symmetric decays
823845
ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3;
@@ -893,6 +915,10 @@ struct TaskPi0FlowEMC {
893915
void processEMCal(CollsWithQvecs const& collisions, EMCalPhotons const& clusters)
894916
{
895917
int nColl = 1;
918+
float energyCorrectionFactor = 1.f;
919+
if (cfgDoReverseScaling.value) {
920+
energyCorrectionFactor = 1.0505f;
921+
}
896922
for (const auto& collision : collisions) {
897923
auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex());
898924

@@ -971,18 +997,14 @@ struct TaskPi0FlowEMC {
971997
continue;
972998
}
973999
}
974-
975-
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
976-
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
977-
if (cfgDoReverseScaling.value) {
978-
// Convert to PxPyPzEVector to modify energy
979-
ROOT::Math::PxPyPzEVector v1Mod(v1);
980-
v1Mod.SetE(v1Mod.E() * 1.0505);
981-
v1 = ROOT::Math::PtEtaPhiMVector(v1Mod.Pt(), v1Mod.Eta(), v1Mod.Phi(), 0.);
982-
ROOT::Math::PxPyPzEVector v2Mod(v2);
983-
v2Mod.SetE(v2Mod.E() * 1.0505);
984-
v2 = ROOT::Math::PtEtaPhiMVector(v2Mod.Pt(), v2Mod.Eta(), v2Mod.Phi(), 0.);
1000+
if (correctionConfig.cfgEnableNonLin.value) {
1001+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy);
9851002
}
1003+
ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.);
1004+
if (correctionConfig.cfgEnableNonLin.value) {
1005+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy);
1006+
}
1007+
ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.);
9861008
ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2;
9871009
float dTheta = v1.Theta() - v2.Theta();
9881010
float dPhi = v1.Phi() - v2.Phi();
@@ -1031,6 +1053,10 @@ struct TaskPi0FlowEMC {
10311053
// Pi0 from EMCal
10321054
void processEMCalMixed(FilteredCollsWithQvecs const& collisions, FilteredEMCalPhotons const& clusters)
10331055
{
1056+
float energyCorrectionFactor = 1.f;
1057+
if (cfgDoReverseScaling.value) {
1058+
energyCorrectionFactor = 1.0505f;
1059+
}
10341060
auto getClustersSize =
10351061
[&clusters, this](FilteredCollsWithQvecs::iterator const& col) {
10361062
auto associatedClusters = clusters.sliceByCached(emccluster::emeventId, col.globalIndex(), this->cache); // it's cached, so slicing/grouping happens only once
@@ -1079,18 +1105,14 @@ struct TaskPi0FlowEMC {
10791105
continue;
10801106
}
10811107
}
1082-
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
1083-
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
1084-
1085-
if (cfgDoReverseScaling.value) {
1086-
// Convert to PxPyPzEVector to modify energy
1087-
ROOT::Math::PxPyPzEVector v1Mod(v1);
1088-
v1Mod.SetE(v1Mod.E() * 1.0505);
1089-
v1 = ROOT::Math::PtEtaPhiMVector(v1Mod.Pt(), v1Mod.Eta(), v1Mod.Phi(), 0.);
1090-
ROOT::Math::PxPyPzEVector v2Mod(v2);
1091-
v2Mod.SetE(v2Mod.E() * 1.0505);
1092-
v2 = ROOT::Math::PtEtaPhiMVector(v2Mod.Pt(), v2Mod.Eta(), v2Mod.Phi(), 0.);
1108+
if (correctionConfig.cfgEnableNonLin.value) {
1109+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy);
1110+
}
1111+
ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.);
1112+
if (correctionConfig.cfgEnableNonLin.value) {
1113+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy);
10931114
}
1115+
ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.);
10941116
ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2;
10951117

10961118
float dTheta = v1.Theta() - v2.Theta();
@@ -1247,8 +1269,9 @@ struct TaskPi0FlowEMC {
12471269
PROCESS_SWITCH(TaskPi0FlowEMC, processResolution, "Process resolution", false);
12481270

12491271
// EMCal calibration
1250-
void processEMCalCalib(CollsWithQvecs const& collisions, EMCalPhotons const& clusters)
1272+
void processEMCalCalib(Colls const& collisions, EMCalPhotons const& clusters)
12511273
{
1274+
float energyCorrectionFactor = 1.f;
12521275
if (!correctionConfig.doEMCalCalib) {
12531276
return;
12541277
}
@@ -1269,10 +1292,6 @@ struct TaskPi0FlowEMC {
12691292
// event selection
12701293
continue;
12711294
}
1272-
if (!isQvecGood(getAllQvec(collision))) {
1273-
// selection based on QVector
1274-
continue;
1275-
}
12761295
runNow = collision.runNumber();
12771296
if (runNow != runBefore) {
12781297
initCCDB(collision);
@@ -1306,18 +1325,16 @@ struct TaskPi0FlowEMC {
13061325
}
13071326
}
13081327

1309-
ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
1310-
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
1311-
if (cfgDoReverseScaling.value) {
1312-
// Convert to PxPyPzEVector to modify energy
1313-
ROOT::Math::PxPyPzEVector v1Mod(v1);
1314-
v1Mod.SetE(v1Mod.E() * 1.0505);
1315-
v1 = ROOT::Math::PtEtaPhiMVector(v1Mod.Pt(), v1Mod.Eta(), v1Mod.Phi(), 0.);
1316-
ROOT::Math::PxPyPzEVector v2Mod(v2);
1317-
v2Mod.SetE(v2Mod.E() * 1.0505);
1318-
v2 = ROOT::Math::PtEtaPhiMVector(v2Mod.Pt(), v2Mod.Eta(), v2Mod.Phi(), 0.);
1328+
if (correctionConfig.cfgEnableNonLin.value) {
1329+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy);
13191330
}
1331+
ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.);
1332+
if (correctionConfig.cfgEnableNonLin.value) {
1333+
energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy);
1334+
}
1335+
ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.);
13201336
ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2;
1337+
13211338
float dTheta = v1.Theta() - v2.Theta();
13221339
float dPhi = v1.Phi() - v2.Phi();
13231340
float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P()));

0 commit comments

Comments
 (0)