Skip to content

Commit b837e96

Browse files
joachimckhjoachimckh
committed
light ion flow event cuts (def OO pass1)
1 parent c585c36 commit b837e96

File tree

1 file changed

+216
-6
lines changed

1 file changed

+216
-6
lines changed

‎Common/TableProducer/eseTableProducer.cxx‎

Lines changed: 216 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
#include <Framework/runDataProcessing.h>
3838

3939
#include <TMath.h>
40+
#include <TF1.h>
4041

4142
#include <chrono>
4243
#include <cstddef>
@@ -110,15 +111,68 @@ struct EseTableProducer {
110111
bool correctionsLoaded = false;
111112
} cfg;
112113

114+
TF1* fMultPVCutLow = nullptr;
115+
TF1* fMultPVCutHigh = nullptr;
116+
TF1* fMultCutLow = nullptr;
117+
TF1* fMultCutHigh = nullptr;
118+
TF1* fMultPVGlobalCutHigh = nullptr;
119+
TF1* fMultGlobalV0ACutLow = nullptr;
120+
TF1* fMultGlobalV0ACutHigh = nullptr;
121+
TF1* fMultGlobalT0ACutLow = nullptr;
122+
TF1* fMultGlobalT0ACutHigh = nullptr;
123+
124+
TF1* fPtDepDCAxy = nullptr;
125+
126+
float vtxZlow = -10.0, vtxZup = 10.0;
127+
int nchbins = 300;
128+
float nchlow = 0;
129+
float nchup = 3000;
130+
std::vector<double> multGlobalCorrCutPars;
131+
std::vector<double> multPVCorrCutPars;
132+
std::vector<double> multGlobalPVCorrCutPars;
133+
std::vector<double> multGlobalV0ACutPars;
134+
std::vector<double> multGlobalT0ACutPars;
135+
113136
Configurable<float> cfgVtxZ{"cfgVtxZ", 10.0f, "max z vertex position"};
114137
Configurable<float> cfgEta{"cfgEta", 0.8f, "max eta"};
115138
Configurable<float> cfgPtmin{"cfgPtmin", 0.2f, "min pt"};
116-
Configurable<float> cfgPtmax{"cfgPtmax", 5.0f, "max pt"};
139+
Configurable<float> cfgPtmax{"cfgPtmax", 3.0f, "max pt"};
117140
Configurable<float> cfgChi2PrITSCls{"cfgChi2PrITSCls", 4.0f, "max chi2 per ITS cluster"};
118141
Configurable<float> cfgChi2PrTPCCls{"cfgChi2PrTPCCls", 2.5f, "max chi2 per TPC cluster"};
119142
Configurable<float> cfgDCAz{"cfgDCAz", 2.0f, "max DCAz cut"};
143+
Configurable<bool> cfgDoOccupancySel{"cfgDoOccupancySel", true, "enable occupancy selection"};
144+
Configurable<int> cfgOccupancySelection{"cfgOccupancySelection", 2000, "Max occupancy selection, -999 to disable"};
145+
Configurable<bool> cfgUseAdditionalEventCut{"cfgUseAdditionalEventCut", false, "Use additional event cut on mult correlations"};
146+
Configurable<bool> cfgMultCut{"cfgMultCut", true, "Use additional event cut on mult correlations"};
147+
Configurable<bool> cfgTVXinTRD{"cfgTVXinTRD", true, "Use kTVXinTRD (reject TRD triggered events)"};
148+
Configurable<bool> cfgNoSameBunchPileupCut{"cfgNoSameBunchPileupCut", true, "kNoSameBunchPileupCut"};
149+
Configurable<bool> cfgIsGoodZvtxFT0vsPV{"cfgIsGoodZvtxFT0vsPV", true, "kIsGoodZvtxFT0vsPV"};
150+
Configurable<bool> cfgNoCollInTimeRangeStandard{"cfgNoCollInTimeRangeStandard", true, "kNoCollInTimeRangeStandard"};
151+
Configurable<bool> cfgIsVertexITSTPC{"cfgIsVertexITSTPC", true, "Selects collisions with at least one ITS-TPC track"};
152+
Configurable<bool> cfgIsGoodITSLayersAll{"cfgIsGoodITSLayersAll", true, "kIsGoodITSLayersAll"};
120153
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency object"};
121154

