Skip to content

Commit 2a7633d

Browse files
jokonigjokonig
authored
PWGEM/PhotonMeson: Add Background calculation to pi0QC task (#1597)
- Add mixed events and rotation background calculation - Fix loop over collisions instead of whole timeframes - Added opening angle cut
1 parent f158ee1 commit 2a7633d

1 file changed

Lines changed: 152 additions & 28 deletions

File tree

‎PWGEM/PhotonMeson/Tasks/emcalPi0QC.cxx‎

Lines changed: 152 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
#include <sstream>
1717
#include <string>
1818
#include <vector>
19-
#include <math.h>
19+
#include <cmath>
2020

2121
#include "Framework/runDataProcessing.h"
2222
#include "Framework/AnalysisTask.h"
@@ -51,7 +51,7 @@
5151

5252
using namespace o2::framework;
5353
using namespace o2::framework::expressions;
54-
using collisionEvSelIt = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>::iterator;
54+
using collisionEvSelIt = o2::aod::Collision;
5555
using selectedClusters = o2::soa::Filtered<o2::aod::EMCALClusters>;
5656
using selectedCluster = o2::soa::Filtered<o2::aod::EMCALCluster>;
5757
using selectedAmbiguousClusters = o2::soa::Filtered<o2::aod::EMCALAmbiguousClusters>;
@@ -94,9 +94,33 @@ struct Meson {
9494
Photon pgamma2;
9595
TLorentzVector pMeson;
9696

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()); };
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 EventMixVec {
103+
104+
void AddEvent(std::vector<Photon> vecGamma)
105+
{
106+
if (vecEvtMix.size() < nEVtMixSize) {
107+
vecEvtMix.push_back(vecGamma);
108+
} else {
109+
vecEvtMix.erase(vecEvtMix.begin() + nEVtMixSize - 1);
110+
vecEvtMix.push_back(vecGamma);
111+
}
112+
}
113+
Photon* getPhoton(unsigned int iEvt, unsigned int iGamma)
114+
{
115+
if (vecEvtMix.size() >= iEvt)
116+
return nullptr;
117+
if (vecEvtMix[iEvt].size() >= iGamma)
118+
return nullptr;
119+
return &vecEvtMix[iEvt][iGamma];
120+
}
121+
122+
std::vector<std::vector<Photon>> vecEvtMix;
123+
unsigned int nEVtMixSize = 20;
100124
};
101125

102126
struct Pi0QCTask {
@@ -115,6 +139,7 @@ struct Pi0QCTask {
115139
Configurable<float> mClusterMaxM02Cut{"MaxM02Cut", 0.7, "apply max M02 cut"};
116140
Configurable<float> mMinEnergyCut{"MinEnergyCut", 0.7, "apply min cluster energy cut"};
117141
Configurable<int> mMinNCellsCut{"MinNCellsCut", 1, "apply min cluster number of cell cut"};
142+
Configurable<float> mMinOpenAngleCut{"OpeningAngleCut", 0.0202, "apply min opening angle cut"};
118143
Configurable<std::string> mClusterDefinition{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"};
119144
std::vector<int> mVetoBCIDs;
120145
std::vector<int> mSelectBCIDs;
@@ -127,6 +152,11 @@ struct Pi0QCTask {
127152

128153
// define container for photons
129154
std::vector<Photon> mPhotons;
155+
// define container for photons for each collision
156+
std::map<int, std::vector<Photon>> mapPhotons;
157+
158+
// event mixing class
159+
EventMixVec evtMix;
130160

131161
/// \brief Create output histograms and initialize geometry
132162
void init(InitContext const&)
@@ -145,6 +175,7 @@ struct Pi0QCTask {
145175

146176
// event properties
147177
mHistManager.add("eventsAll", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}});
178+
mHistManager.add("eventsEMCTrigg", "Number of EMC triggered events", o2HistType::kTH1F, {{1, 0.5, 1.5}});
148179
mHistManager.add("eventsSelected", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}});
149180
mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", o2HistType::kTH1F, {bcAxis});
150181
mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", o2HistType::kTH1F, {bcAxis});
@@ -154,6 +185,7 @@ struct Pi0QCTask {
154185
// cluster properties
155186
mHistManager.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis});
156187
mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{400, 0, 100}});
188+
mHistManager.add("clusterTime", "Time of cluster", o2HistType::kTH1F, {{500, -250, 250}});
157189
mHistManager.add("clusterEtaPhi", "Eta and phi of cluster", o2HistType::kTH2F, {{100, -1, 1}, {100, 0, 2 * TMath::Pi()}});
158190
mHistManager.add("clusterM02", "M02 of cluster", o2HistType::kTH1F, {{400, 0, 5}});
159191
mHistManager.add("clusterM20", "M20 of cluster", o2HistType::kTH1F, {{400, 0, 2.5}});
@@ -164,6 +196,7 @@ struct Pi0QCTask {
164196
// meson related histograms
165197
mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}});
166198
mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}});
199+
mHistManager.add("invMassVsPtMixedBackground", "invariant mass and pT of mixed background meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}});
167200

168201
if (mVetoBCID->length()) {
169202
std::stringstream parser(mVetoBCID.value);
@@ -187,36 +220,38 @@ struct Pi0QCTask {
187220
}
188221
}
189222
/// \brief Process EMCAL clusters that are matched to a collisions
190-
void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::BCs const& bcs)
223+
224+
// void processCollisions(collisionEvSelIt const& collision, selectedClusters const& clusters)
225+
void processCollisions(o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>::iterator const& collision, selectedClusters const& clusters)
191226
{
227+
// for(const auto & collision : theCollisions){
192228
mHistManager.fill(HIST("eventsAll"), 1);
193-
229+
LOG(debug) << "processCollisions";
194230
// do event selection if mDoEventSel is specified
195231
// currently the event selection is hard coded to kINT7
196232
// but other selections are possible that are defined in TriggerAliases.h
197-
if (mDoEventSel && (!theCollision.alias()[kINT7])) {
233+
if (mDoEventSel && (!collision.alias()[kINT7])) {
198234
LOG(debug) << "Event not selected becaus it is not kINT7, skipping";
199235
return;
200236
}
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";
237+
mHistManager.fill(HIST("eventVertexZAll"), collision.posZ());
238+
if (mVertexCut > 0 && std::abs(collision.posZ()) > mVertexCut) {
239+
LOG(debug) << "Event not selected because of z-vertex cut z= " << collision.posZ() << " > " << mVertexCut << " cm, skipping";
204240
return;
205241
}
206242
mHistManager.fill(HIST("eventsSelected"), 1);
207-
mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ());
243+
mHistManager.fill(HIST("eventVertexZSelected"), collision.posZ());
208244

209-
ProcessClusters(theCollision, clusters, bcs);
210-
ProcessMesons(theCollision, clusters, bcs);
245+
ProcessClusters(clusters);
246+
ProcessMesons(clusters);
211247
}
212248
PROCESS_SWITCH(Pi0QCTask, processCollisions, "Process clusters from collision", false);
213249

