Skip to content

Commit 2407b81

Browse files
sgarettisgaretti
authored
Merge pull request #2 from sgaretti/dev
[PWGDQ] Adding Mixed events in the table reader with association
2 parents f28f5c4 + 4b157e4 commit 2407b81

File tree

1 file changed

+336
-0
lines changed

1 file changed

+336
-0
lines changed

‎PWGDQ/Tasks/tableReader_withAssoc.cxx‎

Lines changed: 336 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1205,6 +1205,342 @@ struct AnalysisPrefilterSelection {
12051205
PROCESS_SWITCH(AnalysisPrefilterSelection, processDummy, "Do nothing", true);
12061206
};
12071207

1208+
struct AnalysisEventMixing {
1209+
OutputObj<THashList> fOutputList{"output"};
1210+
// Here one should provide the list of electron and muon candidate cuts in the same order as specified in the above
1211+
// single particle selection tasks to preserve the correspondence between the track cut name and its
1212+
// bit position in the cuts bitmap
1213+
// TODO: Create a configurable to specify exactly on which of the bits one should run the event mixing
1214+
Configurable<std::string> fConfigTrackCuts{"cfgTrackCuts", "", "Comma separated list of barrel track cuts"};
1215+
Configurable<std::string> fConfigMuonCuts{"cfgMuonCuts", "", "Comma separated list of muon cuts"};
1216+
Configurable<int> fConfigMixingDepth{"cfgMixingDepth", 100, "Number of Events stored for event mixing"};
1217+
Configurable<std::string> fConfigAddEventMixingHistogram{"cfgAddEventMixingHistogram", "", "Comma separated list of histograms"};
1218+
Configurable<std::string> ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
1219+
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
1220+
Configurable<bool> fConfigAmbiguousHist{"cfgAmbiHist", false, "Enable Ambiguous histograms for time association studies"};
1221+
Configurable<std::string> ccdbPathFlow{"ccdb-path-flow", "Users/c/chizh/FlowResolution", "path to the ccdb object for flow resolution factors"};
1222+
Configurable<bool> fConfigFlowReso{"cfgFlowReso", false, "Enable loading of flow resolution factors from CCDB"};
1223+
Configurable<bool> fConfigSingleMuCumulants{"cfgSingleMuCumulants", false, "Enable loading of flow resolution factors from CCDB"};
1224+
Configurable<std::string> fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"};
1225+
1226+
Service<o2::ccdb::BasicCCDBManager> ccdb;
1227+
1228+
o2::parameters::GRPMagField* grpmag = nullptr;
1229+
TH1D* ResoFlowSP = nullptr; // Resolution factors for flow analysis, this will be loaded from CCDB
1230+
TH1D* ResoFlowEP = nullptr; // Resolution factors for flow analysis, this will be loaded from CCDB
1231+
TH2D* SingleMuv22m = nullptr; // Single muon v22, loaded from CCDB
1232+
TH2D* SingleMuv24m = nullptr; // Single muon v24, loaded from CCDB
1233+
TH2D* SingleMuv22p = nullptr; // Single antimuon v22, loaded from CCDB
1234+
TH2D* SingleMuv24p = nullptr; // Single antimuon v24, loaded from CCDB
1235+
int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc.
1236+
1237+
Filter filterEventSelected = aod::dqanalysisflags::isEventSelected == 1;
1238+
Filter filterTrackSelected = aod::dqanalysisflags::isBarrelSelected > 0;
1239+
Filter filterMuonTrackSelected = aod::dqanalysisflags::isMuonSelected > 0;
1240+
1241+
HistogramManager* fHistMan;
1242+
// NOTE: The bit mask is required to run pairing just based on the desired electron/muon candidate cuts
1243+
uint32_t fTwoTrackFilterMask = 0;
1244+
uint32_t fTwoMuonFilterMask = 0;
1245+
std::vector<std::vector<TString>> fTrackHistNames;
1246+
std::vector<std::vector<TString>> fMuonHistNames;
1247+
std::vector<std::vector<TString>> fTrackMuonHistNames;
1248+
1249+
NoBinningPolicy<aod::dqanalysisflags::MixingHash> hashBin;
1250+
1251+
void init(o2::framework::InitContext& context)
1252+
{
1253+
if (context.mOptions.get<bool>("processDummy")) {
1254+
return;
1255+
}
1256+
1257+
fCurrentRun = 0;
1258+
1259+
ccdb->setURL(ccdburl.value);
1260+
ccdb->setCaching(true);
1261+
ccdb->setLocalObjectValidityChecking();
1262+
1263+
VarManager::SetDefaultVarNames();
1264+
fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars);
1265+
fHistMan->SetUseDefaultVariableNames(kTRUE);
1266+
fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits);
1267+
1268+
// Keep track of all the histogram class names to avoid composing strings in the event mixing pairing
1269+
TString histNames = "";
1270+
if (context.mOptions.get<bool>("processBarrelSkimmed") || context.mOptions.get<bool>("processBarrelVnSkimmed")) {
1271+
TString cutNames = fConfigTrackCuts.value;
1272+
if (!cutNames.IsNull()) {
1273+
std::unique_ptr<TObjArray> objArray(cutNames.Tokenize(","));
1274+
for (int icut = 0; icut < objArray->GetEntries(); ++icut) {
1275+
std::vector<TString> names = {
1276+
Form("PairsBarrelMEPM_%s", objArray->At(icut)->GetName()),
1277+
Form("PairsBarrelMEPP_%s", objArray->At(icut)->GetName()),
1278+
Form("PairsBarrelMEMM_%s", objArray->At(icut)->GetName())};
1279+
histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data());
1280+
fTrackHistNames.push_back(names);
1281+
fTwoTrackFilterMask |= (static_cast<uint32_t>(1) << icut);
1282+
}
1283+
}
1284+
}
1285+
if (context.mOptions.get<bool>("processMuonSkimmed") || context.mOptions.get<bool>("processMuonVnSkimmed") || context.mOptions.get<bool>("processMuonVnCentrSkimmed") || context.mOptions.get<bool>("processMuonVnExtraSkimmed")) {
1286+
TString cutNames = fConfigMuonCuts.value;
1287+
if (!cutNames.IsNull()) {
1288+
std::unique_ptr<TObjArray> objArray(cutNames.Tokenize(","));
1289+
for (int icut = 0; icut < objArray->GetEntries(); ++icut) {
1290+
std::vector<TString> names = {
1291+
Form("PairsMuonMEPM_%s", objArray->At(icut)->GetName()),
1292+
Form("PairsMuonMEPP_%s", objArray->At(icut)->GetName()),
1293+
Form("PairsMuonMEMM_%s", objArray->At(icut)->GetName())};
1294+
if (fConfigAmbiguousHist) {
1295+
histNames += Form("%s;%s;%s;%s_unambiguous;%s_unambiguous;%s_unambiguous;", names[0].Data(), names[1].Data(), names[2].Data(), names[0].Data(), names[1].Data(), names[2].Data());
1296+
} else {
1297+
histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data());
1298+
}
1299+
fMuonHistNames.push_back(names);
1300+
fTwoMuonFilterMask |= (static_cast<uint32_t>(1) << icut);
1301+
}
1302+
}
1303+
}
1304+
if (context.mOptions.get<bool>("processBarrelMuonSkimmed")) {
1305+
TString cutNamesBarrel = fConfigTrackCuts.value;
1306+
TString cutNamesMuon = fConfigMuonCuts.value;
1307+
if (!cutNamesBarrel.IsNull() && !cutNamesMuon.IsNull()) {
1308+
std::unique_ptr<TObjArray> objArrayBarrel(cutNamesBarrel.Tokenize(","));
1309+
std::unique_ptr<TObjArray> objArrayMuon(cutNamesMuon.Tokenize(","));
1310+
if (objArrayBarrel->GetEntries() == objArrayMuon->GetEntries()) { // one must specify equal number of barrel and muon cuts
1311+
for (int icut = 0; icut < objArrayBarrel->GetEntries(); ++icut) {
1312+
std::vector<TString> names = {
1313+
Form("PairsEleMuMEPM_%s_%s", objArrayBarrel->At(icut)->GetName(), objArrayMuon->At(icut)->GetName()),
1314+
Form("PairsEleMuMEPP_%s_%s", objArrayBarrel->At(icut)->GetName(), objArrayMuon->At(icut)->GetName()),
1315+
Form("PairsEleMuMEMM_%s_%s", objArrayBarrel->At(icut)->GetName(), objArrayMuon->At(icut)->GetName())};
1316+
histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data());
1317+
fTrackMuonHistNames.push_back(names);
1318+
fTwoTrackFilterMask |= (static_cast<uint32_t>(1) << icut);
1319+
fTwoMuonFilterMask |= (static_cast<uint32_t>(1) << icut);
1320+
}
1321+
}
1322+
}
1323+
}
1324+
1325+
DefineHistograms(fHistMan, histNames.Data(), fConfigAddEventMixingHistogram); // define all histograms
1326+
// Additional histograms via JSON
1327+
dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigAddJSONHistograms.value.c_str());
1328+
VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill
1329+
fOutputList.setObject(fHistMan->GetMainHistogramList());
1330+
}
1331+
1332+
template <uint32_t TEventFillMap, int TPairType, typename TTracks1, typename TTracks2>
1333+
void runMixedPairing(TTracks1 const& tracks1, TTracks2 const& tracks2)
1334+
{
1335+
1336+
unsigned int ncuts = fTrackHistNames.size();
1337+
std::vector<std::vector<TString>> histNames = fTrackHistNames;
1338+
if constexpr (TPairType == pairTypeMuMu) {
1339+
ncuts = fMuonHistNames.size();
1340+
histNames = fMuonHistNames;
1341+
}
1342+
if constexpr (TPairType == pairTypeEMu) {
1343+
ncuts = fTrackMuonHistNames.size();
1344+
histNames = fTrackMuonHistNames;
1345+
}
1346+
1347+
uint32_t twoTrackFilter = 0;
1348+
uint32_t mult_dimuons = 0;
1349+
for (auto& track1 : tracks1) {
1350+
for (auto& track2 : tracks2) {
1351+
if constexpr (TPairType == VarManager::kDecayToMuMu) {
1352+
twoTrackFilter = static_cast<uint32_t>(track1.isMuonSelected()) & static_cast<uint32_t>(track2.isMuonSelected()) & fTwoMuonFilterMask;
1353+
}
1354+
if (twoTrackFilter && track1.sign() * track2.sign() < 0) {
1355+
mult_dimuons++;
1356+
}
1357+
} // end for (track2)
1358+
} // end for (track1)
1359+
VarManager::fgValues[VarManager::kMultDimuonsME] = mult_dimuons;
1360+
1361+
twoTrackFilter = 0;
1362+
for (auto& track1 : tracks1) {
1363+
for (auto& track2 : tracks2) {
1364+
if constexpr (TPairType == VarManager::kDecayToEE) {
1365+
twoTrackFilter = static_cast<uint32_t>(track1.isBarrelSelected()) & static_cast<uint32_t>(track2.isBarrelSelected()) & fTwoTrackFilterMask;
1366+
}
1367+
if constexpr (TPairType == VarManager::kDecayToMuMu) {
1368+
twoTrackFilter = static_cast<uint32_t>(track1.isMuonSelected()) & static_cast<uint32_t>(track2.isMuonSelected()) & fTwoMuonFilterMask;
1369+
if (fConfigSingleMuCumulants) {
1370+
VarManager::FillTwoMixEventsCumulants(SingleMuv22m, SingleMuv24m, SingleMuv22p, SingleMuv24p, track1, track2);
1371+
}
1372+
}
1373+
if constexpr (TPairType == VarManager::kElectronMuon) {
1374+
twoTrackFilter = static_cast<uint32_t>(track1.isBarrelSelected()) & static_cast<uint32_t>(track2.isMuonSelected()) & fTwoTrackFilterMask;
1375+
}
1376+
1377+
if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue
1378+
continue;
1379+
}
1380+
VarManager::FillPairME<TEventFillMap, TPairType>(track1, track2);
1381+
1382+
for (unsigned int icut = 0; icut < ncuts; icut++) {
1383+
if (twoTrackFilter & (static_cast<uint32_t>(1) << icut)) {
1384+
if (track1.sign() * track2.sign() < 0) {
1385+
fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues);
1386+
if (fConfigAmbiguousHist && !(track1.isAmbiguous() || track2.isAmbiguous())) {
1387+
fHistMan->FillHistClass(Form("%s_unambiguous", histNames[icut][0].Data()), VarManager::fgValues);
1388+
}
1389+
} else {
1390+
if (track1.sign() > 0) {
1391+
fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues);
1392+
if (fConfigAmbiguousHist && !(track1.isAmbiguous() || track2.isAmbiguous())) {
1393+
fHistMan->FillHistClass(Form("%s_unambiguous", histNames[icut][1].Data()), VarManager::fgValues);
1394+
}
1395+
} else {
1396+
fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues);
1397+
if (fConfigAmbiguousHist && !(track1.isAmbiguous() || track2.isAmbiguous())) {
1398+
fHistMan->FillHistClass(Form("%s_unambiguous", histNames[icut][2].Data()), VarManager::fgValues);
1399+
}
1400+
}
1401+
}
1402+
} // end if (filter bits)
1403+
} // end for (cuts)
1404+
} // end for (track2)
1405+
} // end for (track1)
1406+
}
1407+
1408+
// barrel-barrel and muon-muon event mixing
1409+
template <int TPairType, uint32_t TEventFillMap, typename TEvents, typename TTracks>
1410+
void runSameSide(TEvents& events, TTracks const& tracks, Preslice<TTracks>& preSlice)
1411+
{
1412+
if (events.size() > 0 && fCurrentRun != events.begin().runNumber()) {
1413+
grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, events.begin().timestamp());
1414+
if (grpmag != nullptr) {
1415+
VarManager::SetMagneticField(grpmag->getNominalL3Field());
1416+
} else {
1417+
LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", events.begin().timestamp());
1418+
}
1419+
if (fConfigFlowReso) {
1420+
TString PathFlow = ccdbPathFlow.value;
1421+
TString ccdbPathFlowSP = Form("%s/ScalarProduct", PathFlow.Data());
1422+
TString ccdbPathFlowEP = Form("%s/EventPlane", PathFlow.Data());
1423+
ResoFlowSP = ccdb->getForTimeStamp<TH1D>(ccdbPathFlowSP.Data(), events.begin().timestamp());
1424+
ResoFlowEP = ccdb->getForTimeStamp<TH1D>(ccdbPathFlowEP.Data(), events.begin().timestamp());
1425+
if (ResoFlowSP == nullptr || ResoFlowEP == nullptr) {
1426+
LOGF(fatal, "Resolution factor is not available in CCDB at timestamp=%llu", events.begin().timestamp());
1427+
}
1428+
}
1429+
if (fConfigSingleMuCumulants) {
1430+
TString PathFlow = ccdbPathFlow.value;
1431+
TString ccdbPathMuv22m = Form("%s/SingleMuv22m", PathFlow.Data());
1432+
TString ccdbPathMuv24m = Form("%s/SingleMuv24m", PathFlow.Data());
1433+
TString ccdbPathMuv22p = Form("%s/SingleMuv22p", PathFlow.Data());
1434+
TString ccdbPathMuv24p = Form("%s/SingleMuv24p", PathFlow.Data());
1435+
SingleMuv22m = ccdb->getForTimeStamp<TH2D>(ccdbPathMuv22m.Data(), events.begin().timestamp());
1436+
SingleMuv24m = ccdb->getForTimeStamp<TH2D>(ccdbPathMuv24m.Data(), events.begin().timestamp());
1437+
SingleMuv22p = ccdb->getForTimeStamp<TH2D>(ccdbPathMuv22p.Data(), events.begin().timestamp());
1438+
SingleMuv24p = ccdb->getForTimeStamp<TH2D>(ccdbPathMuv24p.Data(), events.begin().timestamp());
1439+
if (SingleMuv22m == nullptr || SingleMuv24m == nullptr || SingleMuv22p == nullptr || SingleMuv24p == nullptr) {
1440+
LOGF(fatal, "Single muon cumulants are not available in CCDB at timestamp=%llu", events.begin().timestamp());
1441+
}
1442+
}
1443+
fCurrentRun = events.begin().runNumber();
1444+
}
1445+
1446+
events.bindExternalIndices(&tracks);
1447+
int mixingDepth = fConfigMixingDepth.value;
1448+
for (auto& [event1, event2] : selfCombinations(hashBin, mixingDepth, -1, events, events)) {
1449+
VarManager::ResetValues(0, VarManager::kNVars);
1450+
VarManager::FillEvent<TEventFillMap>(event1, VarManager::fgValues);
1451+
1452+
auto tracks1 = tracks.sliceBy(preSlice, event1.globalIndex());
1453+
tracks1.bindExternalIndices(&events);
1454+
1455+
auto tracks2 = tracks.sliceBy(preSlice, event2.globalIndex());
1456+
tracks2.bindExternalIndices(&events);
1457+
1458+
VarManager::FillTwoMixEvents<TEventFillMap>(event1, event2, tracks1, tracks2);
1459+
if (fConfigFlowReso) {
1460+
VarManager::FillTwoMixEventsFlowResoFactor(ResoFlowSP, ResoFlowEP);
1461+
}
1462+
runMixedPairing<TEventFillMap, TPairType>(tracks1, tracks2);
1463+
} // end event loop
1464+
}
1465+
1466+
// barrel-muon event mixing
1467+
template <uint32_t TEventFillMap, typename TEvents, typename TTracks, typename TMuons>
1468+
void runBarrelMuon(TEvents& events, TTracks const& tracks, TMuons const& muons)
1469+
{
1470+
if (events.size() > 0 && fCurrentRun != events.begin().runNumber()) {
1471+
grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, events.begin().timestamp());
1472+
if (grpmag != nullptr) {
1473+
VarManager::SetMagneticField(grpmag->getNominalL3Field());
1474+
} else {
1475+
LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", events.begin().timestamp());
1476+
}
1477+
fCurrentRun = events.begin().runNumber();
1478+
}
1479+
1480+
events.bindExternalIndices(&muons);
1481+
1482+
for (auto& [event1, event2] : selfCombinations(hashBin, 100, -1, events, events)) {
1483+
VarManager::ResetValues(0, VarManager::kNVars);
1484+
VarManager::FillEvent<TEventFillMap>(event1, VarManager::fgValues);
1485+
1486+
auto tracks1 = tracks.sliceBy(perEventsSelectedT, event1.globalIndex());
1487+
tracks1.bindExternalIndices(&events);
1488+
1489+
auto muons2 = muons.sliceBy(perEventsSelectedM, event2.globalIndex());
1490+
muons2.bindExternalIndices(&events);
1491+
1492+
runMixedPairing<TEventFillMap, pairTypeEMu>(tracks1, muons2);
1493+
} // end event loop
1494+
}
1495+
1496+
Preslice<soa::Filtered<MyBarrelTracksSelected>> perEventsSelectedT = aod::reducedtrack::reducedeventId;
1497+
Preslice<soa::Filtered<MyMuonTracksSelected>> perEventsSelectedM = aod::reducedmuon::reducedeventId;
1498+
1499+
void processBarrelSkimmed(soa::Filtered<MyEventsHashSelected>& events, soa::Filtered<MyBarrelTracksSelected> const& tracks)
1500+
{
1501+
runSameSide<pairTypeEE, gkEventFillMap>(events, tracks, perEventsSelectedT);
1502+
}
1503+
void processMuonSkimmed(soa::Filtered<MyEventsHashSelected>& events, soa::Filtered<MyMuonTracksSelected> const& muons)
1504+
{
1505+
runSameSide<pairTypeMuMu, gkEventFillMap>(events, muons, perEventsSelectedM);
1506+
}
1507+
void processBarrelMuonSkimmed(soa::Filtered<MyEventsHashSelected>& events, soa::Filtered<MyBarrelTracksSelected> const& tracks, soa::Filtered<MyMuonTracksSelected> const& muons)
1508+
{
1509+
runBarrelMuon<gkEventFillMap>(events, tracks, muons);
1510+
}
1511+
void processBarrelVnSkimmed(soa::Filtered<MyEventsHashSelectedQvector>& events, soa::Filtered<MyBarrelTracksSelected> const& tracks)
1512+
{
1513+
runSameSide<pairTypeEE, gkEventFillMapWithQvector>(events, tracks, perEventsSelectedT);
1514+
}
1515+
void processMuonVnSkimmed(soa::Filtered<MyEventsHashSelectedQvector>& events, soa::Filtered<MyMuonTracksSelected> const& muons)
1516+
{
1517+
runSameSide<pairTypeMuMu, gkEventFillMapWithQvector>(events, muons, perEventsSelectedM);
1518+
}
1519+
void processMuonVnCentrSkimmed(soa::Filtered<MyEventsHashSelectedQvectorCentr>& events, soa::Filtered<MyMuonTracksSelected> const& muons)
1520+
{
1521+
runSameSide<pairTypeMuMu, gkEventFillMapWithQvectorCentr>(events, muons, perEventsSelectedM);
1522+
}
1523+
void processMuonVnExtraSkimmed(soa::Filtered<MyEventsHashSelectedQvectorExtra>& events, soa::Filtered<MyMuonTracksSelected> const& muons)
1524+
{
1525+
runSameSide<pairTypeMuMu, gkEventFillMapWithCovQvectorExtraWithRefFlow>(events, muons, perEventsSelectedM);
1526+
}
1527+
// TODO: This is a dummy process function for the case when the user does not want to run any of the process functions (no event mixing)
1528+
// If there is no process function enabled, the workflow hangs
1529+
void processDummy(MyEvents&)
1530+
{
1531+
// do nothing
1532+
}
1533+
1534+
PROCESS_SWITCH(AnalysisEventMixing, processBarrelSkimmed, "Run barrel-barrel mixing on skimmed tracks", false);
1535+
PROCESS_SWITCH(AnalysisEventMixing, processMuonSkimmed, "Run muon-muon mixing on skimmed muons", false);
1536+
PROCESS_SWITCH(AnalysisEventMixing, processBarrelMuonSkimmed, "Run barrel-muon mixing on skimmed tracks/muons", false);
1537+
PROCESS_SWITCH(AnalysisEventMixing, processBarrelVnSkimmed, "Run barrel-barrel vn mixing on skimmed tracks", false);
1538+
PROCESS_SWITCH(AnalysisEventMixing, processMuonVnSkimmed, "Run muon-muon vn mixing on skimmed tracks", false);
1539+
PROCESS_SWITCH(AnalysisEventMixing, processMuonVnCentrSkimmed, "Run muon-muon vn mixing on skimmed tracks from central framework", false);
1540+
PROCESS_SWITCH(AnalysisEventMixing, processMuonVnExtraSkimmed, "Run muon-muon vn mixing on skimmed tracks from GFW", false);
1541+
PROCESS_SWITCH(AnalysisEventMixing, processDummy, "Dummy function", false);
1542+
};
1543+
12081544
// Run the same-event pairing
12091545
// This task assumes that both legs of the resonance fulfill the same cuts (symmetric decay channel)
12101546
// Runs combinatorics for barrel-barrel, muon-muon and barrel-muon combinations

0 commit comments

Comments
 (0)