Skip to content

Commit 640f9cd

Browse files
sangwoo184sangwoo
sangwoo184
and
sangwoo
authored
[PWGLF] Add Rotational background method (AliceO2Group#10159)
Co-authored-by: sangwoo <[email protected]>
1 parent a50cd6b commit 640f9cd

File tree

1 file changed

+78
-51
lines changed

1 file changed

+78
-51
lines changed

‎PWGLF/Tasks/Resonances/f0980pbpbanalysis.cxx‎

Lines changed: 78 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,16 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12+
/// \file f0980pbpbanalysis.cxx
13+
/// \brief f0980 resonance analysis in PbPb collisions
1214
/// \author Junlee Kim ([email protected])
1315

1416
#include <Framework/Configurable.h>
1517
#include <cmath>
1618
#include <array>
1719
#include <cstdlib>
1820
#include <chrono>
19-
#include <iostream>
21+
// #include <iostream>
2022
#include <string>
2123

2224
#include "TLorentzVector.h"
@@ -73,11 +75,11 @@ struct f0980pbpbanalysis {
7375
o2::ccdb::CcdbApi ccdbApi;
7476

7577
Configurable<std::string> cfgURL{"cfgURL", "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"};
76-
Configurable<int64_t> nolaterthan{"ccdb-no-later-than", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "Latest acceptable timestamp of creation for the object"};
78+
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "Latest acceptable timestamp of creation for the object"};
7779

7880
Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0, "PV selection"};
7981
Configurable<bool> cfgQvecSel{"cfgQvecSel", true, "Reject events when no QVector"};
80-
Configurable<bool> cfgOccupancySel{"cfgOccupancySe", false, "Occupancy selection"};
82+
Configurable<bool> cfgOccupancySel{"cfgOccupancySel", false, "Occupancy selection"};
8183
Configurable<int> cfgMaxOccupancy{"cfgMaxOccupancy", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
8284
Configurable<int> cfgMinOccupancy{"cfgMinOccupancy", -100, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
8385
Configurable<bool> cfgNCollinTR{"cfgNCollinTR", false, "Additional selection for the number of coll in time range"};
@@ -105,7 +107,7 @@ struct f0980pbpbanalysis {
105107
Configurable<double> cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 5.0, "TPC nSigma cut for Pion"}; // TPC
106108
Configurable<double> cMaxTPCnSigmaPionS{"cMaxTPCnSigmaPionS", 3.0, "TPC nSigma cut for Pion as a standalone"};
107109
Configurable<bool> cfgUSETOF{"cfgUSETOF", false, "TPC usage"};
108-
Configurable<int> SelectType{"SelectType", 0, "PID selection type"};
110+
Configurable<int> cfgSelectType{"cfgSelectType", 0, "PID selection type"};
109111

110112
Configurable<int> cfgnMods{"cfgnMods", 1, "The number of modulations of interest starting from 2"};
111113
Configurable<int> cfgNQvec{"cfgNQvec", 7, "The number of total Qvectors for looping over the task"};
@@ -114,37 +116,44 @@ struct f0980pbpbanalysis {
114116
Configurable<std::string> cfgQvecRefAName{"cfgQvecRefAName", "TPCpos", "The name of detector for reference A"};
115117
Configurable<std::string> cfgQvecRefBName{"cfgQvecRefBName", "TPCneg", "The name of detector for reference B"};
116118

119+
Configurable<bool> cfgRotBkg{"cfgRotBkg", true, "flag to construct rotational backgrounds"};
120+
Configurable<int> cfgNRotBkg{"cfgNRotBkg", 10, "the number of rotational backgrounds"};
121+
117122
ConfigurableAxis massAxis{"massAxis", {400, 0.2, 2.2}, "Invariant mass axis"};
118123
ConfigurableAxis ptAxis{"ptAxis", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.0, 13.0, 20.0}, "Transverse momentum Binning"};
119124
ConfigurableAxis centAxis{"centAxis", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100}, "Centrality interval"};
120125

121126
TF1* fMultPVCutLow = nullptr;
122127
TF1* fMultPVCutHigh = nullptr;
123128

124-
int DetId;
125-
int RefAId;
126-
int RefBId;
129+
int detId;
130+
int refAId;
131+
int refBId;
127132

128-
int QvecDetInd;
129-
int QvecRefAInd;
130-
int QvecRefBInd;
133+
int qVecDetInd;
134+
int qVecRefAInd;
135+
int qVecRefBInd;
131136

132137
float centrality;
133138

134139
double angle;
135140
double relPhi;
141+
double relPhiRot;
136142

137143
double massPi = o2::constants::physics::MassPionCharged;
138144

145+
TRandom* rn = new TRandom();
146+
// float theta2;
147+
139148
Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
140149
Filter acceptanceFilter = (nabs(aod::track::eta) < cfgMaxEta && nabs(aod::track::pt) > cfgMinPt);
141-
Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgMaxDCArToPVcut) && (nabs(aod::track::dcaZ) < cfgMaxDCAzToPVcut);
150+
Filter cutDCAFilter = (nabs(aod::track::dcaXY) < cfgMaxDCArToPVcut) && (nabs(aod::track::dcaZ) < cfgMaxDCAzToPVcut);
142151

