Skip to content

Commit a9d8a1e

Browse files
DanielSamitzDanielSamitz
authored
add efficiency to smearing task (#5967)
1 parent acab90c commit a9d8a1e

3 files changed

Lines changed: 136 additions & 9 deletions

File tree

‎PWGDQ/DataModel/ReducedInfoTables.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -458,10 +458,11 @@ namespace smearedtrack
458458
DECLARE_SOA_COLUMN(PtSmeared, ptSmeared, float);
459459
DECLARE_SOA_COLUMN(EtaSmeared, etaSmeared, float);
460460
DECLARE_SOA_COLUMN(PhiSmeared, phiSmeared, float);
461+
DECLARE_SOA_COLUMN(Efficiency, efficiency, float);
461462
} // namespace smearedtrack
462463

463464
DECLARE_SOA_TABLE(SmearedTracks, "AOD", "SMEAREDTRACK", // use like this Join<ReducedMCTracks, SmearedTracks>
464-
smearedtrack::PtSmeared, smearedtrack::EtaSmeared, smearedtrack::PhiSmeared);
465+
smearedtrack::PtSmeared, smearedtrack::EtaSmeared, smearedtrack::PhiSmeared, smearedtrack::Efficiency);
465466
using SmearedTrack = SmearedTracks::iterator;
466467

467468
namespace dilepton_track_index

‎PWGEM/Dilepton/Tasks/smearing.cxx‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,8 @@ struct ApplySmearing {
3939
Configurable<std::string> fConfigResEtaHistName{"cfgResEtaHistName", "EtaResArr", "histogram name for eta in resolution file"};
4040
Configurable<std::string> fConfigResPhiPosHistName{"cfgResPhiPosHistName", "PhiPosResArr", "histogram name for phi pos in resolution file"};
4141
Configurable<std::string> fConfigResPhiNegHistName{"cfgResPhiNegHistName", "PhiEleResArr", "hisogram for phi neg in resolution file"};
42+
Configurable<std::string> fConfigEffFileName{"cfgEffFileName", "", "name of efficiency file"};
43+
Configurable<std::string> fConfigEffHistName{"cfgEffHistName", "fhwEffpT", "name of efficiency histogram"};
4244

4345
MomentumSmearer smearer;
4446

@@ -49,6 +51,8 @@ struct ApplySmearing {
4951
smearer.setResEtaHistName(TString(fConfigResEtaHistName));
5052
smearer.setResPhiPosHistName(TString(fConfigResPhiPosHistName));
5153
smearer.setResPhiNegHistName(TString(fConfigResPhiNegHistName));
54+
smearer.setEffFileName(TString(fConfigEffFileName));
55+
smearer.setEffHistName(TString(fConfigEffHistName));
5256
smearer.init();
5357
}
5458

@@ -59,6 +63,7 @@ struct ApplySmearing {
5963
float ptgen = mctrack.pt();
6064
float etagen = mctrack.eta();
6165
float phigen = mctrack.phi();
66+
float efficiency = 1.;
6267

6368
int pdgCode = mctrack.pdgCode();
6469
if (abs(pdgCode) == fPdgCode) {
@@ -69,10 +74,12 @@ struct ApplySmearing {
6974
// apply smearing for electrons or muons.
7075
float ptsmeared, etasmeared, phismeared;
7176
smearer.applySmearing(ch, ptgen, etagen, phigen, ptsmeared, etasmeared, phismeared);
72-
smearedtrack(ptsmeared, etasmeared, phismeared);
77+
// get the efficiency
78+
efficiency = smearer.getEfficiency(ptgen, etagen, phigen);
79+
smearedtrack(ptsmeared, etasmeared, phismeared, efficiency);
7380
} else {
7481
// don't apply smearing
75-
smearedtrack(ptgen, etagen, phigen);
82+
smearedtrack(ptgen, etagen, phigen, efficiency);
7683
}
7784
}
7885
}

‎PWGEM/Dilepton/Utils/MomentumSmearer.h‎

Lines changed: 125 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,21 @@ class MomentumSmearer
3636
setResEtaHistName(resEtaHistName);
3737
setResPhiPosHistName(resPhiPosHistName);
3838
setResPhiNegHistName(resPhiNegHistName);
39+
setResPhiNegHistName("");
40+
setResPhiNegHistName("");
41+
init();
42+
}
43+
44+
/// Constructor with resolution histograms and efficiency
45+
MomentumSmearer(TString resFileName, TString resPtHistName, TString resEtaHistName, TString resPhiPosHistName, TString resPhiNegHistName, TString effFileName, TString effHistName)
46+
{
47+
setResFileName(resFileName);
48+
setResPtHistName(resPtHistName);
49+
setResEtaHistName(resEtaHistName);
50+
setResPhiPosHistName(resPhiPosHistName);
51+
setResPhiNegHistName(resPhiNegHistName);
52+
setResPhiNegHistName(resPhiNegHistName);
53+
setResPhiNegHistName(resPhiNegHistName);
3954
init();
4055
}
4156

