Skip to content

Commit 22e2a82

Browse files
YazhenLinYazhenLin
authored
Add files via upload
1 parent 14b28ea commit 22e2a82

File tree

1 file changed

+122
-21
lines changed

1 file changed

+122
-21
lines changed

‎PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx‎

Lines changed: 122 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -2065,6 +2065,7 @@ struct AnalysisDileptonTrack {
20652065
Configurable<std::string> fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"};
20662066
Configurable<int> fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"};
20672067
Configurable<bool> fConfigPublishTripletTable{"cfgPublishTripletTable", false, "Publish the triplet tables, BmesonCandidates"};
2068+
Configurable<bool> fConfigApplyMassEC{"cfgApplyMassEC", false, "Apply fit mass for sideband for the energy correlator study"};
20682069
} fConfigOptions;
20692070

20702071
struct : ConfigurableGroup {
@@ -2137,6 +2138,7 @@ struct AnalysisDileptonTrack {
21372138
bool isDummy = context.mOptions.get<bool>("processDummy");
21382139
bool isMCGen_energycorrelators = context.mOptions.get<bool>("processMCGenEnergyCorrelators") || context.mOptions.get<bool>("processMCGenEnergyCorrelatorsPion");
21392140
bool isMCGen_energycorrelatorsME = context.mOptions.get<bool>("processMCGenEnergyCorrelatorsME") || context.mOptions.get<bool>("processMCGenEnergyCorrelatorsPionME");
2141+
bool isMC_energycorrelatorsUnfolding = context.mOptions.get<bool>("processMCEnergyCorrelatorsUnfolding");
21402142

21412143
if (isDummy) {
21422144
if (isBarrel || isMCGen /*|| isBarrelAsymmetric*/ /*|| isMuon*/) {
@@ -2393,6 +2395,9 @@ struct AnalysisDileptonTrack {
23932395
DefineHistograms(fHistMan, Form("DileptonTrack_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data()), fConfigOptions.fConfigHistogramSubgroups.value.data());
23942396
for (auto& sig : fRecMCSignals) {
23952397
DefineHistograms(fHistMan, Form("DileptonTrackMCMatched_%s_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data(), sig->GetName()), fConfigOptions.fConfigHistogramSubgroups.value.data());
2398+
if (isMC_energycorrelatorsUnfolding) {
2399+
DefineHistograms(fHistMan, Form("DileptonTrackMCUnfold_%s_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data(), sig->GetName()), "");
2400+
}
23962401
}
23972402

23982403
if (!cfgPairing_strCommonTrackCuts.IsNull()) {
@@ -2580,7 +2585,8 @@ struct AnalysisDileptonTrack {
25802585
// compute needed quantities
25812586
VarManager::FillDileptonHadron(dilepton, track, fValuesHadron);
25822587
VarManager::FillDileptonTrackVertexing<TCandidateType, TEventFillMap, TTrackFillMap>(event, lepton1, lepton2, track, fValuesHadron);
2583-
2588+
// for the energy correlator analysis
2589+
VarManager::FillEnergyCorrelator(dilepton, track, fValuesHadron, fConfigOptions.fConfigApplyMassEC);
25842590
if (!track.has_mcParticle()) {
25852591
continue;
25862592
}
@@ -2938,29 +2944,28 @@ struct AnalysisDileptonTrack {
29382944
{
29392945
auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.mcCollisionId());
29402946
groupedMCTracks.bindInternalIndicesTo(&mcTracks);
2941-
for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) {
2947+
for (auto& t1 : groupedMCTracks) {
29422948
auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex());
29432949
// apply kinematic cuts for signal
2944-
if ((t1_raw.pt() < fConfigOptions.fConfigDileptonLowpTCut || t1_raw.pt() > fConfigOptions.fConfigDileptonHighpTCut)) {
2945-
continue;
2946-
}
2947-
if (abs(t1_raw.y()) > fConfigOptions.fConfigDileptonRapCutAbs) {
2950+
if ((t1_raw.pt() < fConfigOptions.fConfigDileptonLowpTCut || t1_raw.pt() > fConfigOptions.fConfigDileptonHighpTCut))
29482951
continue;
2949-
}
2950-
auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex());
2951-
if (TMath::Abs(t2_raw.pdgCode()) == 443 || TMath::Abs(t2_raw.pdgCode()) == 11 || TMath::Abs(t2_raw.pdgCode()) == 22) {
2952-
continue;
2953-
}
2954-
if (t2_raw.pt() < fConfigMCOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigMCOptions.fConfigMCGenHadronEtaAbs.value) {
2955-
continue;
2956-
}
2957-
if (t2_raw.getGenStatusCode() <= 0) {
2952+
if (abs(t1_raw.y()) > fConfigOptions.fConfigDileptonRapCutAbs)
29582953
continue;
2959-
}
2960-
for (auto& sig : fGenMCSignals) {
2961-
if (sig->CheckSignal(true, t1_raw)) {
2962-
VarManager::FillEnergyCorrelatorsMC<THadronMassType>(t1_raw, t2_raw, VarManager::fgValues);
2963-
fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues);
2954+
// for the energy correlators
2955+
for (auto& t2 : groupedMCTracks) {
2956+
auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex());
2957+
if (TMath::Abs(t2_raw.pdgCode()) == 443 || TMath::Abs(t2_raw.pdgCode()) == 11 || TMath::Abs(t2_raw.pdgCode()) == 22)
2958+
continue;
2959+
if (t2_raw.pt() < fConfigMCOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigMCOptions.fConfigMCGenHadronEtaAbs.value) {
2960+
continue;
2961+
}
2962+
if (t2_raw.getGenStatusCode() <= 0)
2963+
continue;
2964+
VarManager::FillEnergyCorrelatorsMC<THadronMassType>(t1_raw, t2_raw, VarManager::fgValues);
2965+
for (auto& sig : fGenMCSignals) {
2966+
if (sig->CheckSignal(true, t1_raw)) {
2967+
fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues);
2968+
}
29642969
}
29652970
}
29662971
}
@@ -3078,6 +3083,99 @@ struct AnalysisDileptonTrack {
30783083
}
30793084
}
30803085

3086+
void processMCEnergyCorrelatorsUnfolding(soa::Filtered<MyEventsSelected> const& events, BCsWithTimestamps const& bcs,
3087+
soa::Join<aod::TrackAssoc, aod::BarrelTrackCuts> const& assocs,
3088+
MyBarrelTracksWithCov const& tracks, soa::Filtered<MyDielectronCandidates> const& dileptons,
3089+
McCollisions const& /* mcEvents*/, McParticles const& mcTracks)
3090+
{
3091+
// set up KF or DCAfitter
3092+
if (events.size() == 0) {
3093+
return;
3094+
}
3095+
if (fCurrentRun != bcs.begin().runNumber()) { // start: runNumber
3096+
initParamsFromCCDB(bcs.begin().timestamp());
3097+
fCurrentRun = bcs.begin().runNumber();
3098+
} // end: runNumber
3099+
for (auto& event : events) {
3100+
if (!event.isEventSelected_bit(0)) {
3101+
continue;
3102+
}
3103+
auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex());
3104+
auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex());
3105+
3106+
uint32_t mcDecision = static_cast<uint32_t>(0);
3107+
size_t isig = 0;
3108+
3109+
for (auto dilepton : groupedDielectrons) {
3110+
// get full track info of tracks based on the index
3111+
auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id());
3112+
auto lepton2 = tracks.rawIteratorAt(dilepton.index1Id());
3113+
if (!lepton1.has_mcParticle() || !lepton2.has_mcParticle()) {
3114+
continue;
3115+
}
3116+
auto lepton1MC = lepton1.mcParticle();
3117+
auto lepton2MC = lepton2.mcParticle();
3118+
if (!lepton1MC.has_mothers())
3119+
continue;
3120+
const auto& motherParticle = lepton1MC.template mothers_first_as<McParticles>();
3121+
// Check that the dilepton has zero charge
3122+
if (dilepton.sign() != 0) {
3123+
continue;
3124+
}
3125+
// loop over track associations
3126+
for (auto& assoc : groupedBarrelAssocs) {
3127+
uint32_t trackSelection = 0;
3128+
// check the cuts fulfilled by this candidate track; if none just continue
3129+
trackSelection = (assoc.isBarrelSelected_raw() & fTrackCutBitMap);
3130+
if (!trackSelection) {
3131+
continue;
3132+
}
3133+
// get the track from this association
3134+
auto track = tracks.rawIteratorAt(assoc.trackId());
3135+
// check that this track is not included in the current dilepton
3136+
if (track.globalIndex() == dilepton.index0Id() || track.globalIndex() == dilepton.index1Id()) {
3137+
continue;
3138+
}
3139+
3140+
if (!track.has_mcParticle()) {
3141+
continue;
3142+
}
3143+
auto trackMC = track.mcParticle();
3144+
VarManager::FillEnergyCorrelatorsMCUnfolding<VarManager::kJpsiHadronMass>(dilepton, track, motherParticle, trackMC, VarManager::fgValues);
3145+
3146+
mcDecision = 0;
3147+
isig = 0;
3148+
for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) {
3149+
if ((*sig)->CheckSignal(true, lepton1MC, lepton2MC, trackMC)) {
3150+
mcDecision |= (static_cast<uint32_t>(1) << isig);
3151+
}
3152+
}
3153+
// Fill histograms for the triplets
3154+
// loop over dilepton / ditrack cuts and MC signals
3155+
for (int icut = 0; icut < fNCuts; icut++) {
3156+
3157+
if (!dilepton.filterMap_bit(icut)) {
3158+
continue;
3159+
}
3160+
3161+
// loop over specified track cuts (the tracks to be combined with the dileptons)
3162+
for (int iTrackCut = 0; iTrackCut < fNCuts; iTrackCut++) {
3163+
3164+
if (!(trackSelection & (static_cast<uint32_t>(1) << iTrackCut))) {
3165+
continue;
3166+
}
3167+
for (uint32_t isig = 0; isig < fRecMCSignals.size(); isig++) {
3168+
if (mcDecision & (static_cast<uint32_t>(1) << isig)) {
3169+
fHistMan->FillHistClass(Form("DileptonTrackMCUnfold_%s_%s_%s", fLegCutNames[icut].Data(), fTrackCutNames[iTrackCut].Data(), fRecMCSignals[isig]->GetName()), VarManager::fgValues);
3170+
}
3171+
}
3172+
}
3173+
} // end loop over track cuts
3174+
} // end loop over dilepton cuts
3175+
} // end loop over dileptons
3176+
} // end loop over events
3177+
}
3178+
30813179
void processDummy(MyEvents&)
30823180
{
30833181
// do nothing
@@ -3091,6 +3189,7 @@ struct AnalysisDileptonTrack {
30913189
PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsPion, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false);
30923190
PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false);
30933191
PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsPionME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false);
3192+
PROCESS_SWITCH(AnalysisDileptonTrack, processMCEnergyCorrelatorsUnfolding, "Loop over MC particle stack and fill response matrix(energy correlators)", false);
30943193
PROCESS_SWITCH(AnalysisDileptonTrack, processDummy, "Dummy function", true);
30953194
};
30963195

@@ -3205,6 +3304,8 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char
32053304
if (classStr.Contains("MCTruthEenergyCorrelators")) {
32063305
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-gen");
32073306
}
3208-
3307+
if (classStr.Contains("DileptonTrackMCUnfold")) {
3308+
dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-unfolding");
3309+
}
32093310
} // end loop over histogram classes
32103311
}

0 commit comments

Comments
 (0)