214250
/// \brief Process EMCAL clusters that are not matched to a collision
215251
/// This is not needed for most users
216-
void processAmbiguous(o2::aod::BC const& bc, selectedAmbiguousClusters const& clusters)
252+
void processAmbiguous(o2::aod::BCs::iterator const& bc, selectedAmbiguousClusters const& clusters)
217253
{
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)
254+
LOG(debug) << "processAmbiguous";
220255
// TODO: remove this loop and put it in separate process function that only takes care of ambiguous clusters
221256
o2::InteractionRecord eventIR;
222257
eventIR.setFromLong(bc.globalBC());
@@ -230,22 +265,35 @@ struct Pi0QCTask {
230265
}
231266
mHistManager.fill(HIST("eventBCSelected"), eventIR.bc);
232267

233-
// ToDo: Add mode if collision is not found
234-
// ProcessClusters(theCollision, clusters, bcs);
235-
// ProcessMesons(theCollision, clusters, bcs);
268+
ProcessAmbigousClusters(clusters);
269+
ProcessMesons(clusters);
236270
}
237271
PROCESS_SWITCH(Pi0QCTask, processAmbiguous, "Process Ambiguous clusters", false);
238272

239273
/// \brief Process EMCAL clusters that are matched to a collisions
240274
template <typename Clusters>
241-
void ProcessClusters(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs)
275+
void ProcessClusters(Clusters const& clusters)
242276
{
277+
LOG(debug) << "ProcessClusters";
243278
// clear photon vector
244279
mPhotons.clear();
280+
mapPhotons.clear();
281+
282+
int globalCollID = -1000;
245283

246284
// loop over all clusters from accepted collision
247285
// auto eventClusters = clusters.select(o2::aod::emcalcluster::bcId == theCollision.bc().globalBC());
248286
for (const auto& cluster : clusters) {
287+
288+
// o2::InteractionRecord eventIR;
289+
auto collID = cluster.collisionId();
290+
if (globalCollID == -1000)
291+
globalCollID = collID;
292+
293+
if (globalCollID != collID) {
294+
LOG(info) << "Something went wrong with the collision ID";
295+
}
296+
249297
// fill histograms of cluster properties
250298
// in this implementation the cluster properties are directly
251299
// loaded from the flat table, in the future one should
@@ -255,6 +303,59 @@ struct Pi0QCTask {
255303
LOG(debug) << "Cluster time: " << cluster.time();
256304
LOG(debug) << "Cluster M02: " << cluster.m02();
257305
mHistManager.fill(HIST("clusterE"), cluster.energy());
306+
mHistManager.fill(HIST("clusterTime"), cluster.time());
307+
mHistManager.fill(HIST("clusterE_SimpleBinning"), cluster.energy());
308+
mHistManager.fill(HIST("clusterEtaPhi"), cluster.eta(), cluster.phi());
309+
mHistManager.fill(HIST("clusterM02"), cluster.m02());
310+
mHistManager.fill(HIST("clusterM20"), cluster.m20());
311+
mHistManager.fill(HIST("clusterNLM"), cluster.nlm());
312+
mHistManager.fill(HIST("clusterNCells"), cluster.nCells());
313+
mHistManager.fill(HIST("clusterDistanceToBadChannel"), cluster.distanceToBadChannel());
314+
315+
// apply basic cluster cuts
316+
if (cluster.energy() < mMinEnergyCut) {
317+
LOG(debug) << "Cluster rejected because of energy cut";
318+
continue;
319+
}
320+
if (cluster.nCells() <= mMinNCellsCut) {
321+
LOG(debug) << "Cluster rejected because of nCells cut";
322+
continue;
323+
}
324+
if (cluster.m02() < mClusterMinM02Cut || cluster.m02() > mClusterMaxM02Cut) {
325+
LOG(debug) << "Cluster rejected because of m02 cut";
326+
continue;
327+
}
328+
if (cluster.time() < mTimeMin || cluster.time() > mTimeMax) {
329+
LOG(debug) << "Cluster rejected because of time cut";
330+
continue;
331+
}
332+
333+
// put clusters in photon vector
334+
mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id()));
335+
}
336+
}
337+
338+
/// \brief Process EMCAL clusters that are matched to a collisions
339+
template <typename Clusters>
340+
void ProcessAmbigousClusters(Clusters const& clusters)
341+
{
342+
LOG(debug) << "ProcessClusters";
343+
// clear photon vector
344+
mPhotons.clear();
345+
346+
// loop over all clusters from accepted collision
347+
for (const auto& cluster : clusters) {
348+
349+
// fill histograms of cluster properties
350+
// in this implementation the cluster properties are directly
351+
// loaded from the flat table, in the future one should
352+
// consider using the AnalysisCluster object to work with
353+
// after loading.
354+
LOG(debug) << "Cluster energy: " << cluster.energy();
355+
LOG(debug) << "Cluster time: " << cluster.time();
356+
LOG(debug) << "Cluster M02: " << cluster.m02();
357+
mHistManager.fill(HIST("clusterE"), cluster.energy());
358+
mHistManager.fill(HIST("clusterTime"), cluster.time());
258359
mHistManager.fill(HIST("clusterE_SimpleBinning"), cluster.energy());
259360
mHistManager.fill(HIST("clusterEtaPhi"), cluster.eta(), cluster.phi());
260361
mHistManager.fill(HIST("clusterM02"), cluster.m02());
@@ -282,32 +383,39 @@ struct Pi0QCTask {
282383
}
283384

284385
// 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
286386
mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id()));
287387
}
288388
}
289389