@@ -47,7 +62,9 @@ class MomentumSmearer
4762
if (fInitialized)
4863
return;
4964

50-
if (fResFileName.BeginsWith("alien://")) {
65+
// ToDo: make it possible to access .root files from CCDB
66+
67+
if (fResFileName.BeginsWith("alien://") || fEffFileName.BeginsWith("alien://")) {
5168
TGrid::Connect("alien://");
5269
}
5370

@@ -56,35 +73,35 @@ class MomentumSmearer
5673
// Get Resolution map
5774
TFile* fFile = TFile::Open(fResFileName);
5875
if (!fFile) {
59-
LOGP(error, "Could not open Resolution file {}", fResFileName.Data());
76+
LOGP(fatal, "Could not open Resolution file {}", fResFileName.Data());
6077
return;
6178
}
6279
TObjArray* ArrResoPt = nullptr;
6380
if (fFile->GetListOfKeys()->Contains(fResPtHistName)) {
6481
ArrResoPt = reinterpret_cast<TObjArray*>(fFile->Get(fResPtHistName));
6582
} else {
66-
LOGP(error, "Could not open {} from file {}", fResPtHistName.Data(), fResFileName.Data());
83+
LOGP(fatal, "Could not open {} from file {}", fResPtHistName.Data(), fResFileName.Data());
6784
}
6885

6986
TObjArray* ArrResoEta = nullptr;
7087
if (fFile->GetListOfKeys()->Contains(fResEtaHistName)) {
7188
ArrResoEta = reinterpret_cast<TObjArray*>(fFile->Get(fResEtaHistName));
7289
} else {
73-
LOGP(error, "Could not open {} from file {}", fResEtaHistName.Data(), fResFileName.Data());
90+
LOGP(fatal, "Could not open {} from file {}", fResEtaHistName.Data(), fResFileName.Data());
7491
}
7592

7693
TObjArray* ArrResoPhi_Pos = nullptr;
7794
if (fFile->GetListOfKeys()->Contains(TString(fResPhiPosHistName))) {
7895
ArrResoPhi_Pos = reinterpret_cast<TObjArray*>(fFile->Get(fResPhiPosHistName));
7996
} else {
80-
LOGP(error, "Could not open {} from file {}", fResPhiPosHistName.Data(), fResFileName.Data());
97+
LOGP(fatal, "Could not open {} from file {}", fResPhiPosHistName.Data(), fResFileName.Data());
8198
}
8299

83100
TObjArray* ArrResoPhi_Neg = nullptr;
84101
if (fFile->GetListOfKeys()->Contains(TString(fResPhiNegHistName))) {
85102
ArrResoPhi_Neg = reinterpret_cast<TObjArray*>(fFile->Get(fResPhiNegHistName));
86103
} else {
87-
LOGP(error, "Could not open {} from file {}", fResPhiNegHistName.Data(), fResFileName.Data());
104+
LOGP(fatal, "Could not open {} from file {}", fResPhiNegHistName.Data(), fResFileName.Data());
88105
}
89106

90107
fArrResoPt = ArrResoPt;
@@ -93,6 +110,35 @@ class MomentumSmearer
93110
fArrResoPhi_Neg = ArrResoPhi_Neg;
94111
fFile->Close();
95112

113+
// get efficiency histo
114+
fEffType = 0;
115+
if (fEffFileName.CompareTo("") != 0) {
116+
LOGP(info, "Set Resolution histo");
117+
TFile* fEffFile = TFile::Open(fEffFileName);
118+
if (!fEffFile) {
119+
LOGP(fatal, "Could not open efficiency file {}", fEffFileName.Data());
120+
return;
121+
}
122+
if (fEffFile->GetListOfKeys()->Contains(fEffHistName.Data())) {
123+
fArrEff = reinterpret_cast<TObject*>(fEffFile->Get(fEffHistName.Data()));
124+
// check which type is used
125+
if (dynamic_cast<TH3*>(fArrEff)) {
126+
fEffType = 3;
127+
LOGP(info, "Use 3d efficiency histo (pt, eta, phi)");
128+
} else if (dynamic_cast<TH2*>(fArrEff)) {
129+
fEffType = 2;
130+
LOGP(info, "Use 2d efficiency histo (pt, eta)");
131+
} else if (dynamic_cast<TH1*>(fArrEff)) {
132+
fEffType = 1;
133+
LOGP(info, "Use 1d efficiency histo (pt)");
134+
} else {
135+
LOGP(fatal, "Could not identify type of histogram {}", fEffHistName.Data());
136+
}
137+
} else {
138+
LOGP(fatal, "Could not find histogram {} in file {}", fEffHistName.Data(), fEffFileName.Data());
139+
}
140+
}
141+
96142
fInitialized = true;
97143
}
98144

@@ -148,23 +194,92 @@ class MomentumSmearer
148194
phismeared = phigen - smearing;
149195
}
150196

