Skip to content

Commit 0daa29d

Browse files
JosueMtzGarJosueMtzGar
authored
[PWGUD] Adding correlation study (#11333)
1 parent e1e8b42 commit 0daa29d

File tree

3 files changed

+615
-115
lines changed

3 files changed

+615
-115
lines changed

‎PWGUD/Core/UPCPairCuts.h‎

Lines changed: 346 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,346 @@
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+
/// \brief Functions which cut on particle pairs (decays, conversions, two-track cuts) adapted for data from UD tables
13+
/// Based on the code "PWGCF/Core/PairCuts.h" made by Jan Fiete Grosse-Oetringhaus
14+
/// Author:
15+
16+
#ifndef PWGUD_CORE_UPCPAIRCUTS_H_
17+
#define PWGUD_CORE_UPCPAIRCUTS_H_
18+
19+
#include <cmath>
20+
21+
#include "Framework/Logger.h"
22+
#include "Framework/HistogramRegistry.h"
23+
#include "CommonConstants/MathConstants.h"
24+
#include "CommonConstants/PhysicsConstants.h"
25+
26+
#include "PWGUD/Core/UPCTauCentralBarrelHelperRL.h"
27+
28+
using namespace o2;
29+
using namespace o2::framework;
30+
using namespace constants::math;
31+
32+
class UPCPairCuts
33+
{
34+
public:
35+
enum Particle { Photon = 0,
36+
K0,
37+
Lambda,
38+
Phi,
39+
Rho,
40+
ParticlesLastEntry };
41+
42+
void setHistogramRegistry(HistogramRegistry* registry) { histogramRegistry = registry; }
43+
44+
void setPairCut(Particle particle, float cut)
45+
{
46+
LOGF(info, "Enabled pair cut for %d with value %f", static_cast<int>(particle), cut);
47+
mCuts[particle] = cut;
48+
if (histogramRegistry != nullptr && histogramRegistry->contains(HIST("ControlConvResonances")) == false) {
49+
histogramRegistry->add("ControlConvResonances", "", {HistType::kTH2F, {{6, -0.5, 5.5, "id"}, {500, -0.5, 0.5, "delta mass"}}});
50+
}
51+
}
52+
53+
void setTwoTrackCuts(float distance = 0.02f, float radius = 0.8f)
54+
{
55+
LOGF(info, "Enabled two-track cut with distance %f and radius %f", distance, radius);
56+
mTwoTrackDistance = distance;
57+
mTwoTrackRadius = radius;
58+
59+
if (histogramRegistry != nullptr && histogramRegistry->contains(HIST("TwoTrackDistancePt_0")) == false) {
60+
histogramRegistry->add("TwoTrackDistancePt_0", "", {HistType::kTH3F, {{100, -0.15, 0.15, "#Delta#eta"}, {100, -0.05, 0.05, "#Delta#varphi^{*}_{min}"}, {20, 0, 10, "#Delta p_{T}"}}});
61+
histogramRegistry->addClone("TwoTrackDistancePt_0", "TwoTrackDistancePt_1");
62+
}
63+
}
64+
65+
template <typename T>
66+
bool conversionCuts(T const& track1, T const& track2);
67+
68+
template <typename T>
69+
bool twoTrackCut(T const& track1, T const& track2);
70+
71+
protected:
72+
float mCuts[ParticlesLastEntry] = {-1};
73+
float mTwoTrackDistance = -1; // distance below which the pair is flagged as to be removed
74+
float mTwoTrackRadius = 0.8f; // radius at which the two track cuts are applied
75+
int magField = 5; // magField: B field in kG
76+
77+
HistogramRegistry* histogramRegistry = nullptr; // if set, control histograms are stored here
78+
79+
template <typename T>
80+
bool conversionCut(T const& track1, T const& track2, Particle conv, double cut);
81+
82+
template <typename T>
83+
double getInvMassSquared(T const& track1, double m0_1, T const& track2, double m0_2);
84+
85+
template <typename T>
86+
double getInvMassSquaredFast(T const& track1, double m0_1, T const& track2, double m0_2);
87+
88+
template <typename T>
89+
float getDPhiStar(T const& track1, T const& track2, float radius, int magField);
90+
};
91+
92+
template <typename T>
93+
bool UPCPairCuts::conversionCuts(T const& track1, T const& track2)
94+
{
95+
// skip if like sign
96+
if (track1.sign() * track2.sign() > 0) {
97+
return false;
98+
}
99+
100+
for (int i = 0; i < static_cast<int>(ParticlesLastEntry); i++) {
101+
Particle particle = static_cast<Particle>(i);
102+
if (mCuts[i] > 0) {
103+
if (conversionCut(track1, track2, particle, mCuts[i])) {
104+
return true;
105+
}
106+
if (particle == Lambda) {
107+
if (conversionCut(track2, track1, particle, mCuts[i])) {
108+
return true;
109+
}
110+
}
111+
}
112+
}
113+
114+
return false;
115+
}
116+
117+
template <typename T>
118+
bool UPCPairCuts::twoTrackCut(T const& track1, T const& track2)
119+
{
120+
// the variables & cut have been developed in Run 1 by the CF - HBT group
121+
//
122+
// Parameters:
123+
// magField: B field in kG
124+
125+
auto deta = eta(track1.px(), track1.py(), track1.pz()) - eta(track2.px(), track2.py(), track2.pz());
126+
127+
// optimization
128+
if (std::fabs(deta) < mTwoTrackDistance * 2.5 * 3) {
129+
// check first boundaries to see if is worth to loop and find the minimum
130+
float dphistar1 = getDPhiStar(track1, track2, mTwoTrackRadius, magField);
131+
float dphistar2 = getDPhiStar(track1, track2, 2.5, magField);
132+
133+
const float kLimit = mTwoTrackDistance * 3;
134+
135+
if (std::fabs(dphistar1) < kLimit || std::fabs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0) {
136+
float dphistarminabs = 1e5;
137+
float dphistarmin = 1e5;
138+
for (Double_t rad = mTwoTrackRadius; rad < 2.51; rad += 0.01) {
139+
float dphistar = getDPhiStar(track1, track2, rad, magField);
140+
141+
float dphistarabs = std::fabs(dphistar);
142+
143+
if (dphistarabs < dphistarminabs) {
144+
dphistarmin = dphistar;
145+
dphistarminabs = dphistarabs;
146+
}
147+
}
148+
149+
if (histogramRegistry != nullptr) {
150+
histogramRegistry->fill(HIST("TwoTrackDistancePt_0"), deta, dphistarmin, std::fabs(track1.pt() - track2.pt()));
151+
}
152+
153+
if (dphistarminabs < mTwoTrackDistance && std::fabs(deta) < mTwoTrackDistance) {
154+
// LOGF(debug, "Removed track pair %ld %ld with %f %f %f %f %d %f %f %d %d", track1.index(), track2.index(), deta, dphistarminabs, track1.phi2(), track1.pt(), track1.sign(), track2.phi2(), track2.pt(), track2.sign(), magField);
155+
return true;
156+
}
157+
158+
if (histogramRegistry != nullptr) {
159+
histogramRegistry->fill(HIST("TwoTrackDistancePt_1"), deta, dphistarmin, std::fabs(track1.pt() - track2.pt()));
160+
}
161+
}
162+
}
163+
164+
return false;
165+
}
166+
167+
template <typename T>
168+
bool UPCPairCuts::conversionCut(T const& track1, T const& track2, Particle conv, double cut)
169+
{
170+
// LOGF(info, "pt is %f %f", track1.pt(), track2.pt());
171+
172+
if (cut < 0) {
173+
return false;
174+
}
175+
176+
double massD1, massD2, massM;
177+
178+
switch (conv) {
179+
case Photon:
180+
massD1 = o2::constants::physics::MassElectron;
181+
massD2 = o2::constants::physics::MassElectron;
182+
massM = 0;
183+
break;
184+
case K0:
185+
massD1 = o2::constants::physics::MassPiPlus;
186+
massD2 = o2::constants::physics::MassPiPlus;
187+
massM = o2::constants::physics::MassK0;
188+
break;
189+
case Lambda:
190+
massD1 = o2::constants::physics::MassProton;
191+
massD2 = o2::constants::physics::MassPiPlus;
192+
massM = o2::constants::physics::MassLambda0;
193+
break;
194+
case Phi:
195+
massD1 = o2::constants::physics::MassKPlus;
196+
massD2 = o2::constants::physics::MassKPlus;
197+
massM = o2::constants::physics::MassPhi;
198+
break;
199+
case Rho:
200+
massD1 = o2::constants::physics::MassPiPlus;
201+
massD2 = o2::constants::physics::MassPiPlus;
202+
massM = 0.770;
203+
break;
204+
default:
205+
LOGF(fatal, "Particle now known");
206+
return false;
207+
break;
208+
}
209+
210+
auto massC = getInvMassSquaredFast(track1, massD1, track2, massD2);
211+
212+
if (std::fabs(massC - massM * massM) > cut * 5) {
213+
return false;
214+
}
215+
216+
massC = getInvMassSquared(track1, massD1, track2, massD2);
217+
218+
if (histogramRegistry != nullptr) {
219+
histogramRegistry->fill(HIST("ControlConvResonances"), static_cast<int>(conv), massC - massM * massM);
220+
}
221+
222+
if (massC > (massM - cut) * (massM - cut) && massC < (massM + cut) * (massM + cut)) {
223+
return true;
224+
}
225+
226+
return false;
227+
}
228+
229+
template <typename T>
230+
double UPCPairCuts::getInvMassSquared(T const& track1, double m0_1, T const& track2, double m0_2)
231+
{
232+
// calculate inv mass squared
233+
// same can be achieved, but with more computing time with
234+
/*TLorentzVector photon, p1, p2;
235+
p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
236+
p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
237+
photon = p1+p2;
238+
photon.M()*/
239+
240+
float tantheta1 = 1e10;
241+
242+
if (eta(track1.px(), track1.py(), track1.pz()) < -1e-10 || eta(track1.px(), track1.py(), track1.pz()) > 1e-10) {
243+
float expTmp = std::exp(-eta(track1.px(), track1.py(), track1.pz()));
244+
tantheta1 = 2.0 * expTmp / (1.0 - expTmp * expTmp);
245+
}
246+
247+
float tantheta2 = 1e10;
248+
if (eta(track2.px(), track2.py(), track2.pz()) < -1e-10 || eta(track2.px(), track2.py(), track2.pz()) > 1e-10) {
249+
float expTmp = std::exp(-eta(track2.px(), track2.py(), track2.pz()));
250+
tantheta2 = 2.0 * expTmp / (1.0 - expTmp * expTmp);
251+
}
252+
253+
float e1squ = m0_1 * m0_1 + track1.pt() * track1.pt() * (1.0 + 1.0 / tantheta1 / tantheta1);
254+
float e2squ = m0_2 * m0_2 + track2.pt() * track2.pt() * (1.0 + 1.0 / tantheta2 / tantheta2);
255+
256+
float mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * (std::sqrt(e1squ * e2squ) - (track1.pt() * track2.pt() * (std::cos(phi(track1.px(), track1.py()) - phi(track2.px(), track2.py())) + 1.0 / tantheta1 / tantheta2)));
257+
258+
// LOGF(debug, "%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2);
259+
260+
return mass2;
261+
}
262+
263+
template <typename T>
264+
double UPCPairCuts::getInvMassSquaredFast(T const& track1, double m0_1, T const& track2, double m0_2)
265+
{
266+
// calculate inv mass squared approximately
267+
268+
const float eta1 = eta(track1.px(), track1.py(), track1.pz());
269+
const float eta2 = eta(track2.px(), track2.py(), track2.pz());
270+
const float phi1 = phi(track1.px(), track1.py());
271+
const float phi2 = phi(track2.px(), track2.py());
272+
const float pt1 = track1.pt();
273+
const float pt2 = track2.pt();
274+
275+
float tantheta1 = 1e10f;
276+
277+
if (eta1 < -1e-10f || eta1 > 1e-10f) {
278+
float expTmp = 1.0f - eta1 + eta1 * eta1 / 2.0f - eta1 * eta1 * eta1 / 6.0f + eta1 * eta1 * eta1 * eta1 / 24.0f;
279+
tantheta1 = 2.0f * expTmp / (1.0f - expTmp * expTmp);
280+
}
281+
282+
float tantheta2 = 1e10f;
283+
if (eta2 < -1e-10f || eta2 > 1e-10f) {
284+
float expTmp = 1.0f - eta2 + eta2 * eta2 / 2.0f - eta2 * eta2 * eta2 / 6.0f + eta2 * eta2 * eta2 * eta2 / 24.0f;
285+
tantheta2 = 2.0f * expTmp / (1.0f - expTmp * expTmp);
286+
}
287+
288+
float e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0f + 1.0f / tantheta1 / tantheta1);
289+
float e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0f + 1.0f / tantheta2 / tantheta2);
290+
291+
// fold onto 0...pi
292+
float deltaPhi = std::fabs(phi1 - phi2);
293+
while (deltaPhi > TwoPI) {
294+
deltaPhi -= TwoPI;
295+
}
296+
if (deltaPhi > PI) {
297+
deltaPhi = TwoPI - deltaPhi;
298+
}
299+
300+
float cosDeltaPhi = 0;
301+
if (deltaPhi < PI / 3.0f) {
302+
cosDeltaPhi = 1.0 - deltaPhi * deltaPhi / 2 + deltaPhi * deltaPhi * deltaPhi * deltaPhi / 24;
303+
} else if (deltaPhi < 2.0f * PI / 3.0f) {
304+
cosDeltaPhi = -(deltaPhi - PI / 2) + 1.0 / 6 * std::pow((deltaPhi - PI / 2), 3);
305+
} else {
306+
cosDeltaPhi = -1.0f + 1.0f / 2.0f * (deltaPhi - PI) * (deltaPhi - PI) - 1.0f / 24.0f * std::pow(deltaPhi - PI, 4.0f);
307+
}
308+
309+
double mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2.0f * (std::sqrt(e1squ * e2squ) - (pt1 * pt2 * (cosDeltaPhi + 1.0f / tantheta1 / tantheta2)));
310+
311+
// LOGF(debug, "%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2);
312+
313+
return mass2;
314+
}
315+
316+
template <typename T>
317+
float UPCPairCuts::getDPhiStar(T const& track1, T const& track2, float radius, int magField)
318+
{
319+
//
320+
// calculates dphistar
321+
//
322+
323+
auto phi1 = phi(track1.px(), track1.py());
324+
auto pt1 = track1.pt();
325+
auto charge1 = track1.sign();
326+
327+
auto phi2 = phi(track2.px(), track2.py());
328+
auto pt2 = track2.pt();
329+
auto charge2 = track2.sign();
330+
331+
float dphistar = phi1 - phi2 - charge1 * std::asin(0.015 * magField * radius / pt1) + charge2 * std::asin(0.015 * magField * radius / pt2);
332+
333+
if (dphistar > PI) {
334+
dphistar = TwoPI - dphistar;
335+
}
336+
if (dphistar < -PI) {
337+
dphistar = -TwoPI - dphistar;
338+
}
339+
if (dphistar > PI) { // might look funny but is needed
340+
dphistar = TwoPI - dphistar;
341+
}
342+
343+
return dphistar;
344+
}
345+
346+
#endif // PWGUD_CORE_UPCPAIRCUTS_H_

‎PWGUD/Tasks/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -146,7 +146,7 @@ o2physics_add_dpl_workflow(exclusive-phi
146146

147147
o2physics_add_dpl_workflow(upc-photonuclear-jmg
148148
SOURCES upcPhotonuclearAnalysisJMG.cxx
149-
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsBase O2::DetectorsCommonDataFormats
149+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore O2::ReconstructionDataFormats O2::DetectorsBase O2::DetectorsCommonDataFormats
150150
COMPONENT_NAME Analysis)
151151

152152
o2physics_add_dpl_workflow(exclusive-two-protons

0 commit comments

Comments
 (0)