290390
/// \brief Process meson candidates, calculate invariant mass and pT and fill histograms
291391
template <typename Clusters>
292-
void ProcessMesons(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs)
392+
void ProcessMesons(Clusters const& clusters)
293393
{
394+
LOG(debug) << "ProcessMesons " << mPhotons.size();
395+
396+
mHistManager.fill(HIST("eventsEMCTrigg"), 1);
397+
294398
// if less then 2 clusters are found, skip event
295-
if (mPhotons.size() < 2) {
399+
if (mPhotons.size() < 2)
296400
return;
297-
}
298401

299402
// loop over all photon combinations and build meson candidates
300403
for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) {
301404
for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) {
302405

303406
// build meson from photons
304407
Meson meson(mPhotons[ig1], mPhotons[ig2]);
305-
mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt());
408+
if (meson.getOpeningAngle() > mMinOpenAngleCut) {
409+
mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt());
410+
}
306411

307412
// calculate background candidates (rotation background)
308413
CalculateBackground(meson, ig1, ig2);
309414
}
415+
CalculateMixedBack(mPhotons[ig1]);
310416
}
417+
418+
evtMix.AddEvent(mPhotons);
311419
}
312420

313421
/// \brief Calculate background (using rotation background method)
@@ -347,8 +455,24 @@ struct Pi0QCTask {
347455
Meson mesonRotated2(rotPhoton2, mPhotons[ig3]);
348456

349457
// Fill histograms
350-
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt());
351-
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt());
458+
if (mesonRotated1.getOpeningAngle() > mMinOpenAngleCut) {
459+
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt());
460+
}
461+
if (mesonRotated2.getOpeningAngle() > mMinOpenAngleCut) {
462+
mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt());
463+
}
464+
}
465+
}
466+
467+
void CalculateMixedBack(Photon gamma)
468+
{
469+
for (unsigned int i = 0; i < evtMix.vecEvtMix.size(); ++i) {
470+
for (unsigned int ig1 = 0; ig1 < evtMix.vecEvtMix[i].size(); ++ig1) {
471+
Meson meson(gamma, evtMix.vecEvtMix[i][ig1]);
472+
if (meson.getOpeningAngle() > mMinOpenAngleCut) {
473+
mHistManager.fill(HIST("invMassVsPtMixedBackground"), meson.getMass(), meson.getPt());
474+
}
475+
}
352476
}
353477
}
354478

0 commit comments

Comments
 (0)