197+
float getEfficiency(float pt, float eta, float phi)
198+
{
199+
200+
if (fEffType == 0) {
201+
return 1.;
202+
}
203+
204+
if (fEffType == 1) {
205+
TH1F* hist = reinterpret_cast<TH1F*>(fArrEff);
206+
int ptbin = hist->GetXaxis()->FindBin(pt);
207+
int ptbin_max = hist->GetXaxis()->GetNbins();
208+
// make sure that no underflow or overflow bins are used
209+
if (ptbin < 1)
210+
ptbin = 1;
211+
else if (ptbin > ptbin_max)
212+
ptbin = ptbin_max;
213+
return hist->GetBinContent(ptbin);
214+
}
215+
216+
if (fEffType == 2) {
217+
TH2F* hist = reinterpret_cast<TH2F*>(fArrEff);
218+
int ptbin = hist->GetXaxis()->FindBin(pt);
219+
int ptbin_max = hist->GetXaxis()->GetNbins();
220+
int etabin = hist->GetYaxis()->FindBin(eta);
221+
int etabin_max = hist->GetYaxis()->GetNbins();
222+
// make sure that no underflow or overflow bins are used
223+
if (ptbin < 1)
224+
ptbin = 1;
225+
else if (ptbin > ptbin_max)
226+
ptbin = ptbin_max;
227+
if (etabin < 1)
228+
etabin = 1;
229+
else if (etabin > etabin_max)
230+
etabin = etabin_max;
231+
return hist->GetBinContent(ptbin, etabin);
232+
}
233+
234+
if (fEffType == 3) {
235+
TH3F* hist = reinterpret_cast<TH3F*>(fArrEff);
236+
int ptbin = hist->GetXaxis()->FindBin(pt);
237+
int ptbin_max = hist->GetXaxis()->GetNbins();
238+
int etabin = hist->GetYaxis()->FindBin(eta);
239+
int etabin_max = hist->GetYaxis()->GetNbins();
240+
int phibin = hist->GetZaxis()->FindBin(phi);
241+
int phibin_max = hist->GetZaxis()->GetNbins();
242+
// make sure that no underflow or overflow bins are used
243+
if (ptbin < 1)
244+
ptbin = 1;
245+
else if (ptbin > ptbin_max)
246+
ptbin = ptbin_max;
247+
if (etabin < 1)
248+
etabin = 1;
249+
else if (etabin > etabin_max)
250+
etabin = etabin_max;
251+
if (phibin < 1)
252+
phibin = 1;
253+
else if (phibin > phibin_max)
254+
phibin = phibin_max;
255+
return hist->GetBinContent(ptbin, etabin, phibin);
256+
}
257+
258+
return 1.;
259+
}
260+
151261
// setters
152262
void setResFileName(TString resFileName) { fResFileName = resFileName; }
153263
void setResPtHistName(TString resPtHistName) { fResPtHistName = resPtHistName; }
154264
void setResEtaHistName(TString resEtaHistName) { fResEtaHistName = resEtaHistName; }
155265
void setResPhiPosHistName(TString resPhiPosHistName) { fResPhiPosHistName = resPhiPosHistName; }
156266
void setResPhiNegHistName(TString resPhiNegHistName) { fResPhiNegHistName = resPhiNegHistName; }
267+
void setEffFileName(TString effFileName) { fEffFileName = effFileName; }
268+
void setEffHistName(TString effHistName) { fEffHistName = effHistName; }
157269

158270
// getters
159271
TString getResFileName() { return fResFileName; }
160272
TString getResPtHistName() { return fResPtHistName; }
161273
TString getResEtaHistName() { return fResEtaHistName; }
162274
TString getResPhiPosHistName() { return fResPhiPosHistName; }
163275
TString getResPhiNegHistName() { return fResPhiNegHistName; }
276+
TString getEffFileName() { return fEffFileName; }
277+
TString getEffHistName() { return fEffHistName; }
164278
TObjArray* getArrResoPt() { return fArrResoPt; }
165279
TObjArray* getArrResoEta() { return fArrResoEta; }
166280
TObjArray* getArrResoPhiPos() { return fArrResoPhi_Pos; }
167281
TObjArray* getArrResoPhiNeg() { return fArrResoPhi_Neg; }
282+
TObject* getArrEff() { return fArrEff; }
168283

169284
private:
170285
bool fInitialized = false;
@@ -173,10 +288,14 @@ class MomentumSmearer
173288
TString fResEtaHistName;
174289
TString fResPhiPosHistName;
175290
TString fResPhiNegHistName;
291+
TString fEffFileName;
292+
TString fEffHistName;
293+
int fEffType = 0;
176294
TObjArray* fArrResoPt;
177295
TObjArray* fArrResoEta;
178296
TObjArray* fArrResoPhi_Pos;
179297
TObjArray* fArrResoPhi_Neg;
298+
TObject* fArrEff;
180299
};
181300

182301
#endif // PWGEM_DILEPTON_UTILS_MOMENTUMSMEARER_H_

0 commit comments

Comments
 (0)