Skip to content

Commit 932dbb9

Browse files
jokonigjokonig
authored
[EMCAL-737] Adding EMCal pi0 qc task (AliceO2Group#826)
- Task to obtain invariant mass vs. pt histograms for basic pi0 analysis - Goal to validate mass position and width of the pi0 using O2 and compare to AliPhysics - Basic classes for photons and mesons added - Task is NOT to be used for physics analyses and only serves the purpose of QC! - Task based on clustermonitor task in PWGJE
1 parent f1b5cce commit 932dbb9

2 files changed

Lines changed: 390 additions & 0 deletions

File tree

‎PWGEM/PhotonMeson/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,3 +18,8 @@ o2physics_add_dpl_workflow(gammaconversionstruthonlymc
1818
SOURCES gammaConversionsTruthOnlyMc.cxx
1919
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::DetectorsVertexing
2020
COMPONENT_NAME Analysis)
21+
22+
o2physics_add_dpl_workflow(emc-pi0-qc
23+
SOURCES emcalPi0QC.cxx
24+
PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
25+
COMPONENT_NAME Analysis)

‎PWGEM/PhotonMeson/Tasks/emcalPi0QC.cxx‎

Lines changed: 385 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,385 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
#include <climits>
13+
#include <cstdlib>
14+
#include <map>
15+
#include <memory>
16+
#include <sstream>
17+
#include <string>
18+
#include <vector>
19+
#include <math.h>
20+
21+
#include "Framework/runDataProcessing.h"
22+
#include "Framework/AnalysisTask.h"
23+
#include "Framework/AnalysisDataModel.h"
24+
#include "Framework/ASoA.h"
25+
#include "Framework/HistogramRegistry.h"
26+
27+
#include "Common/DataModel/EventSelection.h"
28+
#include "Common/DataModel/Centrality.h"
29+
30+
#include "EMCALBase/Geometry.h"
31+
#include "PWGJE/DataModel/EMCALClusters.h"
32+
#include "DataFormatsEMCAL/Cell.h"
33+
#include "DataFormatsEMCAL/Constants.h"
34+
#include "DataFormatsEMCAL/AnalysisCluster.h"
35+
36+
#include "CommonDataFormat/InteractionRecord.h"
37+
38+
#include "TLorentzVector.h"
39+
#include "TVector3.h"
40+
41+
// \struct Pi0QCTask
42+
/// \brief Simple monitoring task for EMCal clusters
43+
/// \author Joshua Koenig <joshua.konig@cern.ch>, Goethe University Frankfurt
44+
/// \since 25.05.2022
45+
///
46+
/// This task is meant to be used for QC for the emcal using properties of the pi0
47+
/// - ...
48+
/// Simple event selection using the flag doEventSel is provided, which selects INT7 events if set to 1
49+
/// For pilot beam data, instead of relying on the event selection, one can veto specific BC IDS using the flag
50+
/// fDoVetoBCID.
51+
52+
using namespace o2::framework;
53+
using namespace o2::framework::expressions;
54+
using collisionEvSelIt = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>::iterator;
55+
using selectedClusters = o2::soa::Filtered<o2::aod::EMCALClusters>;
56+
using selectedCluster = o2::soa::Filtered<o2::aod::EMCALCluster>;
57+
using selectedAmbiguousClusters = o2::soa::Filtered<o2::aod::EMCALAmbiguousClusters>;
58+
using selectedAmbiguousCluster = o2::soa::Filtered<o2::aod::EMCALAmbiguousCluster>;
59+
60+
struct Photon {
61+
Photon(float eta_tmp, float phi_tmp, float energy_tmp, int clusid = 0)
62+
{
63+
eta = eta_tmp;
64+
phi = phi_tmp;
65+
energy = energy_tmp;
66+
theta = 2 * std::atan2(std::exp(-eta), 1);
67+
px = energy * std::sin(theta) * std::cos(phi);
68+
py = energy * std::sin(theta) * std::sin(phi);
69+
pz = energy * std::cos(theta);
70+
pt = std::sqrt(px * px + py * py);
71+
photon.SetPxPyPzE(px, py, pz, energy);
72+
id = clusid;
73+
}
74+
75+
TLorentzVector photon;
76+
float pt;
77+
float px;
78+
float py;
79+
float pz;
80+
float eta;
81+
float phi;
82+
float energy;
83+
float theta;
84+
int id;
85+
};
86+
87+
struct Meson {
88+
Meson(Photon p1, Photon p2) : pgamma1(p1),
89+
pgamma2(p2)
90+
{
91+
pMeson = p1.photon + p2.photon;
92+
}
93+
Photon pgamma1;
94+
Photon pgamma2;
95+
TLorentzVector pMeson;
96+
97+
float getMass() const { return pMeson.M(); };
98+
float getPt() const { return pMeson.Pt(); };
99+
float getOpeningAngle() const { return pgamma1.photon.Angle(pgamma2.photon.Vect()); };
100+
};
101+
102+
struct Pi0QCTask {
103+
HistogramRegistry mHistManager{"NeutralMesonHistograms"};
104+
o2::emcal::Geometry* mGeometry = nullptr;
105+
106+
// configurable parameters
107+
// TODO adapt mDoEventSel switch to also allow selection of other triggers (e.g. EMC7)
108+
Configurable<bool> mDoEventSel{"doEventSel", 0, "demand kINT7"};
109+
Configurable<std::string> mVetoBCID{"vetoBCID", "", "BC ID(s) to be excluded, this should be used as an alternative to the event selection"};
110+
Configurable<std::string> mSelectBCID{"selectBCID", "all", "BC ID(s) to be included, this should be used as an alternative to the event selection"};
111+
Configurable<double> mVertexCut{"vertexCut", -1, "apply z-vertex cut with value in cm"};
112+
Configurable<int> mTimeMin{"TimeMinCut", -600, "apply min timing cut (in ns)"};
113+
Configurable<int> mTimeMax{"TimeMaxCut", 900, "apply min timing cut (in ns)"};
114+
Configurable<float> mClusterMinM02Cut{"MinM02Cut", 0.1, "apply min M02 cut"};
115+
Configurable<float> mClusterMaxM02Cut{"MaxM02Cut", 0.7, "apply max M02 cut"};
116+
Configurable<float> mMinEnergyCut{"MinEnergyCut", 0.7, "apply min cluster energy cut"};
117+
Configurable<int> mMinNCellsCut{"MinNCellsCut", 1, "apply min cluster number of cell cut"};
118+
Configurable<std::string> mClusterDefinition{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"};
119+
std::vector<int> mVetoBCIDs;
120+
std::vector<int> mSelectBCIDs;
121+
122+
// define cluster filter. It selects only those clusters which are of the type
123+
// specified in the string mClusterDefinition,e.g. kV3Default, which is V3 clusterizer with default
124+
// clusterization parameters
125+
o2::aod::EMCALClusterDefinition clusDef = o2::aod::emcalcluster::getClusterDefinitionFromString(mClusterDefinition.value);
126+
Filter clusterDefinitionSelection = o2::aod::emcalcluster::definition == static_cast<int>(clusDef);
127+
128+
// define container for photons
129+
std::vector<Photon> mPhotons;
130+
131+
/// \brief Create output histograms and initialize geometry
132+
void init(InitContext const&)
133+
{
134+
// create histograms
135+
using o2HistType = HistType;
136+
using o2Axis = AxisSpec;
137+
138+
// load geometry just in case we need it
139+
mGeometry = o2::emcal::Geometry::GetInstanceFromRunNumber(300000);
140+
141+
// create common axes
142+
LOG(info) << "Creating histograms";
143+
const o2Axis bcAxis{3501, -0.5, 3500.5};
144+
const o2Axis energyAxis{makeClusterBinning(), "#it{p}_{T} (GeV)"};
145+
146+
// event properties
147+
mHistManager.add("eventsAll", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}});
148+
mHistManager.add("eventsSelected", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}});
149+
mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", o2HistType::kTH1F, {bcAxis});
150+
mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", o2HistType::kTH1F, {bcAxis});
151+
mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", o2HistType::kTH1F, {{200, -20, 20}});
152+
mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", o2HistType::kTH1F, {{200, -20, 20}});
153+
154+
// cluster properties
155+
mHistManager.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis});
156+
mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{400, 0, 100}});
157+
mHistManager.add("clusterEtaPhi", "Eta and phi of cluster", o2HistType::kTH2F, {{100, -1, 1}, {100, 0, 2 * TMath::Pi()}});
158+
mHistManager.add("clusterM02", "M02 of cluster", o2HistType::kTH1F, {{400, 0, 5}});
159+
mHistManager.add("clusterM20", "M20 of cluster", o2HistType::kTH1F, {{400, 0, 2.5}});
160+
mHistManager.add("clusterNLM", "Number of local maxima of cluster", o2HistType::kTH1I, {{10, 0, 10}});
161+
mHistManager.add("clusterNCells", "Number of cells in cluster", o2HistType::kTH1I, {{50, 0, 50}});
162+
mHistManager.add("clusterDistanceToBadChannel", "Distance to bad channel", o2HistType::kTH1F, {{100, 0, 100}});
163+
164+
// meson related histograms
165+
mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}});
166+
mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}});
167+
168+
if (mVetoBCID->length()) {
169+
std::stringstream parser(mVetoBCID.value);
170+
std::string token;
171+
int bcid;
172+
while (std::getline(parser, token, ',')) {
173+
bcid = std::stoi(token);
174+
LOG(info) << "Veto BCID " << bcid;
175+
mVetoBCIDs.push_back(bcid);
176+
}
177+
}
178+
if (mSelectBCID.value != "all") {
179+
std::stringstream parser(mSelectBCID.value);
180+
std::string token;
181+
int bcid;
182+
while (std::getline(parser, token, ',')) {
183+
bcid = std::stoi(token);
184+
LOG(info) << "Select BCID " << bcid;
185+
mSelectBCIDs.push_back(bcid);
186+
}
187+
}
188+
}
189+
/// \brief Process EMCAL clusters that are matched to a collisions
190+
void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::BCs const& bcs)
191+
{
192+
mHistManager.fill(HIST("eventsAll"), 1);
193+
194+
// do event selection if mDoEventSel is specified
195+
// currently the event selection is hard coded to kINT7
196+
// but other selections are possible that are defined in TriggerAliases.h
197+
if (mDoEventSel && (!theCollision.alias()[kINT7])) {
198+
LOG(debug) << "Event not selected becaus it is not kINT7, skipping";
199+
return;
200+
}
201+
mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ());
202+
if (mVertexCut > 0 && TMath::Abs(theCollision.posZ()) > mVertexCut) {
203+
LOG(debug) << "Event not selected because of z-vertex cut z= " << theCollision.posZ() << " > " << mVertexCut << " cm, skipping";
204+
return;
205+
}
206+
mHistManager.fill(HIST("eventsSelected"), 1);
207+
mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ());
208+
209+
ProcessClusters(theCollision, clusters, bcs);
210+
ProcessMesons(theCollision, clusters, bcs);
211+
}
212+
PROCESS_SWITCH(Pi0QCTask, processCollisions, "Process clusters from collision", false);
213+
214+
/// \brief Process EMCAL clusters that are not matched to a collision
215+
/// This is not needed for most users
216+
void processAmbiguous(o2::aod::BC const bc, selectedAmbiguousClusters const& clusters)
217+
{
218+
// loop over bc , if requested (mVetoBCID >= 0), reject everything from a certain BC
219+
// this can be used as alternative to event selection (e.g. for pilot beam data)
220+
// TODO: remove this loop and put it in separate process function that only takes care of ambiguous clusters
221+
o2::InteractionRecord eventIR;
222+
eventIR.setFromLong(bc.globalBC());
223+
mHistManager.fill(HIST("eventBCAll"), eventIR.bc);
224+
if (std::find(mVetoBCIDs.begin(), mVetoBCIDs.end(), eventIR.bc) != mVetoBCIDs.end()) {
225+
LOG(info) << "Event rejected because of veto BCID " << eventIR.bc;
226+
return;
227+
}
228+
if (mSelectBCIDs.size() && (std::find(mSelectBCIDs.begin(), mSelectBCIDs.end(), eventIR.bc) == mSelectBCIDs.end())) {
229+
return;
230+
}
231+
mHistManager.fill(HIST("eventBCSelected"), eventIR.bc);
232+
233+
// ToDo: Add mode if collision is not found
234+
// ProcessClusters(theCollision, clusters, bcs);
235+
// ProcessMesons(theCollision, clusters, bcs);
236+
}
237+
PROCESS_SWITCH(Pi0QCTask, processAmbiguous, "Process Ambiguous clusters", false);
238+
239+
/// \brief Process EMCAL clusters that are matched to a collisions
240+
template <typename Clusters>
241+
void ProcessClusters(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs)
242+
{
243+
// clear photon vector
244+
mPhotons.clear();
245+
246+
// loop over all clusters from accepted collision
247+
// auto eventClusters = clusters.select(o2::aod::emcalcluster::bcId == theCollision.bc().globalBC());
248+
for (const auto& cluster : clusters) {
249+
// fill histograms of cluster properties
250+
// in this implementation the cluster properties are directly
251+
// loaded from the flat table, in the future one should
252+
// consider using the AnalysisCluster object to work with
253+
// after loading.
254+
LOG(debug) << "Cluster energy: " << cluster.energy();
255+
LOG(debug) << "Cluster time: " << cluster.time();
256+
LOG(debug) << "Cluster M02: " << cluster.m02();
257+
mHistManager.fill(HIST("clusterE"), cluster.energy());
258+
mHistManager.fill(HIST("clusterE_SimpleBinning"), cluster.energy());
259+
mHistManager.fill(HIST("clusterEtaPhi"), cluster.eta(), cluster.phi());
260+
mHistManager.fill(HIST("clusterM02"), cluster.m02());
261+
mHistManager.fill(HIST("clusterM20"), cluster.m20());
262+
mHistManager.fill(HIST("clusterNLM"), cluster.nlm());
263+
mHistManager.fill(HIST("clusterNCells"), cluster.nCells());
264+
mHistManager.fill(HIST("clusterDistanceToBadChannel"), cluster.distanceToBadChannel());
265+
266+
// apply basic cluster cuts
267+
if (cluster.energy() < mMinEnergyCut) {
268+
LOG(debug) << "Cluster rejected because of energy cut";
269+
continue;
270+
}
271+
if (cluster.nCells() <= mMinNCellsCut) {
272+
LOG(debug) << "Cluster rejected because of nCells cut";
273+
continue;
274+
}
275+
if (cluster.m02() < mClusterMinM02Cut || cluster.m02() > mClusterMaxM02Cut) {
276+
LOG(debug) << "Cluster rejected because of m02 cut";
277+
continue;
278+
}
279+
if (cluster.time() < mTimeMin || cluster.time() > mTimeMax) {
280+
LOG(debug) << "Cluster rejected because of time cut";
281+
continue;
282+
}
283+
284+
// put clusters in photon vector
285+
// ToDo: At the moment, the eta and phi values are not corrected for a shift of the primary vertex! Should only be a small effect but has to be corrected
286+
mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id()));
287+
}
288+
}
289+
290+
/// \brief Process meson candidates, caluclate invariant mass and pT and fill histograms
291+
template <typename Clusters>
292+
void ProcessMesons(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs)
293+
{
294+
// if less then 2 clusters are found, skip event
295+
if (mPhotons.size() < 2) {
296+
return;
297+
}
298+
299+
// loop over all photon combinations and build meson candidates
300+
for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) {
301+
for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) {
302+
303+
// build meson from photons
304+
Meson meson(mPhotons[ig1], mPhotons[ig2]);
305+
mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt());
306+
307+
// calculate background candidates (rotation background)
308+
CalculateBackground(meson, ig1, ig2);
309+
}
310+
}
311+
}
312+
313+
/// \brief Calculate background (using rotation background method)
314+
void CalculateBackground(const Meson& meson, unsigned int ig1, unsigned int ig2)
315+
{
316+
// if less than 3 clusters are present, skip event
317+
if (mPhotons.size() < 3) {
318+
return;
319+
}
320+
const double rotationAngle = M_PI / 2.0; //0.78539816339; // rotaion angle 90°
321+
322+
TLorentzVector lvRotationPhoton1; // photon candidates which get rotated
323+
TLorentzVector lvRotationPhoton2; // photon candidates which get rotated
324+
TVector3 lvRotationPion; // rotation axis
325+
for (unsigned int ig3 = 0; ig3 < mPhotons.size(); ++ig3) {
326+
// continue if photons are identical
327+
if (ig3 == ig1 || ig3 == ig2) {
328+
continue;
329+
}
330+
// calculate rotation axis
331+
lvRotationPion = (meson.pMeson).Vect();
332+
333+
// initialize photons for rotation
334+
lvRotationPhoton1.SetPxPyPzE(mPhotons[ig1].px, mPhotons[ig1].py, mPhotons[ig1].pz, mPhotons[ig1].energy);
335+
lvRotationPhoton2.SetPxPyPzE(mPhotons[ig2].px, mPhotons[ig2].py, mPhotons[ig2].pz, mPhotons[ig2].energy);
336+
337+
// rotate photons around rotation axis
338+
lvRotationPhoton1.Rotate(rotationAngle, lvRotationPion);
339+
lvRotationPhoton2.Rotate(rotationAngle, lvRotationPion);
340+
341+
// initialize Photon objects for rotated photons
342+
Photon rotPhoton1(lvRotationPhoton1.Eta(), lvRotationPhoton1.Phi(), lvRotationPhoton1.E(), mPhotons[ig1].id);
343+
Photon rotPhoton2(lvRotationPhoton2.Eta(), lvRotationPhoton2.Phi(), lvRotationPhoton2.E(), mPhotons[ig2].id);
344+
345+
// build meson from rotated photons
346+
Meson mesonRotated1(rotPhoton1, mPhotons[ig3]);
347+
Meson mesonRotated2(rotPhoton2, mPhotons[ig3]);
348+
349+
// Fill histograms
350+
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt());
351+
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt());
352+
}
353+
}
354+
355+
/// \brief Create binning for cluster energy/pT axis (variable bin size)
356+
/// direct port from binning often used in AliPhysics for debugging
357+
/// \return vector with bin limits
358+
std::vector<double> makeClusterBinning() const
359+
{
360+
361+
std::vector<double> result;
362+
int nBinsPt = 179;
363+
double maxPt = 60;
364+
for (Int_t i = 0; i < nBinsPt + 1; i++) {
365+
if (i < 100) {
366+
result.emplace_back(0.10 * i);
367+
} else if (i < 140) {
368+
result.emplace_back(10. + 0.25 * (i - 100));
369+
} else if (i < 180) {
370+
result.emplace_back(20. + 1.00 * (i - 140));
371+
} else {
372+
result.emplace_back(maxPt);
373+
}
374+
}
375+
return result;
376+
}
377+
};
378+
379+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
380+
{
381+
WorkflowSpec workflow{
382+
adaptAnalysisTask<Pi0QCTask>(cfgc, TaskName{"EMCPi0QCTask"}, SetDefaultProcesses{{{"processCollisions", true}, {"processAmbiguous", false}}}),
383+
adaptAnalysisTask<Pi0QCTask>(cfgc, TaskName{"EMCPi0QCTaskAmbiguous"}, SetDefaultProcesses{{{"processCollisions", false}, {"processAmbiguous", true}}})};
384+
return workflow;
385+
}

0 commit comments

Comments
 (0)