Skip to content

Commit d3935f2

Browse files
YazhenLinYazhenLin
authored
Add files via upload
1 parent 6635ecc commit d3935f2

File tree

1 file changed

+124
-38
lines changed

1 file changed

+124
-38
lines changed

‎PWGDQ/Tasks/dqEfficiency_withAssoc.cxx‎

Lines changed: 124 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1322,6 +1322,7 @@ struct AnalysisSameEventPairing {
13221322
Configurable<std::string> track{"cfgTrackCuts", "jpsiO2MCdebugCuts2", "Comma separated list of barrel track cuts"};
13231323
Configurable<std::string> muon{"cfgMuonCuts", "", "Comma separated list of muon cuts"};
13241324
Configurable<std::string> pair{"cfgPairCuts", "", "Comma separated list of pair cuts, !!! Use only if you know what you are doing, otherwise leave empty"};
1325+
Configurable<std::string> MCgenAcc{"cfgMCGenAccCut", "", "cut for MC generated particles acceptance"};
13251326
// TODO: Add pair cuts via JSON
13261327
} fConfigCuts;
13271328

@@ -1355,7 +1356,6 @@ struct AnalysisSameEventPairing {
13551356
Configurable<std::string> recSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"};
13561357
Configurable<std::string> recSignalsJSON{"cfgMCRecSignalsJSON", "", "Comma separated list of MC signals (reconstructed) via JSON"};
13571358
Configurable<bool> skimSignalOnly{"cfgSkimSignalOnly", false, "Configurable to select only matched candidates"};
1358-
// Configurable<bool> runMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"};
13591359
} fConfigMC;
13601360

