Skip to content

Commit 22630c8

Browse files
joachimckhjoachimckh
authored
[PWGJE] Updates to process task for MC: bkg RM + occupancy cuts + var (#9481)
1 parent 674f29d commit 22630c8

File tree

1 file changed

+124
-54
lines changed

1 file changed

+124
-54
lines changed

‎PWGJE/Tasks/jetSpectraEseTask.cxx‎

Lines changed: 124 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ using namespace o2::framework::expressions;
4545

4646
struct JetSpectraEseTask {
4747
ConfigurableAxis binJetPt{"binJetPt", {200, 0., 200.}, ""};
48-
ConfigurableAxis bindPhi{"bindPhi", {100, -o2::constants::math::PI - 1, o2::constants::math::PI + 1}, ""};
48+
ConfigurableAxis bindPhi{"bindPhi", {180, -o2::constants::math::PI, o2::constants::math::PI}, ""};
4949
ConfigurableAxis binESE{"binESE", {100, 0, 100}, ""};
5050
ConfigurableAxis binCos{"binCos", {100, -1.05, 1.05}, ""};
5151
ConfigurableAxis binOccupancy{"binOccupancy", {5000, 0, 25000}, ""};
@@ -57,13 +57,13 @@ struct JetSpectraEseTask {
5757
Configurable<std::vector<float>> centRange{"centRange", {30, 50}, "centrality region of interest"};
5858
Configurable<double> leadingJetPtCut{"leadingJetPtCut", 5.0, "leading jet pT cut"};
5959

60-
Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
60+
Configurable<std::string> eventSelections{"eventSelections", "sel8FullPbPb", "choose event selection"};
6161
Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};
6262

6363
Configurable<int> fColSwitch{"fColSwitch", 0, "collision switch"};
6464

6565
Configurable<bool> cfgEvSelOccupancy{"cfgEvSelOccupancy", false, "Flag for occupancy cut"};
66-
Configurable<std::vector<int>> cfgCutOccupancy{"cfgCutOccupancy", {0, 500}, "Occupancy cut"};
66+
Configurable<std::vector<int>> cfgCutOccupancy{"cfgCutOccupancy", {0, 1000}, "Occupancy cut"};
6767
Configurable<std::vector<float>> cfgOccupancyPtCut{"cfgOccupancyPtCut", {0, 100}, "pT cut"};
6868

6969
Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "total qvector number // look in Qvector table for this number"};
@@ -82,8 +82,11 @@ struct JetSpectraEseTask {
8282

8383
HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
8484

85-
int eventSelection = -1;
86-
int trackSelection = -1;
85+
int eventSelection{-1};
86+
int trackSelection{-1};
87+
88+
using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>>;
89+
Preslice<ChargedMCDJets> mcdjetsPerJCollision = o2::aod::jet::collisionId;
8790

8891
enum class DetID { FT0C,
8992
FT0A,
@@ -104,6 +107,7 @@ struct JetSpectraEseTask {
104107
case 0:
105108
LOGF(info, "JetSpectraEseTask::init() - using data");
106109
registry.add("hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
110+
registry.add("hCentralityMult", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
107111
registry.add("hJetPt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
108112
registry.add("hJetPt_bkgsub", "jet pT background sub;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
109113
registry.add("hJetEta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
@@ -130,23 +134,34 @@ struct JetSpectraEseTask {
130134
registry.addClone("hPsi2FT0C", "hEPUncorV2");
131135
registry.addClone("hPsi2FT0C", "hEPRectrV2");
132136
registry.addClone("hPsi2FT0C", "hEPTwistV2");
137+
133138
break;
134139
case 1:
135140
LOGF(info, "JetSpectraEseTask::init() - using MC");
136141
registry.add("hMCEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
142+
registry.add("hMCPureEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
143+
registry.add("hMCDMatchedEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
144+
registry.add("hCentralityMult", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
137145
registry.add("hPartJetPt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
146+
registry.add("hPartJetPtSubBkg", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
138147
registry.add("hPartJetEta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
139148
registry.add("hPartJetPhi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});
140149
registry.add("hPartJetPtMatch", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
150+
registry.add("hPartJetPtMatchSubBkg", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
141151
registry.add("hPartJetEtaMatch", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
142152
registry.add("hPartJetPhiMatch", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});
143153
registry.add("hDetectorJetPt", "detector level jet pT;#it{p}_{T,jet det} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
154+
registry.add("hDetectorJetPtSubBkg", "detector level jet pT;#it{p}_{T,jet det} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
144155
registry.add("hDetectorJetEta", "detector level jet #eta;#eta_{jet det};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
145156
registry.add("hDetectorJetPhi", "detector level jet #phi;#phi_{jet det};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});
146157
registry.add("hMatchedJetsPtDelta", "#it{p}_{T,jet part}; det - part", {HistType::kTH2F, {{jetPtAxis}, {200, -20., 20.0}}});
147-
registry.add("hMatchedJetsEtaDelta", "#eta_{jet part}; det - part", {HistType::kTH2F, {{100, -1.0, 1.0}, {200, -20.0, 20.0}}});
148-
registry.add("hMatchedJetsPhiDelta", "#phi_{jet part}; det - part", {HistType::kTH2F, {{80, -1.0, 7.}, {200, -20.0, 20.}}});
149-
registry.add("hResponseMatrixMatch", "#it{p}_{T, jet det}; #it{p}_{T, jet part}", HistType::kTH2F, {jetPtAxis, jetPtAxis});
158+
registry.add("hMatchedJetsPtDeltaSubBkg", "#it{p}_{T,jet part}; det - part", {HistType::kTH2F, {{jetPtAxis}, {200, -20., 20.0}}});
159+
registry.add("hGenMatchedJetsPtDeltadPt", "#it{p}_{T,jet part}; det - part / part", {HistType::kTHnSparseF, {{100, 0, 100}, {jetPtAxis}, {200, -20., 20.0}}});
160+
registry.add("hRecoMatchedJetsPtDeltadPt", "#it{p}_{T,jet det}; det - part / det", {HistType::kTHnSparseF, {{100, 0, 100}, {jetPtAxis}, {200, -20., 20.0}}});
161+
registry.add("hMatchedJetsEtaDelta", "#eta_{jet part}; det - part", {HistType::kTH2F, {{100, -1.0, 1.0}, {200, -0.8, 0.8}}});
162+
registry.add("hMatchedJetsPhiDelta", "#phi_{jet part}; det - part", {HistType::kTH2F, {{80, -1.0, 7.}, {200, -10.0, 10.}}});
163+
registry.add("hRespMcDMcPMatch", ";Centrality,#it{p}_{T, jet det}; #it{p}_{T, jet part}", HistType::kTHnSparseF, {{100, 0, 100}, jetPtAxis, jetPtAxis});
164+
registry.add("hRespMcDMcPMatchSubBkg", ";Centrality,#it{p}_{T, jet det}; #it{p}_{T, jet part}", HistType::kTHnSparseF, {{100, 0, 100}, jetPtAxis, jetPtAxis});
150165
break;
151166
case 2:
152167
LOGF(info, "JetSpectraEseTask::init() - using Occupancy processing");
@@ -170,31 +185,26 @@ struct JetSpectraEseTask {
170185
{
171186
float counter{0.5f};
172187
registry.fill(HIST("hEventCounter"), counter++);
188+
if (!jetderiveddatautilities::selectCollision(collision, eventSelection))
189+
return;
173190
registry.fill(HIST("hEventCounter"), counter++);
174191

175-
const auto vPsi2 = procEP<true>(collision);
176-
const auto qPerc = collision.qPERCFT0C();
177-
if (qPerc[0] < 0)
192+
if (cfgEvSelOccupancy && !isOccupancyWithin(collision))
178193
return;
179194
registry.fill(HIST("hEventCounter"), counter++);
180195

181-
if (!jetderiveddatautilities::selectCollision(collision, eventSelection))
196+
const auto vPsi2{procEP<true>(collision)};
197+
const auto qPerc{collision.qPERCFT0C()};
198+
if (qPerc[0] < 0)
182199
return;
183-
184200
registry.fill(HIST("hEventCounter"), counter++);
185201

186202
if (!isAcceptedLeadTrack(tracks))
187203
return;
188204

189-
if (cfgEvSelOccupancy) {
190-
auto occupancy = collision.trackOccupancyInTimeRange();
191-
if (occupancy < cfgCutOccupancy->at(0) || occupancy > cfgCutOccupancy->at(1))
192-
registry.fill(HIST("hEventCounter"), counter++);
193-
return;
194-
}
195-
196205
registry.fill(HIST("hEventCounter"), counter++);
197206
registry.fill(HIST("hRho"), collision.rho());
207+
registry.fill(HIST("hCentralityMult"), collision.centrality());
198208
for (auto const& jet : jets) {
199209
float jetpTbkgsub = jet.pt() - (collision.rho() * jet.area());
200210
registry.fill(HIST("hJetPt"), jet.pt());
@@ -203,13 +213,13 @@ struct JetSpectraEseTask {
203213
registry.fill(HIST("hJetPhi"), jet.phi());
204214
registry.fill(HIST("hJetArea"), jet.area());
205215

206-
float dPhi = RecoDecay::constrainAngle(jet.phi() - vPsi2, -o2::constants::math::PI);
216+
float dPhi{RecoDecay::constrainAngle(jet.phi() - vPsi2, -o2::constants::math::PI)};
207217
registry.fill(HIST("hdPhi"), dPhi);
208218
registry.fill(HIST("hCentJetPtdPhiq2"), collision.centrality(), jetpTbkgsub, dPhi, qPerc[0]);
209219
}
210220
registry.fill(HIST("hEventCounter"), counter++);
211221

212-
if (collision.centrality() < centRange->at(0) || collision.centrality() > centRange->at(1)) /* for counting */
222+
if (collision.centrality() < centRange->at(0) || collision.centrality() > centRange->at(1))
213223
return;
214224
registry.fill(HIST("hEventCounter"), counter++);
215225
}
@@ -220,10 +230,10 @@ struct JetSpectraEseTask {
220230
{
221231
float count{0.5};
222232
registry.fill(HIST("hEventCounterOcc"), count++);
223-
const auto vPsi2 = procEP<false>(collision);
224-
const auto qPerc = collision.qPERCFT0C();
233+
const auto vPsi2{procEP<false>(collision)};
234+
const auto qPerc{collision.qPERCFT0C()};
225235

226-
auto occupancy = collision.trackOccupancyInTimeRange();
236+
auto occupancy{collision.trackOccupancyInTimeRange()};
227237
registry.fill(HIST("hPsiOccupancy"), collision.centrality(), vPsi2, occupancy);
228238
registry.fill(HIST("hOccupancy"), occupancy);
229239

@@ -244,47 +254,97 @@ struct JetSpectraEseTask {
244254
}
245255
PROCESS_SWITCH(JetSpectraEseTask, processESEOccupancy, "process occupancy", false);
246256

247-
void processMCParticleLevel(soa::Filtered<aod::ChargedMCParticleLevelJets>::iterator const& jet)
257+
void processMCParticleLevel(soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>::iterator const& col,
258+
soa::Filtered<aod::ChargedMCParticleLevelJets> const& jets)
248259
{
249-
registry.fill(HIST("hPartJetPt"), jet.pt());
250-
registry.fill(HIST("hPartJetEta"), jet.eta());
251-
registry.fill(HIST("hPartJetPhi"), jet.phi());
260+
float counter{0.5f};
261+
registry.fill(HIST("hMCPureEventCounter"), counter++);
262+
if (col.size() < 1) {
263+
return;
264+
}
265+
registry.fill(HIST("hMCPureEventCounter"), counter++);
266+
if (!(std::abs(col.posZ()) < vertexZCut)) {
267+
return;
268+
}
269+
registry.fill(HIST("hMCPureEventCounter"), counter++);
270+
271+
for (const auto& jet : jets) {
272+
registry.fill(HIST("hPartJetPt"), jet.pt());
273+
registry.fill(HIST("hPartJetPtSubBkg"), jet.pt() - (col.rho() * jet.area()));
274+
registry.fill(HIST("hPartJetEta"), jet.eta());
275+
registry.fill(HIST("hPartJetPhi"), jet.phi());
276+
}
252277
}
253278
PROCESS_SWITCH(JetSpectraEseTask, processMCParticleLevel, "jets on particle level MC", false);
254279

255280
using JetMCPTable = soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>>;
256-
void processMCChargedMatched(soa::Filtered<aod::JetCollisionsMCD>::iterator const& collision,
257-
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
281+
void processMCChargedMatched(/*soa::Filtered<aod::JetCollisionsMCD>::iterator const& collision*/
282+
soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>::iterator const& mcCol,
283+
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>> const& collisions,
284+
ChargedMCDJets const& mcdjets,
258285
JetMCPTable const&,
259286
aod::JetTracks const&,
260287
aod::JetParticles const&)
261288
{
262-
if (!jetderiveddatautilities::selectCollision(collision, eventSelection))
263-
return;
264-
265289
float counter{0.5f};
266290
registry.fill(HIST("hMCEventCounter"), counter++);
267-
for (const auto& mcdjet : mcdjets) {
268-
269-
registry.fill(HIST("hDetectorJetPt"), mcdjet.pt());
270-
registry.fill(HIST("hDetectorJetEta"), mcdjet.eta());
271-
registry.fill(HIST("hDetectorJetPhi"), mcdjet.phi());
272-
for (const auto& mcpjet : mcdjet.template matchedJetGeo_as<JetMCPTable>()) {
273-
274-
registry.fill(HIST("hPartJetPtMatch"), mcpjet.pt());
275-
registry.fill(HIST("hPartJetEtaMatch"), mcpjet.eta());
276-
registry.fill(HIST("hPartJetPhiMatch"), mcpjet.phi());
291+
if (mcCol.size() < 1) {
292+
return;
293+
}
294+
registry.fill(HIST("hMCEventCounter"), counter++);
295+
if (!(std::abs(mcCol.posZ()) < vertexZCut)) {
296+
return;
297+
}
298+
registry.fill(HIST("hMCEventCounter"), counter++);
277299

278-
registry.fill(HIST("hMatchedJetsPtDelta"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt());
279-
registry.fill(HIST("hMatchedJetsPhiDelta"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi());
280-
registry.fill(HIST("hMatchedJetsEtaDelta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta());
300+
for (const auto& collision : collisions) {
301+
float secCount{0.5f};
302+
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);
303+
if (!(std::abs(collision.posZ()) < vertexZCut)) {
304+
return;
305+
}
281306

282-
registry.fill(HIST("hResponseMatrixMatch"), mcdjet.pt(), mcpjet.pt());
307+
if (collision.centrality() < centRange->at(0) || collision.centrality() > centRange->at(1))
308+
return;
309+
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);
310+
311+
if (!jetderiveddatautilities::selectCollision(collision, eventSelection))
312+
return;
313+
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);
314+
315+
if (cfgEvSelOccupancy && !isOccupancyWithin(collision))
316+
return;
317+
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);
318+
319+
registry.fill(HIST("hCentralityMult"), collision.centrality());
320+
auto collmcdJets{mcdjets.sliceBy(mcdjetsPerJCollision, collision.globalIndex())};
321+
for (const auto& mcdjet : collmcdJets) {
322+
auto mcdPt{mcdjet.pt() - (collision.rho() * mcdjet.area())};
323+
registry.fill(HIST("hDetectorJetPt"), mcdjet.pt());
324+
registry.fill(HIST("hDetectorJetPtSubBkg"), mcdPt);
325+
registry.fill(HIST("hDetectorJetEta"), mcdjet.eta());
326+
registry.fill(HIST("hDetectorJetPhi"), mcdjet.phi());
327+
for (const auto& mcpjet : mcdjet.template matchedJetGeo_as<JetMCPTable>()) {
328+
auto mcpPt{mcpjet.pt() - (mcCol.rho() * mcpjet.area())};
329+
registry.fill(HIST("hPartJetPtMatch"), mcpjet.pt());
330+
registry.fill(HIST("hPartJetPtMatchSubBkg"), mcpPt);
331+
registry.fill(HIST("hPartJetEtaMatch"), mcpjet.eta());
332+
registry.fill(HIST("hPartJetPhiMatch"), mcpjet.phi());
333+
registry.fill(HIST("hMatchedJetsPtDelta"), mcpPt, mcdPt - mcpPt);
334+
registry.fill(HIST("hGenMatchedJetsPtDeltadPt"), collision.centrality(), mcpPt, (mcdPt - mcpPt) / mcpPt);
335+
registry.fill(HIST("hRecoMatchedJetsPtDeltadPt"), collision.centrality(), mcdPt, (mcdPt - mcpPt) / mcdPt);
336+
registry.fill(HIST("hMatchedJetsPhiDelta"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi());
337+
registry.fill(HIST("hMatchedJetsEtaDelta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta());
338+
339+
registry.fill(HIST("hRespMcDMcPMatch"), collision.centrality(), mcdjet.pt(), mcpjet.pt());
340+
registry.fill(HIST("hRespMcDMcPMatchSubBkg"), collision.centrality(), mcdPt, mcpPt);
341+
}
283342
}
343+
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);
284344
}
285345
registry.fill(HIST("hMCEventCounter"), counter++);
286346
}
287-
PROCESS_SWITCH(JetSpectraEseTask, processMCChargedMatched, "jet matched mcp and mcd", false);
347+
PROCESS_SWITCH(JetSpectraEseTask, processMCChargedMatched, "jet geto matched mcp and mcd", false);
288348

289349
template <typename T>
290350
bool isAcceptedLeadTrack(T const& tracks)
@@ -322,7 +382,7 @@ struct JetSpectraEseTask {
322382
epCorrContainer[2] = cosPsi(epMap[cfgEPRefB], epMap[cfgEPRefC]);
323383
fillEPCos(/*std::make_index_sequence<3>{},*/ vec, epCorrContainer);
324384
}
325-
return epMap["FT0A"];
385+
return epMap[cfgEPRefA];
326386
}
327387
template </*std::size_t... Idx,*/ typename collision>
328388
void fillEPCos(/*const std::index_sequence<Idx...>&,*/ const collision& col, const std::array<float, 3>& Corr)
@@ -347,7 +407,7 @@ struct JetSpectraEseTask {
347407
}
348408
constexpr std::string_view RemovePrefix(std::string_view str)
349409
{
350-
constexpr std::string_view Prefix = "hPsi2";
410+
constexpr std::string_view Prefix{"hPsi2"};
351411
return str.substr(Prefix.size());
352412
}
353413

@@ -375,9 +435,9 @@ struct JetSpectraEseTask {
375435
template <DetID id, bool fill, typename Col>
376436
std::vector<float> qVecNoESE(Col collision)
377437
{
378-
const int nmode = 2;
379-
int detId = detIDN(id);
380-
int detInd = detId * 4 + cfgnTotalSystem * 4 * (nmode - 2);
438+
const int nmode{2};
439+
int detId{detIDN(id)};
440+
int detInd{detId * 4 + cfgnTotalSystem * 4 * (nmode - 2)};
381441
if constexpr (fill) {
382442
if (collision.qvecAmp()[detInd] > 1e-8) {
383443
registry.fill(HIST("hQvecUncorV2"), collision.centrality(), collision.qvecRe()[detInd], collision.qvecIm()[detInd]);
@@ -395,6 +455,16 @@ struct JetSpectraEseTask {
395455
qVec.push_back(collision.qvecAmp()[detId]);
396456
return qVec;
397457
}
458+
459+
template <typename col>
460+
bool isOccupancyWithin(const col& collision)
461+
{
462+
auto occupancy{collision.trackOccupancyInTimeRange()};
463+
if (occupancy < cfgCutOccupancy->at(0) || occupancy > cfgCutOccupancy->at(1))
464+
return false;
465+
else
466+
return true;
467+
}
398468
};
399469

400470
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetSpectraEseTask>(cfgc)}; }

0 commit comments

Comments
 (0)