155+
// Cut Configurables
156+
Configurable<std::vector<double>> cfgMultGlobalCutPars{"cfgMultGlobalCutPars", std::vector<double>{2272.16, -76.6932, 1.01204, -0.00631545, 1.59868e-05, 136.336, -4.97006, 0.121199, -0.0015921, 7.66197e-06}, "Global vs FT0C multiplicity cut parameter values"};
157+
Configurable<std::vector<double>> cfgMultPVCutPars{"cfgMultPVCutPars", std::vector<double>{3074.43, -106.192, 1.46176, -0.00968364, 2.61923e-05, 182.128, -7.43492, 0.193901, -0.00256715, 1.22594e-05}, "PV vs FT0C multiplicity cut parameter values"};
158+
Configurable<std::vector<double>> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector<double>{-0.223013, 0.715849, 0.664242, 0.0829653, -0.000503733, 1.21185e-06}, "Global vs PV multiplicity cut parameter values"};
159+
Configurable<std::string> cfgMultCorrHighCutFunction{"cfgMultCorrHighCutFunction", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"};
160+
Configurable<std::string> cfgMultCorrLowCutFunction{"cfgMultCorrLowCutFunction", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"};
161+
Configurable<std::string> cfgMultGlobalPVCorrCutFunction{"cfgMultGlobalPVCorrCutFunction", "[0] + [1]*x + 3*([2] + [3]*x + [4]*x*x + [5]*x*x*x)", "Functional for global vs pv multiplicity correlation cut"};
162+
struct : ConfigurableGroup {
163+
Configurable<std::string> cfgMultGlobalASideCorrCutFunction{"cfgMultGlobalASideCorrCutFunction","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x + [10]*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", "Functional for global vs V0A multiplicity low correlation cut"};
164+
Configurable<std::vector<double>> cfgMultGlobalV0ACutPars{"cfgMultGlobalV0ACutPars", std::vector<double>{567.785, 172.715, 0.77888, -0.00693466, 1.40564e-05, 679.853, 66.8068, -0.444332, 0.00115002, -4.92064e-07}, "Global vs FV0A multiplicity cut parameter values"};
165+
Configurable<std::vector<double>> cfgMultGlobalT0ACutPars{"cfgMultGlobalT0ACutPars", std::vector<double>{241.618, 61.8402, 0.348049, -0.00306078, 6.20357e-06, 315.235, 29.1491, -0.188639, 0.00044528, -9.08912e-08}, "Global vs FT0A multiplicity cut parameter values"};
166+
Configurable<float> cfgGlobalV0ALowSigma{"cfgGlobalV0ALowSigma", -3.0f, "Number of sigma deviations below expected value in global vs V0A correlation"};
167+
Configurable<float> cfgGlobalV0AHighSigma{"cfgGlobalV0AHighSigma", 4.0f, "Number of sigma deviations above expected value in global vs V0A correlation"};
168+
Configurable<float> cfgGlobalT0ALowSigma{"cfgGlobalT0ALowSigma", -3.0f, "Number of sigma deviations below expected value in global vs T0A correlation"};
169+
Configurable<float> cfgGlobalT0AHighSigma{"cfgGlobalT0AHighSigma", 4.0f, "Number of sigma deviations above expected value in global vs T0A correlation"};
170+
} cfgGlobalAsideCorrCuts;
171+
Configurable<std::string> cfgDCAxy{"cfgDCAxy", "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"};
172+
Configurable<float> cfgDCAxyNSigma{"cfgDCAxyNSigma", 7.0f, "Cut on number of sigma deviations from expected DCA in the transverse direction"};
173+
174+
//
175+
122176
// o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ;
123177
o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == static_cast<uint8_t>(true))) && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz;
124178