13611361
struct : ConfigurableGroup {
@@ -1383,6 +1383,8 @@ struct AnalysisSameEventPairing {
13831383
std::vector<MCSignal*> fGenMCSignals;
13841384

13851385
std::vector<AnalysisCompositeCut> fPairCuts;
1386+
AnalysisCompositeCut fMCGenAccCut;
1387+
bool fUseMCGenAccCut = false;
13861388

13871389
uint32_t fTrackFilterMask; // mask for the track cuts required in this task to be applied on the barrel cuts produced upstream
13881390
uint32_t fMuonFilterMask; // mask for the muon cuts required in this task to be applied on the muon cuts produced upstream
@@ -1402,7 +1404,7 @@ struct AnalysisSameEventPairing {
14021404
if (context.mOptions.get<bool>("processDummy")) {
14031405
return;
14041406
}
1405-
bool isMCGen = context.mOptions.get<bool>("processMCGen") || context.mOptions.get<bool>("processMCGenWithGrouping");
1407+
bool isMCGen = context.mOptions.get<bool>("processMCGen") || context.mOptions.get<bool>("processMCGenWithGrouping") || context.mOptions.get<bool>("processBarrelOnlySkimmed");
14061408
VarManager::SetDefaultVarNames();
14071409

14081410
fEnableBarrelHistos = context.mOptions.get<bool>("processAllSkimmed") || context.mOptions.get<bool>("processBarrelOnlySkimmed") || context.mOptions.get<bool>("processBarrelOnlyWithCollSkimmed");
@@ -1470,6 +1472,16 @@ struct AnalysisSameEventPairing {
14701472
}
14711473
}
14721474

1475+
// get the mc generated acceptance cut
1476+
TString mcGenAccCutStr = fConfigCuts.MCgenAcc.value;
1477+
if (mcGenAccCutStr != "") {
1478+
AnalysisCut* cut = dqcuts::GetAnalysisCut(mcGenAccCutStr.Data());
1479+
if (cut != nullptr) {
1480+
fMCGenAccCut.AddCut(cut);
1481+
}
1482+
fUseMCGenAccCut = true;
1483+
}
1484+
14731485
// check that the barrel track cuts array required in this task is not empty
14741486
if (!trackCutsStr.IsNull()) {
14751487
// tokenize and loop over the barrel cuts produced by the barrel track selection task
@@ -2112,44 +2124,69 @@ struct AnalysisSameEventPairing {
21122124
} // end loop over events
21132125
}
21142126

2115-
// Preslice<ReducedMCTracks> perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId;
21162127
PresliceUnsorted<ReducedMCTracks> perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId;
21172128

21182129
template <int TPairType>
2119-
void runMCGen(ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
2130+
void runMCGenWithGrouping(MyEventsVtxCovSelected const& events, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
21202131
{
2121-
// loop over mc stack and fill histograms for pure MC truth signals
2122-
// group all the MC tracks which belong to the MC event corresponding to the current reconstructed event
2123-
// auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex());
21242132
uint32_t mcDecision = 0;
21252133
int isig = 0;
2134+
21262135
for (auto& mctrack : mcTracks) {
21272136
VarManager::FillTrackMC(mcTracks, mctrack);
2137+
// if we have a mc generated acceptance cut, apply it here
2138+
if (fUseMCGenAccCut) {
2139+
if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) {
2140+
continue;
2141+
}
2142+
}
21282143
// NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete.
21292144
// NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member.
21302145
// TODO: Use the mcReducedFlags to select signals
21312146
for (auto& sig : fGenMCSignals) {
2132-
if (sig->GetNProngs() != 1) { // NOTE: 1-prong signals required here
2147+
if (sig->CheckSignal(true, mctrack)) {
2148+
fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues);
2149+
}
2150+
}
2151+
}
2152+
// Fill Generated histograms taking into account selected collisions
2153+
for (auto& event : events) {
2154+
if (!event.isEventSelected_bit(0)) {
2155+
continue;
2156+
}
2157+
if (!event.has_reducedMCevent()) {
2158+
continue;
2159+
}
2160+
2161+
for (auto& track : mcTracks) {
2162+
if (track.reducedMCeventId() != event.reducedMCeventId()) {
21332163
continue;
21342164
}
2135-
bool checked = false;
2136-
/*if constexpr (soa::is_soa_filtered_v<TTracksMC>) {
2137-
auto mctrack_raw = groupedMCTracks.rawIteratorAt(mctrack.globalIndex());
2138-
checked = sig.CheckSignal(true, mctrack_raw);
2139-
} else {*/
2140-
checked = sig->CheckSignal(true, mctrack);
2141-
//}
2142-
if (checked) {
2143-
mcDecision |= (static_cast<uint32_t>(1) << isig);
2144-
fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues);
2145-
if (useMiniTree.fConfigMiniTree) {
2146-
auto mcEvent = mcEvents.rawIteratorAt(mctrack.reducedMCeventId());
2147-
dileptonMiniTreeGen(mcDecision, mcEvent.impactParameter(), mctrack.pt(), mctrack.eta(), mctrack.phi(), -999, -999, -999);
2165+
VarManager::FillTrackMC(mcTracks, track);
2166+
// if we have a mc generated acceptance cut, apply it here
2167+
if (fUseMCGenAccCut) {
2168+
if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) {
2169+
continue;
21482170
}
21492171
}
2172+
auto track_raw = mcTracks.rawIteratorAt(track.globalIndex());
2173+
mcDecision = 0;
2174+
isig = 0;
2175+
for (auto& sig : fGenMCSignals) {
2176+
if (sig->CheckSignal(true, track_raw)) {
2177+
mcDecision |= (static_cast<uint32_t>(1) << isig);
2178+
fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues);
2179+
MCTruthTableEffi(VarManager::fgValues[VarManager::kMCPt], VarManager::fgValues[VarManager::kMCEta], VarManager::fgValues[VarManager::kMCY], VarManager::fgValues[VarManager::kMCPhi], VarManager::fgValues[VarManager::kMCVz], VarManager::fgValues[VarManager::kMCVtxZ], VarManager::fgValues[VarManager::kMultFT0A], VarManager::fgValues[VarManager::kMultFT0C], VarManager::fgValues[VarManager::kCentFT0M], VarManager::fgValues[VarManager::kVtxNcontribReal]);
2180+
2181+
if (useMiniTree.fConfigMiniTree) {
2182+
auto mcEvent = mcEvents.rawIteratorAt(track_raw.reducedMCeventId());
2183+
dileptonMiniTreeGen(mcDecision, mcEvent.impactParameter(), track_raw.pt(), track_raw.eta(), track_raw.phi(), -999, -999, -999);
2184+
}
2185+
}
2186+
isig++;
2187+
}
21502188
}
2151-
isig++;
2152-
}
2189+
} // end loop over reconstructed events
21532190

21542191
if (fHasTwoProngGenMCsignals) {
21552192
for (auto& [t1, t2] : combinations(mcTracks, mcTracks)) {
@@ -2161,20 +2198,60 @@ struct AnalysisSameEventPairing {
21612198
continue;
21622199
}
21632200
if (sig->CheckSignal(true, t1_raw, t2_raw)) {
2164-
mcDecision |= (static_cast<uint32_t>(1) << isig);
21652201
VarManager::FillPairMC<TPairType>(t1, t2);
2166-
fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues);
2167-
if (useMiniTree.fConfigMiniTree) {
2168-
// WARNING! To be checked
2169-
dileptonMiniTreeGen(mcDecision, -999, t1.pt(), t1.eta(), t1.phi(), t2.pt(), t2.eta(), t2.phi());
2202+
if (fUseMCGenAccCut) {
2203+
if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) {
2204+
continue;
2205+
}
21702206
}
2207+
fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues);
21712208
}
2172-
isig++;
21732209
}
21742210
}
21752211
}
21762212
}
2177-
} // end runMCGen
2213+
for (auto& event : events) {
2214+
if (!event.isEventSelected_bit(0)) {
2215+
continue;
2216+
}
2217+
if (!event.has_reducedMCevent()) {
2218+
continue;
2219+
}
2220+
// CURRENTLY ONLY FOR 1-GENERATION 2-PRONG SIGNALS
2221+
if (fHasTwoProngGenMCsignals) {
2222+
auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.reducedMCeventId());
2223+
groupedMCTracks.bindInternalIndicesTo(&mcTracks);
2224+
for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) {
2225+
auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex());
2226+
auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex());
2227+
if (t1_raw.reducedMCeventId() == t2_raw.reducedMCeventId()) {
2228+
mcDecision = 0;
2229+
isig = 0;
2230+
for (auto& sig : fGenMCSignals) {
2231+
if (sig->GetNProngs() != 2) { // NOTE: 2-prong signals required here
2232+
continue;
2233+
}
2234+
if (sig->CheckSignal(true, t1_raw, t2_raw)) {
2235+
mcDecision |= (static_cast<uint32_t>(1) << isig);
2236+
VarManager::FillPairMC<TPairType>(t1, t2);
2237+
if (fUseMCGenAccCut) {
2238+
if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) {
2239+
continue;
2240+
}
2241+
}
2242+
fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), VarManager::fgValues);
2243+
if (useMiniTree.fConfigMiniTree) {
2244+
// WARNING! To be checked
2245+
dileptonMiniTreeGen(mcDecision, -999, t1.pt(), t1.eta(), t1.phi(), t2.pt(), t2.eta(), t2.phi());
2246+
}
2247+
}
2248+
isig++;
2249+
}
2250+
}
2251+
}
2252+
} // end loop over reconstructed events
2253+
}
2254+
}
21782255