143152
using EventCandidates = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::TPCMults, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0Cs, aod::CentFT0As, aod::Mults, aod::Qvectors>>;
144153
using TrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullPi, aod::pidTOFFullPi>>;
145154

146155
template <typename T>
147-
int GetDetId(const T& name)
156+
int getDetId(const T& name)
148157
{
149158
if (name.value == "FT0C") {
150159
return 0;
@@ -188,7 +197,7 @@ struct f0980pbpbanalysis {
188197
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
189198
return 0;
190199
}
191-
if (cfgQvecSel && (collision.qvecAmp()[DetId] < 1e-4 || collision.qvecAmp()[RefAId] < 1e-4 || collision.qvecAmp()[RefAId] < 1e-4)) {
200+
if (cfgQvecSel && (collision.qvecAmp()[detId] < 1e-4 || collision.qvecAmp()[refAId] < 1e-4 || collision.qvecAmp()[refAId] < 1e-4)) {
192201
return 0;
193202
}
194203
if (cfgOccupancySel && (collision.trackOccupancyInTimeRange() > cfgMaxOccupancy || collision.trackOccupancyInTimeRange() < cfgMinOccupancy)) {
@@ -239,9 +248,9 @@ struct f0980pbpbanalysis {
239248
}
240249

241250
template <typename TrackType>
242-
bool PIDSelected(const TrackType track)
251+
bool selectionPID(const TrackType track)
243252
{
244-
if (SelectType == 0) {
253+
if (cfgSelectType == 0) {
245254
if (cfgUSETOF) {
246255
if (std::fabs(track.tofNSigmaPi()) > cMaxTOFnSigmaPion) {
247256
return 0;
@@ -254,7 +263,7 @@ struct f0980pbpbanalysis {
254263
return 0;
255264
}
256265
}
257-
if (SelectType == 1) {
266+
if (cfgSelectType == 1) {
258267
if (cfgUSETOF) {
259268
if (track.hasTOF()) {
260269
if (std::fabs(track.tofNSigmaPi()) > cMaxTOFnSigmaPion) {
@@ -279,37 +288,43 @@ struct f0980pbpbanalysis {
279288
}
280289

281290
template <bool IsMC, typename CollisionType, typename TracksType>
282-
void FillHistograms(const CollisionType& collision,
291+
void fillHistograms(const CollisionType& collision,
283292
const TracksType& dTracks, int nmode)
284293
{
285-
QvecDetInd = DetId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
286-
QvecRefAInd = RefAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
287-
QvecRefBInd = RefBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
294+
qVecDetInd = detId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
295+
qVecRefAInd = refAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
296+
qVecRefBInd = refBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
288297

289-
double eventPlaneDet = TMath::ATan2(collision.qvecIm()[QvecDetInd], collision.qvecRe()[QvecDetInd]) / static_cast<float>(nmode);
290-
double eventPlaneRefA = TMath::ATan2(collision.qvecIm()[QvecRefAInd], collision.qvecRe()[QvecRefAInd]) / static_cast<float>(nmode);
291-
double eventPlaneRefB = TMath::ATan2(collision.qvecIm()[QvecRefBInd], collision.qvecRe()[QvecRefBInd]) / static_cast<float>(nmode);
298+
double eventPlaneDet = std::atan2(collision.qvecIm()[qVecDetInd], collision.qvecRe()[qVecDetInd]) / static_cast<float>(nmode);
299+
double eventPlaneRefA = std::atan2(collision.qvecIm()[qVecRefAInd], collision.qvecRe()[qVecRefAInd]) / static_cast<float>(nmode);
300+
double eventPlaneRefB = std::atan2(collision.qvecIm()[qVecRefBInd], collision.qvecRe()[qVecRefBInd]) / static_cast<float>(nmode);
292301

293302
histos.fill(HIST("QA/EPhist"), centrality, eventPlaneDet);
294-
histos.fill(HIST("QA/EPResAB"), centrality, TMath::Cos(static_cast<float>(nmode) * (eventPlaneDet - eventPlaneRefA)));
295-
histos.fill(HIST("QA/EPResAC"), centrality, TMath::Cos(static_cast<float>(nmode) * (eventPlaneDet - eventPlaneRefB)));
296-
histos.fill(HIST("QA/EPResBC"), centrality, TMath::Cos(static_cast<float>(nmode) * (eventPlaneRefA - eventPlaneRefB)));
303+
histos.fill(HIST("QA/EPResAB"), centrality, std::cos(static_cast<float>(nmode) * (eventPlaneDet - eventPlaneRefA)));
304+
histos.fill(HIST("QA/EPResAC"), centrality, std::cos(static_cast<float>(nmode) * (eventPlaneDet - eventPlaneRefB)));
305+
histos.fill(HIST("QA/EPResBC"), centrality, std::cos(static_cast<float>(nmode) * (eventPlaneRefA - eventPlaneRefB)));
306+
307+
TLorentzVector pion1, pion2, pion2Rot, reco, recoRot;
308+
for (const auto& trk1 : dTracks) {
309+
if (!trackSelected(trk1)) {
310+
continue;
311+
}
297312

298-
TLorentzVector Pion1, Pion2, Reco;
299-
for (auto& trk1 : dTracks) {
300-
if (!trackSelected(trk1))
313+
if (!selectionPID(trk1)) {
301314
continue;
315+
}
316+
302317
histos.fill(HIST("QA/Nsigma_TPC"), trk1.pt(), trk1.tpcNSigmaPi());
303318
histos.fill(HIST("QA/Nsigma_TOF"), trk1.pt(), trk1.tofNSigmaPi());
304319
histos.fill(HIST("QA/TPC_TOF"), trk1.tpcNSigmaPi(), trk1.tofNSigmaPi());
305320

306-
for (auto& trk2 : dTracks) {
307-
if (!trackSelected(trk1) || !trackSelected(trk2)) {
321+
for (const auto& trk2 : dTracks) {
322+
if (!trackSelected(trk2)) {
308323
continue;
309324
}
310325

311326
// PID
312-
if (!PIDSelected(trk1) || !PIDSelected(trk2)) {
327+
if (!selectionPID(trk2)) {
313328
continue;
314329
}
315330

@@ -319,35 +334,46 @@ struct f0980pbpbanalysis {
319334
histos.fill(HIST("QA/TPC_TOF_selected"), trk1.tpcNSigmaPi(), trk1.tofNSigmaPi());
320335
}
321336

322-
Pion1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massPi);
323-
Pion2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi);
324-
Reco = Pion1 + Pion2;
337+
pion1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massPi);
338+
pion2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi);
339+
reco = pion1 + pion2;
325340

326-
if (Reco.Rapidity() > cfgMaxRap || Reco.Rapidity() < cfgMinRap) {
341+
if (reco.Rapidity() > cfgMaxRap || reco.Rapidity() < cfgMinRap) {
327342
continue;
328343
}
329344

330-
relPhi = TVector2::Phi_0_2pi((Reco.Phi() - eventPlaneDet) * static_cast<float>(nmode));
345+
relPhi = TVector2::Phi_0_2pi((reco.Phi() - eventPlaneDet) * static_cast<float>(nmode));
331346

332347
if (trk1.sign() * trk2.sign() < 0) {
333-
histos.fill(HIST("hInvMass_f0980_US_EPA"), Reco.M(), Reco.Pt(), centrality, relPhi);
348+
histos.fill(HIST("hInvMass_f0980_US_EPA"), reco.M(), reco.Pt(), centrality, relPhi);
334349
} else if (trk1.sign() > 0 && trk2.sign() > 0) {
335-
histos.fill(HIST("hInvMass_f0980_LSpp_EPA"), Reco.M(), Reco.Pt(), centrality, relPhi);
350+
histos.fill(HIST("hInvMass_f0980_LSpp_EPA"), reco.M(), reco.Pt(), centrality, relPhi);
336351
} else if (trk1.sign() < 0 && trk2.sign() < 0) {
337-
histos.fill(HIST("hInvMass_f0980_LSmm_EPA"), Reco.M(), Reco.Pt(), centrality, relPhi);
352+
histos.fill(HIST("hInvMass_f0980_LSmm_EPA"), reco.M(), reco.Pt(), centrality, relPhi);
353+
}
354+
355+
if (cfgRotBkg && trk1.sign() * trk2.sign() < 0) {
356+
for (int nr = 0; nr < cfgNRotBkg; nr++) {
357+
auto randomPhi = rn->Uniform(o2::constants::math::PI * 5.0 / 6.0, o2::constants::math::PI * 7.0 / 6.0);
358+
randomPhi += pion2.Phi();
359+
pion2Rot.SetXYZM(pion2.Pt() * std::cos(randomPhi), pion2.Pt() * std::sin(randomPhi), trk2.pz(), massPi);
360+
recoRot = pion1 + pion2Rot;
361+
relPhiRot = TVector2::Phi_0_2pi((recoRot.Phi() - eventPlaneDet) * static_cast<float>(nmode));
362+
histos.fill(HIST("hInvMass_f0980_USRot_EPA"), recoRot.M(), recoRot.Pt(), centrality, relPhiRot);
363+
}
338364
}
339365
}
340366
}
341367
}
342368

343369
void init(o2::framework::InitContext&)
344370
{
345-
AxisSpec epAxis = {6, 0.0, 2.0 * constants::math::PI};
371+
AxisSpec epAxis = {6, 0.0, 2.0 * o2::constants::math::PI};
346372
AxisSpec centQaAxis = {110, 0, 110};
347373
AxisSpec vzQaAxis = {100, -20, 20};
348374
AxisSpec PIDqaAxis = {100, -10, 10};
349375
AxisSpec pTqaAxis = {200, 0, 20};
350-
AxisSpec epQaAxis = {100, -1.0 * constants::math::PI, constants::math::PI};
376+
AxisSpec epQaAxis = {100, -1.0 * o2::constants::math::PI, o2::constants::math::PI};
351377
AxisSpec epresAxis = {102, -1.02, 1.02};
352378

353379
histos.add("QA/CentDist", "", {HistType::kTH1F, {centQaAxis}});
@@ -372,21 +398,22 @@ struct f0980pbpbanalysis {
372398
{HistType::kTHnSparseF, {massAxis, ptAxis, centAxis, epAxis}});
373399
histos.add("hInvMass_f0980_LSmm_EPA", "-- invariant mass",
374400
{HistType::kTHnSparseF, {massAxis, ptAxis, centAxis, epAxis}});
375-
401+
histos.add("hInvMass_f0980_USRot_EPA", "unlike invariant mass Rotation",
402+
{HistType::kTHnSparseF, {massAxis, ptAxis, centAxis, epAxis}});
376403
// if (doprocessMCLight) {
377404
// histos.add("MCL/hpT_f0980_GEN", "generated f0 signals", HistType::kTH1F, {pTqaAxis});
378405
// histos.add("MCL/hpT_f0980_REC", "reconstructed f0 signals", HistType::kTH3F, {massAxis, pTqaAxis, centAxis});
379406
// }
380407

381-
DetId = GetDetId(cfgQvecDetName);
382-
RefAId = GetDetId(cfgQvecRefAName);
383-
RefBId = GetDetId(cfgQvecRefBName);
408+
detId = getDetId(cfgQvecDetName);
409+
refAId = getDetId(cfgQvecRefAName);
410+
refBId = getDetId(cfgQvecRefBName);
384411

385-
if (DetId == RefAId || DetId == RefBId || RefAId == RefBId) {
412+
if (detId == refAId || detId == refBId || refAId == refBId) {
386413
LOGF(info, "Wrong detector configuration \n The FT0C will be used to get Q-Vector \n The TPCpos and TPCneg will be used as reference systems");
387-
DetId = 0;
388-
RefAId = 4;
389-
RefBId = 5;
414+
detId = 0;
415+
refAId = 4;
416+
refBId = 5;
390417
}
391418

392419
fMultPVCutLow = new TF1("fMultPVCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 2.5*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100);
@@ -415,7 +442,7 @@ struct f0980pbpbanalysis {
415442
histos.fill(HIST("QA/CentDist"), centrality, 1.0);
416443
histos.fill(HIST("QA/Vz"), collision.posZ(), 1.0);
417444

418-
FillHistograms<false>(collision, tracks, 2); // second order
445+
fillHistograms<false>(collision, tracks, 2); // second order
419446
};
420447
PROCESS_SWITCH(f0980pbpbanalysis, processData, "Process Event for data", true);
421448
};

0 commit comments

Comments
 (0)