@@ -127,8 +181,24 @@ struct EseTableProducer {
127181

128182
void init(o2::framework::InitContext&)
129183
{
130-
131184
LOGF(info, "ESETable::init()");
185+
EseTableProducer::vtxZlow = -cfgVtxZ;
186+
EseTableProducer::vtxZup = cfgVtxZ;
187+
EseTableProducer::multPVCorrCutPars = cfgMultPVCutPars;
188+
EseTableProducer::multGlobalCorrCutPars = cfgMultGlobalCutPars;
189+
EseTableProducer::multGlobalV0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalV0ACutPars;
190+
EseTableProducer::multGlobalT0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalT0ACutPars;
191+
EseTableProducer::multGlobalPVCorrCutPars = cfgMultGlobalPVCutPars;
192+
193+
194+
AxisSpec t0cAxis = {1000, 0, 10000, "N_{ch} (T0C)"};
195+
AxisSpec t0aAxis = {300, 0, 30000, "N_{ch} (T0A)"};
196+
AxisSpec v0aAxis = {800, 0, 80000, "N_{ch} (V0A)"};
197+
std::vector<double> nchbinning;
198+
int nchskip = (nchup - nchlow) / nchbins;
199+
for (int i = 0; i <= nchbins; ++i) {
200+
nchbinning.push_back(nchskip * i + nchlow + 0.5);
201+
}
132202

133203
registry.add("hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}});
134204
registry.add("hESEstat", "ese status;ese status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}});
@@ -151,6 +221,52 @@ struct EseTableProducer {
151221
weightsFFit->setResolution(cfgnResolution);
152222
weightsFFit->setQnType(veccfg);
153223
weightsFFit->init();
224+
225+
fPtDepDCAxy = new TF1("ptDepDCAxy", Form("[0]*%s", cfgDCAxy->c_str()), 0.001, 100);
226+
fPtDepDCAxy->SetParameter(0, cfgDCAxyNSigma);
227+
LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", cfgDCAxy->c_str()));
228+
if (cfgUseAdditionalEventCut) {
229+
fMultPVCutLow = new TF1("fMultPVCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100);
230+
fMultPVCutLow->SetParameters(&(EseTableProducer::multPVCorrCutPars[0]));
231+
fMultPVCutHigh = new TF1("fMultPVCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100);
232+
fMultPVCutHigh->SetParameters(&(EseTableProducer::multPVCorrCutPars[0]));
233+
fMultCutLow = new TF1("fMultCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100);
234+
fMultCutLow->SetParameters(&(EseTableProducer::multGlobalCorrCutPars[0]));
235+
fMultCutHigh = new TF1("fMultCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100);
236+
fMultCutHigh->SetParameters(&(EseTableProducer::multGlobalCorrCutPars[0]));
237+
fMultPVGlobalCutHigh = new TF1("fMultPVGlobalCutHigh", cfgMultGlobalPVCorrCutFunction->c_str(), 0, nchbinning.back());
238+
fMultPVGlobalCutHigh->SetParameters(&(EseTableProducer::multGlobalPVCorrCutPars[0]));
239+
240+
LOGF(info, "Global V0A function: %s in range 0-%g", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), v0aAxis.binEdges.back());
241+
fMultGlobalV0ACutLow = new TF1("fMultGlobalV0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back());
242+
for (std::size_t i = 0; i < EseTableProducer::multGlobalV0ACutPars.size(); ++i)
243+
fMultGlobalV0ACutLow->SetParameter(i, EseTableProducer::multGlobalV0ACutPars[i]);
244+
fMultGlobalV0ACutLow->SetParameter(EseTableProducer::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0ALowSigma);
245+
for (int i = 0; i < fMultGlobalV0ACutLow->GetNpar(); ++i)
246+
LOGF(info, "fMultGlobalV0ACutLow par %d = %g", i, fMultGlobalV0ACutLow->GetParameter(i));
247+
248+
fMultGlobalV0ACutHigh = new TF1("fMultGlobalV0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back());
249+
for (std::size_t i = 0; i < EseTableProducer::multGlobalV0ACutPars.size(); ++i)
250+
fMultGlobalV0ACutHigh->SetParameter(i, EseTableProducer::multGlobalV0ACutPars[i]);
251+
fMultGlobalV0ACutHigh->SetParameter(EseTableProducer::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0AHighSigma);
252+
for (int i = 0; i < fMultGlobalV0ACutHigh->GetNpar(); ++i)
253+
LOGF(info, "fMultGlobalV0ACutHigh par %d = %g", i, fMultGlobalV0ACutHigh->GetParameter(i));
254+
255+
LOGF(info, "Global T0A function: %s", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str());
256+
fMultGlobalT0ACutLow = new TF1("fMultGlobalT0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back());
257+
for (std::size_t i = 0; i < EseTableProducer::multGlobalT0ACutPars.size(); ++i)
258+
fMultGlobalT0ACutLow->SetParameter(i, EseTableProducer::multGlobalT0ACutPars[i]);
259+
fMultGlobalT0ACutLow->SetParameter(EseTableProducer::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0ALowSigma);
260+
for (int i = 0; i < fMultGlobalT0ACutLow->GetNpar(); ++i)
261+
LOGF(info, "fMultGlobalT0ACutLow par %d = %g", i, fMultGlobalT0ACutLow->GetParameter(i));
262+
263+
fMultGlobalT0ACutHigh = new TF1("fMultGlobalT0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back());
264+
for (std::size_t i = 0; i < EseTableProducer::multGlobalT0ACutPars.size(); ++i)
265+
fMultGlobalT0ACutHigh->SetParameter(i, EseTableProducer::multGlobalT0ACutPars[i]);
266+
fMultGlobalT0ACutHigh->SetParameter(EseTableProducer::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0AHighSigma);
267+
for (int i = 0; i < fMultGlobalT0ACutHigh->GetNpar(); ++i)
268+
LOGF(info, "fMultGlobalT0ACutHigh par %d = %g", i, fMultGlobalT0ACutHigh->GetParameter(i));
269+
}
154270
}
155271

156272
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
@@ -329,17 +445,112 @@ struct EseTableProducer {
329445
}
330446
PROCESS_SWITCH(EseTableProducer, processESE, "process q vectors to calculate reduced q-vector", true);
331447

332-
void processMeanPt(soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
448+
struct XAxis {
449+
float centrality;
450+
int64_t multiplicity;
451+
};
452+
453+
template <typename TCollision>
454+
bool eventSelected(TCollision collision, const int& multTrk, const float& centrality)
333455
{
456+
if (cfgTVXinTRD) {
457+
if (collision.alias_bit(kTVXinTRD)) {
458+
return 0;
459+
}
460+
}
461+
if (cfgNoSameBunchPileupCut) {
462+
if (!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
463+
return 0;
464+
}
465+
}
466+
if (cfgIsGoodZvtxFT0vsPV) {
467+
if (!collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) {
468+
return 0;
469+
}
470+
}
471+
if (cfgNoCollInTimeRangeStandard) {
472+
if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) {
473+
return 0;
474+
}
475+
}
476+
477+
if (cfgIsVertexITSTPC) {
478+
if (!collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) {
479+
return 0;
480+
}
481+
}
482+
483+
if (cfgIsGoodITSLayersAll) {
484+
if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
485+
return 0;
486+
}
487+
}
334488

335-
std::vector<float> meanPt{-1};
489+
float vtxz = -999;
490+
if (collision.numContrib() > 1) {
491+
vtxz = collision.posZ();
492+
float zRes = std::sqrt(collision.covZZ());
493+
float minZRes = 0.25;
494+
int minNContrib = 20;
495+
if (zRes > minZRes && collision.numContrib() < minNContrib)
496+
vtxz = -999;
497+
}
498+
auto multNTracksPV = collision.multNTracksPV();
499+
500+
if (vtxz > EseTableProducer::vtxZup || vtxz < EseTableProducer::vtxZlow)
501+
return 0;
502+
503+
if (cfgMultCut) {
504+
if (multNTracksPV < fMultPVCutLow->Eval(centrality))
505+
return 0;
506+
if (multNTracksPV > fMultPVCutHigh->Eval(centrality))
507+
return 0;
508+
if (multTrk < fMultCutLow->Eval(centrality))
509+
return 0;
510+
if (multTrk > fMultCutHigh->Eval(centrality))
511+
return 0;
512+
if (multTrk > fMultPVGlobalCutHigh->Eval(collision.multNTracksPV()))
513+
return 0;
514+
515+
if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast<double>(collision.multFV0A()) < fMultGlobalV0ACutLow->Eval(multTrk))
516+
return 0;
517+
if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast<double>(collision.multFV0A()) > fMultGlobalV0ACutHigh->Eval(multTrk))
518+
return 0;
519+
if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast<double>(collision.multFT0A()) < fMultGlobalT0ACutLow->Eval(multTrk))
520+
return 0;
521+
if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast<double>(collision.multFT0A()) > fMultGlobalT0ACutHigh->Eval(multTrk))
522+
return 0;
523+
}
524+
return 1;
525+
}
526+
527+
void processMeanPt(soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
528+
{
529+
std::vector<float> meanPt{-999};
336530
std::vector<float> meanPtShape{-1};
337531
if (collision.posZ() < -cfgVtxZ || collision.posZ() > cfgVtxZ) {
338532
meanPts(meanPt);
339533
meanPtShapes(meanPtShape);
340534
return;
341535
}
342-
536+
if (!collision.sel8()) {
537+
meanPts(meanPt);
538+
meanPtShapes(meanPtShape);
539+
return;
540+
}
541+
if (cfgDoOccupancySel) {
542+
int occupancy = collision.trackOccupancyInTimeRange();
543+
if (occupancy < 0 || occupancy > cfgOccupancySelection) {
544+
meanPts(meanPt);
545+
meanPtShapes(meanPtShape);
546+
return;
547+
}
548+
}
549+
550+
const XAxis xaxis{collision.centFT0C(), tracks.size()};
551+
if (cfgUseAdditionalEventCut && !eventSelected(collision, xaxis.multiplicity, xaxis.centrality))
552+
return;
553+
343554
registry.fill(HIST("hMeanPtStat"), 0.5);
344555
const auto centrality = collision.centFT0C();
345556
const auto mean = calculateMeanPt(tracks, centrality);
@@ -363,7 +574,6 @@ struct EseTableProducer {
363574
}
364575
}
365576
}
366-
367577
meanPts(meanPt);
368578
meanPtShapes(meanPtShape);
369579
}

0 commit comments

Comments
 (0)