21792256
void processAllSkimmed(MyEventsVtxCovSelected const& events,
21802257
soa::Join<aod::ReducedTracksAssoc, aod::BarrelTrackCuts, aod::Prefilter> const& barrelAssocs, MyBarrelTracksWithCovWithAmbiguities const& barrelTracks,
@@ -2195,10 +2272,7 @@ struct AnalysisSameEventPairing {
21952272
MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks)
21962273
{
21972274
runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCov, gkTrackFillMapWithCov>(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks);
2198-
// Feature replaced by processMCGen
2199-
/*if (fConfigMC.runMCGenPair) {
2200-
runMCGen<VarManager::kDecayToEE>(mcEvents, mcTracks);
2201-
}*/
2275+
runMCGenWithGrouping<VarManager::kDecayToEE>(events, mcEvents, mcTracks);
22022276
}
22032277

22042278
void processBarrelOnlyWithCollSkimmed(MyEventsVtxCovSelected const& events,
@@ -4261,16 +4335,22 @@ struct AnalysisDileptonTrack {
42614335
for (auto t1 : mcTrackIndices) {
42624336
auto track1 = mcTracks.rawIteratorAt(*(&t1));
42634337
for (auto t2 : mcTrackIndices) {
4264-
if (t1 == t2 || t2 < t1)
4338+
if (t1 == t2)
42654339
continue;
4340+
// if (t2 < t1) continue;
42664341
auto track2 = mcTracks.rawIteratorAt(*(&t2));
42674342
for (auto t3 : mcTrackIndices) {
4268-
if (t3 == t1 || t3 == t2)
4343+
if (t3 == t1)
4344+
continue;
4345+
if (t3 == t2)
42694346
continue;
42704347
auto track3 = mcTracks.rawIteratorAt(*(&t3));
42714348

42724349
for (auto& sig : fRecMCSignals) {
42734350
if (sig->CheckSignal(true, track1, track2, track3)) {
4351+
4352+
VarManager::FillTripleMC<VarManager::kBtoJpsiEEK>(track1, track2, track3, VarManager::fgValues); // nb! hardcoded for jpsiK
4353+
42744354
fHistMan->FillHistClass(Form("MCTruthGenSelBR_%s", sig->GetName()), VarManager::fgValues);
42754355

42764356
// apply kinematic cuts
@@ -4514,13 +4594,19 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char
45144594
}
45154595

45164596
if (classStr.Contains("MCTruthGenPair")) {
4517-
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair");
4597+
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair", histName);
45184598
}
45194599

4520-
if (classStr.Contains("MCTruthGen")) {
4600+
if (classStr.Contains("MCTruthGenSelBR")) {
4601+
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_triple");
4602+
} else if (classStr.Contains("MCTruthGen")) {
45214603
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_track");
45224604
}
45234605

4606+
// if (classStr.Contains("MCTruthGen")) {
4607+
// dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_track");
4608+
// }
4609+
45244610
if (classStr.Contains("DileptonsSelected")) {
45254611
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair", "barrel,vertexing");
45264612
}

0 commit comments

Comments
 (0)