Skip to content

Commit 7847e60

Browse files
mhemmer-cernmhemmer-cern
authored
[Common] Run2 PCM QC Analysis (#11957)
1 parent 71e1a65 commit 7847e60

File tree

2 files changed

+319
-0
lines changed

2 files changed

+319
-0
lines changed

‎Common/LegacyDataQA/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,3 +28,8 @@ o2physics_add_dpl_workflow(pmd-qa
2828
SOURCES pmdQa.cxx
2929
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
3030
COMPONENT_NAME Analysis)
31+
32+
o2physics_add_dpl_workflow(pcm-run2
33+
SOURCES pcmRun2.cxx
34+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
35+
COMPONENT_NAME Analysis)

‎Common/LegacyDataQA/pcmRun2.cxx‎

Lines changed: 314 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,314 @@
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+
/// \file pcmRun2.cxx
13+
/// \brief Analysis using PCM photons from Run2
14+
/// \author M. Hemmer, marvin.hemmer@cern.ch
15+
16+
#include "Common/DataModel/PIDResponseTPC.h"
17+
18+
#include <Framework/ASoAHelpers.h>
19+
#include <Framework/AnalysisDataModel.h>
20+
#include <Framework/AnalysisTask.h>
21+
#include <Framework/Configurable.h>
22+
#include <Framework/HistogramRegistry.h>
23+
#include <Framework/HistogramSpec.h>
24+
#include <Framework/InitContext.h>
25+
#include <Framework/OutputObjHeader.h>
26+
#include <Framework/runDataProcessing.h>
27+
28+
#include <Math/Vector4D.h>
29+
30+
#include <array>
31+
#include <cmath>
32+
#include <string>
33+
#include <tuple>
34+
#include <utility>
35+
#include <vector>
36+
37+
using namespace o2;
38+
using namespace o2::aod;
39+
using namespace o2::framework;
40+
using namespace o2::framework::expressions;
41+
using namespace o2::soa;
42+
struct PcmRun2 {
43+
struct : ConfigurableGroup {
44+
std::string prefix = "trackcuts";
45+
Configurable<float> tpcFindOverFoundable{"tpcFindOverFoundable", 0.6, "Ratio of number of found TPC clusters to the number of foundable TPC clusters a track must have at least."};
46+
Configurable<float> minPt{"minPt", 0.05, "Minimum pT cut for the tracks."};
47+
} trackcuts;
48+
49+
struct : ConfigurableGroup {
50+
std::string prefix = "pidcuts";
51+
Configurable<float> minNSigmaEl{"minNSigmaEl", -3., "Minimum NSgimal electron allowed for V0 leg."};
52+
Configurable<float> maxNSigmaEl{"maxNSigmaEl", +4., "Maximum NSgimal electron allowed for V0 leg."};
53+
Configurable<bool> doPionRejection{"doPionRejection", true, "Flag to enable pion rejection based on TPC PID for V0 legs."};
54+
Configurable<float> minNSigmaPiLowP{"minNSigmaPiLowP", 1., "Minimum NSgimal pion to reject V0 legs for low 0.4 < pT/(GeV/c) < 3.5."};
55+
Configurable<float> minNSigmaPiHighP{"minNSigmaPiHighP", +0.5, "Minimum NSgimal pion to reject V0 legs for pT > 3.5 GeV/c."};
56+
} pidcuts;
57+
58+
struct : ConfigurableGroup {
59+
std::string prefix = "v0cuts";
60+
Configurable<float> minPt{"minPt", 0.02, "Minimum pT cut for the V0."};
61+
Configurable<float> maxEta{"maxEta", 0.8, "Maximum absolut pseudorapidity cut for the V0."};
62+
Configurable<float> maxZconv{"maxZconv", 0.8, "Maximum absolut z conversion position cut for the V0."};
63+
Configurable<float> qTFactor{"qTFactor", 0.125, "qT < this * pT"};
64+
Configurable<float> cosP{"cosP", 0.85, "cos of pointing angle > this value."};
65+
Configurable<bool> rejectTooCloseV0{"rejectTooCloseV0", true, "."};
66+
} v0cuts;
67+
68+
using TracksWithPID = soa::Join<o2::aod::FullTracks, o2::aod::pidTPCPi, o2::aod::pidTPCEl>;
69+
70+
SliceCache cache;
71+
Preslice<o2::aod::Run2OTFV0s> perCollisionV0 = aod::run2::oftv0::collisionId;
72+
73+
HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
74+
75+
void init(InitContext&)
76+
{
77+
78+
const AxisSpec thnAxisInvMass{400, 0.0, 0.8, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"};
79+
const AxisSpec thnAxisPt{100, 0., 20., "#it{p}_{T} (GeV/#it{c})"};
80+
const AxisSpec thnAxisCent{20, 0., 100., "Centrality (%)"};
81+
const AxisSpec thnAxisCuts{5, -0.5, 4.5, "cuts"};
82+
const AxisSpec thnAxisNSgimaE{254, -6.35f, 6.35f, "n#sigma_{e}"};
83+
84+
registry.add("hSameEvent2D", "hSameEvent2D", HistType::kTH2D, {thnAxisInvMass, thnAxisPt});
85+
registry.add("hNSigmaEPt", "hNSigmaEPt", HistType::kTH2D, {thnAxisNSgimaE, thnAxisPt});
86+
87+
registry.add("hCuts", "hCuts", HistType::kTH1D, {thnAxisCuts});
88+
registry.add("hCutsPair", "hCutsPair", HistType::kTH1D, {thnAxisCuts});
89+
90+
}; // end init
91+
92+
template <typename CollisionIter, typename V0Iter>
93+
float getCosP(CollisionIter const& coll, V0Iter const& v0)
94+
{
95+
// V0 momentum
96+
float px = v0.px();
97+
float py = v0.py();
98+
float pz = v0.pz();
99+
100+
// V0 decay position
101+
float xV0 = v0.x();
102+
float yV0 = v0.y();
103+
float zV0 = v0.z();
104+
105+
// Primary vertex
106+
float xPV = coll.posX();
107+
float yPV = coll.posY();
108+
float zPV = coll.posZ();
109+
110+
// Displacement vector r = V0 - PV
111+
float rx = xV0 - xPV;
112+
float ry = yV0 - yPV;
113+
float rz = zV0 - zPV;
114+
115+
// Dot product p·r
116+
float dot = px * rx + py * ry + pz * rz;
117+
118+
// Magnitudes |p| and |r|
119+
float magP = std::sqrt(px * px + py * py + pz * pz);
120+
float magR = std::sqrt(rx * rx + ry * ry + rz * rz);
121+
122+
// Cosine of pointing angle
123+
return dot / (magP * magR);
124+
}
125+
126+
template <typename V0Iter>
127+
bool checkLegs(V0Iter const& v0)
128+
{
129+
130+
auto posLeg = v0.template posTrack_as<TracksWithPID>();
131+
auto negLeg = v0.template negTrack_as<TracksWithPID>();
132+
133+
registry.fill(HIST("hNSigmaEPt"), posLeg.tpcNSigmaEl(), posLeg.pt());
134+
registry.fill(HIST("hNSigmaEPt"), negLeg.tpcNSigmaEl(), negLeg.pt());
135+
136+
// first let's check the positive leg
137+
if (posLeg.pt() <= trackcuts.minPt.value) {
138+
registry.fill(HIST("hCuts"), 1);
139+
return false;
140+
}
141+
if (posLeg.tpcFoundOverFindableCls() <= trackcuts.tpcFindOverFoundable.value) {
142+
registry.fill(HIST("hCuts"), 1);
143+
return false;
144+
}
145+
if (posLeg.tpcNSigmaEl() <= pidcuts.minNSigmaEl.value || posLeg.tpcNSigmaEl() >= pidcuts.maxNSigmaEl.value) {
146+
registry.fill(HIST("hCuts"), 2);
147+
return false;
148+
}
149+
150+
float minP = 0.4f; // minimum momentum for legs
151+
float midP = 3.5f; // momentum where we change NSigma window
152+
if (pidcuts.doPionRejection.value && posLeg.p() > minP) {
153+
if (posLeg.p() < midP && posLeg.tpcNSigmaPi() <= pidcuts.minNSigmaPiLowP.value) {
154+
registry.fill(HIST("hCuts"), 2);
155+
return false;
156+
} else if (posLeg.p() >= midP && posLeg.tpcNSigmaPi() <= pidcuts.minNSigmaPiHighP.value) {
157+
registry.fill(HIST("hCuts"), 2);
158+
return false;
159+
}
160+
}
161+
162+
// second let's check the negative leg
163+
if (negLeg.pt() <= trackcuts.minPt.value) {
164+
registry.fill(HIST("hCuts"), 1);
165+
return false;
166+
}
167+
if (negLeg.tpcFoundOverFindableCls() <= trackcuts.tpcFindOverFoundable.value) {
168+
registry.fill(HIST("hCuts"), 1);
169+
return false;
170+
}
171+
if (negLeg.tpcNSigmaEl() <= pidcuts.minNSigmaEl.value || negLeg.tpcNSigmaEl() >= pidcuts.maxNSigmaEl.value) {
172+
registry.fill(HIST("hCuts"), 2);
173+
return false;
174+
}
175+
176+
if (pidcuts.doPionRejection.value && negLeg.p() > minP) {
177+
if (negLeg.p() < midP && negLeg.tpcNSigmaPi() <= pidcuts.minNSigmaPiLowP.value) {
178+
registry.fill(HIST("hCuts"), 2);
179+
return false;
180+
} else if (negLeg.p() >= midP && negLeg.tpcNSigmaPi() <= pidcuts.minNSigmaPiHighP.value) {
181+
registry.fill(HIST("hCuts"), 2);
182+
return false;
183+
}
184+
}
185+
return true;
186+
}
187+
188+
// Pi0 from EMCal
189+
void process(o2::aod::Collisions const& collisions, o2::aod::Run2OTFV0s const& v0s, TracksWithPID const& /*tracks*/)
190+
{
191+
// TODO:
192+
// - Add TooCloseToV0 cut: if two V0s are within dAngle < 0.02 and dR < 6. then remove the V0 with higher Chi2
193+
// - Everything here!
194+
// nothing yet
195+
196+
const float zRslope = std::tan(2.f * std::atan(std::exp(-1.f * v0cuts.maxEta.value)));
197+
const float z0 = 7.f; // 7 cm
198+
const float maxZ = 10.f;
199+
const float minR = 5.f;
200+
const float maxR = 280.f;
201+
const float minRExclude = 55.f;
202+
const float maxRExclude = 72.f;
203+
204+
for (const auto& collision : collisions) {
205+
if (std::fabs(collision.posZ()) > maxZ) {
206+
continue;
207+
}
208+
auto photonsPerCollision = v0s.sliceBy(perCollisionV0, collision.globalIndex());
209+
210+
for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) {
211+
registry.fill(HIST("hCutsPair"), 0);
212+
// next V0 cuts
213+
float cX1 = g1.x();
214+
float cY1 = g1.y();
215+
float cZ1 = g1.z();
216+
float cR1 = std::sqrt(std::pow(cX1, 2.) + std::pow(cY1, 2.));
217+
218+
float cX2 = g2.x();
219+
float cY2 = g2.y();
220+
float cZ2 = g2.z();
221+
float cR2 = std::sqrt(std::pow(cX2, 2.) + std::pow(cY2, 2.));
222+
223+
ROOT::Math::PxPyPzEVector v4Photon2(g2.px(), g2.py(), g2.pz(), g2.e());
224+
ROOT::Math::PxPyPzEVector v4Photon1(g1.px(), g1.py(), g1.pz(), g1.e());
225+
226+
if (!checkLegs(g1)) {
227+
registry.fill(HIST("hCutsPair"), 1);
228+
continue;
229+
}
230+
231+
if (!checkLegs(g2)) {
232+
registry.fill(HIST("hCutsPair"), 1);
233+
continue;
234+
}
235+
236+
if (cR1 < minR || (cR1 > minRExclude && cR1 < maxRExclude) || cR1 > maxR) {
237+
registry.fill(HIST("hCutsPair"), 2);
238+
continue;
239+
}
240+
if (v4Photon1.Pt() < v0cuts.minPt.value) {
241+
registry.fill(HIST("hCutsPair"), 2);
242+
continue;
243+
}
244+
if (std::fabs(v4Photon1.Eta()) >= v0cuts.maxEta.value) {
245+
registry.fill(HIST("hCutsPair"), 2);
246+
continue;
247+
}
248+
if (std::fabs(cZ1) > v0cuts.maxZconv.value) {
249+
registry.fill(HIST("hCutsPair"), 2);
250+
continue;
251+
}
252+
if (cR1 <= ((std::fabs(cZ1) * zRslope) - z0)) {
253+
registry.fill(HIST("hCutsPair"), 2);
254+
continue;
255+
}
256+
if (g1.psiPair() >= (0.18f * std::exp(-0.055f * g1.chi2NDF()))) {
257+
registry.fill(HIST("hCutsPair"), 2);
258+
continue;
259+
}
260+
if (g1.qt() >= (v0cuts.qTFactor.value * v4Photon1.Pt())) {
261+
registry.fill(HIST("hCutsPair"), 2);
262+
continue;
263+
}
264+
if (getCosP(collision, g1) <= v0cuts.cosP.value) {
265+
registry.fill(HIST("hCutsPair"), 2);
266+
continue;
267+
}
268+
269+
if (cR2 < minR || (cR2 > minRExclude && cR2 < maxRExclude) || cR2 > maxR) {
270+
registry.fill(HIST("hCutsPair"), 2);
271+
continue;
272+
}
273+
if (v4Photon2.Pt() < v0cuts.minPt.value) {
274+
registry.fill(HIST("hCutsPair"), 2);
275+
continue;
276+
}
277+
if (std::fabs(v4Photon2.Eta()) >= v0cuts.maxEta.value) {
278+
registry.fill(HIST("hCutsPair"), 2);
279+
continue;
280+
}
281+
if (std::fabs(cZ2) > v0cuts.maxZconv.value) {
282+
registry.fill(HIST("hCutsPair"), 2);
283+
continue;
284+
}
285+
if (cR2 <= ((std::fabs(cZ2) * zRslope) - z0)) {
286+
registry.fill(HIST("hCutsPair"), 2);
287+
continue;
288+
}
289+
if (g2.psiPair() >= (0.18f * std::exp(-0.055f * g2.chi2NDF()))) {
290+
registry.fill(HIST("hCutsPair"), 2);
291+
continue;
292+
}
293+
if (g2.qt() >= (v0cuts.qTFactor.value * v4Photon2.Pt())) {
294+
registry.fill(HIST("hCutsPair"), 2);
295+
continue;
296+
}
297+
if (getCosP(collision, g2) <= v0cuts.cosP.value) {
298+
registry.fill(HIST("hCutsPair"), 2);
299+
continue;
300+
}
301+
302+
ROOT::Math::PxPyPzEVector vMeson = v4Photon1 + v4Photon2;
303+
registry.fill(HIST("hCutsPair"), 3);
304+
registry.fill(HIST("hSameEvent2D"), vMeson.M(), vMeson.Pt());
305+
} // end of loop over photon pairs
306+
} // end of loop over collision
307+
}
308+
PROCESS_SWITCH(PcmRun2, process, "Default process function", true);
309+
}; // End struct PcmRun2
310+
311+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
312+
{
313+
return WorkflowSpec{adaptAnalysisTask<PcmRun2>(cfgc)};
314+
}

0 commit comments

Comments
 (0)