Skip to content

Commit d29b02a

Browse files
ThorkjThorkj
authored
[PWGCF] added two track cut and MFT switch (#10499)
1 parent b5dad3c commit d29b02a

File tree

1 file changed

+168
-28
lines changed

1 file changed

+168
-28
lines changed

‎PWGCF/TwoParticleCorrelations/Tasks/corrSparse.cxx‎

Lines changed: 168 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -13,27 +13,33 @@
1313
/// \brief Provides a sparse with usefull two particle correlation info
1414
/// \author Thor Jensen (thor.kjaersgaard.jensen@cern.ch) and Debojit Sarkar (debojit.sarkar@cern.ch)
1515

16+
#include <CCDB/BasicCCDBManager.h>
17+
#include "TRandom3.h"
1618
#include <vector>
1719

1820
#include "Framework/runDataProcessing.h"
1921
#include "Framework/AnalysisTask.h"
2022
#include "Framework/AnalysisDataModel.h"
2123
#include "Framework/ASoAHelpers.h"
22-
#include "Framework/ASoA.h"
24+
#include "Framework/StepTHn.h"
2325
#include "Framework/HistogramRegistry.h"
2426
#include "Framework/RunningWorkflowInfo.h"
2527
#include "CommonConstants/MathConstants.h"
26-
#include "CCDB/BasicCCDBManager.h"
2728
#include "Common/Core/RecoDecay.h"
2829

2930
#include "Common/DataModel/EventSelection.h"
3031
#include "Common/DataModel/TrackSelectionTables.h"
3132
#include "Common/DataModel/Centrality.h"
32-
#include "Common/DataModel/Multiplicity.h"
3333
#include "PWGCF/DataModel/CorrelationsDerived.h"
34+
#include "Common/DataModel/CollisionAssociationTables.h"
3435
#include "PWGCF/Core/CorrelationContainer.h"
3536
#include "PWGCF/Core/PairCuts.h"
37+
#include "DataFormatsParameters/GRPObject.h"
38+
#include "DataFormatsParameters/GRPMagField.h"
3639

40+
using namespace o2;
41+
using namespace o2::framework;
42+
using namespace o2::framework::expressions;
3743
namespace o2::aod
3844
{
3945
namespace corrsparse
@@ -45,18 +51,14 @@ DECLARE_SOA_TABLE(Multiplicity, "AOD", "MULTIPLICITY",
4551

4652
} // namespace o2::aod
4753

48-
using namespace o2;
49-
using namespace o2::framework;
50-
using namespace o2::framework::expressions;
51-
5254
// define the filtered collisions and tracks
5355
#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
5456

5557
struct CalcNch {
5658
O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range")
5759
O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT")
5860
O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT")
59-
O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut")
61+
O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 10.0f, "Eta cut")
6062
O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix")
6163

6264
Filter trackFilter = (nabs(aod::track::eta) < cfgEtaCut) && (aod::track::pt > cfgPtCutMin) && (aod::track::pt < cfgPtCutMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true));
@@ -79,32 +81,45 @@ struct CalcNch {
7981

8082
void process(AodCollisions::iterator const& collision, AodTracks const& tracks)
8183
{
84+
// LOGF(info, "multiplicity from tracks %d, multiplicity from collision %d", tracks.size(), collision.numContrib());
8285
multiplicityNch(tracks.size());
8386
registry.fill(HIST("Ncharge"), tracks.size());
8487
registry.fill(HIST("zVtx_all"), collision.posZ());
8588
}
8689
};
8790

8891
struct CorrSparse {
92+
Service<ccdb::BasicCCDBManager> ccdb;
93+
8994
O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range")
9095
O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT")
9196
O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT")
9297
O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut")
9398
O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix")
9499
O2_DEFINE_CONFIGURABLE(cfgMinMult, int, 0, "Minimum multiplicity for collision")
95100
O2_DEFINE_CONFIGURABLE(cfgMaxMult, int, 10, "Maximum multiplicity for collision")
101+
O2_DEFINE_CONFIGURABLE(cfgMergingCut, float, 0.0, "Merging cut on track merge")
102+
O2_DEFINE_CONFIGURABLE(cfgRadiusLow, float, 0.8, "Low radius for merging cut")
103+
O2_DEFINE_CONFIGURABLE(cfgRadiusHigh, float, 2.5, "High radius for merging cut")
104+
O2_DEFINE_CONFIGURABLE(etaMftTrackMin, float, -5.0, "Minimum eta for MFT track")
105+
O2_DEFINE_CONFIGURABLE(etaMftTrackMax, float, 0.0, "Maximum eta for MFT track")
106+
O2_DEFINE_CONFIGURABLE(nClustersMftTrack, int, 5, "Minimum number of clusters for MFT track")
107+
O2_DEFINE_CONFIGURABLE(cfgSampleSize, double, 10, "Sample size for mixed event")
108+
109+
Configurable<bool> processMFT{"processMFT", true, "Associate particle from MFT"};
96110

97111
ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"};
98112
ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"};
99113
ConfigurableAxis axisPhi{"axisPhi", {72, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
100114
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt axis for histograms"};
101115
ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"};
102-
ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -2, 2}, "delta eta axis for histograms"};
116+
ConfigurableAxis axisDeltaEta{"axisDeltaEta", {48, -2.4, 2.4}, "delta eta axis for histograms"};
103117
ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"};
104118
ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt associated axis for histograms"};
105119
ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for histograms"};
106120
ConfigurableAxis vtxMix{"vtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"};
107121
ConfigurableAxis multMix{"multMix", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for mixed event histograms"};
122+
ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"};
108123

109124
ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"};
110125
ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"};
@@ -114,18 +129,26 @@ struct CorrSparse {
114129
Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut) && (aod::corrsparse::multiplicity) > cfgMinMult && (aod::corrsparse::multiplicity) < cfgMaxMult && (aod::evsel::sel8) == true;
115130
Filter trackFilter = (nabs(aod::track::eta) < cfgEtaCut) && (aod::track::pt > cfgPtCutMin) && (aod::track::pt < cfgPtCutMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true));
116131

117-
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSel, o2::aod::Multiplicity>>; // aod::CentFT0Cs
118-
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra>>;
119-
120132
// Define the outputs
121133
OutputObj<CorrelationContainer> same{Form("sameEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult))};
122134
OutputObj<CorrelationContainer> mixed{Form("mixedEvent_%i_%i", static_cast<int>(cfgMinMult), static_cast<int>(cfgMaxMult))};
123135

124136
HistogramRegistry registry{"registry"};
125137

138+
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSel, o2::aod::Multiplicity>>; // aod::CentFT0Cs
139+
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra>>;
140+
126141
void init(InitContext&)
127142
{
143+
ccdb->setURL("http://alice-ccdb.cern.ch");
144+
ccdb->setCaching(true);
145+
ccdb->setLocalObjectValidityChecking();
146+
147+
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
148+
ccdb->setCreatedNotAfter(now);
149+
128150
LOGF(info, "Starting init");
151+
129152
// Make histograms to check the distributions after cuts
130153
registry.add("deltaEta_deltaPhi_same", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); // check to see the delta eta and delta phi distribution
131154
registry.add("deltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}});
@@ -139,7 +162,7 @@ struct CorrSparse {
139162

140163
registry.add("eventcount", "bin", {HistType::kTH1F, {{3, 0, 3, "bin"}}}); // histogram to see how many events are in the same and mixed event
141164

142-
std::vector<AxisSpec> corrAxis = {{axisMultiplicity, "Nch"},
165+
std::vector<AxisSpec> corrAxis = {{axisSample, "Sample"},
143166
{axisVertex, "z-vtx (cm)"},
144167
{axisPtTrigger, "p_{T} (GeV/c)"},
145168
{axisPtAssoc, "p_{T} (GeV/c)"},
@@ -160,6 +183,37 @@ struct CorrSparse {
160183
MixedEvent = 3
161184
};
162185

186+
template <typename TTrackAssoc>
187+
bool isAcceptedMftTrack(TTrackAssoc const& mftTrack)
188+
{
189+
// cut on the eta of MFT tracks
190+
if (mftTrack.eta() > etaMftTrackMax || mftTrack.eta() < etaMftTrackMin) {
191+
return false;
192+
}
193+
194+
// cut on the number of clusters of the reconstructed MFT track
195+
if (mftTrack.nClusters() < nClustersMftTrack) {
196+
return false;
197+
}
198+
199+
return true;
200+
}
201+
202+
int getMagneticField(uint64_t timestamp)
203+
{
204+
// Get the magnetic field
205+
static o2::parameters::GRPMagField* grpo = nullptr;
206+
if (grpo == nullptr) {
207+
grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>("/GLO/Config/GRPMagField", timestamp);
208+
if (grpo == nullptr) {
209+
LOGF(fatal, "GRP object not found for timestamp %llu", timestamp);
210+
return 0;
211+
}
212+
LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field());
213+
}
214+
return grpo->getNominalL3Field();
215+
}
216+
163217
// fill multiple histograms
164218
template <typename TCollision, typename TTracks>
165219
void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms.
@@ -174,9 +228,38 @@ struct CorrSparse {
174228
}
175229
}
176230

177-
template <CorrelationContainer::CFStep step, typename TTracks>
178-
void fillCorrelations(TTracks tracks1, TTracks tracks2, float posZ, int system, float Nch) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms
231+
template <typename TTrack, typename TTrackAssoc>
232+
float getDPhiStar(TTrack const& track1, TTrackAssoc const& track2, float radius, int magField)
233+
{
234+
float charge1 = track1.sign();
235+
float charge2 = track2.sign();
236+
237+
float phi1 = track1.phi();
238+
float phi2 = track2.phi();
239+
240+
float pt1 = track1.pt();
241+
float pt2 = track2.pt();
242+
243+
int fbSign = (magField > 0) ? 1 : -1;
244+
245+
float dPhiStar = phi1 - phi2 - charge1 * fbSign * std::asin(0.075 * radius / pt1) + charge2 * fbSign * std::asin(0.075 * radius / pt2);
246+
247+
if (dPhiStar > constants::math::PI)
248+
dPhiStar = constants::math::TwoPI - dPhiStar;
249+
if (dPhiStar < -constants::math::PI)
250+
dPhiStar = -constants::math::TwoPI - dPhiStar;
251+
252+
return dPhiStar;
253+
}
254+
255+
//
256+
template <CorrelationContainer::CFStep step, typename TTracks, typename TTracksAssoc>
257+
void fillCorrelations(TTracks tracks1, TTracksAssoc tracks2, float posZ, int system, float Nch, int magneticField) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms
179258
{
259+
TRandom3* gRandom = new TRandom3();
260+
261+
int fSampleIndex = gRandom->Uniform(0, cfgSampleSize);
262+
180263
// loop over all tracks
181264
for (auto const& track1 : tracks1) {
182265

@@ -187,47 +270,104 @@ struct CorrSparse {
187270
for (auto const& track2 : tracks2) {
188271

189272
if (track1.pt() <= track2.pt())
190-
continue; // skip if the trigger pt is less than the associate p
273+
continue; // skip if the trigger pt is less than the associate pt
274+
275+
if (processMFT) {
276+
if constexpr (std::is_same_v<aod::MFTTracks, TTracksAssoc>) {
277+
if (!isAcceptedMftTrack(track2)) {
278+
continue;
279+
}
280+
}
281+
}
191282

192283
float deltaPhi = RecoDecay::constrainAngle(track1.phi() - track2.phi(), -PIHalf);
193284
float deltaEta = track1.eta() - track2.eta();
194285

286+
if (std::abs(deltaEta) < cfgMergingCut) {
287+
288+
double dPhiStarHigh = getDPhiStar(track1, track2, cfgRadiusHigh, magneticField);
289+
double dPhiStarLow = getDPhiStar(track1, track2, cfgRadiusLow, magneticField);
290+
291+
const double kLimit = 3.0 * cfgMergingCut;
292+
293+
bool bIsBelow = kFALSE;
294+
295+
if (std::abs(dPhiStarLow) < kLimit || std::abs(dPhiStarHigh) < kLimit || dPhiStarLow * dPhiStarHigh < 0) {
296+
for (double rad(cfgRadiusLow); rad < cfgRadiusHigh; rad += 0.01) {
297+
double dPhiStar = getDPhiStar(track1, track2, rad, magneticField);
298+
if (std::abs(dPhiStar) < kLimit) {
299+
bIsBelow = kTRUE;
300+
break;
301+
}
302+
}
303+
if (bIsBelow)
304+
continue;
305+
}
306+
}
307+
195308
// fill the right sparse and histograms
196309
if (system == SameEvent) {
197-
same->getPairHist()->Fill(step, Nch, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
310+
311+
same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
198312
registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta);
199313
} else if (system == MixedEvent) {
200-
mixed->getPairHist()->Fill(step, Nch, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
314+
315+
mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta);
201316
registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta);
202317
}
203318
}
204319
}
205320
}
206321

207-
void processSame(AodCollisions::iterator const& collision, AodTracks const& tracks)
322+
void processSame(AodCollisions::iterator const& collision, AodTracks const& tracks, aod::MFTTracks const& mfts, aod::BCsWithTimestamps const&)
208323
{
209324

325+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
326+
210327
registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin
211328
fillYield(collision, tracks);
212-
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks, tracks, collision.posZ(), SameEvent, tracks.size()); // fill the SE histogram and Sparse
329+
330+
if (processMFT) {
331+
332+
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks, mfts, collision.posZ(), SameEvent, tracks.size(), getMagneticField(bc.timestamp()));
333+
334+
} else {
335+
336+
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks, tracks, collision.posZ(), SameEvent, tracks.size(), getMagneticField(bc.timestamp()));
337+
}
213338
}
214339
PROCESS_SWITCH(CorrSparse, processSame, "Process same event", true);
215340

216341
// event mixing
217-
218342
SliceCache cache;
219343
using MixedBinning = ColumnBinningPolicy<aod::collision::PosZ, aod::corrsparse::Multiplicity>;
220344

221345
// the process for filling the mixed events
222-
void processMixed(AodCollisions const& collisions, AodTracks const& tracks)
346+
void processMixed(AodCollisions const& collisions, AodTracks const& tracks, aod::MFTTracks const& MFTtracks, aod::BCsWithTimestamps const&)
223347
{
224-
MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default)
225-
auto tracksTuple = std::make_tuple(tracks);
226-
SameKindPair<AodCollisions, AodTracks, MixedBinning> pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip
227348

228-
for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) {
229-
registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin
230-
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks1, tracks2, collision1.posZ(), MixedEvent, tracks1.size());
349+
if (processMFT) {
350+
MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default)
351+
auto tracksTuple = std::make_tuple(tracks, MFTtracks);
352+
SameKindPair<AodCollisions, AodTracks, MixedBinning> pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip
353+
354+
for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) {
355+
registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin
356+
auto bc = collision1.bc_as<aod::BCsWithTimestamps>();
357+
358+
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks1, tracks2, collision1.posZ(), MixedEvent, tracks1.size(), getMagneticField(bc.timestamp()));
359+
}
360+
} else {
361+
MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default)
362+
auto tracksTuple = std::make_tuple(tracks, tracks);
363+
SameKindPair<AodCollisions, AodTracks, MixedBinning> pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip
364+
365+
for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) {
366+
registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin
367+
auto bc = collision1.bc_as<aod::BCsWithTimestamps>();
368+
369+
fillCorrelations<CorrelationContainer::kCFStepReconstructed>(tracks1, tracks2, collision1.posZ(), MixedEvent, tracks1.size(), getMagneticField(bc.timestamp()));
370+
}
231371
}
232372
}
233373
PROCESS_SWITCH(CorrSparse, processMixed, "Process mixed events", true);

0 commit comments

Comments
 (0)