Skip to content

Commit 18a3aa3

Browse files
cnkosteralibuild
cnkoster
and
alibuild
authored
[PWGCF] FlowSP: add rct flags for zdc and shift correction (#12956)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent e0481d5 commit 18a3aa3

File tree

2 files changed

+544
-372
lines changed

2 files changed

+544
-372
lines changed

‎PWGCF/Flow/TableProducer/zdcQVectors.cxx‎

Lines changed: 136 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@ using namespace o2::framework;
6363
using namespace o2::framework::expressions;
6464
using namespace o2::aod::track;
6565
using namespace o2::aod::evsel;
66+
using namespace o2::aod::rctsel;
6667

6768
namespace o2::analysis::qvectortask
6869
{
@@ -103,6 +104,15 @@ struct ZdcQVectors {
103104

104105
Produces<aod::SPTableZDC> spTableZDC;
105106

107+
struct : ConfigurableGroup {
108+
Configurable<bool> cfgEvtUseRCTFlagChecker{"cfgEvtUseRCTFlagChecker", false, "Evt sel: use RCT flag checker"};
109+
Configurable<std::string> cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label (CBT, CBT_hadronPID)"}; // all Labels can be found in Common/CCDB/RCTSelectionFlags.h
110+
Configurable<bool> cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"};
111+
Configurable<bool> cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", false, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"};
112+
} rctFlags;
113+
114+
RCTFlagsChecker rctChecker;
115+
106116
ConfigurableAxis axisCent{"axisCent", {90, 0, 90}, "Centrality axis in 1% bins"};
107117
ConfigurableAxis axisCent10{"axisCent10", {9, 0, 90}, "Centrality axis in 10% bins"};
108118
ConfigurableAxis axisQ{"axisQ", {100, -2, 2}, "Q vector (xy) in ZDC"};
@@ -121,10 +131,10 @@ struct ZdcQVectors {
121131

122132
O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10.0f, "Accepted z-vertex range")
123133
O2_DEFINE_CONFIGURABLE(cfgMagField, float, 99999, "Configurable magnetic field; default CCDB will be queried")
124-
O2_DEFINE_CONFIGURABLE(cfgEnergyCal, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4/Energy", "ccdb path for energy calibration histos")
125-
O2_DEFINE_CONFIGURABLE(cfgMeanv, std::string, "Users/c/ckoster/ZDC/LHC23_zzh_pass4/vmean", "ccdb path for mean v histos")
134+
O2_DEFINE_CONFIGURABLE(cfgEnergyCal, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Energy", "ccdb path for energy calibration histos")
135+
O2_DEFINE_CONFIGURABLE(cfgMeanv, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/vmean", "ccdb path for mean v histos")
126136
O2_DEFINE_CONFIGURABLE(cfgMinEntriesSparseBin, int, 100, "Minimal number of entries allowed in 4D recentering histogram to use for recentering.")
127-
O2_DEFINE_CONFIGURABLE(cfgRec, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass4", "ccdb path for recentering histos");
137+
O2_DEFINE_CONFIGURABLE(cfgRec, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5", "ccdb path for recentering histos");
128138
O2_DEFINE_CONFIGURABLE(cfgFillCommonRegistry, bool, true, "Fill common registry with histograms");
129139

130140
// Additional event selections
@@ -138,13 +148,17 @@ struct ZdcQVectors {
138148
O2_DEFINE_CONFIGURABLE(cfgEvSelsCentMin, float, 0, "Minimum cenrality for selected events");
139149
O2_DEFINE_CONFIGURABLE(cfgEvSelsCentMax, float, 90, "Maximum cenrality for selected events");
140150

151+
O2_DEFINE_CONFIGURABLE(cfgUseShift, bool, false, "Use shift for PsiA and PsiC ZDC");
152+
O2_DEFINE_CONFIGURABLE(cfgCCDBdir_Shift, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Shift", "CCDB directory for Shift ZDC");
153+
141154
// define my.....
142155
// Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ;
143156
using UsedCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNGlobals>;
144157
using BCsRun3 = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels, aod::Run3MatchedToBCSparse>;
145158

146159
enum SelectionCriteria {
147160
evSel_FilteredEvent,
161+
evSel_RCTFlagsZDC,
148162
evSel_Zvtx,
149163
evSel_sel8,
150164
evSel_occupancy,
@@ -175,6 +189,12 @@ struct ZdcQVectors {
175189
std::vector<bool> calibfilesLoaded = std::vector<bool>(3, false);
176190
int atStep = 0;
177191
int atIteration = 0;
192+
193+
TProfile3D* shiftprofileC = nullptr;
194+
TProfile3D* shiftprofileA = nullptr;
195+
196+
bool isShiftProfileFound = false;
197+
178198
} cal;
179199

180200
enum FillType {
@@ -191,6 +211,8 @@ struct ZdcQVectors {
191211
int64_t now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
192212
ccdb->setCreatedNotAfter(now);
193213

214+
rctChecker.init(rctFlags.cfgEvtRCTFlagCheckerLabel, rctFlags.cfgEvtRCTFlagCheckerZDCCheck, rctFlags.cfgEvtRCTFlagCheckerLimitAcceptAsBad);
215+
194216
std::vector<const char*> sides = {"A", "C"};
195217
std::vector<const char*> capCOORDS = {"X", "Y"};
196218

@@ -236,6 +258,13 @@ struct ZdcQVectors {
236258
registry.add<TProfile>("QA/ZNA_Energy", "ZNA_Energy", kTProfile, {{8, 0, 8}});
237259
registry.add<TProfile>("QA/ZNC_Energy", "ZNC_Energy", kTProfile, {{8, 0, 8}});
238260

261+
registry.add<TProfile3D>("QA/ShiftZDCC", "ShiftZDCC", kTProfile3D, {{100, 0, 100}, {2, 0, 2}, {10, 0, 10}});
262+
registry.add<TProfile3D>("QA/ShiftZDCA", "ShiftZDCA", kTProfile3D, {{100, 0, 100}, {2, 0, 2}, {10, 0, 10}});
263+
registry.add<TH1>("QA/psiZDCA", "psiZDCA", kTH1D, {{100, -4, 4}});
264+
registry.add<TH1>("QA/psiZDCA_shift", "psiZDCA_shift", kTH1D, {{100, -4, 4}});
265+
registry.add<TH1>("QA/psiZDCC", "psiZDCC", kTH1D, {{100, -4, 4}});
266+
registry.add<TH1>("QA/psiZDCC_shift", "psiZDCC_shift", kTH1D, {{100, -4, 4}});
267+
239268
registry.add<TProfile>("QA/before/ZNA_pmC", "ZNA_pmC", kTProfile, {{1, 0, 1.}});
240269
registry.add<TProfile>("QA/before/ZNA_pm1", "ZNA_pm1", kTProfile, {{1, 0, 1.}});
241270
registry.add<TProfile>("QA/before/ZNA_pm2", "ZNA_pm2", kTProfile, {{1, 0, 1.}});
@@ -295,6 +324,7 @@ struct ZdcQVectors {
295324

296325
registry.add("hEventCount", "Number of Event; Cut; #Events Passed Cut", {HistType::kTH1D, {{nEventSelections, 0, nEventSelections}}});
297326
registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_FilteredEvent + 1, "Filtered events");
327+
registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_RCTFlagsZDC + 1, "RCT Flags ZDC");
298328
registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_Zvtx + 1, "Z vertex cut event");
299329
registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_sel8 + 1, "Sel8");
300330
registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(evSel_occupancy + 1, "kOccupancy");
@@ -560,34 +590,46 @@ struct ZdcQVectors {
560590
if (cfgNGlobal)
561591
cent = collision.centNGlobal();
562592

593+
v[0] = collision.posX();
594+
v[1] = collision.posY();
595+
v[2] = collision.posZ();
596+
centrality = cent;
597+
598+
const auto& foundBC = collision.foundBC_as<BCsRun3>();
599+
runnumber = foundBC.runNumber();
600+
563601
if (cfgFillCommonRegistry)
564602
registry.fill(HIST("QA/centrality_before"), cent);
565603

566604
registry.fill(HIST("hEventCount"), evSel_FilteredEvent);
567605

568-
if (!eventSelected(collision, cent)) {
606+
if (rctFlags.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) {
569607
// event not selected
570608
isSelected = false;
571609
spTableZDC(runnumber, cent, v[0], v[1], v[2], 0, 0, 0, 0, isSelected, 0, 0);
572610
counter++;
611+
lastRunNumber = runnumber;
573612
return;
574613
}
614+
registry.fill(HIST("hEventCount"), evSel_RCTFlagsZDC);
575615

576-
const auto& foundBC = collision.foundBC_as<BCsRun3>();
616+
if (!eventSelected(collision, cent)) {
617+
// event not selected
618+
isSelected = false;
619+
spTableZDC(runnumber, cent, v[0], v[1], v[2], 0, 0, 0, 0, isSelected, 0, 0);
620+
counter++;
621+
lastRunNumber = runnumber;
622+
return;
623+
}
577624

578625
if (!foundBC.has_zdc()) {
579626
isSelected = false;
580627
spTableZDC(runnumber, cent, v[0], v[1], v[2], 0, 0, 0, 0, isSelected, 0, 0);
581628
counter++;
629+
lastRunNumber = runnumber;
582630
return;
583631
}
584632

585-
v[0] = collision.posX();
586-
v[1] = collision.posY();
587-
v[2] = collision.posZ();
588-
centrality = cent;
589-
runnumber = foundBC.runNumber();
590-
591633
// load new calibrations for new runs only
592634
if (runnumber != lastRunNumber) {
593635
cal.calibfilesLoaded[0] = false;
@@ -598,8 +640,6 @@ struct ZdcQVectors {
598640

599641
cal.calibfilesLoaded[2] = false;
600642
cal.calibList[2] = nullptr;
601-
602-
lastRunNumber = runnumber;
603643
}
604644

605645
const auto& zdcCol = foundBC.zdc();
@@ -663,6 +703,7 @@ struct ZdcQVectors {
663703
counter++;
664704
isSelected = false;
665705
spTableZDC(runnumber, centrality, v[0], v[1], v[2], 0, 0, 0, 0, isSelected, 0, 0);
706+
lastRunNumber = runnumber;
666707
return;
667708
}
668709

@@ -673,6 +714,7 @@ struct ZdcQVectors {
673714
counter++;
674715
isSelected = false;
675716
spTableZDC(runnumber, centrality, v[0], v[1], v[2], 0, 0, 0, 0, isSelected, 0, 0);
717+
lastRunNumber = runnumber;
676718
return;
677719
}
678720

@@ -814,6 +856,7 @@ struct ZdcQVectors {
814856

815857
spTableZDC(runnumber, centrality, v[0], v[1], v[2], q[0], q[1], q[2], q[3], isSelected, 0, 0);
816858
counter++;
859+
lastRunNumber = runnumber;
817860
return;
818861
} else {
819862
if (cfgFillCommonRegistry)
@@ -867,14 +910,93 @@ struct ZdcQVectors {
867910
registry.get<TProfile>(HIST("QA/after/ZNC_Qy"))->Fill(Form("%d", runnumber), qRec[3]);
868911
}
869912

870-
spTableZDC(runnumber, centrality, v[0], v[1], v[2], qRec[0], qRec[1], qRec[2], qRec[3], isSelected, cal.atIteration, cal.atStep);
913+
// do shift for psi.
914+
double psiZDCA = 1.0 * std::atan2(qRec[1], qRec[0]);
915+
double psiZDCC = 1.0 * std::atan2(qRec[3], qRec[2]);
916+
917+
int nshift = 10; // no. of iterations
918+
919+
double psiZDCAshift = psiZDCA;
920+
double psiZDCCshift = psiZDCC;
921+
922+
double deltaPsiZDCA = 0;
923+
double deltaPsiZDCC = 0;
924+
925+
if (cfgUseShift && !cfgCCDBdir_Shift.value.empty()) {
926+
if (lastRunNumber != runnumber) {
927+
cal.isShiftProfileFound = false;
928+
LOGF(info, "Getting shift profile from CCDB for runnumber: %d", runnumber);
929+
TList* hcorrList = ccdb->getForTimeStamp<TList>(cfgCCDBdir_Shift.value, foundBC.timestamp());
930+
cal.shiftprofileC = reinterpret_cast<TProfile3D*>(hcorrList->FindObject("ShiftZDCC"));
931+
cal.shiftprofileA = reinterpret_cast<TProfile3D*>(hcorrList->FindObject("ShiftZDCA"));
932+
if (!cal.shiftprofileC || !cal.shiftprofileA) {
933+
LOGF(error, "Shift profile not found in CCDB for runnumber: %d", runnumber);
934+
} else {
935+
LOGF(error, "Shift profile found in CCDB for runnumber: %d", runnumber);
936+
cal.isShiftProfileFound = true;
937+
}
938+
}
939+
}
940+
941+
float coeffshiftxZDCC = 0.0;
942+
float coeffshiftyZDCC = 0.0;
943+
float coeffshiftxZDCA = 0.0;
944+
float coeffshiftyZDCA = 0.0;
945+
946+
for (int ishift = 1; ishift <= nshift; ishift++) {
947+
if (cfgFillCommonRegistry) {
948+
registry.fill(HIST("QA/ShiftZDCC"), centrality, 0.5, ishift - 0.5, std::sin(ishift * 1.0 * psiZDCC));
949+
registry.fill(HIST("QA/ShiftZDCC"), centrality, 1.5, ishift - 0.5, std::cos(ishift * 1.0 * psiZDCC));
950+
registry.fill(HIST("QA/ShiftZDCA"), centrality, 0.5, ishift - 0.5, std::sin(ishift * 1.0 * psiZDCA));
951+
registry.fill(HIST("QA/ShiftZDCA"), centrality, 1.5, ishift - 0.5, std::cos(ishift * 1.0 * psiZDCA));
952+
}
953+
954+
if (cal.isShiftProfileFound) {
955+
int binshiftxZDCC = cal.shiftprofileC->FindBin(centrality, 0.5, ishift - 0.5);
956+
int binshiftyZDCC = cal.shiftprofileC->FindBin(centrality, 1.5, ishift - 0.5);
957+
int binshiftxZDCA = cal.shiftprofileA->FindBin(centrality, 0.5, ishift - 0.5);
958+
int binshiftyZDCA = cal.shiftprofileA->FindBin(centrality, 1.5, ishift - 0.5);
959+
960+
if (binshiftxZDCC > 0 && binshiftyZDCC > 0 && binshiftxZDCA > 0 && binshiftyZDCA > 0) {
961+
coeffshiftxZDCC = cal.shiftprofileC->GetBinContent(binshiftxZDCC);
962+
coeffshiftyZDCC = cal.shiftprofileC->GetBinContent(binshiftyZDCC);
963+
coeffshiftxZDCA = cal.shiftprofileA->GetBinContent(binshiftxZDCA);
964+
coeffshiftyZDCA = cal.shiftprofileA->GetBinContent(binshiftyZDCA);
965+
}
966+
deltaPsiZDCC += deltaPsiZDCC + ((2 / (1.0 * ishift)) * (-coeffshiftxZDCC * std::cos(ishift * 1.0 * psiZDCC) + coeffshiftyZDCC * std::sin(ishift * 1.0 * psiZDCC)));
967+
deltaPsiZDCA += deltaPsiZDCA + ((2 / (1.0 * ishift)) * (-coeffshiftxZDCA * std::cos(ishift * 1.0 * psiZDCA) + coeffshiftyZDCA * std::sin(ishift * 1.0 * psiZDCA)));
968+
}
969+
}
970+
971+
psiZDCCshift += deltaPsiZDCC;
972+
psiZDCAshift += deltaPsiZDCA;
973+
974+
// Normalize angles to [-pi, pi]
975+
psiZDCCshift = std::atan2(std::sin(psiZDCCshift), std::cos(psiZDCCshift));
976+
psiZDCAshift = std::atan2(std::sin(psiZDCAshift), std::cos(psiZDCAshift));
977+
978+
if (cfgFillCommonRegistry) {
979+
registry.fill(HIST("QA/psiZDCA"), psiZDCA);
980+
registry.fill(HIST("QA/psiZDCC"), psiZDCC);
981+
registry.fill(HIST("QA/psiZDCA_shift"), psiZDCAshift);
982+
registry.fill(HIST("QA/psiZDCC_shift"), psiZDCCshift);
983+
}
984+
985+
double qXaShift = std::hypot(qRec[1], qRec[0]) * std::cos(psiZDCAshift);
986+
double qYaShift = std::hypot(qRec[1], qRec[0]) * std::sin(psiZDCAshift);
987+
double qXcShift = std::hypot(qRec[2], qRec[3]) * std::cos(psiZDCCshift);
988+
double qYcShift = std::hypot(qRec[2], qRec[3]) * std::sin(psiZDCCshift);
989+
990+
spTableZDC(runnumber, centrality, v[0], v[1], v[2], qXaShift, qYaShift, qXcShift, qYcShift, isSelected, cal.atIteration, cal.atStep);
871991

872992
qRec.clear();
873993

874994
counter++;
995+
lastRunNumber = runnumber;
875996
return;
876997
}
877998
LOGF(warning, "We return without saving table... -> THis is a problem");
999+
lastRunNumber = runnumber;
8781000
} // end of process
8791001
};
8801002

0 commit comments

Comments
 (0)