diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 9673dd32a9e..68e65dd9cfe 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -2253,6 +2253,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMCCosChi"] = kMCCosChi; fgVarNamesMap["kMCHadronPt"] = kMCHadronPt; fgVarNamesMap["kMCWeight_before"] = kMCWeight_before; + fgVarNamesMap["kMCEWeight_before"] = kMCEWeight_before; fgVarNamesMap["kMCdeltaeta"] = kMCdeltaeta; fgVarNamesMap["kMCHadronPt"] = kMCHadronPt; fgVarNamesMap["kMCHadronEta"] = kMCHadronEta; @@ -2267,6 +2268,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMCdeltaphi_randomPhi_toward"] = kMCdeltaphi_randomPhi_toward; fgVarNamesMap["kMCdeltaphi_randomPhi_away"] = kMCdeltaphi_randomPhi_away; fgVarNamesMap["kMCdeltaphi_randomPhi_trans"] = kMCdeltaphi_randomPhi_trans; + fgVarNamesMap["kMCHadronpt_randomPhi_trans"] = kMCHadronpt_randomPhi_trans; fgVarNamesMap["kMCCosChi_gen"] = kMCCosChi_gen; fgVarNamesMap["kMCWeight_gen"] = kMCWeight_gen; fgVarNamesMap["kMCdeltaeta_gen"] = kMCdeltaeta_gen; @@ -2504,6 +2506,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kWeight"] = kWeight; fgVarNamesMap["kECWeight"] = kECWeight; fgVarNamesMap["kEWeight_before"] = kEWeight_before; + fgVarNamesMap["kWeight_before"] = kWeight_before; fgVarNamesMap["kPtDau"] = kPtDau; fgVarNamesMap["kEtaDau"] = kEtaDau; fgVarNamesMap["kPhiDau"] = kPhiDau; @@ -2516,6 +2519,7 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kdeltaphi_randomPhi_trans"] = kdeltaphi_randomPhi_trans; fgVarNamesMap["kdeltaphi_randomPhi_toward"] = kdeltaphi_randomPhi_toward; fgVarNamesMap["kdeltaphi_randomPhi_away"] = kdeltaphi_randomPhi_away; + fgVarNamesMap["kPtDau_randomPhi_trans"] = kPtDau_randomPhi_trans; fgVarNamesMap["kdileptonmass"] = kdileptonmass; fgVarNamesMap["kNCorrelationVariables"] = kNCorrelationVariables; fgVarNamesMap["kQuadMass"] = kQuadMass; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index c23ef7a1903..e39a7ae2aaf 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -692,7 +692,9 @@ class VarManager : public TObject kMCdeltaphi_randomPhi_toward, kMCdeltaphi_randomPhi_away, kMCdeltaphi_randomPhi_trans, + kMCHadronpt_randomPhi_trans, kMCWeight_before, + kMCEWeight_before, kMCCosChi_gen, kMCWeight_gen, kMCdeltaeta_gen, @@ -945,6 +947,7 @@ class VarManager : public TObject kPtDau, kCosTheta, kEWeight_before, + kWeight_before, kCosChi_randomPhi_trans, kCosChi_randomPhi_toward, kCosChi_randomPhi_away, @@ -955,6 +958,7 @@ class VarManager : public TObject kdeltaphi_randomPhi_toward, kdeltaphi_randomPhi_away, kdileptonmass, + kPtDau_randomPhi_trans, // Dilepton-track-track variables kQuadMass, @@ -3370,6 +3374,7 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va values[kMCAccweight] = Accweight; values[kMCCosChi] = CosChi; values[kMCWeight_before] = t1.pt() / o2::constants::physics::MassJPsi * Accweight; + values[kMCEWeight_before] = t1.e() / o2::constants::physics::MassJPsi * Accweight; values[kMCCosTheta] = CosTheta; values[kMCdeltaphi] = deltaphi; values[kMCdeltaeta] = deltaeta; @@ -3397,6 +3402,7 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); + values[kMCHadronpt_randomPhi_trans] = v2.pt(); ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, MassHadron); values[kMCCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans); values[kMCWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * Accweight; @@ -5877,7 +5883,8 @@ void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2 values[kWeight] = weight; values[kECWeight] = E_boost / v1.M() * weight; values[kCosTheta] = LorentzTransformJpsihadroncosChi("costheta", v1, v2); - values[kEWeight_before] = v2.Pt() / v1.M() * weight; + values[kEWeight_before] = v2.E() / v1.M() * weight; + values[kWeight_before] = v2.Pt() / v1.M() * weight; values[kPtDau] = v2.pt(); values[kEtaDau] = v2.eta(); values[kPhiDau] = RecoDecay::constrainAngle(v2.phi(), -o2::constants::math::PIHalf); @@ -5901,7 +5908,7 @@ void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2 randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); - + values[kPtDau_randomPhi_trans] = v2.pt(); ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, o2::constants::physics::MassPionCharged); values[kCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans); values[kWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * weight; @@ -5969,8 +5976,7 @@ void VarManager::FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 cons ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_rec(v2_rec.pt(), v2_rec.eta(), randomPhi_trans_rec, o2::constants::physics::MassPionCharged); values[kMCCosChi_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_randomPhi_trans_rec); values[kMCWeight_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("weight_boost", v1_rec, v2_randomPhi_trans_rec) / v1_rec.M() * Effweight_rec; - float randomPhi_trans_gen = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); - ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_gen, MassHadron); + ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_rec, MassHadron); values[kMCCosChi_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("coschi", v1_gen, v2_randomPhi_trans_gen); values[kMCWeight_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_randomPhi_trans_gen) / v1_gen.M() * Accweight_gen; } diff --git a/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx b/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx index 611227d3585..b7b64010f6a 100644 --- a/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx +++ b/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx @@ -104,6 +104,8 @@ struct AnalysisEnergyCorrelator { Configurable fConfigTrackCuts{"cfgTrackCuts", "electronSelection1_ionut", "Comma separated list of barrel track cuts for electrons"}; Configurable fConfigTrackCutsJSON{"cfgTrackCutsJSON", "", "Additional track cuts in JSON"}; Configurable fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Track histograms"}; + Configurable fConfigMCRecTrackSignals{"cfgMCRecTrackSignals", "", "Comma separated list of MC signals (reconstructed)"}; + Configurable fConfigMCRecTrackSignalsJSON{"cfgMCRecTrackSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"}; Configurable fConfigTrackQA{"cfgTrackQA", false, "If true, fill Track QA histograms"}; } fConfigTrackOptions; @@ -155,12 +157,14 @@ struct AnalysisEnergyCorrelator { std::vector fTrackCutNames; std::vector fHadronCutNames; std::vector fHistNamesReco; + std::vector fHistNamesMCMatched; std::map> fTrackHistNames; std::map> fBarrelHistNamesMCmatched; std::map fSelMap; - std::vector fRecMCSignals; // MC signals for reconstructed pairs + std::vector fRecMCTrackSignals; // MC signals for reconstructed tracks + std::vector fRecMCSignals; // MC signals for reconstructed pairs std::vector fGenMCSignals; std::vector fRecMCTripleSignals; // MC signals for reconstructed triples @@ -175,8 +179,8 @@ struct AnalysisEnergyCorrelator { TH2F* hAcceptance_rec; TH2F* hAcceptance_gen; - TH1F* hEfficiency_dilepton; - TH1F* hEfficiency_hadron; + TH2F* hEfficiency_dilepton; + TH2F* hEfficiency_hadron; TH1F* hMasswindow; void init(o2::framework::InitContext& context) @@ -332,6 +336,28 @@ struct AnalysisEnergyCorrelator { } } + // TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function + TString sigRecTrackNamesStr = fConfigTrackOptions.fConfigMCRecTrackSignals.value; + std::unique_ptr objRecTrackSigArray(sigRecTrackNamesStr.Tokenize(",")); + for (int isig = 0; isig < objRecTrackSigArray->GetEntries(); isig++) { + MCSignal* sig_RecTrack = o2::aod::dqmcsignals::GetMCSignal(objRecTrackSigArray->At(isig)->GetName()); + if (sig_RecTrack) { + fRecMCTrackSignals.push_back(sig_RecTrack); + } + } + + // Add the MCSignals from the JSON config + TString addRecTrackSignalsGenStr = fConfigTrackOptions.fConfigMCRecTrackSignalsJSON.value; + if (addRecTrackSignalsGenStr != "") { + std::vector addMCRecTrackSignals = dqmcsignals::GetMCSignalsFromJSON(addRecTrackSignalsGenStr.Data()); + for (auto& mcIt : addMCRecTrackSignals) { + if (mcIt->GetNProngs() > 2) { // NOTE: only 2 prong signals + continue; + } + fRecMCTrackSignals.push_back(mcIt); + } + } + VarManager::SetUseVars(AnalysisCut::fgUsedVars); fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars); @@ -352,6 +378,11 @@ struct AnalysisEnergyCorrelator { TString nameStr = Form("AssocsBarrel_%s", cut->GetName()); fHistNamesReco.push_back(nameStr); histClasses += Form("%s;", nameStr.Data()); + for (auto& sig : fRecMCTrackSignals) { + TString nameStr2 = Form("AssocsBarrelMatched_%s_%s", cut->GetName(), sig->GetName()); + fHistNamesMCMatched.push_back(nameStr2); + histClasses += Form("%s;", nameStr2.Data()); + } } DefineHistograms(fHistMan, histClasses.Data(), fConfigTrackOptions.fConfigAddTrackHistogram.value.data()); } @@ -435,8 +466,8 @@ struct AnalysisEnergyCorrelator { if (!listAccs) { LOG(fatal) << "Problem getting TList object with efficiencies!"; } - hEfficiency_dilepton = static_cast(listAccs->FindObject("hEfficiency_dilepton")); - hEfficiency_hadron = static_cast(listAccs->FindObject("hEfficiency_hadron")); + hEfficiency_dilepton = static_cast(listAccs->FindObject("hEfficiency_dilepton")); + hEfficiency_hadron = static_cast(listAccs->FindObject("hEfficiency_hadron")); hAcceptance_rec = static_cast(listAccs->FindObject("hAcceptance_rec")); hAcceptance_gen = static_cast(listAccs->FindObject("hAcceptance_gen")); hMasswindow = static_cast(listAccs->FindObject("hMasswindow")); @@ -479,9 +510,9 @@ struct AnalysisEnergyCorrelator { float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI); Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi); Accweight_gen = hAcceptance_gen->Interpolate(dilepton_eta - hadron_eta, deltaphi); - float Effdilepton = hEfficiency_dilepton->Interpolate(VarManager::fgValues[VarManager::kPt]); + float Effdilepton = hEfficiency_dilepton->Interpolate(VarManager::fgValues[VarManager::kPt], dilepton_eta); float Masswindow = hMasswindow->Interpolate(VarManager::fgValues[VarManager::kPt]); - float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt()); + float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt(), hadron.eta()); Accweight_gen = Accweight_gen * Effdilepton * Effhadron; if (fConfigDileptonHadronOptions.fConfigApplyEfficiencyME) { Effweight_rec = Effdilepton * Effhadron * Masswindow; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs @@ -647,6 +678,17 @@ struct AnalysisEnergyCorrelator { fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", VarManager::fgValues); } + uint32_t mcDecision = static_cast(0); + // run MC matching for this pair + int isig = 0; + mcDecision = 0; + for (auto sig = fRecMCTrackSignals.begin(); sig != fRecMCTrackSignals.end(); sig++, isig++) { + if (t1.has_mcParticle()) { + if ((*sig)->CheckSignal(true, t1.mcParticle())) { + mcDecision |= (static_cast(1) << isig); + } + } + } // Apply electron cuts and fill histograms int iCut1 = 0; for (auto cut1 = fTrackCuts.begin(); cut1 != fTrackCuts.end(); cut1++, iCut1++) { @@ -654,10 +696,14 @@ struct AnalysisEnergyCorrelator { filter1 |= (static_cast(1) << iCut1); if (fConfigTrackOptions.fConfigTrackQA) { fHistMan->FillHistClass(fHistNamesReco[iCut1], VarManager::fgValues); + for (size_t isig = 0; isig < fRecMCTrackSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(fHistNamesMCMatched[iCut1 * fRecMCTrackSignals.size() + isig], VarManager::fgValues); // matched signal + } + } } } } - // Check opposite charge with t2 for (auto& a2 : groupedAssocs) { auto t2 = a2.template track_as(); diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 40d15f3dad3..34dd7fc7314 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -3493,8 +3493,8 @@ struct AnalysisDileptonTrack { TH2F* hAcceptance_rec; TH2F* hAcceptance_gen; - TH1F* hEfficiency_dilepton; - TH1F* hEfficiency_hadron; + TH2F* hEfficiency_dilepton; + TH2F* hEfficiency_hadron; TH1F* hMasswindow; void init(o2::framework::InitContext& context) @@ -3760,8 +3760,8 @@ struct AnalysisDileptonTrack { if (!listAccs) { LOG(fatal) << "Problem getting TList object with efficiencies!"; } - hEfficiency_dilepton = static_cast(listAccs->FindObject("hEfficiency_dilepton")); - hEfficiency_hadron = static_cast(listAccs->FindObject("hEfficiency_hadron")); + hEfficiency_dilepton = static_cast(listAccs->FindObject("hEfficiency_dilepton")); + hEfficiency_hadron = static_cast(listAccs->FindObject("hEfficiency_hadron")); hAcceptance_rec = static_cast(listAccs->FindObject("hAcceptance_rec")); hAcceptance_gen = static_cast(listAccs->FindObject("hAcceptance_gen")); hMasswindow = static_cast(listAccs->FindObject("hMasswindow")); @@ -3854,9 +3854,9 @@ struct AnalysisDileptonTrack { float hadron_phi = track.phi(); float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI); Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi); - float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt()); + float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta); float Masswindow = hMasswindow->Interpolate(dilepton.pt()); - float Effhadron = hEfficiency_hadron->Interpolate(track.pt()); + float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta); Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; } std::vector fTransRange = fConfigTransRange; @@ -4088,9 +4088,9 @@ struct AnalysisDileptonTrack { float hadron_phi = track.phi(); float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI); Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi); - float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt()); + float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta); float Masswindow = hMasswindow->Interpolate(dilepton.pt()); - float Effhadron = hEfficiency_hadron->Interpolate(track.pt()); + float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta); if (fConfigApplyEfficiencyME) { Effweight_rec = Effdilepton * Effhadron * Masswindow; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs } else { @@ -4098,7 +4098,7 @@ struct AnalysisDileptonTrack { } } std::vector fTransRange = fConfigTransRange; - VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec); + VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, VarManager::fgValues, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec); // loop over dilepton leg cuts and track cuts and fill histograms separately for each combination for (int icut = 0; icut < fNCuts; icut++) {