Skip to content

Commit 6f434d2

Browse files
skundu692skundu692
authored
[PWGLF] Add event plane dependent efficiency from MC as weight in phi meson reconstruction for flow calculation (#9711)
1 parent 3f98ec8 commit 6f434d2

File tree

1 file changed

+98
-29
lines changed

1 file changed

+98
-29
lines changed

‎PWGLF/Tasks/Resonances/phipbpb.cxx‎

Lines changed: 98 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,7 @@
5454
#include "DataFormatsParameters/GRPObject.h"
5555
#include "DataFormatsParameters/GRPMagField.h"
5656
#include "CCDB/BasicCCDBManager.h"
57+
#include "CCDB/CcdbApi.h"
5758
#include "Common/DataModel/PIDResponseITS.h"
5859
#include "PWGMM/Mult/DataModel/Index.h" // for Particles2Tracks table
5960

@@ -66,7 +67,15 @@ struct phipbpb {
6667
int mRunNumber;
6768
int multEstimator;
6869
float d_bz;
70+
71+
struct : ConfigurableGroup {
72+
Configurable<std::string> cfgURL{"cfgURL", "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"};
73+
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"};
74+
} cfgCcdbParam;
75+
76+
// Enable access to the CCDB for the offset and correction constants and save them in dedicated variables.
6977
Service<o2::ccdb::BasicCCDBManager> ccdb;
78+
o2::ccdb::CcdbApi ccdbApi;
7079
Service<o2::framework::O2DatabasePDG> pdg;
7180

7281
// CCDB options
@@ -109,7 +118,7 @@ struct phipbpb {
109118
ConfigurableAxis configThnAxisCosThetaStar{"configThnAxisCosThetaStar", {10, -1.0, 1.}, "cos(#vartheta)"};
110119
ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {8, 0., 80}, "Centrality"};
111120
ConfigurableAxis configThnAxisPhiminusPsi{"configThnAxisPhiminusPsi", {6, 0.0, TMath::Pi()}, "#phi - #psi"};
112-
ConfigurableAxis configThnAxisV2{"configThnAxisV2", {200, -1, 1}, "V2"};
121+
ConfigurableAxis configThnAxisV2{"configThnAxisV2", {200, -6, 6}, "V2"};
113122
ConfigurableAxis configThnAxisRapidity{"configThnAxisRapidity", {8, 0, 0.8}, "Rapidity"};
114123
ConfigurableAxis configThnAxisSA{"configThnAxisSA", {200, -1, 1}, "SA"};
115124
ConfigurableAxis configThnAxiscosthetaSA{"configThnAxiscosthetaSA", {200, 0, 1}, "costhetaSA"};
@@ -118,6 +127,9 @@ struct phipbpb {
118127
Configurable<bool> genacceptancecut{"genacceptancecut", true, "use acceptance cut for generated"};
119128
Configurable<bool> avoidsplitrackMC{"avoidsplitrackMC", false, "avoid split track in MC"};
120129
Configurable<bool> islike{"islike", false, "use like"};
130+
Configurable<bool> useWeight{"useWeight", true, "use EP dep effi weight"};
131+
Configurable<std::string> ConfWeightPath{"ConfWeightPath", "Users/s/skundu/My/Object/mcweight", "Path to gain calibration"};
132+
121133
Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
122134
Filter centralityFilter = nabs(aod::cent::centFT0C) < cfgCutCentrality;
123135
Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPT);
@@ -232,10 +244,12 @@ struct phipbpb {
232244

233245
// weight
234246
histos.add("hImpactParameter", "Impact parameter", kTH1F, {{200, 0.0f, 20.0f}});
235-
histos.add("hEventPlaneAngle", "hEventPlaneAngle", kTH1F, {{100, 0.0f, 2.0f * TMath::Pi()}});
247+
histos.add("hEventPlaneAngle", "hEventPlaneAngle", kTH1F, {{200, -2.0f * TMath::Pi(), 2.0f * TMath::Pi()}});
248+
histos.add("hEventPlaneAngleRec", "hEventPlaneAngleRec", kTH1F, {{200, -2.0f * TMath::Pi(), 2.0f * TMath::Pi()}});
236249
histos.add("hNchVsImpactParameter", "hNchVsImpactParameter", kTH2F, {{200, 0.0f, 20.0f}, {500, -0.5f, 5000.5f}});
237-
histos.add("hSparseMCGenWeight", "hSparseMCGenWeight", HistType::kTHnSparseF, {thnAxisCentrality, {36, 0.0f, 2.0f * TMath::Pi()}, axisPtKaonWeight, {8, -0.8, 0.8}});
238-
histos.add("hSparseMCRecWeight", "hSparseMCRecWeight", HistType::kTHnSparseF, {thnAxisCentrality, {36, 0.0f, 2.0f * TMath::Pi()}, axisPtKaonWeight, {8, -0.8, 0.8}});
250+
histos.add("hSparseMCGenWeight", "hSparseMCGenWeight", HistType::kTHnSparseF, {thnAxisCentrality, {36, 0.0f, TMath::Pi()}, axisPtKaonWeight, {8, -0.8, 0.8}});
251+
histos.add("hSparseMCRecWeight", "hSparseMCRecWeight", HistType::kTHnSparseF, {thnAxisCentrality, {36, 0.0f, TMath::Pi()}, axisPtKaonWeight, {8, -0.8, 0.8}});
252+
histos.add("hSparseMCRecAllTrackWeight", "hSparseMCRecAllTrackWeight", HistType::kTHnSparseF, {thnAxisCentrality, {36, 0.0, TMath::Pi()}, axisPtKaonWeight, {8, -0.8, 0.8}});
239253
}
240254
// Event selection cut additional - Alex
241255
if (additionalEvsel) {
@@ -250,6 +264,12 @@ struct phipbpb {
250264
fMultMultPVCut = new TF1("fMultMultPVCut", "[0]+[1]*x+[2]*x*x", 0, 5000);
251265
fMultMultPVCut->SetParameters(-0.1, 0.785, -4.7e-05);
252266
}
267+
268+
ccdb->setURL(cfgCcdbParam.cfgURL);
269+
ccdbApi.init("http://alice-ccdb.cern.ch");
270+
ccdb->setCaching(true);
271+
ccdb->setLocalObjectValidityChecking();
272+
ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
253273
}
254274

255275
double massKa = o2::constants::physics::MassKPlus;
@@ -372,8 +392,10 @@ struct phipbpb {
372392
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::cent::CentFT0C, o2::aod::evsel::NumTracksInTimeRange>;
373393
ROOT::Math::PxPyPzMVector PhiMesonMother, KaonPlus, KaonMinus, fourVecDauCM;
374394
ROOT::Math::XYZVector threeVecDauCM, threeVecDauCMXY, eventplaneVec, eventplaneVecNorm, beamvector;
375-
376-
void processSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& /*tracks*/, aod::BCs const&)
395+
int currentRunNumber = -999;
396+
int lastRunNumber = -999;
397+
TH3D* hweight;
398+
void processSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& /*tracks, aod::BCs const&*/, aod::BCsWithTimestamps const&)
377399
{
378400
if (!collision.sel8() || !collision.triggereventep() || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || !collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) {
379401
return;
@@ -428,7 +450,16 @@ struct phipbpb {
428450

429451
histos.fill(HIST("hCentrality"), centrality);
430452
histos.fill(HIST("hVtxZ"), collision.posZ());
453+
454+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
455+
currentRunNumber = collision.bc_as<aod::BCsWithTimestamps>().runNumber();
456+
if (useWeight && (currentRunNumber != lastRunNumber)) {
457+
hweight = ccdb->getForTimeStamp<TH3D>(ConfWeightPath.value, bc.timestamp());
458+
}
459+
lastRunNumber = currentRunNumber;
431460
int Npostrack = 0;
461+
float weight1 = 1.0;
462+
float weight2 = 1.0;
432463
for (auto track1 : posThisColl) {
433464
// track selection
434465
if (!selectionTrack(track1)) {
@@ -451,6 +482,13 @@ struct phipbpb {
451482
histos.fill(HIST("hNsigmaKaonTPC"), track1.tpcNSigmaKa());
452483
histos.fill(HIST("hNsigmaKaonTOF"), track1.tofNSigmaKa());
453484
auto track1ID = track1.globalIndex();
485+
if (useWeight) {
486+
if (track1.pt() < 10.0 && track1.pt() > 0.15) {
487+
weight1 = hweight->GetBinContent(hweight->FindBin(centrality, GetPhiInRange(track1.phi() - psiFT0C), track1.pt() + 0.000005));
488+
} else {
489+
weight1 = 1;
490+
}
491+
}
454492
for (auto track2 : negThisColl) {
455493
// track selection
456494
if (!selectionTrack(track2)) {
@@ -482,6 +520,13 @@ struct phipbpb {
482520
if (useGlobalTrack && track2.p() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Kaon>(track2) > -2.5 && itsResponse.nSigmaITS<o2::track::PID::Kaon>(track2) < 2.5)) {
483521
continue;
484522
}
523+
if (useWeight) {
524+
if (track2.pt() < 10.0 && track2.pt() > 0.15) {
525+
weight2 = hweight->GetBinContent(hweight->FindBin(centrality, GetPhiInRange(track2.phi() - psiFT0C), track2.pt() + 0.000005));
526+
} else {
527+
weight2 = 1;
528+
}
529+
}
485530
KaonPlus = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa);
486531
KaonMinus = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massKa);
487532
PhiMesonMother = KaonPlus + KaonMinus;
@@ -490,15 +535,31 @@ struct phipbpb {
490535
auto v2sin = TMath::Sin(2.0 * phiminuspsi);
491536
auto phimother = PhiMesonMother.Phi();
492537
histos.fill(HIST("hpTvsRapidity"), PhiMesonMother.Pt(), PhiMesonMother.Rapidity());
538+
auto totalweight = weight1 * weight2;
539+
if (totalweight <= 0.0005) {
540+
totalweight = 1.0;
541+
}
542+
// LOGF(info, Form("weight %f %f",weight1, weight2));
493543
if (TMath::Abs(PhiMesonMother.Rapidity()) < confRapidity) {
494-
histos.fill(HIST("hSparseV2SameEventCosDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * QFT0C, centrality);
495-
histos.fill(HIST("hSparseV2SameEventCosDeltaPhiSquare"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * v2, centrality);
496-
histos.fill(HIST("hSparseV2SameEventSinDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2sin * QFT0C, centrality);
497-
498-
histos.fill(HIST("hSparseV2SameEventCosPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * phimother), centrality);
499-
histos.fill(HIST("hSparseV2SameEventSinPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * phimother), centrality);
500-
histos.fill(HIST("hSparseV2SameEventCosPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * psiFT0C), centrality);
501-
histos.fill(HIST("hSparseV2SameEventSinPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * psiFT0C), centrality);
544+
if (useWeight) {
545+
histos.fill(HIST("hSparseV2SameEventCosDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * QFT0C, centrality, 1 / totalweight);
546+
histos.fill(HIST("hSparseV2SameEventCosDeltaPhiSquare"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * v2, centrality, 1 / totalweight);
547+
histos.fill(HIST("hSparseV2SameEventSinDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2sin * QFT0C, centrality, 1 / totalweight);
548+
549+
histos.fill(HIST("hSparseV2SameEventCosPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * phimother), centrality, 1 / totalweight);
550+
histos.fill(HIST("hSparseV2SameEventSinPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * phimother), centrality, 1 / totalweight);
551+
histos.fill(HIST("hSparseV2SameEventCosPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * psiFT0C), centrality, 1 / totalweight);
552+
histos.fill(HIST("hSparseV2SameEventSinPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * psiFT0C), centrality, 1 / totalweight);
553+
} else {
554+
histos.fill(HIST("hSparseV2SameEventCosDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * QFT0C, centrality);
555+
histos.fill(HIST("hSparseV2SameEventCosDeltaPhiSquare"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2 * v2, centrality);
556+
histos.fill(HIST("hSparseV2SameEventSinDeltaPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), v2sin * QFT0C, centrality);
557+
558+
histos.fill(HIST("hSparseV2SameEventCosPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * phimother), centrality);
559+
histos.fill(HIST("hSparseV2SameEventSinPhi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * phimother), centrality);
560+
histos.fill(HIST("hSparseV2SameEventCosPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Cos(2.0 * psiFT0C), centrality);
561+
histos.fill(HIST("hSparseV2SameEventSinPsi"), PhiMesonMother.M(), PhiMesonMother.Pt(), TMath::Sin(2.0 * psiFT0C), centrality);
562+
}
502563
}
503564
if (fillSA) {
504565
ROOT::Math::Boost boost{PhiMesonMother.BoostToCM()};
@@ -510,8 +571,13 @@ struct phipbpb {
510571
auto cosPhistarminuspsi = GetPhiInRange(fourVecDauCM.Phi() - psiFT0C);
511572
auto SA = TMath::Cos(2.0 * cosPhistarminuspsi);
512573
auto cosThetaStar = eventplaneVecNorm.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(eventplaneVecNorm.Mag2());
513-
histos.fill(HIST("hSparseV2SameEventSA"), PhiMesonMother.M(), PhiMesonMother.Pt(), SA, TMath::Abs(PhiMesonMother.Rapidity()), centrality);
514-
histos.fill(HIST("hSparseV2SameEventCosThetaStar"), PhiMesonMother.M(), PhiMesonMother.Pt(), cosThetaStar, TMath::Abs(PhiMesonMother.Rapidity()), centrality);
574+
if (useWeight) {
575+
histos.fill(HIST("hSparseV2SameEventSA"), PhiMesonMother.M(), PhiMesonMother.Pt(), SA, TMath::Abs(PhiMesonMother.Rapidity()), centrality, 1 / totalweight);
576+
histos.fill(HIST("hSparseV2SameEventCosThetaStar"), PhiMesonMother.M(), PhiMesonMother.Pt(), cosThetaStar, TMath::Abs(PhiMesonMother.Rapidity()), centrality, 1 / totalweight);
577+
} else {
578+
histos.fill(HIST("hSparseV2SameEventSA"), PhiMesonMother.M(), PhiMesonMother.Pt(), SA, TMath::Abs(PhiMesonMother.Rapidity()), centrality);
579+
histos.fill(HIST("hSparseV2SameEventCosThetaStar"), PhiMesonMother.M(), PhiMesonMother.Pt(), cosThetaStar, TMath::Abs(PhiMesonMother.Rapidity()));
580+
}
515581
}
516582
}
517583
Npostrack = Npostrack + 1;
@@ -807,12 +873,11 @@ struct phipbpb {
807873

808874
} // process MC
809875
PROCESS_SWITCH(phipbpb, processMC, "Process MC", false);
810-
811876
using recoTracks = soa::Join<aod::TracksIU, aod::TracksExtra>;
812877
void processMCweight(aod::McCollision const& mcCollision, soa::Join<aod::McParticles, aod::ParticlesToTracks> const& mcParticles, recoTracks const&)
813878
{
814879
float imp = mcCollision.impactParameter();
815-
float evPhi = mcCollision.eventPlaneAngle();
880+
float evPhi = mcCollision.eventPlaneAngle() / 2.0;
816881
float centclass = -999;
817882
if (imp >= 0 && imp < 3.49) {
818883
centclass = 2.5;
@@ -841,8 +906,8 @@ struct phipbpb {
841906
if (imp >= 13.1 && imp < 14) {
842907
centclass = 75.0;
843908
}
844-
if (evPhi < 0)
845-
evPhi += 2. * TMath::Pi();
909+
// if (evPhi < 0)
910+
// evPhi += 2. * TMath::Pi();
846911

847912
int nCh = 0;
848913

@@ -851,36 +916,40 @@ struct phipbpb {
851916
histos.fill(HIST("hImpactParameter"), imp);
852917
histos.fill(HIST("hEventPlaneAngle"), evPhi);
853918
for (auto const& mcParticle : mcParticles) {
919+
920+
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
854921
// focus on bulk: e, mu, pi, k, p
855922
int pdgCode = TMath::Abs(mcParticle.pdgCode());
856923
if (checkAllCharge && pdgCode != 11 && pdgCode != 13 && pdgCode != 211 && pdgCode != 321 && pdgCode != 2212)
857924
continue;
858-
if (pdgCode != 321)
925+
if (!checkAllCharge && pdgCode != 321)
859926
continue;
860927
if (!mcParticle.isPhysicalPrimary())
861928
continue;
862929
if (TMath::Abs(mcParticle.eta()) > 0.8) // main acceptance
863930
continue;
864-
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
865-
if (deltaPhi < 0)
866-
deltaPhi += 2. * TMath::Pi();
867-
if (deltaPhi > 2. * TMath::Pi())
868-
deltaPhi -= 2. * TMath::Pi();
869-
870-
histos.fill(HIST("hSparseMCGenWeight"), centclass, deltaPhi, mcParticle.pt(), mcParticle.eta());
931+
histos.fill(HIST("hSparseMCGenWeight"), centclass, GetPhiInRange(deltaPhi), mcParticle.pt(), mcParticle.eta());
871932
nCh++;
872933
bool validGlobal = false;
934+
bool validAny = false;
873935
if (mcParticle.has_tracks()) {
874936
auto const& tracks = mcParticle.tracks_as<recoTracks>();
875937
for (auto const& track : tracks) {
876938
if (track.hasTPC() && track.hasITS()) {
877939
validGlobal = true;
878940
}
941+
if (track.hasTPC() || track.hasITS()) {
942+
validAny = true;
943+
}
879944
}
880945
}
881946
// if valid global, fill
882947
if (validGlobal) {
883-
histos.fill(HIST("hSparseMCRecWeight"), centclass, deltaPhi, mcParticle.pt(), mcParticle.eta());
948+
histos.fill(HIST("hSparseMCRecWeight"), centclass, GetPhiInRange(deltaPhi), mcParticle.pt(), mcParticle.eta());
949+
}
950+
if (validAny) {
951+
histos.fill(HIST("hSparseMCRecAllTrackWeight"), centclass, GetPhiInRange(deltaPhi), mcParticle.pt(), mcParticle.eta());
952+
histos.fill(HIST("hEventPlaneAngleRec"), GetPhiInRange(deltaPhi));
884953
}
885954
// if any track present, fill
886955
}

0 commit comments

Comments
 (0)