diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 5cbc6a50c62..27f6c17d839 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -3333,6 +3333,10 @@ struct AnalysisDileptonTrack { Configurable> fConfigFitmassEC{"cfgTFitmassEC", std::vector{-0.541438, 2.8, 3.2}, "parameter from the fit fuction and fit range"}; Configurable> fConfigTransRange{"cfgTransRange", std::vector{0.333333, 0.666667}, "Transverse region for the energy correlstor analysis"}; + Configurable fConfigApplyEfficiency{"cfgApplyEfficiency", false, "If true, apply efficiency correction for the energy correlator study"}; + Configurable fConfigApplyEfficiencyME{"cfgApplyEfficiencyME", false, "If true, apply efficiency correction for the energy correlator study"}; + Configurable fConfigAccCCDBPath{"AccCCDBPath", "Users/y/yalin/pptest/test2", "Path of the efficiency corrections"}; + int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. int fNCuts; // number of dilepton leg cuts int fNLegCuts; @@ -3366,6 +3370,12 @@ struct AnalysisDileptonTrack { TF1* fMassBkg = nullptr; + TH2F* hAcceptance_rec; + TH2F* hAcceptance_gen; + TH1F* hEfficiency_dilepton; + TH1F* hEfficiency_hadron; + TH1F* hMasswindow; + void init(o2::framework::InitContext& context) { bool isBarrel = context.mOptions.get("processBarrelSkimmed"); @@ -3623,6 +3633,22 @@ struct AnalysisDileptonTrack { } } + void initAccFromCCDB(uint64_t timestamp) + { + TList* listAccs = fCCDB->getForTimeStamp(fConfigAccCCDBPath, timestamp); + 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")); + hAcceptance_rec = static_cast(listAccs->FindObject("hAcceptance_rec")); + hAcceptance_gen = static_cast(listAccs->FindObject("hAcceptance_gen")); + hMasswindow = static_cast(listAccs->FindObject("hMasswindow")); + if (!hAcceptance_rec || !hAcceptance_gen || !hEfficiency_dilepton || !hEfficiency_hadron || !hMasswindow) { + LOG(fatal) << "Problem getting histograms from the TList object with efficiencies!"; + } + } + // Template function to run pair - hadron combinations template void runDileptonHadron(TEvent const& event, TTrackAssocs const& assocs, TTracks const& tracks, TDileptons const& dileptons) @@ -3699,8 +3725,21 @@ struct AnalysisDileptonTrack { VarManager::FillDileptonTrackVertexing(event, lepton1, lepton2, track, fValuesHadron); // for the energy correlator analysis + float Effweight_rec = 1.0f; + if (fConfigApplyEfficiency) { + float dilepton_eta = dilepton.eta(); + float dilepton_phi = dilepton.phi(); + float hadron_eta = track.eta(); + 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 Masswindow = hMasswindow->Interpolate(dilepton.pt()); + float Effhadron = hEfficiency_hadron->Interpolate(track.pt()); + Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; + } std::vector fTransRange = fConfigTransRange; - VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom()); + VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec); // table to be written out for ML analysis BmesonsTable(event.runNumber(), event.globalIndex(), event.timestamp(), fValuesHadron[VarManager::kPairMass], dilepton.mass(), fValuesHadron[VarManager::kDeltaMass], fValuesHadron[VarManager::kPairPt], fValuesHadron[VarManager::kPairEta], fValuesHadron[VarManager::kPairPhi], fValuesHadron[VarManager::kPairRap], @@ -3805,6 +3844,9 @@ struct AnalysisDileptonTrack { } if (fCurrentRun != events.begin().runNumber()) { // start: runNumber initParamsFromCCDB(events.begin().timestamp()); + if (fConfigApplyEfficiency) { + initAccFromCCDB(events.begin().timestamp()); + } fCurrentRun = events.begin().runNumber(); } // end: runNumber for (auto& event : events) { @@ -3862,6 +3904,14 @@ struct AnalysisDileptonTrack { if (events.size() == 0) { return; } + + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + if (fConfigApplyEfficiency) { + initAccFromCCDB(events.begin().timestamp()); + } + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + events.bindExternalIndices(&dileptons); events.bindExternalIndices(&assocs); @@ -3909,8 +3959,25 @@ struct AnalysisDileptonTrack { VarManager::FillDileptonHadron(dilepton, track, VarManager::fgValues); // for the energy correlator analysis + float Effweight_rec = 1.0f; + if (fConfigApplyEfficiency) { + float dilepton_eta = dilepton.eta(); + float dilepton_phi = dilepton.phi(); + float hadron_eta = track.eta(); + 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 Masswindow = hMasswindow->Interpolate(dilepton.pt()); + float Effhadron = hEfficiency_hadron->Interpolate(track.pt()); + 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 { + Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; // apply acceptance and efficiency correction for the real pairs + } + } std::vector fTransRange = fConfigTransRange; - VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom()); + VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, 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++) { diff --git a/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx b/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx index ae6192a9dbb..ffd5c2b3939 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx @@ -63,6 +63,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; using MyCollisions = o2::soa::Join; +using MyMCCollisions = o2::soa::Join; using MyBCs = o2::soa::Join; struct Photon { @@ -96,6 +97,11 @@ struct Photon { bool onDCal; // Checks whether photon is in phi region of the DCal, otherwise: EMCal }; +enum SubGeneratorId { + none = -1, + mbGap = 0 +}; + struct Meson { Meson(Photon p1, Photon p2) : pgamma1(p1), pgamma2(p2) @@ -156,6 +162,7 @@ struct EmcalPi0Qc { Configurable mMinOpenAngleCut{"mMinOpenAngleCut", 0.0202, "apply min opening angle cut"}; Configurable mClusterDefinition{"mClusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"}; Configurable mSplitEMCalDCal{"mSplitEMCalDCal", 0, "Create and fill inv mass histograms for photons on EMCal and DCal individually"}; + Configurable mDoSumw2{"mDoSumw2", 1, "enable Sumw2 for all histograms"}; std::vector mVetoBCIDs; std::vector mSelectBCIDs; @@ -173,6 +180,7 @@ struct EmcalPi0Qc { // event mixing class EventMixVec evtMix; + float mWeight = 1.0f; // current event weight, by default 1.0 o2::ccdb::CcdbApi ccdbApi; int lastRunNumber = -1; // get the runnumber to obtain the SOR of the run to get t - SOR in (s) later @@ -197,8 +205,27 @@ struct EmcalPi0Qc { const AxisSpec invmassAxis{invmassBinning, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; const AxisSpec ptAxis{pTBinning, "#it{p}_{T} (GeV/#it{c})"}; - if (doprocessCollision) { - mHistManager.add("events", "events;;#it{count}", HistType::kTH1F, {{7, 0.5, 7.5}}); + if (doprocessCollisionMC) { + mHistManager.add("eventsWithoutWeight", "events without weight;;#it{count}", HistType::kTH1F, {{8, 0.5, 8.5}}, mDoSumw2.value); + auto heventWithoutWeight = mHistManager.get(HIST("eventsWithoutWeight")); + heventWithoutWeight->GetXaxis()->SetBinLabel(1, "All events"); + heventWithoutWeight->GetXaxis()->SetBinLabel(2, "Has MC collision"); + heventWithoutWeight->GetXaxis()->SetBinLabel(3, "sel8"); + heventWithoutWeight->GetXaxis()->SetBinLabel(4, "EMCal readout"); + heventWithoutWeight->GetXaxis()->SetBinLabel(5, "1+ Contributor"); + heventWithoutWeight->GetXaxis()->SetBinLabel(6, "z<10cm"); + heventWithoutWeight->GetXaxis()->SetBinLabel(7, "unique col"); + heventWithoutWeight->GetXaxis()->SetBinLabel(8, "EMCal cell>0"); + + // histogram the number of gap events and signal events (2 bins, bin 1 gap bin 2 signal) + mHistManager.add("signalGapEvents", "number of signal and gap events;;#it{count}", HistType::kTH1F, {{2, 0.5, 2.5}}, mDoSumw2.value); + auto hsignalGapEvents = mHistManager.get(HIST("signalGapEvents")); + hsignalGapEvents->GetXaxis()->SetBinLabel(1, "Gap events"); + hsignalGapEvents->GetXaxis()->SetBinLabel(2, "Signal events"); + } + + if (doprocessCollision || doprocessCollisionMC) { + mHistManager.add("events", "events;;#it{count}", HistType::kTH1F, {{7, 0.5, 7.5}}, mDoSumw2.value); auto heventType = mHistManager.get(HIST("events")); heventType->GetXaxis()->SetBinLabel(1, "All events"); heventType->GetXaxis()->SetBinLabel(2, "sel8"); @@ -207,49 +234,65 @@ struct EmcalPi0Qc { heventType->GetXaxis()->SetBinLabel(5, "z<10cm"); heventType->GetXaxis()->SetBinLabel(6, "unique col"); heventType->GetXaxis()->SetBinLabel(7, "EMCal cell>0"); - mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", HistType::kTH1F, {{200, -20, 20}}); - mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", HistType::kTH1F, {{200, -20, 20}}); - mHistManager.add("hEventPerTime", "number of events per time", HistType::kTH1F, {collisionTimeAxis}); + mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", HistType::kTH1F, {{200, -20, 20}}, mDoSumw2.value); + mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", HistType::kTH1F, {{200, -20, 20}}, mDoSumw2.value); + mHistManager.add("hEventPerTime", "number of events per time", HistType::kTH1F, {collisionTimeAxis}, mDoSumw2.value); + + // emcal hardware triggers + mHistManager.add("eventsEMCALHardwareTriggers", "events with EMCal hardware triggers;;#it{count}", HistType::kTH1F, {{12, 0.5, 12.5}}, mDoSumw2.value); + auto heventsEMCALHardwareTriggers = mHistManager.get(HIST("eventsEMCALHardwareTriggers")); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(1, "kTVXinEMC"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(2, "kEMC7"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(3, "kDMC7"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(4, "kEG1"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(5, "kEG2"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(6, "kDG1"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(7, "kDG2"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(8, "kEJ1"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(9, "kEJ2"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(10, "kDJ1"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(11, "kDJ2"); + heventsEMCALHardwareTriggers->GetXaxis()->SetBinLabel(12, "All"); } if (doprocessAmbiguous) { - mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", HistType::kTH1F, {bcAxis}); - mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", HistType::kTH1F, {bcAxis}); + mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", HistType::kTH1F, {bcAxis}, mDoSumw2.value); + mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", HistType::kTH1F, {bcAxis}, mDoSumw2.value); } // cluster properties for (const bool& iBeforeCuts : {false, true}) { const char* clusterDirectory = iBeforeCuts ? "ClustersBeforeCuts" : "ClustersAfterCuts"; - mHistManager.add(Form("%s/clusterE", clusterDirectory), "Energy of cluster", HistType::kTH1F, {energyAxis}); - mHistManager.add(Form("%s/clusterE_SimpleBinning", clusterDirectory), "Energy of cluster", HistType::kTH1F, {{400, 0, 100, "#it{E} (GeV)"}}); - mHistManager.add(Form("%s/clusterTime", clusterDirectory), "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}); - mHistManager.add(Form("%s/clusterEtaPhi", clusterDirectory), "Eta and phi of cluster", HistType::kTH2F, {{100, -1, 1, "#eta"}, {100, 0, o2::constants::math::TwoPI, "#phi"}}); - mHistManager.add(Form("%s/clusterM02", clusterDirectory), "M02 of cluster", HistType::kTH1F, {{400, 0, 5, "#it{M}_{02}"}}); - mHistManager.add(Form("%s/clusterM20", clusterDirectory), "M20 of cluster", HistType::kTH1F, {{400, 0, 2.5, "#it{M}_{20}"}}); - mHistManager.add(Form("%s/clusterNLM", clusterDirectory), "Number of local maxima of cluster", HistType::kTH1I, {{10, 0, 10, "#it{N}_{local maxima}"}}); - mHistManager.add(Form("%s/clusterNCells", clusterDirectory), "Number of cells in cluster", HistType::kTH1I, {{50, 0, 50, "#it{N}_{cells}"}}); - mHistManager.add(Form("%s/clusterDistanceToBadChannel", clusterDirectory), "Distance to bad channel", HistType::kTH1F, {{100, 0, 100, "#it{d}"}}); + mHistManager.add(Form("%s/clusterE", clusterDirectory), "Energy of cluster", HistType::kTH1F, {energyAxis}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterE_SimpleBinning", clusterDirectory), "Energy of cluster", HistType::kTH1F, {{400, 0, 100, "#it{E} (GeV)"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterTime", clusterDirectory), "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterEtaPhi", clusterDirectory), "Eta and phi of cluster", HistType::kTH2F, {{100, -1, 1, "#eta"}, {100, 0, o2::constants::math::TwoPI, "#phi"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterM02", clusterDirectory), "M02 of cluster", HistType::kTH1F, {{400, 0, 5, "#it{M}_{02}"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterM20", clusterDirectory), "M20 of cluster", HistType::kTH1F, {{400, 0, 2.5, "#it{M}_{20}"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterNLM", clusterDirectory), "Number of local maxima of cluster", HistType::kTH1I, {{10, 0, 10, "#it{N}_{local maxima}"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterNCells", clusterDirectory), "Number of cells in cluster", HistType::kTH1I, {{50, 0, 50, "#it{N}_{cells}"}}, mDoSumw2.value); + mHistManager.add(Form("%s/clusterDistanceToBadChannel", clusterDirectory), "Distance to bad channel", HistType::kTH1F, {{100, 0, 100, "#it{d}"}}, mDoSumw2.value); } // meson related histograms - mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtMixedBackground", "invariant mass and pT of mixed background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtMixedBackground", "invariant mass and pT of mixed background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); if (mSplitEMCalDCal) { - mHistManager.add("invMassVsPt_EMCal", "invariant mass and pT of meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtBackground_EMCal", "invariant mass and pT of background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtMixedBackground_EMCal", "invariant mass and pT of mixed background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPt_DCal", "invariant mass and pT of meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtBackground_DCal", "invariant mass and pT of background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); - mHistManager.add("invMassVsPtMixedBackground_DCal", "invariant mass and pT of mixed background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPt_EMCal", "invariant mass and pT of meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtBackground_EMCal", "invariant mass and pT of background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtMixedBackground_EMCal", "invariant mass and pT of mixed background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPt_DCal", "invariant mass and pT of meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtBackground_DCal", "invariant mass and pT of background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); + mHistManager.add("invMassVsPtMixedBackground_DCal", "invariant mass and pT of mixed background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}, mDoSumw2.value); } // add histograms per supermodule for (int ism = 0; ism < 20; ++ism) { - mHistManager.add(Form("clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM%d", ism), Form("Cluster time vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {clusterTimeAxis, collisionTimeAxis}); - mHistManager.add(Form("clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM%d", ism), Form("Cluster number of cells vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{50, 0, 50}, collisionTimeAxis}); - mHistManager.add(Form("clusterM02VsTimeStamp/clusterM02VsTimeStampSM%d", ism), Form("Cluster M02 vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{400, 0, 5}, collisionTimeAxis}); - mHistManager.add(Form("mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM%d", ism), Form("invariant mass vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {invmassAxis, collisionTimeAxis}); + mHistManager.add(Form("clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM%d", ism), Form("Cluster time vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {clusterTimeAxis, collisionTimeAxis}, mDoSumw2.value); + mHistManager.add(Form("clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM%d", ism), Form("Cluster number of cells vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{50, 0, 50}, collisionTimeAxis}, mDoSumw2.value); + mHistManager.add(Form("clusterM02VsTimeStamp/clusterM02VsTimeStampSM%d", ism), Form("Cluster M02 vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{400, 0, 5}, collisionTimeAxis}, mDoSumw2.value); + mHistManager.add(Form("mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM%d", ism), Form("invariant mass vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {invmassAxis, collisionTimeAxis}, mDoSumw2.value); } if (mVetoBCID->length()) { @@ -277,6 +320,7 @@ struct EmcalPi0Qc { LOG(info) << "mRequireCaloReadout = " << mRequireCaloReadout.value; LOG(info) << "mRequireEMCalCells = " << mRequireEMCalCells.value; LOG(info) << "mSplitEMCalDCal = " << mSplitEMCalDCal.value; + LOG(info) << "mDoSumw2 = " << mDoSumw2.value; } template @@ -285,16 +329,16 @@ struct EmcalPi0Qc { static constexpr std::string_view ClusterTimeHistSM[20] = {"clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM0", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM1", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM2", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM3", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM4", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM5", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM6", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM7", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM8", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM9", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM10", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM11", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM12", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM13", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM14", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM15", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM16", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM17", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM18", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM19"}; static constexpr std::string_view ClusterNcellHistSM[20] = {"clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM0", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM1", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM2", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM3", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM4", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM5", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM6", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM7", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM8", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM9", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM10", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM11", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM12", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM13", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM14", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM15", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM16", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM17", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM18", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM19"}; static constexpr std::string_view ClusterM02HistSM[20] = {"clusterM02VsTimeStamp/clusterM02VsTimeStampSM0", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM1", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM2", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM3", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM4", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM5", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM6", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM7", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM8", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM9", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM10", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM11", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM12", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM13", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM14", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM15", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM16", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM17", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM18", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM19"}; - mHistManager.fill(HIST(ClusterTimeHistSM[supermoduleID]), time, timeSinceSOR); - mHistManager.fill(HIST(ClusterNcellHistSM[supermoduleID]), NCell, timeSinceSOR); - mHistManager.fill(HIST(ClusterM02HistSM[supermoduleID]), m02, timeSinceSOR); + mHistManager.fill(HIST(ClusterTimeHistSM[supermoduleID]), time, timeSinceSOR, mWeight); + mHistManager.fill(HIST(ClusterNcellHistSM[supermoduleID]), NCell, timeSinceSOR, mWeight); + mHistManager.fill(HIST(ClusterM02HistSM[supermoduleID]), m02, timeSinceSOR, mWeight); } template void supermoduleHistHelperMeson(float minv, float timeSinceSOR) { static constexpr std::string_view MesonInvMassHistSM[20] = {"mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM0", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM1", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM2", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM3", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM4", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM5", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM6", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM7", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM8", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM9", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM10", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM11", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM12", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM13", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM14", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM15", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM16", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM17", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM18", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM19"}; - mHistManager.fill(HIST(MesonInvMassHistSM[supermoduleID]), minv, timeSinceSOR); + mHistManager.fill(HIST(MesonInvMassHistSM[supermoduleID]), minv, timeSinceSOR, mWeight); } void fillSupermoduleHistogramsPhoton(int supermoduleID, float time, float m02, int NCell, float timeSinceSOR) @@ -455,13 +499,50 @@ struct EmcalPi0Qc { } for (const auto& collision : collisions) { + mHistManager.fill(HIST("events"), 1); // Fill "All events" bin of event histogram + // emcal hardware triggers + if (collision.alias_bit(kTVXinEMC)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 1); + } + if (collision.alias_bit(kEMC7)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 2); + } + if (collision.alias_bit(kDMC7)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 3); + } + if (collision.alias_bit(kEG1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 4); + } + if (collision.alias_bit(kEG2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 5); + } + if (collision.alias_bit(kDG1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 6); + } + if (collision.alias_bit(kDG2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 7); + } + if (collision.alias_bit(kEJ1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 8); + } + if (collision.alias_bit(kEJ2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 9); + } + if (collision.alias_bit(kDJ1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 10); + } + if (collision.alias_bit(kDJ2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 11); + } + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 12); + if (mDoEventSel.value && (!collision.sel8())) { // Check sel8 continue; } - mHistManager.fill(HIST("events"), 2); // Fill sel8 - if (mRequireCaloReadout.value && !collision.alias_bit(kTVXinEMC)) { // Check whether EMC was read out + mHistManager.fill(HIST("events"), 2); + if (mRequireCaloReadout.value && !(collision.alias_bit(kTVXinEMC) || collision.alias_bit(kEMC7) || collision.alias_bit(kDMC7) || collision.alias_bit(kEG1) || collision.alias_bit(kEG2) || collision.alias_bit(kDG1) || collision.alias_bit(kDG2) || collision.alias_bit(kEJ1) || collision.alias_bit(kEJ2) || collision.alias_bit(kDJ1) || collision.alias_bit(kDJ2))) { // Check whether EMC was read out continue; } mHistManager.fill(HIST("events"), 3); // Fill readout @@ -521,6 +602,138 @@ struct EmcalPi0Qc { } PROCESS_SWITCH(EmcalPi0Qc, processCollision, "Process clusters from collision", false); + /// \brief Process EMCAL clusters from MC-tagged collisions and set event weight + void processCollisionMC(MyBCs const& bcs, MyMCCollisions const& collisions, o2::aod::McCollisions const&, o2::aod::EMCALClusters const& clusters, o2::soa::Filtered const& cells, o2::aod::EMCALClusterCells const& clusterCells) + { + auto cellIter = cells.begin(); + auto bcIter = bcs.begin(); + int runNumber = bcIter.runNumber(); + std::unordered_map cellGlobalBCs; + // Build map of number of cells for corrected BCs using global BCs + // used later in the determination whether a BC has EMC cell content (for speed reason) + for (const auto& cell : cells) { + cellGlobalBCs[cell.bc_as().globalBC()]++; + } + + for (const auto& collision : collisions) { + mWeight = 1.0f; + // emcal hardware triggers + if (collision.alias_bit(kTVXinEMC)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 1); + } + if (collision.alias_bit(kEMC7)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 2); + } + if (collision.alias_bit(kDMC7)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 3); + } + if (collision.alias_bit(kEG1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 4); + } + if (collision.alias_bit(kEG2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 5); + } + if (collision.alias_bit(kDG1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 6); + } + if (collision.alias_bit(kDG2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 7); + } + if (collision.alias_bit(kEJ1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 8); + } + if (collision.alias_bit(kEJ2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 9); + } + if (collision.alias_bit(kDJ1)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 10); + } + if (collision.alias_bit(kDJ2)) { + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 11); + } + mHistManager.fill(HIST("eventsEMCALHardwareTriggers"), 12); + + if (collision.has_mcCollision()) { + mWeight = collision.mcCollision().weight(); + if (collision.mcCollision().getSubGeneratorId() == SubGeneratorId::mbGap) { + mHistManager.fill(HIST("signalGapEvents"), 1); // Fill gap events bin of signalGapEvents histogram + } else { + mHistManager.fill(HIST("signalGapEvents"), 2); // Fill signal events bin of signalGapEvents histogram + } + } + + mHistManager.fill(HIST("events"), 1, mWeight); // Fill "All events" bin of event histogram + mHistManager.fill(HIST("eventsWithoutWeight"), 1); // Fill "All events" bin of event histogram without weight + if (collision.has_mcCollision()) { + mHistManager.fill(HIST("eventsWithoutWeight"), 2); // Fill "Has MC collision" bin of event histogram without weight + } + if (mDoEventSel.value && (!collision.sel8())) { // Check sel8 + continue; + } + + mHistManager.fill(HIST("events"), 2, mWeight); // Fill sel8 + mHistManager.fill(HIST("eventsWithoutWeight"), 3); // Fill sel8 bin of event histogram without weight + if (mRequireCaloReadout.value && !(collision.alias_bit(kTVXinEMC) || collision.alias_bit(kEMC7) || collision.alias_bit(kDMC7) || collision.alias_bit(kEG1) || collision.alias_bit(kEG2) || collision.alias_bit(kDG1) || collision.alias_bit(kDG2) || collision.alias_bit(kEJ1) || collision.alias_bit(kEJ2) || collision.alias_bit(kDJ1) || collision.alias_bit(kDJ2))) { // Check whether EMC was read out + continue; + } + mHistManager.fill(HIST("events"), 3, mWeight); // Fill readout + mHistManager.fill(HIST("eventsWithoutWeight"), 4); // Fill readout bin of event histogram without weight + if (mDoEventSel.value && collision.numContrib() < 0.5) { // Skip collisions without contributors + continue; + } + mHistManager.fill(HIST("events"), 4, mWeight); // Fill >1 vtx contr. bin of event histogram + mHistManager.fill(HIST("eventsWithoutWeight"), 5); // Fill >1 vtx contr. bin of event histogram without weight + mHistManager.fill(HIST("eventVertexZAll"), collision.posZ(), mWeight); + if (mVertexCut > 0 && std::abs(collision.posZ()) > mVertexCut) { + continue; + } + mHistManager.fill(HIST("events"), 5, mWeight); // Fill z-Vertex selected bin of event histogram + mHistManager.fill(HIST("eventsWithoutWeight"), 6); // Fill z-Vertex selected bin of event histogram without weight + mHistManager.fill(HIST("eventVertexZSelected"), collision.posZ(), mWeight); + + if (mDoEventSel.value && collision.ambiguous()) { // Skip ambiguous collisions (those that are in BCs including multiple collisions) + continue; + } + mHistManager.fill(HIST("events"), 6, mWeight); // Fill "One collision in BC" bin of event histogram + mHistManager.fill(HIST("eventsWithoutWeight"), 7); // Fill "One collision in BC" bin of event histogram without weight + if (mDoEventSel.value) { + auto found = cellGlobalBCs.find(collision.foundBC_as().globalBC()); + if (mRequireEMCalCells.value && (found == cellGlobalBCs.end() || found->second == 0)) { // Skip collisions without any readout EMCal cells + continue; + } + } + mHistManager.fill(HIST("events"), 7, mWeight); // Fill at least one non0 cell in EMCal of event histogram (Selected) + mHistManager.fill(HIST("eventsWithoutWeight"), 8); // Fill at least one non0 cell in EMCal of event histogram (Selected) without weight + // Get BC and run number + int64_t foundBCId = collision.foundBCId(); + if (foundBCId >= 0) { + bcIter.setCursor(foundBCId); + } + runNumber = bcIter.runNumber(); + + // Fetch SOR only when run changes + if (runNumber != lastRunNumber) { + std::map headers, metadata; + headers = ccdbApi.retrieveHeaders(Form("RCT/Info/RunInformation/%i", runNumber), metadata, -1); + tsSOR = atol(headers["SOR"].c_str()); + // LOGP(info, "Run {} | SOR = {} ms", runNumber, tsSOR); + lastRunNumber = runNumber; + } + + // Time since SOR in minutes (bc.timestamp() is in ms) + float timeSinceSORMin = (bcIter.timestamp() - tsSOR) / 1000.0f / 60.f; + mHistManager.fill(HIST("hEventPerTime"), timeSinceSORMin, mWeight); + + auto clustersPerColl = clusters.sliceBy(perCollision, collision.globalIndex()); + if (clustersPerColl.size() == 0) { + continue; + } + processClusters(clustersPerColl, clusterCells, cellIter, timeSinceSORMin); + processMesons(timeSinceSORMin); + } + } + PROCESS_SWITCH(EmcalPi0Qc, processCollisionMC, "Process clusters from MC-tagged collisions", false); + /// \brief Process EMCAL clusters that are not matched to a collision /// This is not needed for most users void processAmbiguous(o2::aod::BCs::iterator const& bc, o2::aod::EMCALAmbiguousClusters const& clusters) @@ -614,15 +827,15 @@ struct EmcalPi0Qc { static constexpr std::string_view ClusterQAHistNLM[2] = {"ClustersBeforeCuts/clusterNLM", "ClustersAfterCuts/clusterNLM"}; static constexpr std::string_view ClusterQAHistNCells[2] = {"ClustersBeforeCuts/clusterNCells", "ClustersAfterCuts/clusterNCells"}; static constexpr std::string_view ClusterQAHistDistanceToBadChannel[2] = {"ClustersBeforeCuts/clusterDistanceToBadChannel", "ClustersAfterCuts/clusterDistanceToBadChannel"}; - mHistManager.fill(HIST(ClusterQAHistEnergy[BeforeCuts]), cluster.energy()); - mHistManager.fill(HIST(ClusterQAHistEnergySimpleBinning[BeforeCuts]), cluster.energy()); - mHistManager.fill(HIST(ClusterQAHistTime[BeforeCuts]), cluster.time()); - mHistManager.fill(HIST(ClusterQAHistEtaPhi[BeforeCuts]), cluster.eta(), cluster.phi()); - mHistManager.fill(HIST(ClusterQAHistM02[BeforeCuts]), cluster.m02()); - mHistManager.fill(HIST(ClusterQAHistM20[BeforeCuts]), cluster.m20()); - mHistManager.fill(HIST(ClusterQAHistNLM[BeforeCuts]), cluster.nlm()); - mHistManager.fill(HIST(ClusterQAHistNCells[BeforeCuts]), cluster.nCells()); - mHistManager.fill(HIST(ClusterQAHistDistanceToBadChannel[BeforeCuts]), cluster.distanceToBadChannel()); + mHistManager.fill(HIST(ClusterQAHistEnergy[BeforeCuts]), cluster.energy(), mWeight); + mHistManager.fill(HIST(ClusterQAHistEnergySimpleBinning[BeforeCuts]), cluster.energy(), mWeight); + mHistManager.fill(HIST(ClusterQAHistTime[BeforeCuts]), cluster.time(), mWeight); + mHistManager.fill(HIST(ClusterQAHistEtaPhi[BeforeCuts]), cluster.eta(), cluster.phi(), mWeight); + mHistManager.fill(HIST(ClusterQAHistM02[BeforeCuts]), cluster.m02(), mWeight); + mHistManager.fill(HIST(ClusterQAHistM20[BeforeCuts]), cluster.m20(), mWeight); + mHistManager.fill(HIST(ClusterQAHistNLM[BeforeCuts]), cluster.nlm(), mWeight); + mHistManager.fill(HIST(ClusterQAHistNCells[BeforeCuts]), cluster.nCells(), mWeight); + mHistManager.fill(HIST(ClusterQAHistDistanceToBadChannel[BeforeCuts]), cluster.distanceToBadChannel(), mWeight); } /// \brief Return a boolean that states, whether a cluster should be rejected by the applied cluster cuts @@ -668,7 +881,7 @@ struct EmcalPi0Qc { // build meson from photons Meson meson(mPhotons[ig1], mPhotons[ig2]); if (meson.getOpeningAngle() > mMinOpenAngleCut) { - mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt()); + mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt(), mWeight); uint8_t sm1 = mPhotons[ig1].sm; uint8_t sm2 = mPhotons[ig2].sm; @@ -678,9 +891,9 @@ struct EmcalPi0Qc { if (mSplitEMCalDCal) { if (!mPhotons[ig1].onDCal && !mPhotons[ig2].onDCal) { - mHistManager.fill(HIST("invMassVsPt_EMCal"), meson.getMass(), meson.getPt()); + mHistManager.fill(HIST("invMassVsPt_EMCal"), meson.getMass(), meson.getPt(), mWeight); } else if (mPhotons[ig1].onDCal && mPhotons[ig2].onDCal) { - mHistManager.fill(HIST("invMassVsPt_DCal"), meson.getMass(), meson.getPt()); + mHistManager.fill(HIST("invMassVsPt_DCal"), meson.getMass(), meson.getPt(), mWeight); } } } @@ -688,6 +901,7 @@ struct EmcalPi0Qc { // calculate background candidates (rotation background) calculateBackground(meson, ig1, ig2); } + // TODO: for now this part makes no sense when running JJ MC due to weights calculateMixedBack(mPhotons[ig1]); } @@ -736,22 +950,22 @@ struct EmcalPi0Qc { // Fill histograms if (mesonRotated1.getOpeningAngle() > mMinOpenAngleCut) { - mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt(), mWeight); if (mSplitEMCalDCal) { if (!mPhotons[ig1].onDCal && !mPhotons[ig2].onDCal && !mPhotons[ig3].onDCal) { - mHistManager.fill(HIST("invMassVsPtBackground_EMCal"), mesonRotated1.getMass(), mesonRotated1.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground_EMCal"), mesonRotated1.getMass(), mesonRotated1.getPt(), mWeight); } else if (mPhotons[ig1].onDCal && mPhotons[ig2].onDCal && mPhotons[ig3].onDCal) { - mHistManager.fill(HIST("invMassVsPtBackground_DCal"), mesonRotated1.getMass(), mesonRotated1.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground_DCal"), mesonRotated1.getMass(), mesonRotated1.getPt(), mWeight); } } } if (mesonRotated2.getOpeningAngle() > mMinOpenAngleCut) { - mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt(), mWeight); if (mSplitEMCalDCal) { if (!mPhotons[ig1].onDCal && !mPhotons[ig2].onDCal && !mPhotons[ig3].onDCal) { - mHistManager.fill(HIST("invMassVsPtBackground_EMCal"), mesonRotated2.getMass(), mesonRotated2.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground_EMCal"), mesonRotated2.getMass(), mesonRotated2.getPt(), mWeight); } else if (mPhotons[ig1].onDCal && mPhotons[ig2].onDCal && mPhotons[ig3].onDCal) { - mHistManager.fill(HIST("invMassVsPtBackground_DCal"), mesonRotated2.getMass(), mesonRotated2.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground_DCal"), mesonRotated2.getMass(), mesonRotated2.getPt(), mWeight); } } } @@ -803,7 +1017,7 @@ struct EmcalPi0Qc { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{ - adaptAnalysisTask(cfgc, TaskName{"EmcalPi0QcAssociate"}, SetDefaultProcesses{{{"processCollision", true}, {"processAmbiguous", false}}}), // o2-linter: disable=name/o2-task (adapted multiple times) - adaptAnalysisTask(cfgc, TaskName{"EmcalPi0QcAmbiguous"}, SetDefaultProcesses{{{"processCollision", false}, {"processAmbiguous", true}}})}; // o2-linter: disable=name/o2-task (adapted multiple times) + adaptAnalysisTask(cfgc, TaskName{"EmcalPi0QcAssociate"}, SetDefaultProcesses{{{"processCollision", true}, {"processCollisionMC", false}, {"processAmbiguous", false}}}), // o2-linter: disable=name/o2-task (adapted multiple times) + adaptAnalysisTask(cfgc, TaskName{"EmcalPi0QcAmbiguous"}, SetDefaultProcesses{{{"processCollision", false}, {"processCollisionMC", false}, {"processAmbiguous", true}}})}; // o2-linter: disable=name/o2-task (adapted multiple times) return workflow; } diff --git a/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx b/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx index aea4ee26bea..e2d209966dc 100644 --- a/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx +++ b/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx @@ -8,10 +8,10 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// \file taskSingleMuonSource.cxx -// \brief Task used to seperate single muons source in Monte Carlo simulation. -// \author Maolin Zhang , CCNU +/// +/// \file taskSingleMuonSource.cxx +/// \brief Task used to seperate single muons source in Monte Carlo simulation. +/// \author Maolin Zhang , CCNU #include "Common/Core/RecoDecay.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -24,11 +24,11 @@ #include #include #include +#include #include #include #include -#include #include #include @@ -78,6 +78,8 @@ struct HfTaskSingleMuonSource { Configurable charge{"charge", 0, "Muon track charge, validated values are 0, 1 and -1, 0 represents both 1 and -1"}; Configurable pairSource{"pairSource", true, "check also the source of like-sign muon pairs"}; + Service pdgDB; + double pDcaMax = 594.0; // p*DCA maximum value for small Rab double pDcaMax2 = 324.0; // p*DCA maximum value for large Rabs double rAbsMid = 26.5; // R at absorber end minimum value @@ -116,7 +118,7 @@ struct HfTaskSingleMuonSource { HistogramConfigSpec const h1ColNumber{HistType::kTH1F, {axisColNumber}}; HistogramConfigSpec const h1Pt{HistType::kTH1F, {axisPt}}; - HistogramConfigSpec h1Mass{HistType::kTH1F, {axisMass}}; + HistogramConfigSpec const h1Mass{HistType::kTH1F, {axisMass}}; HistogramConfigSpec const h2PtDCA{HistType::kTH2F, {axisPt, axisDCA}}; HistogramConfigSpec const h2PtChi2{HistType::kTH2F, {axisPt, axisChi2}}; HistogramConfigSpec const h2PtDeltaPt{HistType::kTH2F, {axisPt, axisDeltaPt}}; @@ -125,9 +127,11 @@ struct HfTaskSingleMuonSource { registry.add("h1MuBeforeCuts", "", h1Pt); registry.add("h1MuonMass", "", h1Mass); registry.add("h1BeautyMass", "", h1Mass); + registry.add("h1CorrBeautyMass", "", h1Mass); registry.add("h1OtherMass", "", h1Mass); registry.add("h1MuonMassGen", "", h1Mass); registry.add("h1BeautyMassGen", "", h1Mass); + registry.add("h1CorrBeautyMassGen", "", h1Mass); registry.add("h1OtherMassGen", "", h1Mass); for (const auto& src : muonSources) { registry.add(Form("h1%sPt", src.Data()), "", h1Pt); @@ -141,6 +145,8 @@ struct HfTaskSingleMuonSource { uint8_t getMask(const McMuons::iterator& muon) { uint8_t mask(0); + const int diquarkEdge = 1000; + const int hadronEdge = 10000; if (muon.has_mcParticle()) { SETBIT(mask, IsIdentified); } else { @@ -159,7 +165,7 @@ struct HfTaskSingleMuonSource { mcPart = *(mcPart.mothers_first_as()); const auto pdgAbs(std::abs(mcPart.pdgCode())); - if (pdgAbs < 10 || pdgAbs == 21) { + if (pdgAbs < kElectron || pdgAbs == kGluon) { break; // Quark and gluon } @@ -168,40 +174,39 @@ struct HfTaskSingleMuonSource { continue; } - if (pdgAbs == kTauMinus) { - // Tau + if (pdgAbs == kTauMinus) { // Tau SETBIT(mask, HasTauParent); continue; } const int pdgRem(pdgAbs % 100000); + const int pdgRemRem(pdgRem % 100); if (pdgRem == kProton) { continue; } // Beam particle - if ((pdgRem < 100) || (pdgRem >= 10000)) { + if ((pdgRem < kPi0) || (pdgRem >= hadronEdge)) { continue; } - if ((pdgRem % 100 == 1 || pdgRem % 100 == 3) && pdgRem > 1000) { // diquarks + if ((pdgRemRem == kDown || pdgRemRem == kStrange) && pdgRem > diquarkEdge) { // diquarks continue; } // compute the flavor of constituent quark const int flv(pdgRem / std::pow(10, static_cast(std::log10(pdgRem)))); - if (flv > 6) { + if (flv > kTop) { // no more than 6 flavors continue; } - if (flv < 4) { + if (flv < kCharm) { // light flavor SETBIT(mask, HasLightParent); continue; } - - auto* pdgData(TDatabasePDG::Instance()->GetParticle(mcPart.pdgCode())); + auto* pdgData = pdgDB->GetParticle(mcPart.pdgCode()); if ((pdgData != nullptr) && (pdgData->AntiParticle() == nullptr)) { SETBIT(mask, HasQuarkoniumParent); - } else if (flv == 4) { + } else if (flv == kCharm) { SETBIT(mask, HasCharmParent); } else { SETBIT(mask, HasBeautyParent); @@ -274,11 +279,13 @@ struct HfTaskSingleMuonSource { // fill the histograms of each particle types void fillHistograms(const McMuons::iterator& muon) { + const int type0 = 0; + const int type2 = 2; const auto mask(getMask(muon)); const auto pt(muon.pt()), chi2(muon.chi2MatchMCHMFT()); const auto dca(RecoDecay::sqrtSumOfSquares(muon.fwdDcaX(), muon.fwdDcaY())); - if (trackType == 0 || trackType == 2) { + if (trackType == type0 || trackType == type2) { if (!muon.has_matchMCHTrack()) { return; } @@ -344,6 +351,8 @@ struct HfTaskSingleMuonSource { int traceAncestor(const McMuons::iterator& muon, aod::McParticles const& mctracks) { int mcNum = 0; + const int hadronStatus = 80; + const int diquarkEdge = 1000; if (!muon.has_mcParticle()) { return 0; } @@ -353,36 +362,38 @@ struct HfTaskSingleMuonSource { } while (mcPart.has_mothers()) { // the first hadron after hadronization auto mother = mcPart.mothers_first_as(); - if (std::abs(mother.getGenStatusCode()) < 80) { + if (std::abs(mother.getGenStatusCode()) < hadronStatus) { break; } mcPart = mother; } int flv = mcPart.pdgCode() / std::pow(10, static_cast(std::log10(std::abs(mcPart.pdgCode())))); - if (abs(flv) == 5 && mcPart.pdgCode() < 1000) + if (std::abs(flv) == kBottom && mcPart.pdgCode() < diquarkEdge) { flv = -flv; + } for (int i = (mcPart.mothers_first_as()).globalIndex(); i <= (mcPart.mothers_last_as()).globalIndex(); i++) { // loop over the lund string - for (auto mctrack : mctracks) { + for (const auto& mctrack : mctracks) { if (mctrack.globalIndex() != i) { continue; } - if ((mctrack.pdgCode() != flv) && (abs(mctrack.pdgCode()) < abs(flv) * 1000)) { + if ((mctrack.pdgCode() != flv) && (std::abs(mctrack.pdgCode()) < std::abs(flv) * 1000)) { continue; } - while (mctrack.has_mothers()) { - int motherflv = (mctrack.mothers_first_as()).pdgCode() / std::pow(10, static_cast(std::log10(abs((mctrack.mothers_first_as()).pdgCode())))); // find the mother with same flavor - auto mother = (abs(motherflv) == abs(flv)) ? (mctrack.mothers_first_as()) : (mctrack.mothers_last_as()); - if ((mother.pdgCode() != mctrack.pdgCode()) && (abs(mctrack.pdgCode()) < 10)) { // both mother is not the the quark with same flavor - mcNum = mctrack.globalIndex(); + auto currentTrk = mctrack; + while (currentTrk.has_mothers()) { + int motherflv = (currentTrk.mothers_first_as()).pdgCode() / std::pow(10, static_cast(std::log10(std::abs((currentTrk.mothers_first_as()).pdgCode())))); // find the mother with same flavor + auto mother = (std::abs(motherflv) == std::abs(flv)) ? (currentTrk.mothers_first_as()) : (currentTrk.mothers_last_as()); + if ((mother.pdgCode() != currentTrk.pdgCode()) && (std::abs(currentTrk.pdgCode()) < kElectron)) { // both mother is not the the quark with same flavor + mcNum = currentTrk.globalIndex(); return mcNum; } - mctrack = mother; + currentTrk = mother; } } } return 0; } - bool Corr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts) + bool isCorr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts) { int moth11(0), moth12(0), moth21(1), moth22(1); @@ -391,7 +402,7 @@ struct HfTaskSingleMuonSource { if (anc1 == 0 || anc2 == 0) { return false; } - for (auto mcPart : mcParts) { + for (const auto& mcPart : mcParts) { if (mcPart.globalIndex() == anc1) { moth11 = (mcPart.mothers_first_as()).globalIndex(); moth12 = (mcPart.mothers_last_as()).globalIndex(); @@ -408,7 +419,8 @@ struct HfTaskSingleMuonSource { } void fillPairs(const McMuons::iterator& muon, const McMuons::iterator& muon2, aod::McParticles const& mcParts) { - if (trackType != 3) { + const int type3 = 3; + if (trackType != type3) { return; } float mm = o2::constants::physics::MassMuon; @@ -419,7 +431,7 @@ struct HfTaskSingleMuonSource { ROOT::Math::PtEtaPhiMVector mu1Vec(muon.pt(), muon.eta(), muon.phi(), mm); ROOT::Math::PtEtaPhiMVector mu2Vec(muon2.pt(), muon2.eta(), muon2.phi(), mm); ROOT::Math::PtEtaPhiMVector dimuVec = mu1Vec + mu2Vec; - auto InvM = dimuVec.M(); + auto invMass = dimuVec.M(); if (!muon.has_mcParticle() || !muon2.has_mcParticle()) { return; @@ -430,18 +442,22 @@ struct HfTaskSingleMuonSource { ROOT::Math::PtEtaPhiMVector mu1VecGen(mcPart1.pt(), mcPart1.eta(), mcPart1.phi(), mm); ROOT::Math::PtEtaPhiMVector mu2VecGen(mcPart2.pt(), mcPart2.eta(), mcPart2.phi(), mm); ROOT::Math::PtEtaPhiMVector dimuVecGen = mu1VecGen + mu2VecGen; - auto InvMGen = dimuVecGen.M(); + auto invMassGen = dimuVecGen.M(); if (isMuon(mask1) && isMuon(mask2)) { - registry.fill(HIST("h1MuonMass"), InvM); - registry.fill(HIST("h1MuonMassGen"), InvMGen); + registry.fill(HIST("h1MuonMass"), invMass); + registry.fill(HIST("h1MuonMassGen"), invMassGen); } - if (Corr(muon, muon2, mcParts) && isBeautyMu(mask1) && isBeautyMu(mask2)) { - registry.fill(HIST("h1BeautyMass"), InvM); - registry.fill(HIST("h1BeautyMassGen"), InvMGen); + if (isBeautyMu(mask1) && isBeautyMu(mask2)) { + registry.fill(HIST("h1BeautyMass"), invMass); + registry.fill(HIST("h1BeautyMassGen"), invMassGen); + if (isCorr(muon, muon2, mcParts)) { + registry.fill(HIST("h1CorrBeautyMass"), invMass); + registry.fill(HIST("h1CorrBeautyMassGen"), invMassGen); + } } else { - registry.fill(HIST("h1OtherMass"), InvM); - registry.fill(HIST("h1OtherMassGen"), InvMGen); + registry.fill(HIST("h1OtherMass"), invMass); + registry.fill(HIST("h1OtherMassGen"), invMassGen); } } @@ -477,7 +493,10 @@ struct HfTaskSingleMuonSource { continue; } } - if ((muon.chi2() >= 1e6) || (muon.chi2() < 0)) { + if (muon.chi2() < 0) { + continue; + } + if (muon.chi2MatchMCHMID() < 0) { continue; } if (charge != 0 && muon.sign() != charge) { @@ -517,10 +536,10 @@ struct HfTaskSingleMuonSource { } } - if ((muon2.chi2() >= 1e6) || (muon2.chi2() < 0)) { + if (muon2.chi2() < 0) { continue; } - if ((muon2.chi2MatchMCHMID() >= 1e6) || (muon2.chi2MatchMCHMID() < 0)) { + if (muon2.chi2MatchMCHMID() < 0) { continue; } fillPairs(muon, muon2, mcParts); diff --git a/PWGLF/Tasks/GlobalEventProperties/PseudorapidityDensityMFT.cxx b/PWGLF/Tasks/GlobalEventProperties/PseudorapidityDensityMFT.cxx index 370f7da75a9..fc7bc57cf19 100644 --- a/PWGLF/Tasks/GlobalEventProperties/PseudorapidityDensityMFT.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/PseudorapidityDensityMFT.cxx @@ -84,6 +84,9 @@ AxisSpec dcaYAxis = {6000, -30, 30}; // previous AxisSpec dcaYAxis = {2000, -10 // AxisSpec dcaZAxis = {600, -0.15f, 0.15f}; // AxisSpec dcaXAxis = {600, -0.15f, 0.15f}; // AxisSpec dcaYAxis = {600, -0.15f, 0.15f}; +// bin width 0.0005 cm: range [-30, 30] cm => 60/0.0005 = 120000 bins +// Keep bin width = 0.0005 cm (5 um): range [-1, 1] cm => 2.0/0.0005 = 4000 bins +// AxisSpec axisBinsDCA = {600, -0.15f, 0.15f, "#it{dca}_{xy} (cm)"}; AxisSpec centAxis = {{0, 10, 20, 30, 40, 50, 60, 70, 80, 100}}; @@ -113,16 +116,16 @@ static constexpr TrackSelectionFlags::flagtype TrackSelectionDca = // replace your alias with the extension included: using FullBCs = soa::Join; -// using MFTTracksLabeled = -// soa::Join; using MFTTracksLabeled = soa::Join; +using MFTTracksLabeled3d = + soa::Join; + using MFTTracksLabeled2d = - soa::Join; using MFTTracksLabeledOrg = @@ -139,7 +142,6 @@ struct PseudorapidityDensityMFT { Preslice perCol = o2::aod::fwdtrack::collisionId; Preslice perMcCol = aod::mcparticle::mcCollisionId; Preslice perColCentral = aod::track::collisionId; - Service pdg; // --- CCDB magnetic field (needed for propagateToDCAhelix in this device) --- @@ -236,6 +238,14 @@ struct PseudorapidityDensityMFT { Selected, Rejected }; + enum class NeitherReasonBin : int { + NotTrueByLabel = 1, + BestColInvalid, + BestColMissingInRecoToMc, + ClassifiedRight, + ClassifiedWrong + }; + static constexpr float ForwardEtaMax = -2.0f; static constexpr float ForwardEtaMin = -3.9f; @@ -430,7 +440,7 @@ struct PseudorapidityDensityMFT { " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaXyAxis}}}); - if (doprocessGenReco || doprocessGenRecoTimeCom) { + if (doprocessGenReco3d || doprocessGenReco2d || doprocessGenRecoTimeCom) { registry.add({"EventsRecoCuts_GenReco", ";cut;events", {HistType::kTH1F, {{16, 0.5, 16.5}}}}); @@ -763,6 +773,15 @@ struct PseudorapidityDensityMFT { hrw->GetXaxis()->SetBinLabel(static_cast(RightWrongBin::Neither), "neither"); hrw->GetXaxis()->SetBinLabel(static_cast(RightWrongBin::Both), "both"); + registry.add("Purity/NeitherReason", "Purity/NeitherReason", kTH1D, {{5, 0.5, 5.5}}, false); + auto hNeitherReason = registry.get(HIST("Purity/NeitherReason")); + auto* xNeitherReason = hNeitherReason->GetXaxis(); + xNeitherReason->SetBinLabel(static_cast(NeitherReasonBin::NotTrueByLabel), "NotTrueByLabel"); + xNeitherReason->SetBinLabel(static_cast(NeitherReasonBin::BestColInvalid), "BestColInvalid"); + xNeitherReason->SetBinLabel(static_cast(NeitherReasonBin::BestColMissingInRecoToMc), "BestColMissingInRecoToMc"); + xNeitherReason->SetBinLabel(static_cast(NeitherReasonBin::ClassifiedRight), "ClassifiedRight"); + xNeitherReason->SetBinLabel(static_cast(NeitherReasonBin::ClassifiedWrong), "ClassifiedWrong"); + registry.add({"Purity/RightWrongLater", ";category;counts", {HistType::kTH1F, {{4, 0.5, 4.5}}}}); @@ -973,6 +992,9 @@ struct PseudorapidityDensityMFT { registry.add({"Purity/BestRecoColNotFound", ";events;counts", {HistType::kTH1F, {{1, 0.5, 1.5}}}}); + registry.add({"Purity/TrueColNotFound", + ";events;counts", + {HistType::kTH1F, {{1, 0.5, 1.5}}}}); } if (doprocessGen) { @@ -1643,7 +1665,6 @@ struct PseudorapidityDensityMFT { if (midtracks.size() > 0 && retrack.ambDegree() == 1 && retrack.ambDegree() != 0) { uniqueCollisions.insert(collision.globalIndex()); } - if (retrack.ambDegree() != 0) { registry.fill(HIST("Tracks/Control/woOrp/woOrpTracksEtaZvtx"), track.eta(), z); @@ -2027,26 +2048,70 @@ struct PseudorapidityDensityMFT { registry.fill(HIST("EventsRecoCuts_GenReco"), static_cast(bin)); }; fillGenRecoCut(GenRecoCutBin::AllRecoCollisions); + std::unordered_map recoToMc; - std::unordered_map> mcToReco; // MC collision id -> list of reco collision globalIndex + std::unordered_map> mcToReco; + std::unordered_set acceptedRecoCols; + std::unordered_set recoCollisionIds; + std::unordered_set trueMCCollisionIds; - for (const auto& collision : collisions) { - int nSavedRows = 0; - std::unordered_set uniqueRecoColsSaved; - int recoCol = collision.globalIndex(); // reconstructed vertex index - int mcCol = collision.mcCollisionId(); // true MC collision index + recoToMc.reserve(collisions.size()); + mcToReco.reserve(collisions.size()); + acceptedRecoCols.reserve(collisions.size()); + recoCollisionIds.reserve(collisions.size()); + trueMCCollisionIds.reserve(collisions.size()); + + const auto countAndPassEvSelGenReco = [&](auto const& collision) { + struct EvSelStep { + bool enabled; + decltype(aod::evsel::kIsTriggerTVX) bit; + GenRecoCutBin bin; + }; - if (mcCol >= 0) { - recoToMc[recoCol] = mcCol; - mcToReco[mcCol].push_back(recoCol); + const std::array steps = {{ + {true, aod::evsel::kIsTriggerTVX, GenRecoCutBin::IsTriggerTVX}, + {true, aod::evsel::kNoTimeFrameBorder, GenRecoCutBin::NoTimeFrameBorder}, + {true, aod::evsel::kNoITSROFrameBorder, GenRecoCutBin::NoITSROFrameBorder}, + {useNoSameBunchPileup, aod::evsel::kNoSameBunchPileup, GenRecoCutBin::NoSameBunchPileup}, + {useGoodZvtxFT0vsPV, aod::evsel::kIsGoodZvtxFT0vsPV, GenRecoCutBin::GoodZvtxFT0vsPV}, + {useNoCollInRofStandard, aod::evsel::kNoCollInRofStandard, GenRecoCutBin::NoCollInRofStandard}, + {useNoCollInRofStrict, aod::evsel::kNoCollInRofStrict, GenRecoCutBin::NoCollInRofStrict}, + {useNoCollInTimeRangeStandard, aod::evsel::kNoCollInTimeRangeStandard, GenRecoCutBin::NoCollInTimeRangeStandard}, + {useNoCollInTimeRangeStrict, aod::evsel::kNoCollInTimeRangeStrict, GenRecoCutBin::NoCollInTimeRangeStrict}, + {useNoHighMultCollInPrevRof, aod::evsel::kNoHighMultCollInPrevRof, GenRecoCutBin::NoHighMultCollInPrevRof}, + }}; + + if (!useEvSel) { + for (const auto& step : steps) { + fillGenRecoCut(step.bin); + } + fillGenRecoCut(GenRecoCutBin::RctMFT); + return true; + } - ++nSavedRows; - uniqueRecoColsSaved.insert(recoCol); + for (const auto& step : steps) { + if (!step.enabled) { + fillGenRecoCut(step.bin); + continue; + } + + if (!collision.selection_bit(step.bit)) { + return false; + } + fillGenRecoCut(step.bin); } - registry.fill(HIST("Purity/HashTableRowCounts"), - static_cast(HashTableRowCountsBin::RowsSaved), nSavedRows); - registry.fill(HIST("Purity/HashTableRowCounts"), - static_cast(HashTableRowCountsBin::UniqueRecoColsSaved), uniqueRecoColsSaved.size()); + + if (useRctMFT && !myChecker(collision)) { + return false; + } + fillGenRecoCut(GenRecoCutBin::RctMFT); + + return true; + }; + + for (const auto& collision : collisions) { + int nSavedRows = 0; + std::unordered_set uniqueRecoColsSaved; registry.fill(HIST("Purity/reco/CollisionNumContrib"), collision.numContrib()); @@ -2056,59 +2121,11 @@ struct PseudorapidityDensityMFT { fillGenRecoCut(GenRecoCutBin::UseContBestCollisionIndex); if (!collision.has_mcCollision()) { - LOGF(warning, "No MC collision found..."); - return; + LOGP(warning, "Reco collision {} has no MC collision label, skipping", collision.globalIndex()); + continue; } fillGenRecoCut(GenRecoCutBin::HasMcCollision); - auto countAndPassEvSelGenReco = [&](auto const& collision) { - struct EvSelStep { - bool enabled; - decltype(aod::evsel::kIsTriggerTVX) bit; - GenRecoCutBin bin; - }; - - const std::array steps = {{ - {true, aod::evsel::kIsTriggerTVX, GenRecoCutBin::IsTriggerTVX}, - {true, aod::evsel::kNoTimeFrameBorder, GenRecoCutBin::NoTimeFrameBorder}, - {true, aod::evsel::kNoITSROFrameBorder, GenRecoCutBin::NoITSROFrameBorder}, - {useNoSameBunchPileup, aod::evsel::kNoSameBunchPileup, GenRecoCutBin::NoSameBunchPileup}, - {useGoodZvtxFT0vsPV, aod::evsel::kIsGoodZvtxFT0vsPV, GenRecoCutBin::GoodZvtxFT0vsPV}, - {useNoCollInRofStandard, aod::evsel::kNoCollInRofStandard, GenRecoCutBin::NoCollInRofStandard}, - {useNoCollInRofStrict, aod::evsel::kNoCollInRofStrict, GenRecoCutBin::NoCollInRofStrict}, - {useNoCollInTimeRangeStandard, aod::evsel::kNoCollInTimeRangeStandard, GenRecoCutBin::NoCollInTimeRangeStandard}, - {useNoCollInTimeRangeStrict, aod::evsel::kNoCollInTimeRangeStrict, GenRecoCutBin::NoCollInTimeRangeStrict}, - {useNoHighMultCollInPrevRof, aod::evsel::kNoHighMultCollInPrevRof, GenRecoCutBin::NoHighMultCollInPrevRof}, - }}; - - if (!useEvSel) { - for (const auto& step : steps) { - fillGenRecoCut(step.bin); - } - fillGenRecoCut(GenRecoCutBin::RctMFT); - return true; - } - - for (const auto& step : steps) { - if (!step.enabled) { - fillGenRecoCut(step.bin); - continue; - } - - if (!collision.selection_bit(step.bit)) { - return false; - } - fillGenRecoCut(step.bin); - } - - if (useRctMFT && !myChecker(collision)) { - return false; - } - fillGenRecoCut(GenRecoCutBin::RctMFT); - - return true; - }; - if (!countAndPassEvSelGenReco(collision)) { continue; } @@ -2125,146 +2142,163 @@ struct PseudorapidityDensityMFT { } fillGenRecoCut(GenRecoCutBin::InelGt0); - // constexpr uint8_t kFakeMcMask = 1u << 7; - for (const auto& track : tracks) { - float ndf = getTrackNdf(track); - const float chi2ndf = track.chi2() / ndf; - float phi = track.phi(); - // const float dcaXyCut = track.bestDCAXY(); - // const float dcaZCut = track.bestDCAZ(); - const float ptCut = track.pt(); - o2::math_utils::bringTo02Pi(phi); + const int recoCol = collision.globalIndex(); + const int mcCol = collision.mcCollisionId(); - const bool failTrackCuts = - track.nClusters() < cfgnCluster || - track.eta() <= cfgnEta1 || - track.eta() >= cfgnEta2 || - chi2ndf >= cfgChi2NDFMax || - phi <= cfgPhiCut1 || - phi >= cfgPhiCut2 || - (usePhiCut && - ((phi <= PhiVetoLow) || - ((phi >= PhiVetoPiMin) && (phi <= PhiVetoPiMax)) || - (phi >= PhiVetoHigh))) || - // (useDCAxyCut && dcaxyCut > maxDCAxy) || - // (useDCAzCut && std::abs(dcazCut) > maxDCAz) || - (usePtCut && ptCut > cfgnPt); + acceptedRecoCols.insert(recoCol); + recoCollisionIds.insert(recoCol); + trueMCCollisionIds.insert(mcCol); - if (failTrackCuts) { - continue; - } - const bool hasMcLabel = track.has_mcParticle(); - const bool isFakeByLabel = hasMcLabel ? (track.mcMask() != 0) : false; - const bool isTrueByLabel = hasMcLabel && !isFakeByLabel; - const bool hasNoMcLabel = !hasMcLabel; - const bool isPrimaryCharged = hasMcLabel && !isFakeByLabel && track.mcParticle().isPhysicalPrimary(); - const bool isSecondaryCharged = hasMcLabel && !isFakeByLabel && !track.mcParticle().isPhysicalPrimary(); - const float eta = track.eta(); - if (!passGenRecoTrackMode(track)) { // 0-> All nonorphans, 1->Non-Amb, 2->Amb - continue; - } - int bin = static_cast(RightWrongBin::Neither); - bool recoOfTrueExists = false; + if (mcCol >= 0) { + recoToMc[recoCol] = mcCol; + mcToReco[mcCol].push_back(recoCol); + ++nSavedRows; + uniqueRecoColsSaved.insert(recoCol); + } + + registry.fill(HIST("Purity/HashTableRowCounts"), + static_cast(HashTableRowCountsBin::RowsSaved), nSavedRows); + registry.fill(HIST("Purity/HashTableRowCounts"), + static_cast(HashTableRowCountsBin::UniqueRecoColsSaved), uniqueRecoColsSaved.size()); + } - if (isTrueByLabel) { - int recoCol = track.collisionId(); - auto itRecoToMc = recoToMc.find(recoCol); - const int mcOfTrack = track.mcParticle().mcCollisionId(); + for (const auto& track : tracks) { + float ndf = getTrackNdf(track); + const float chi2ndf = track.chi2() / ndf; + float phi = track.phi(); + const float ptCut = track.pt(); + o2::math_utils::bringTo02Pi(phi); + + const bool failTrackCuts = + track.nClusters() < cfgnCluster || + track.eta() <= cfgnEta1 || + track.eta() >= cfgnEta2 || + chi2ndf >= cfgChi2NDFMax || + phi <= cfgPhiCut1 || + phi >= cfgPhiCut2 || + (usePhiCut && + ((phi <= PhiVetoLow) || + ((phi >= PhiVetoPiMin) && (phi <= PhiVetoPiMax)) || + (phi >= PhiVetoHigh))) || + (usePtCut && ptCut > cfgnPt); + + if (failTrackCuts) { + continue; + } - // Check whether any reco vertex exists for the true MC collision of this track - for (const auto& [recoId, mcId] : recoToMc) { - if (mcId == mcOfTrack) { - recoOfTrueExists = true; - break; - } - } + if (!passGenRecoTrackMode(track)) { + continue; + } - if (recoCol >= 0 && itRecoToMc != recoToMc.end()) { - int mcFromReco = itRecoToMc->second; - bin = (mcFromReco == mcOfTrack) - ? static_cast(RightWrongBin::Right) - : static_cast(RightWrongBin::Wrong); - } - } + const int recoCol = track.collisionId(); + if (acceptedRecoCols.find(recoCol) == acceptedRecoCols.end()) { + continue; + } - registry.fill(HIST("RightWrong"), bin); + const bool hasMcLabel = track.has_mcParticle(); + const bool isFakeByLabel = hasMcLabel ? (track.mcMask() != 0) : false; + const bool isTrueByLabel = hasMcLabel && !isFakeByLabel; + const bool hasNoMcLabel = !hasMcLabel; + const bool isPrimaryCharged = hasMcLabel && !isFakeByLabel && track.mcParticle().isPhysicalPrimary(); + const bool isSecondaryCharged = hasMcLabel && !isFakeByLabel && !track.mcParticle().isPhysicalPrimary(); + const float eta = track.eta(); - if (bin == static_cast(RightWrongBin::Wrong)) { - registry.fill(HIST("Purity/WrongVertexRecoExists"), recoOfTrueExists ? static_cast(BoolBin::Yes) : static_cast(BoolBin::No)); - } + int bin = static_cast(RightWrongBin::Neither); + bool recoOfTrueExists = false; - const auto fillTrackLabelSummary = [&](TrackLabelSummaryBin bin) { - registry.fill(HIST("Purity/TrackLabelSummary"), static_cast(bin)); - }; - const auto fillTrackEtaCategory = [&](TrackLabelSummaryBin bin) { - constexpr float EtaSentinel = -999.f; + /// const int bestColID = track.bestCollisionId(); + const int mcOfTrack = isTrueByLabel ? track.mcParticle().mcCollisionId() : InvalidCollisionId; - float etaAll = EtaSentinel; - float etaNoMc = EtaSentinel; - float etaFake = EtaSentinel; - float etaTrue = EtaSentinel; - float etaPrimary = EtaSentinel; - float etaSecondary = EtaSentinel; + if (isTrueByLabel) { + auto itRecoList = mcToReco.find(mcOfTrack); + if (itRecoList != mcToReco.end() && !itRecoList->second.empty()) { + recoOfTrueExists = true; + } - switch (bin) { - case TrackLabelSummaryBin::AllTracks: - etaAll = eta; - break; - case TrackLabelSummaryBin::NoMcLabel: - etaNoMc = eta; - break; - case TrackLabelSummaryBin::FakeTracks: - etaFake = eta; - break; - case TrackLabelSummaryBin::TrueTracks: - etaTrue = eta; - break; - case TrackLabelSummaryBin::PrimaryTracks: - etaPrimary = eta; - break; - case TrackLabelSummaryBin::SecondaryTracks: - etaSecondary = eta; - break; - } + auto itRecoToMc = recoToMc.find(recoCol); + if (recoCol >= 0 && itRecoToMc != recoToMc.end()) { + const int mcFromReco = itRecoToMc->second; + bin = (mcFromReco == mcOfTrack) + ? static_cast(RightWrongBin::Right) + : static_cast(RightWrongBin::Wrong); + } + } - registry.fill(HIST("Purity/TrackEtaCategorySparse"), - etaAll, etaNoMc, etaFake, etaTrue, etaPrimary, etaSecondary); - }; + registry.fill(HIST("RightWrong"), bin); - // registry.fill(HIST("Purity/TrackLabelStatus"), hasMcLabel ? 1.0 : 0.0); - // registry.fill(HIST("Purity/TrackFakeStatus"), isFakeByLabel ? 1.0 : 0.0); + if (bin == static_cast(RightWrongBin::Wrong)) { + registry.fill(HIST("Purity/WrongVertexRecoExists"), + recoOfTrueExists ? static_cast(WrongVertexRecoExistsBin::RecoOfTrueExists) + : static_cast(WrongVertexRecoExistsBin::RecoOfTrueMissing)); + } - fillTrackLabelSummary(TrackLabelSummaryBin::AllTracks); - fillTrackEtaCategory(TrackLabelSummaryBin::AllTracks); + const auto fillTrackLabelSummary = [&](TrackLabelSummaryBin binSummary) { + registry.fill(HIST("Purity/TrackLabelSummary"), static_cast(binSummary)); + }; - if (hasNoMcLabel) { - fillTrackLabelSummary(TrackLabelSummaryBin::NoMcLabel); - fillTrackEtaCategory(TrackLabelSummaryBin::NoMcLabel); - continue; + const auto fillTrackEtaCategory = [&](TrackLabelSummaryBin binSummary) { + constexpr float EtaSentinel = -999.f; + + float etaAll = EtaSentinel; + float etaNoMc = EtaSentinel; + float etaFake = EtaSentinel; + float etaTrue = EtaSentinel; + float etaPrimary = EtaSentinel; + float etaSecondary = EtaSentinel; + + switch (binSummary) { + case TrackLabelSummaryBin::AllTracks: + etaAll = eta; + break; + case TrackLabelSummaryBin::NoMcLabel: + etaNoMc = eta; + break; + case TrackLabelSummaryBin::FakeTracks: + etaFake = eta; + break; + case TrackLabelSummaryBin::TrueTracks: + etaTrue = eta; + break; + case TrackLabelSummaryBin::PrimaryTracks: + etaPrimary = eta; + break; + case TrackLabelSummaryBin::SecondaryTracks: + etaSecondary = eta; + break; } - if (isFakeByLabel) { - fillTrackLabelSummary(TrackLabelSummaryBin::FakeTracks); - fillTrackEtaCategory(TrackLabelSummaryBin::FakeTracks); - continue; - } + registry.fill(HIST("Purity/TrackEtaCategorySparse"), + etaAll, etaNoMc, etaFake, etaTrue, etaPrimary, etaSecondary); + }; - fillTrackLabelSummary(TrackLabelSummaryBin::TrueTracks); - fillTrackEtaCategory(TrackLabelSummaryBin::TrueTracks); + fillTrackLabelSummary(TrackLabelSummaryBin::AllTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::AllTracks); - if (isPrimaryCharged) { - fillTrackLabelSummary(TrackLabelSummaryBin::PrimaryTracks); - fillTrackEtaCategory(TrackLabelSummaryBin::PrimaryTracks); - } + if (hasNoMcLabel) { + fillTrackLabelSummary(TrackLabelSummaryBin::NoMcLabel); + fillTrackEtaCategory(TrackLabelSummaryBin::NoMcLabel); + continue; + } - if (isSecondaryCharged) { - fillTrackLabelSummary(TrackLabelSummaryBin::SecondaryTracks); - fillTrackEtaCategory(TrackLabelSummaryBin::SecondaryTracks); - } - // registry.fill(HIST("RightWrong"), bin); + if (isFakeByLabel) { + fillTrackLabelSummary(TrackLabelSummaryBin::FakeTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::FakeTracks); + continue; + } - } // Track loop 1 - } // Collision + fillTrackLabelSummary(TrackLabelSummaryBin::TrueTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::TrueTracks); + + if (isPrimaryCharged) { + fillTrackLabelSummary(TrackLabelSummaryBin::PrimaryTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::PrimaryTracks); + } + + if (isSecondaryCharged) { + fillTrackLabelSummary(TrackLabelSummaryBin::SecondaryTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::SecondaryTracks); + } + } } PROCESS_SWITCH(PseudorapidityDensityMFT, processGenRecoTimeCom, "Process for MC time compatible", false); @@ -2277,75 +2311,99 @@ struct PseudorapidityDensityMFT { // aod::McCollisions::iterator const& mcCollision // McCollisionsWithExtra::iterator const& mcCollision - // void processMCeff(soa::Join::iterator const& mcCollision, CollisionMCRecTable const& RecCols, TrackMCTrueTable const& GenParticles, FilTrackMCRecTable const& RecTracks) - // soa::Join::iterator const& mcCollision //This worked + template void processGenReco(McCollisionsWithExtra::iterator const& mcCollision, o2::soa::SmallGroups> const& collisions, FullBCs const& bcs, - MFTTracksLabeled const& tracks, + MFTTracksT const& tracks, FiCentralTracks const& midtracks, aod::McParticles const&) - { - const auto fillGenRecoCut = [&](GenRecoCutBin bin) { registry.fill(HIST("EventsRecoCuts_GenReco"), static_cast(bin)); }; fillGenRecoCut(GenRecoCutBin::AllRecoCollisions); - // if (maxGenRecoEvents >= 0 && nProcessedGenReco >= maxGenRecoEvents) { - // return; - // } - // ++nProcessedGenReco; // for HIR std::unordered_map recoToMc; - std::unordered_map> mcToReco; // MC collision id -> list of reco collision globalIndex - std::unordered_map recoToMcVZ; + std::unordered_map> mcToReco; + std::unordered_set acceptedRecoCols; std::unordered_map recoVtxX; std::unordered_map recoVtxY; std::unordered_map recoVtxZ; + std::unordered_set recoCollisionIds; + std::unordered_set trueMCCollisionIds; std::unordered_map> recoVtxByRecoId; std::unordered_map> recoVtxByMcId; + recoToMc.reserve(collisions.size()); + mcToReco.reserve(collisions.size()); + acceptedRecoCols.reserve(collisions.size()); + recoCollisionIds.reserve(collisions.size()); + trueMCCollisionIds.reserve(collisions.size()); recoVtxByRecoId.reserve(collisions.size()); recoVtxByMcId.reserve(collisions.size()); - mcToReco.reserve(collisions.size()); - // --- Make sure magnetic field exists in THIS device before any propagation --- - // IMPORTANT: calling collision.bc_as() requires the BC table to be subscribed. - // We subscribe by taking `FullBCs const& bcs` in the process signature and init once here. bool magInited = false; for (auto const& bc : bcs) { initMagField(bc); magInited = true; - break; // once is enough (initMagField is internally guarded) + break; } if (!magInited) { LOGF(fatal, "BC table is empty: cannot initialize magnetic field"); } - //_______________________________________________________________________________ + const auto countAndPassEvSelGenReco = [&](auto const& collision) { + struct EvSelStep { + bool enabled; + decltype(aod::evsel::kIsTriggerTVX) bit; + GenRecoCutBin bin; + }; - for (const auto& collision : collisions) { - int nSavedRows = 0; - std::unordered_set uniqueRecoColsSaved; - int recoCol = collision.globalIndex(); // reconstructed vertex index - int mcCol = collision.mcCollisionId(); // true MC collision index + const std::array steps = {{ + {true, aod::evsel::kIsTriggerTVX, GenRecoCutBin::IsTriggerTVX}, + {true, aod::evsel::kNoTimeFrameBorder, GenRecoCutBin::NoTimeFrameBorder}, + {true, aod::evsel::kNoITSROFrameBorder, GenRecoCutBin::NoITSROFrameBorder}, + {useNoSameBunchPileup, aod::evsel::kNoSameBunchPileup, GenRecoCutBin::NoSameBunchPileup}, + {useGoodZvtxFT0vsPV, aod::evsel::kIsGoodZvtxFT0vsPV, GenRecoCutBin::GoodZvtxFT0vsPV}, + {useNoCollInRofStandard, aod::evsel::kNoCollInRofStandard, GenRecoCutBin::NoCollInRofStandard}, + {useNoCollInRofStrict, aod::evsel::kNoCollInRofStrict, GenRecoCutBin::NoCollInRofStrict}, + {useNoCollInTimeRangeStandard, aod::evsel::kNoCollInTimeRangeStandard, GenRecoCutBin::NoCollInTimeRangeStandard}, + {useNoCollInTimeRangeStrict, aod::evsel::kNoCollInTimeRangeStrict, GenRecoCutBin::NoCollInTimeRangeStrict}, + {useNoHighMultCollInPrevRof, aod::evsel::kNoHighMultCollInPrevRof, GenRecoCutBin::NoHighMultCollInPrevRof}, + }}; + + if (!useEvSel) { + for (const auto& step : steps) { + fillGenRecoCut(step.bin); + } + fillGenRecoCut(GenRecoCutBin::RctMFT); + return true; + } - if (mcCol >= 0) { - recoToMc[recoCol] = mcCol; - mcToReco[mcCol].push_back(recoCol); + for (const auto& step : steps) { + if (!step.enabled) { + fillGenRecoCut(step.bin); + continue; + } - ++nSavedRows; - uniqueRecoColsSaved.insert(recoCol); + if (!collision.selection_bit(step.bit)) { + return false; + } + fillGenRecoCut(step.bin); } - registry.fill(HIST("Purity/HashTableRowCounts"), - static_cast(HashTableRowCountsBin::RowsSaved), nSavedRows); - registry.fill(HIST("Purity/HashTableRowCounts"), - static_cast(HashTableRowCountsBin::UniqueRecoColsSaved), uniqueRecoColsSaved.size()); - recoVtxX[recoCol] = collision.posX(); - recoVtxY[recoCol] = collision.posY(); - recoVtxZ[recoCol] = collision.posZ(); + if (useRctMFT && !myChecker(collision)) { + return false; + } + fillGenRecoCut(GenRecoCutBin::RctMFT); + + return true; + }; + + for (const auto& collision : collisions) { + int nSavedRows = 0; + std::unordered_set uniqueRecoColsSaved; registry.fill(HIST("Purity/reco/CollisionNumContrib"), collision.numContrib()); @@ -2355,59 +2413,11 @@ struct PseudorapidityDensityMFT { fillGenRecoCut(GenRecoCutBin::UseContBestCollisionIndex); if (!collision.has_mcCollision()) { - LOGF(warning, "No MC collision found..."); - return; + LOGP(warning, "Reco collision {} has no MC collision label, skipping", collision.globalIndex()); + continue; } fillGenRecoCut(GenRecoCutBin::HasMcCollision); - auto countAndPassEvSelGenReco = [&](auto const& collision) { - struct EvSelStep { - bool enabled; - decltype(aod::evsel::kIsTriggerTVX) bit; - GenRecoCutBin bin; - }; - - const std::array steps = {{ - {true, aod::evsel::kIsTriggerTVX, GenRecoCutBin::IsTriggerTVX}, - {true, aod::evsel::kNoTimeFrameBorder, GenRecoCutBin::NoTimeFrameBorder}, - {true, aod::evsel::kNoITSROFrameBorder, GenRecoCutBin::NoITSROFrameBorder}, - {useNoSameBunchPileup, aod::evsel::kNoSameBunchPileup, GenRecoCutBin::NoSameBunchPileup}, - {useGoodZvtxFT0vsPV, aod::evsel::kIsGoodZvtxFT0vsPV, GenRecoCutBin::GoodZvtxFT0vsPV}, - {useNoCollInRofStandard, aod::evsel::kNoCollInRofStandard, GenRecoCutBin::NoCollInRofStandard}, - {useNoCollInRofStrict, aod::evsel::kNoCollInRofStrict, GenRecoCutBin::NoCollInRofStrict}, - {useNoCollInTimeRangeStandard, aod::evsel::kNoCollInTimeRangeStandard, GenRecoCutBin::NoCollInTimeRangeStandard}, - {useNoCollInTimeRangeStrict, aod::evsel::kNoCollInTimeRangeStrict, GenRecoCutBin::NoCollInTimeRangeStrict}, - {useNoHighMultCollInPrevRof, aod::evsel::kNoHighMultCollInPrevRof, GenRecoCutBin::NoHighMultCollInPrevRof}, - }}; - - if (!useEvSel) { - for (const auto& step : steps) { - fillGenRecoCut(step.bin); - } - fillGenRecoCut(GenRecoCutBin::RctMFT); - return true; - } - - for (const auto& step : steps) { - if (!step.enabled) { - fillGenRecoCut(step.bin); - continue; - } - - if (!collision.selection_bit(step.bit)) { - return false; - } - fillGenRecoCut(step.bin); - } - - if (useRctMFT && !myChecker(collision)) { - return false; - } - fillGenRecoCut(GenRecoCutBin::RctMFT); - - return true; - }; - if (!countAndPassEvSelGenReco(collision)) { continue; } @@ -2424,202 +2434,342 @@ struct PseudorapidityDensityMFT { } fillGenRecoCut(GenRecoCutBin::InelGt0); - //________________________________________________________________________________ + const int recoCol = collision.globalIndex(); + const int mcCol = collision.mcCollisionId(); + + acceptedRecoCols.insert(recoCol); + recoCollisionIds.insert(recoCol); + trueMCCollisionIds.insert(mcCol); + + if (mcCol >= 0) { + recoToMc[recoCol] = mcCol; + mcToReco[mcCol].push_back(recoCol); + ++nSavedRows; + uniqueRecoColsSaved.insert(recoCol); + } + + registry.fill(HIST("Purity/HashTableRowCounts"), + static_cast(HashTableRowCountsBin::RowsSaved), nSavedRows); + registry.fill(HIST("Purity/HashTableRowCounts"), + static_cast(HashTableRowCountsBin::UniqueRecoColsSaved), uniqueRecoColsSaved.size()); + + recoVtxX[recoCol] = collision.posX(); + recoVtxY[recoCol] = collision.posY(); + recoVtxZ[recoCol] = collision.posZ(); + recoVtxByRecoId[recoCol] = {collision.posX(), collision.posY(), collision.posZ()}; + recoVtxByMcId[mcCol] = {mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()}; registry.fill(HIST("Purity/xReco"), collision.posX()); registry.fill(HIST("Purity/xTrue"), mcCollision.posX()); registry.fill(HIST("Purity/yReco"), collision.posY()); - registry.fill(HIST("Purity/yTrue"), mcCollision.posY()); registry.fill(HIST("Purity/zReco"), collision.posZ()); registry.fill(HIST("Purity/zTrue"), mcCollision.posZ()); - - // --- Vertex position THnSparse: status axis (1=reco, 2=true) --- registry.fill(HIST("Purity/VtxXYZTruth"), mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); registry.fill(HIST("Purity/VtxXYZReco"), collision.posX(), collision.posY(), collision.posZ()); - - // --- Delta vertex THnSparse (reco - true) --- registry.fill(HIST("Purity/DeltaVtxXYZ"), - collision.posX() - mcCollision.posX(), // Reco - Truth + collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ()); + } - recoVtxByRecoId[collision.globalIndex()] = {collision.posX(), collision.posY(), collision.posZ()}; - recoVtxByMcId[collision.mcCollisionId()] = {mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()}; + int64_t woOrpCount = 0; + bool filledRight = false; + bool filledWrong = false; + int nMftSelectedAfterCuts = 0; + std::unordered_set uniqueBestRecoCols; - int64_t woOrpCount = 0; - bool filledRight = false; - bool filledWrong = false; - int nMftSelectedAfterCuts = 0; - std::unordered_set uniqueBestRecoCols; - // auto particle = templatedTrack.template mcParticle_as(); + if (tracks.size() > 0) { + bool countedPrimary = false; + for (const auto& track : tracks) { + float ndf = getTrackNdf(track); + float chi2ndf = track.chi2() / ndf; + float phi = track.phi(); + float dcaXyCut = track.bestDCAXY(); + float dcaZCut = 0.f; + bool failDCAzCut = false; + float ptCut = track.pt(); + constexpr bool hasBestDCAZ = requires { track.bestDCAZ(); }; + + if constexpr (hasBestDCAZ) { + dcaZCut = track.bestDCAZ(); + failDCAzCut = useDCAzCut && (std::abs(dcaZCut) > maxDCAz); + } - if (tracks.size() > 0) { - bool countedPrimary = false; - for (const auto& track : tracks) { // track loop starts - // All compatible collisions assigned by track-to-collision-associator - if (!(midtracks.size() > 0)) - continue; - // std::cout <<" midtracks.size() " <= cfgnEta2 || + chi2ndf >= cfgChi2NDFMax || + phi <= cfgPhiCut1 || + phi >= cfgPhiCut2 || + (usePhiCut && + ((phi <= PhiVetoLow) || + ((phi >= PhiVetoPiMin) && (phi <= PhiVetoPiMax)) || + (phi >= PhiVetoHigh))) || + (useDCAxyCut && dcaXyCut > maxDCAxy) || + failDCAzCut || + (usePtCut && ptCut > cfgnPt); - o2::math_utils::bringTo02Pi(phi); - const float etaReco = track.eta(); - const float dcaXYReco = dcaXyCut; // track.bestDCAXY() - const float dcaZReco = dcaZCut; // track.bestDCAZ() - const float dcaXReco = dcaXYReco * std::cos(phi); - const float dcaYReco = dcaXYReco * std::sin(phi); - - const bool failTrackCuts = - track.nClusters() < cfgnCluster || - etaReco <= cfgnEta1 || - etaReco >= cfgnEta2 || - chi2ndf >= cfgChi2NDFMax || - phi <= cfgPhiCut1 || - phi >= cfgPhiCut2 || - (usePhiCut && - ((phi <= PhiVetoLow) || - ((phi >= PhiVetoPiMin) && (phi <= PhiVetoPiMax)) || - (phi >= PhiVetoHigh))) || - (useDCAxyCut && dcaXyCut > maxDCAxy) || - (useDCAzCut && std::abs(dcaZCut) > maxDCAz) || - (usePtCut && ptCut > cfgnPt); + if (failTrackCuts) { + continue; + } - if (failTrackCuts) { - continue; + if (!passGenRecoTrackMode(track)) { + continue; + } + + const int recoCol = track.collisionId(); + if (acceptedRecoCols.find(recoCol) == acceptedRecoCols.end()) { + continue; + } + + auto itRecoVz = recoVtxZ.find(recoCol); + if (itRecoVz == recoVtxZ.end()) { + continue; + } + const float z = itRecoVz->second; + + const bool hasMcLabel = track.has_mcParticle(); + const bool isFakeByLabel = hasMcLabel ? (track.mcMask() != 0) : false; + const bool isTrueByLabel = hasMcLabel && !isFakeByLabel; + const bool hasNoMcLabel = !hasMcLabel; + const bool isPrimaryCharged = hasMcLabel && !isFakeByLabel && track.mcParticle().isPhysicalPrimary(); + const bool isSecondaryCharged = hasMcLabel && !isFakeByLabel && !track.mcParticle().isPhysicalPrimary(); + const auto fillTrackLabelSummary = [&](TrackLabelSummaryBin binSummary) { + registry.fill(HIST("Purity/TrackLabelSummary"), static_cast(binSummary)); + }; + + const auto fillTrackEtaCategory = [&](TrackLabelSummaryBin binSummary) { + constexpr float EtaSentinel = -999.f; + + float etaAll = EtaSentinel; + float etaNoMc = EtaSentinel; + float etaFake = EtaSentinel; + float etaTrue = EtaSentinel; + float etaPrimary = EtaSentinel; + float etaSecondary = EtaSentinel; + + switch (binSummary) { + case TrackLabelSummaryBin::AllTracks: + etaAll = etaReco; + break; + case TrackLabelSummaryBin::NoMcLabel: + etaNoMc = etaReco; + break; + case TrackLabelSummaryBin::FakeTracks: + etaFake = etaReco; + break; + case TrackLabelSummaryBin::TrueTracks: + etaTrue = etaReco; + break; + case TrackLabelSummaryBin::PrimaryTracks: + etaPrimary = etaReco; + break; + case TrackLabelSummaryBin::SecondaryTracks: + etaSecondary = etaReco; + break; } - const bool hasMcLabel = track.has_mcParticle(); - const bool isFakeByLabel = hasMcLabel ? (track.mcMask() != 0) : false; - const bool isTrueByLabel = hasMcLabel && !isFakeByLabel; - const bool isPrimaryCharged = hasMcLabel && !isFakeByLabel && track.mcParticle().isPhysicalPrimary(); - const bool isSecondaryCharged = hasMcLabel && !isFakeByLabel && !track.mcParticle().isPhysicalPrimary(); - const auto mcColObj = track.mcParticle().mcCollision_as(); - const auto mcPart = track.mcParticle(); + registry.fill(HIST("Purity/TrackEtaCategorySparse"), + etaAll, etaNoMc, etaFake, etaTrue, etaPrimary, etaSecondary); + }; - if (!passGenRecoTrackMode(track)) { // - continue; + fillTrackLabelSummary(TrackLabelSummaryBin::AllTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::AllTracks); + + if (hasMcLabel) { + const auto mcPartForMother = track.mcParticle(); + if (!isPrimaryCharged && mcPartForMother.has_mothers()) { + auto mcpartMother = mcPartForMother.template mothers_as().front(); + if (mcpartMother.pdgCode() == PDG_t::kK0Short || + std::abs(mcpartMother.pdgCode()) == PDG_t::kLambda0) { + registry.fill(HIST("Purity/reco/weakStrange/SelectedTracksEta"), track.eta()); + registry.fill(HIST("Purity/reco/weakStrange/SelectedTracksEtaZvtx"), track.eta(), z); + } } - int bin = static_cast(RightWrongBin::Neither); - bool recoOfTrueExists = false; - bool recoOfTrueInCompatible = false; + } - const int bestColID = track.bestCollisionId(); // same as track.collisionId(); - if (isTrueByLabel) { - auto itRecoToMc = recoToMc.find(bestColID); - const int mcOfTrack = track.mcParticle().mcCollisionId(); - const auto compatibleIds = track.compatibleCollIds(); - auto itRecoList = mcToReco.find(mcOfTrack); - - // 1) First check whether the correct reco collision of the true MC collision - // is present in the compatible-collision list. - if (!compatibleIds.empty() && itRecoList != mcToReco.end() && !itRecoList->second.empty()) { - for (const auto& trueRecoId : itRecoList->second) { - for (const auto& compatibleId : compatibleIds) { - if (compatibleId == trueRecoId) { - recoOfTrueInCompatible = true; - break; - } - } - if (recoOfTrueInCompatible) { + if (hasNoMcLabel) { + fillTrackLabelSummary(TrackLabelSummaryBin::NoMcLabel); + fillTrackEtaCategory(TrackLabelSummaryBin::NoMcLabel); + } else if (isFakeByLabel) { + fillTrackLabelSummary(TrackLabelSummaryBin::FakeTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::FakeTracks); + } else { + fillTrackLabelSummary(TrackLabelSummaryBin::TrueTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::TrueTracks); + + if (isPrimaryCharged) { + fillTrackLabelSummary(TrackLabelSummaryBin::PrimaryTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::PrimaryTracks); + } + + if (isSecondaryCharged) { + fillTrackLabelSummary(TrackLabelSummaryBin::SecondaryTracks); + fillTrackEtaCategory(TrackLabelSummaryBin::SecondaryTracks); + } + } + + int bin = static_cast(RightWrongBin::Neither); + bool recoOfTrueExists = false; + bool recoOfTrueInCompatible = false; + + const int bestColID = track.bestCollisionId(); + const int mcOfTrack = isTrueByLabel ? track.mcParticle().mcCollisionId() : InvalidCollisionId; + + const bool foundRecoColInRecoList = + recoCollisionIds.find(recoCol) != recoCollisionIds.end(); + const bool foundBestColInRecoList = + recoCollisionIds.find(bestColID) != recoCollisionIds.end(); + const bool foundInMCTrueList = + isTrueByLabel && (trueMCCollisionIds.find(mcOfTrack) != trueMCCollisionIds.end()); + + static constexpr int RecoColMissingBin = 1; + static constexpr int BestRecoColMissingBin = 2; + static constexpr int TrueColMissingBin = 1; + if (!foundRecoColInRecoList) { + registry.fill(HIST("Purity/BestRecoColNotFound"), RecoColMissingBin); + } + if (!foundBestColInRecoList) { + registry.fill(HIST("Purity/BestRecoColNotFound"), BestRecoColMissingBin); + } + if (isTrueByLabel && !foundInMCTrueList) { + registry.fill(HIST("Purity/TrueColNotFound"), TrueColMissingBin); + } + + if (!isTrueByLabel) { + registry.fill(HIST("Purity/NeitherReason"), + static_cast(NeitherReasonBin::NotTrueByLabel)); + } else { + auto itRecoToMc = recoToMc.find(bestColID); + const auto compatibleIds = track.compatibleCollIds(); + auto itRecoList = mcToReco.find(mcOfTrack); + + if (!compatibleIds.empty() && itRecoList != mcToReco.end() && !itRecoList->second.empty()) { + for (const auto& trueRecoId : itRecoList->second) { + for (const auto& compatibleId : compatibleIds) { + if (compatibleId == trueRecoId) { + recoOfTrueInCompatible = true; break; } } + if (recoOfTrueInCompatible) { + break; + } } + } - // 2) Then check whether any reco collision for the true MC collision exists at all. - if (itRecoList != mcToReco.end() && !itRecoList->second.empty()) { - recoOfTrueExists = true; - } + if (itRecoList != mcToReco.end() && !itRecoList->second.empty()) { + recoOfTrueExists = true; + } - // 3) Finally classify the actually chosen reco collision as right/wrong. - if (bestColID >= 0 && itRecoToMc != recoToMc.end()) { - const int mcFromReco = itRecoToMc->second; - bin = (mcFromReco == mcOfTrack) - ? static_cast(RightWrongBin::Right) - : static_cast(RightWrongBin::Wrong); + if (bestColID < 0) { + registry.fill(HIST("Purity/NeitherReason"), + static_cast(NeitherReasonBin::BestColInvalid)); + } else if (itRecoToMc == recoToMc.end()) { + registry.fill(HIST("Purity/NeitherReason"), + static_cast(NeitherReasonBin::BestColMissingInRecoToMc)); + } else { + const int mcFromReco = itRecoToMc->second; + if (mcFromReco == mcOfTrack) { + bin = static_cast(RightWrongBin::Right); + registry.fill(HIST("Purity/NeitherReason"), + static_cast(NeitherReasonBin::ClassifiedRight)); + } else { + bin = static_cast(RightWrongBin::Wrong); + registry.fill(HIST("Purity/NeitherReason"), + static_cast(NeitherReasonBin::ClassifiedWrong)); } } + } - registry.fill(HIST("RightWrong"), bin); - registry.fill(HIST("Purity/RecoOfTrueExists"), + registry.fill(HIST("RightWrong"), bin); + registry.fill(HIST("Purity/RecoOfTrueExists"), + recoOfTrueExists ? static_cast(BoolBin::Yes) + : static_cast(BoolBin::No)); + registry.fill(HIST("Purity/RecoOfTrueInCompatible"), + recoOfTrueInCompatible ? static_cast(BoolBin::Yes) + : static_cast(BoolBin::No)); + + if (bestColID >= 0) { + uniqueBestRecoCols.insert(bestColID); + } + + if (bin == static_cast(RightWrongBin::Wrong)) { + registry.fill(HIST("Purity/WrongVertexRecoExists"), + recoOfTrueExists ? static_cast(WrongVertexRecoExistsBin::RecoOfTrueExists) + : static_cast(WrongVertexRecoExistsBin::RecoOfTrueMissing)); + registry.fill(HIST("Purity/RecoOfTrueExistsW"), recoOfTrueExists ? static_cast(BoolBin::Yes) : static_cast(BoolBin::No)); - registry.fill(HIST("Purity/RecoOfTrueInCompatible"), + registry.fill(HIST("Purity/RecoOfTrueInCompatibleW"), recoOfTrueInCompatible ? static_cast(BoolBin::Yes) : static_cast(BoolBin::No)); + } - if (bin == static_cast(RightWrongBin::Wrong)) { - registry.fill(HIST("Purity/WrongVertexRecoExists"), - recoOfTrueExists ? static_cast(WrongVertexRecoExistsBin::RecoOfTrueExists) - : static_cast(WrongVertexRecoExistsBin::RecoOfTrueMissing)); - registry.fill(HIST("Purity/RecoOfTrueExistsW"), - recoOfTrueExists ? static_cast(BoolBin::Yes) - : static_cast(BoolBin::No)); - registry.fill(HIST("Purity/RecoOfTrueInCompatibleW"), - recoOfTrueInCompatible ? static_cast(BoolBin::Yes) - : static_cast(BoolBin::No)); - } + if (bin == static_cast(RightWrongBin::Right)) { + registry.fill(HIST("Purity/RecoOfTrueExistsR"), + recoOfTrueExists ? static_cast(BoolBin::Yes) + : static_cast(BoolBin::No)); + registry.fill(HIST("Purity/RecoOfTrueInCompatibleR"), + recoOfTrueInCompatible ? static_cast(BoolBin::Yes) + : static_cast(BoolBin::No)); + } - if (bin == static_cast(RightWrongBin::Right)) { - registry.fill(HIST("Purity/RecoOfTrueExistsR"), - recoOfTrueExists ? static_cast(BoolBin::Yes) - : static_cast(BoolBin::No)); - registry.fill(HIST("Purity/RecoOfTrueInCompatibleR"), - recoOfTrueInCompatible ? static_cast(BoolBin::Yes) - : static_cast(BoolBin::No)); - } + registry.fill(HIST("Purity/RecoSparseAll"), + etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); - registry.fill(HIST("Purity/RecoSparseAll"), + if (isPrimaryCharged) { + registry.fill(HIST("Purity/RecoSparsePrimary"), etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); + } else { + registry.fill(HIST("Purity/RecoSparseSecondary"), + etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); + } - if (isPrimaryCharged) { - registry.fill(HIST("Purity/RecoSparsePrimary"), - etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); - } else { - registry.fill(HIST("Purity/RecoSparseSecondary"), - etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); - } + registry.fill(HIST("RecoSparseAllBest"), + etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); - registry.fill(HIST("RecoSparseAllBest"), - etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); - if (bin == static_cast(RightWrongBin::Wrong)) { - float vzBest = 999.; - float vzTrue = 999.; - auto itVzBest = recoVtxZ.find(bestColID); - if (itVzBest != recoVtxZ.end()) { - vzBest = itVzBest->second; - } - auto itVzTrue = recoToMcVZ.find(vzBest); - if (itVzTrue != recoToMcVZ.end()) { - vzTrue = itVzBest->second; - } - double_t vztrueParticle = mcColObj.posZ(); - double_t diff1 = vzBest - vztrueParticle; - double_t diff2 = vzBest - vzTrue; - registry.fill(HIST("deltaVZ_fromReco"), diff1); - registry.fill(HIST("deltaVZ_fromTrue"), diff2); - registry.fill(HIST("RecoSparseAllBestWrong"), - etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); + if (bin == static_cast(RightWrongBin::Wrong)) { + float vzBest = 999.f; + float vzTrue = 999.f; + auto itVzBest = recoVtxZ.find(bestColID); + if (itVzBest != recoVtxZ.end()) { + vzBest = itVzBest->second; + } + auto itVzTrue = recoVtxZ.find(recoCol); + if (itVzTrue != recoVtxZ.end()) { + vzTrue = itVzTrue->second; } + double_t vztrueParticle = track.mcParticle().template mcCollision_as().posZ(); + double_t diff1 = vzBest - vztrueParticle; + double_t diff2 = vzBest - vzTrue; + registry.fill(HIST("deltaVZ_fromReco"), diff1); + registry.fill(HIST("deltaVZ_fromTrue"), diff2); + registry.fill(HIST("RecoSparseAllBestWrong"), + etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); + } - // Truth DCA components w.r.t. MC collision vertex - const auto dcaXtruth = mcPart.vx() - mcColObj.posX(); // here mcColObj -> mcPart.mcCollision() + if (hasMcLabel) { + const auto mcColObj = track.mcParticle().template mcCollision_as(); + const auto mcPart = track.mcParticle(); + + const auto dcaXtruth = mcPart.vx() - mcColObj.posX(); const auto dcaYtruth = mcPart.vy() - mcColObj.posY(); - const auto dcaZtruth = mcPart.vz() - mcColObj.posZ(); - const auto dcaXYtruth = std::sqrt(dcaXtruth * dcaXtruth + - dcaYtruth * dcaYtruth); + const auto dcaZtruth = hasBestDCAZ ? (mcPart.vz() - mcColObj.posZ()) : 0.f; + const auto dcaXYtruth = std::sqrt(dcaXtruth * dcaXtruth + dcaYtruth * dcaYtruth); const float etaTruth = mcPart.eta(); const bool isPrimaryTruth = mcPart.isPhysicalPrimary(); - // Base truth histograms (independent of right/wrong vertex) registry.fill(HIST("Tracks/dca/Truth/THnDCAxyBestGenTruthAll"), etaTruth, dcaXYtruth, dcaZtruth, dcaXtruth, dcaYtruth); if (isPrimaryTruth) { @@ -2632,82 +2782,68 @@ struct PseudorapidityDensityMFT { registry.fill(HIST("Purity/reco/woOrp/woOrpTracksEtaZvtx"), track.eta(), z); registry.fill(HIST("Purity/reco/woOrp/woOrpTracksPtZvtx"), track.pt(), z); - if (perCollisionSampleCentral.size() > 0) { - registry.fill(HIST("Purity/reco/woOrp/woOrpEtaZvtx_gt0"), track.eta(), z); - registry.fill(HIST("Purity/reco/woOrp/woOrpPtZvtx_gt0"), track.pt(), z); - registry.fill(HIST("Purity/reco/woOrp/woOrpTracksDCAxyZvtx_gt0"), dcaXyCut, z); - registry.fill(HIST("Purity/reco/woOrp/woOrpTracksDCAzZvtx_gt0"), dcaZCut, z); - } + registry.fill(HIST("Purity/reco/woOrp/woOrpEtaZvtx_gt0"), track.eta(), z); + registry.fill(HIST("Purity/reco/woOrp/woOrpPtZvtx_gt0"), track.pt(), z); + registry.fill(HIST("Purity/reco/woOrp/woOrpTracksDCAxyZvtx_gt0"), dcaXyCut, z); + registry.fill(HIST("Purity/reco/woOrp/woOrpTracksDCAzZvtx_gt0"), dcaZCut, z); registry.fill(HIST("Purity/reco/woOrp/woOrpTracksPhiEta"), phi, track.eta()); if (isFakeByLabel) { - // std::cout << " track.eta() " < 0) { - registry.fill(HIST("Purity/reco/woOrp_fake/woOrpEtaZvtx_gt0"), track.eta(), z); - registry.fill(HIST("Purity/reco/woOrp_fake/woOrpPtZvtx_gt0"), track.pt(), z); - } + registry.fill(HIST("Purity/reco/woOrp_fake/woOrpEtaZvtx_gt0"), track.eta(), z); + registry.fill(HIST("Purity/reco/woOrp_fake/woOrpPtZvtx_gt0"), track.pt(), z); } if (isTrueByLabel) { - // Has MC particle - // std::cout << " track.eta() has mc particle " < 0) { - registry.fill(HIST("Purity/reco/woOrp_hasMC/woOrpEtaZvtx_gt0"), track.eta(), z); - registry.fill(HIST("Purity/reco/woOrp_hasMC/woOrpPtZvtx_gt0"), track.pt(), z); - } + registry.fill(HIST("Purity/reco/woOrp_hasMC/woOrpEtaZvtx_gt0"), track.eta(), z); + registry.fill(HIST("Purity/reco/woOrp_hasMC/woOrpPtZvtx_gt0"), track.pt(), z); } if (isSecondaryCharged) { registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpTracksEtaZvtx"), track.eta(), z); registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpTracksPtZvtx"), track.pt(), z); registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpTracksPhiEta"), phi, track.eta()); - if (perCollisionSampleCentral.size() > 0) { - registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpEtaZvtx_gt0"), track.eta(), z); - registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpPtZvtx_gt0"), track.pt(), z); - } + registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpEtaZvtx_gt0"), track.eta(), z); + registry.fill(HIST("Purity/reco/woOrp_secondary/woOrpPtZvtx_gt0"), track.pt(), z); } if (isPrimaryCharged) { registry.fill(HIST("Purity/reco/woOrp_primary/woOrpTracksEtaZvtx"), track.eta(), z); registry.fill(HIST("Purity/reco/woOrp_primary/woOrpTracksPtZvtx"), track.pt(), z); registry.fill(HIST("Purity/reco/woOrp_primary/woOrpTracksPhiEta"), phi, track.eta()); - if (perCollisionSampleCentral.size() > 0) { - registry.fill(HIST("Purity/reco/woOrp_primary/woOrpEtaZvtx_gt0"), track.eta(), z); - registry.fill(HIST("Purity/reco/woOrp_primary/woOrpPtZvtx_gt0"), track.pt(), z); - } + registry.fill(HIST("Purity/reco/woOrp_primary/woOrpEtaZvtx_gt0"), track.eta(), z); + registry.fill(HIST("Purity/reco/woOrp_primary/woOrpPtZvtx_gt0"), track.pt(), z); } ++woOrpCount; - // Category 2: all reco tracks after selections (woOrp) - // --- Primary vs Fake accounting --- const float xTrue = mcColObj.posX(); const float yTrue = mcColObj.posY(); const float zTrue = mcColObj.posZ(); std::array dcaInfOrig{999., 999., 999.}; - std::array dcaChosen{999., 999.}; // (DCAxy, DCAz) to chosen reco vertex - std::array dcaRight{999., 999.}; // (DCAxy, DCAz) to truth MC vertex (reference) - std::array dcaChosenXYZ{999., 999., 999.}; // (DCAx, DCAy, DCAz) to chosen reco vertex + std::array dcaChosen{999., 999.}; + std::array dcaRight{999., 999.}; + std::array dcaChosenXYZ{999., 999., 999.}; const double bZ = o2::base::Propagator::Instance()->getNominalBz(); - // Build forward track parameters once per track, then use a fresh copy for each vertex propagation - std::vector v1; // empty -> null cov + std::vector v1; SMatrix55 tcovs(v1.begin(), v1.end()); SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); o2::track::TrackParCovFwd trackPar0{track.z(), tpars, tcovs, track.chi2()}; - auto trackPar = trackPar0; // copy + auto trackPar = trackPar0; dcaInfOrig = {999., 999., 999.}; - const std::array vtxChosen{collision.posX(), collision.posY(), collision.posZ()}; - trackPar.propagateToDCAhelix(bZ, vtxChosen, dcaInfOrig); - dcaChosenXYZ = dcaInfOrig; // store (DCAx, DCAy, DCAz) - dcaChosen[0] = std::sqrt(dcaInfOrig[0] * dcaInfOrig[0] + dcaInfOrig[1] * dcaInfOrig[1]); - dcaChosen[1] = dcaInfOrig[2]; + auto itVtxChosen = recoVtxByRecoId.find(bestColID); + if (itVtxChosen != recoVtxByRecoId.end()) { + trackPar.propagateToDCAhelix(bZ, itVtxChosen->second, dcaInfOrig); + dcaChosenXYZ = dcaInfOrig; + dcaChosen[0] = std::sqrt(dcaInfOrig[0] * dcaInfOrig[0] + dcaInfOrig[1] * dcaInfOrig[1]); + dcaChosen[1] = dcaInfOrig[2]; + } dcaInfOrig = {999., 999., 999.}; const std::array vtxTruth{xTrue, yTrue, zTrue}; @@ -2716,30 +2852,21 @@ struct PseudorapidityDensityMFT { dcaRight[1] = dcaInfOrig[2]; registry.fill(HIST("Purity/DCAyVsDCAx_Right"), dcaChosenXYZ[2], dcaChosenXYZ[1]); - // Fill only for WRONG-vertex tracks (the diagnostic you want) if (bin == static_cast(RightWrongBin::Wrong)) { - // std::cout <<"dcaxy choosen " << dcaXyCut<<" dcaxy calculated "<(RightWrongBin::Right)) { - // std::cout <<"dcaxy choosen " << dcaXyCut<<" dcaxy calculated "<(RightWrongBin::Right)) { - // Right-vertex: all / primary / secondary registry.fill(HIST("Purity/RecoSparseRightAll"), etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); if (!filledRight) { registry.fill(HIST("Purity/RecoSparseRightAll_EventCount"), static_cast(SingleCountBin::Count)); - filledRight = true; } if (isPrimaryCharged) { @@ -2750,7 +2877,6 @@ struct PseudorapidityDensityMFT { etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); } } else if (bin == static_cast(RightWrongBin::Wrong)) { - // Wrong-vertex: all / primary / secondary registry.fill(HIST("Purity/RecoSparseWrongAll"), etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); if (!filledWrong) { @@ -2765,45 +2891,38 @@ struct PseudorapidityDensityMFT { etaReco, dcaXYReco, dcaZReco, dcaXReco, dcaYReco); } } - // Reconstructed vertex position for bestbestColID, if available - auto itVtxX = recoVtxX.find(bestColID); + auto itVtxX = recoVtxX.find(bestColID); if (itVtxX != recoVtxX.end()) { - const float xReco = itVtxX->second; const float yReco = recoVtxY[bestColID]; const float zReco = recoVtxZ[bestColID]; const bool recoVzIn = (zReco >= cfgVzCut1) && (zReco <= cfgVzCut2); const bool trueVzIn = (zTrue >= cfgVzCut1) && (zTrue <= cfgVzCut2); - if (!(recoVzIn && trueVzIn)) { - continue; // skip filling Delta* histos for this track - } - - const float deltaXvtx = xReco - xTrue; - const float deltaYvtx = yReco - yTrue; - const float deltaZvtx = zReco - zTrue; + if (recoVzIn && trueVzIn) { + const float deltaXvtx = xReco - xTrue; + const float deltaYvtx = yReco - yTrue; + const float deltaZvtx = zReco - zTrue; - // We are interested in how far the WRONG collisions are from the true one - if (bin == static_cast(RightWrongBin::Wrong)) { - registry.fill(HIST("Purity/DeltaXWrong"), deltaXvtx); - registry.fill(HIST("Purity/DeltaYWrong"), deltaYvtx); - registry.fill(HIST("Purity/DeltaZWrong"), deltaZvtx); - } - if (bin == static_cast(RightWrongBin::Right)) { - registry.fill(HIST("Purity/DeltaXRight"), deltaXvtx); - registry.fill(HIST("Purity/DeltaYRight"), deltaYvtx); - registry.fill(HIST("Purity/DeltaZRight"), deltaZvtx); + if (bin == static_cast(RightWrongBin::Wrong)) { + registry.fill(HIST("Purity/DeltaXWrong"), deltaXvtx); + registry.fill(HIST("Purity/DeltaYWrong"), deltaYvtx); + registry.fill(HIST("Purity/DeltaZWrong"), deltaZvtx); + } + if (bin == static_cast(RightWrongBin::Right)) { + registry.fill(HIST("Purity/DeltaXRight"), deltaXvtx); + registry.fill(HIST("Purity/DeltaYRight"), deltaYvtx); + registry.fill(HIST("Purity/DeltaZRight"), deltaZvtx); + } } } - // --- Delta-DCA components: truth minus reco --- const float deltaDCAxy = dcaXYtruth - dcaXYReco; const float deltaDCAz = dcaZtruth - dcaZReco; const float deltaDCAx = dcaXtruth - dcaXReco; const float deltaDCAy = dcaYtruth - dcaYReco; if (bin == static_cast(RightWrongBin::Right)) { - // Right-vertex: all / primary / secondary registry.fill(HIST("Tracks/dca/Truth/THnDeltaDCARightAll"), deltaDCAxy, deltaDCAz, deltaDCAx, deltaDCAy); if (isPrimaryCharged) { @@ -2814,7 +2933,6 @@ struct PseudorapidityDensityMFT { deltaDCAxy, deltaDCAz, deltaDCAx, deltaDCAy); } } else { - // Wrong-vertex: all / primary / secondary registry.fill(HIST("Tracks/dca/Truth/THnDeltaDCAWrongAll"), deltaDCAxy, deltaDCAz, deltaDCAx, deltaDCAy); if (isPrimaryCharged) { @@ -2826,42 +2944,58 @@ struct PseudorapidityDensityMFT { } } - // auto mcColObj = track.mcParticle().mcCollision_as(); - // True primary match (purity numerator) registry.fill(HIST("Purity/mc/PrimaryAll"), static_cast(SingleCountBin::Count)); registry.fill(HIST("Purity/mc/PrimaryAllEta"), mcPart.eta()); registry.fill(HIST("Purity/mc/PrimaryTracksEtaZvtx"), mcPart.eta(), mcCollision.posZ()); - if (perCollisionSampleCentral.size() > 0) { - registry.fill(HIST("Purity/mc/PrimaryTracksEtaZvtx_gt0"), mcPart.eta(), mcCollision.posZ()); - registry.fill(HIST("Purity/mc/PrimaryTracksPtZvtx_gt0"), mcPart.pt(), mcCollision.posZ()); - registry.fill(HIST("Purity/mc/PrimaryTracksDCAxyZvtx_gt0"), dcaXyCut, mcCollision.posZ()); - registry.fill(HIST("Purity/mc/PrimaryTracksDCAzZvtx_gt0"), dcaZCut, mcCollision.posZ()); - } + registry.fill(HIST("Purity/mc/PrimaryTracksEtaZvtx_gt0"), mcPart.eta(), mcCollision.posZ()); + registry.fill(HIST("Purity/mc/PrimaryTracksPtZvtx_gt0"), mcPart.pt(), mcCollision.posZ()); + registry.fill(HIST("Purity/mc/PrimaryTracksDCAxyZvtx_gt0"), dcaXyCut, mcCollision.posZ()); + registry.fill(HIST("Purity/mc/PrimaryTracksDCAzZvtx_gt0"), dcaZCut, mcCollision.posZ()); registry.fill(HIST("Purity/mc/PrimaryTracksPhiEta"), mcPart.phi(), mcPart.eta()); registry.fill(HIST("Purity/SelectedAfterDCAxy/PrimaryAll"), static_cast(SingleCountBin::Count)); registry.fill(HIST("Purity/SelectedAfterDCAxy/PrimaryAllEta"), mcPart.eta()); countedPrimary = true; - // --- Purity profiles (mean of indicator gives purity) --- registry.fill(HIST("Purity/PurityOverall"), static_cast(SingleCountBin::Count), countedPrimary ? static_cast(BoolBin::Yes) : static_cast(BoolBin::No)); - registry.fill(HIST("Purity/PurityVsEta"), track.eta(), countedPrimary ? static_cast(BoolBin::Yes) : static_cast(BoolBin::No)); - } // track loop - registry.fill(HIST("Purity/HashTableRowCounts"), - static_cast(HashTableRowCountsBin::UniqueBestRecoCols), uniqueBestRecoCols.size()); + } } - registry.fill(HIST("Purity/reco/woOrp/nTrk"), woOrpCount); - registry.fill(HIST("Purity/reco/PNchMFT_afterCuts"), nMftSelectedAfterCuts); - } // collision + } + + registry.fill(HIST("Purity/HashTableRowCounts"), + static_cast(HashTableRowCountsBin::UniqueBestRecoCols), uniqueBestRecoCols.size()); + registry.fill(HIST("Purity/reco/woOrp/nTrk"), woOrpCount); + registry.fill(HIST("Purity/reco/PNchMFT_afterCuts"), nMftSelectedAfterCuts); + } + void processGenReco3d(McCollisionsWithExtra::iterator const& mcCollision, + o2::soa::SmallGroups> const& collisions, + FullBCs const& bcs, + MFTTracksLabeled3d const& tracks, + FiCentralTracks const& midtracks, + aod::McParticles const& mcParticles) + { + processGenReco(mcCollision, collisions, bcs, tracks, midtracks, mcParticles); } - PROCESS_SWITCH(PseudorapidityDensityMFT, processGenReco, - "Process particle-level info of pt", false); + void processGenReco2d(McCollisionsWithExtra::iterator const& mcCollision, + o2::soa::SmallGroups> const& collisions, + FullBCs const& bcs, + MFTTracksLabeled2d const& tracks, + FiCentralTracks const& midtracks, + aod::McParticles const& mcParticles) + { + processGenReco(mcCollision, collisions, bcs, tracks, midtracks, mcParticles); + } + PROCESS_SWITCH(PseudorapidityDensityMFT, processGenReco3d, + "Process gen-reco info with BestCollisionsFwd3d", true); + + PROCESS_SWITCH(PseudorapidityDensityMFT, processGenReco2d, + "Process gen-reco info with BestCollisionsFwd", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Tasks/GlobalEventProperties/studyPnch.cxx b/PWGLF/Tasks/GlobalEventProperties/studyPnch.cxx index 62f5e93d830..35ad555d467 100644 --- a/PWGLF/Tasks/GlobalEventProperties/studyPnch.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/studyPnch.cxx @@ -76,6 +76,7 @@ AxisSpec axisEvent{10, 0.5, 10.5, "#Event", "EventAxis"}; AxisSpec axisVtxZ{40, -20, 20, "Vertex Z", "VzAxis"}; AxisSpec axisEta{40, -2, 2, "#eta", "EtaAxis"}; AxisSpec axisPhi{629, 0, o2::constants::math::TwoPI, "#phi"}; +AxisSpec axisCollSel{5, 0.5, 5.5, "#Event", "CollSelAxis"}; auto static constexpr kMinCharge = 3.f; struct StudyPnch { @@ -97,6 +98,7 @@ struct StudyPnch { ConfigurableAxis ft0aMultHistBin{"ft0aMultHistBin", {501, -0.5, 500.5}, ""}; ConfigurableAxis ft0cMultHistBin{"ft0cMultHistBin", {501, -0.5, 500.5}, ""}; ConfigurableAxis ptHistBin{"ptHistBin", {200, 0., 20.}, ""}; + ConfigurableAxis binsDCA{"binsDCA", {500, -10.0f, 10.0f}, ""}; ConfigurableAxis countNumberTracks{"countNumberTracks", {10, -0.5, 9.5}, ""}; Configurable isApplyTFcut{"isApplyTFcut", true, "Enable TimeFrameBorder cut"}; @@ -105,6 +107,8 @@ struct StudyPnch { Configurable isApplyInelgt0{"isApplyInelgt0", false, "Enable INEL > 0 condition"}; Configurable isApplyExtraPhiCut{"isApplyExtraPhiCut", false, "Enable extra phi cut"}; Configurable isApplyTVX{"isApplyTVX", false, "Enable TVX trigger sel"}; + Configurable isApplyCheckID{"isApplyCheckID", true, "Select Tracks evaluating Collision ID"}; + Configurable isApplyDuplicatedTrack{"isApplyDuplicatedTrack", true, "Select tracks that are not duplicated"}; void init(InitContext const&) { @@ -115,6 +119,7 @@ struct StudyPnch { AxisSpec axisFt0cMult = {ft0cMultHistBin, "ft0c", "FT0CMultAxis"}; AxisSpec axisPt = {ptHistBin, "pT", "pTAxis"}; AxisSpec axisCountNumberTracks = {countNumberTracks, "Count", "CountAxis"}; + AxisSpec dcaAxis = {binsDCA, "DCA vs PV"}; histos.add("EventHist", "EventHist", kTH1D, {axisEvent}, false); histos.add("VtxZHist", "VtxZHist", kTH1D, {axisVtxZ}, false); @@ -129,8 +134,19 @@ struct StudyPnch { x->SetBinLabel(6, "INEL > 0"); x->SetBinLabel(7, "|vz| < 10"); - if (doprocessData || doprocessCorrelation || doprocessMonteCarlo || doprocessModifiedMonteCarlo) { + histos.add("SelCollsHist", "SelCollsHist", kTH1D, {axisCollSel}, false); + auto hstat_colls = histos.get(HIST("SelCollsHist")); + auto* xColls = hstat_colls->GetXaxis(); + xColls->SetBinLabel(1, "All collisions"); + xColls->SetBinLabel(2, "Best Collision Selection"); + xColls->SetBinLabel(3, "Has MC Collision Selection"); + + if (doprocessData || doprocessCorrelation || doprocessMonteCarlo) { histos.add("PhiVsEtaHist", "PhiVsEtaHist", kTH2F, {axisPhi, axisEta}, false); + histos.add("EtaHist", "EtaHist", kTH1D, {axisEta}, false); + histos.add("PhiHist", "PhiHist", kTH1D, {axisPhi}, false); + histos.add("hdcaxy", "dca to pv in the xy plane", kTH1D, {dcaAxis}, false); + histos.add("hdcaz", "dca to pv in the z axis", kTH1D, {dcaAxis}, false); } if (doprocessData) { histos.add("hMultiplicityData", "hMultiplicityData", kTH1F, {axisMult}, true); @@ -143,17 +159,14 @@ struct StudyPnch { histos.add("NPVtracks_vs_GlobalMult", "NPVtracks_vs_GlobalMult", kTH2F, {axisPV, axisMult}, true); } if (doprocessMonteCarlo) { + histos.add("PhiVsEtaGenHist", "PhiVsEtaGenHist", kTH2F, {axisPhi, axisEta}, false); + histos.add("EtaGenHist", "EtaGenHist", kTH1D, {axisEta}, false); + histos.add("PhiGenHist", "PhiGenHist", kTH1D, {axisPhi}, false); histos.add("hMultiplicityMCrec", "hMultiplicityMCrec", kTH1F, {axisMult}, true); histos.add("hMultiplicityMCgen", "hMultiplicityMCgen", kTH1F, {axisMult}, true); histos.add("hResponseMatrix", "hResponseMatrix", kTH2F, {axisMult, axisMult}, true); histos.add("hCountNTracks", "hCountNTracks", kTH1F, {axisCountNumberTracks}, true); } - if (doprocessModifiedMonteCarlo) { - histos.add("hMultiplicityMCrecMod", "hMultiplicityMCrecMod", kTH1F, {axisMult}, true); - histos.add("hMultiplicityMCgenMod", "hMultiplicityMCgenMod", kTH1F, {axisMult}, true); - histos.add("hResponseMatrixMod", "hResponseMatrixMod", kTH2F, {axisMult, axisMult}, true); - histos.add("hCountNTracksMod", "hCountNTracksMod", kTH1F, {axisCountNumberTracks}, true); - } if (doprocessEvtLossSigLossMC) { histos.add("MCEventHist", "MCEventHist", kTH1F, {axisEvent}, false); auto hstat = histos.get(HIST("MCEventHist")); @@ -242,6 +255,10 @@ struct StudyPnch { if (!isTrackSelected(track)) { continue; } + histos.fill(HIST("hdcaxy"), track.dcaXY()); + histos.fill(HIST("hdcaz"), track.dcaZ()); + histos.fill(HIST("EtaHist"), track.eta()); + histos.fill(HIST("PhiHist"), track.phi()); histos.fill(HIST("PhiVsEtaHist"), track.phi(), track.eta()); nTrk++; } @@ -256,11 +273,12 @@ struct StudyPnch { if (!isGenTrackSelected(track)) { continue; } - if (track.mcCollisionId() != McCol.globalIndex()) { + if (track.mcCollisionId() != McCol.mcCollisionId()) { continue; } - - histos.fill(HIST("PhiVsEtaHist"), track.phi(), track.eta()); + histos.fill(HIST("EtaGenHist"), track.eta()); + histos.fill(HIST("PhiGenHist"), track.phi()); + histos.fill(HIST("PhiVsEtaGenHist"), track.phi(), track.eta()); nTrk++; } return nTrk; @@ -270,30 +288,28 @@ struct StudyPnch { int countNTracksMcCol(countTrk const& tracks, McColType const& McCol) { auto nTrk = 0; - std::unordered_map recoFrequencies; // Map that stores globalIndex and the times it appears + std::vector mcRecIDs; for (const auto& track : tracks) { if (!isTrackSelected(track)) { continue; } if (track.has_mcParticle()) { auto particle = track.mcParticle(); - if (particle.mcCollisionId() != McCol.mcCollisionId()) { + if (isApplyCheckID && particle.mcCollisionId() != McCol.mcCollisionId()) { continue; } - auto globalIndex = particle.globalIndex(); - recoFrequencies[globalIndex]++; // Increment the count for this globalIndex - } - histos.fill(HIST("PhiVsEtaHist"), track.phi(), track.eta()); - } - // Loop to fill the histogram without cloned tracks - for (const auto& [globalIndex, frequency] : recoFrequencies) { - histos.fill(HIST("hCountNTracks"), frequency); - // Fill histogram with not cloned tracks - if (frequency == 1) { + if (isApplyDuplicatedTrack && find(mcRecIDs.begin(), mcRecIDs.end(), particle.globalIndex()) != mcRecIDs.end()) { + continue; + } + mcRecIDs.push_back(particle.globalIndex()); nTrk++; } + histos.fill(HIST("hdcaxy"), track.dcaXY()); + histos.fill(HIST("hdcaz"), track.dcaZ()); + histos.fill(HIST("EtaHist"), track.eta()); + histos.fill(HIST("PhiHist"), track.phi()); + histos.fill(HIST("PhiVsEtaHist"), track.phi(), track.eta()); } - // return recoFrequencies; return nTrk; } @@ -328,47 +344,31 @@ struct StudyPnch { histos.fill(HIST("NPVtracks_vs_GlobalMult"), cols.multNTracksPV(), mult); } - void processMonteCarlo(ColMCTrueTable::iterator const& mcCollision, ColMCRecTable const& RecCols, TrackMCTrueTable const& GenParticles, FilTrackMCRecTable const& RecTracks) - { - for (const auto& RecCol : RecCols) { - if (!isEventSelected(RecCol)) { - continue; - } - auto recTracksPart = RecTracks.sliceBy(perCollision, RecCol.globalIndex()); - auto multrec = countNTracksMcCol(recTracksPart, RecCol); - if (multrec > 0) { - histos.fill(HIST("hMultiplicityMCrec"), multrec); - } - auto multgen = countGenTracks(GenParticles, mcCollision); - if (multgen > 0 && multrec > 0) { - histos.fill(HIST("hMultiplicityMCgen"), multgen); - histos.fill(HIST("hResponseMatrix"), multrec, multgen); - } - } - } - - void processModifiedMonteCarlo(soa::Join::iterator const& mcCollision, ColMCRecTable const& RecCols, TrackMCTrueTable const& GenParticles, FilTrackMCRecTable const& RecTracks) + void processMonteCarlo(soa::Join::iterator const& mcCollision, ColMCRecTable const& RecCols, TrackMCTrueTable const& GenParticles, FilTrackMCRecTable const& RecTracks) { for (const auto& RecCol : RecCols) { if (!isEventSelected(RecCol)) { continue; } + histos.fill(HIST("SelCollsHist"), 1); // Evaluation of reconstructed collisions with more than 1 contributor if (RecCol.globalIndex() != mcCollision.bestCollisionIndex()) { continue; } + histos.fill(HIST("SelCollsHist"), 2); if (!RecCol.has_mcCollision()) { continue; } + histos.fill(HIST("SelCollsHist"), 3); auto recTracksPart = RecTracks.sliceBy(perCollision, RecCol.globalIndex()); auto multrec = countNTracksMcCol(recTracksPart, RecCol); if (multrec > 0) { - histos.fill(HIST("hMultiplicityMCrecMod"), multrec); + histos.fill(HIST("hMultiplicityMCrec"), multrec); } - auto multgen = countGenTracks(GenParticles, mcCollision); + auto multgen = countGenTracks(GenParticles, RecCol); if (multgen > 0 && multrec > 0) { - histos.fill(HIST("hMultiplicityMCgenMod"), multgen); - histos.fill(HIST("hResponseMatrixMod"), multrec, multgen); + histos.fill(HIST("hMultiplicityMCgen"), multgen); + histos.fill(HIST("hResponseMatrix"), multrec, multgen); } } } @@ -386,9 +386,15 @@ struct StudyPnch { } // All generated events histos.fill(HIST("MCEventHist"), 1); - auto multAll = countGenTracks(GenParticles, mcCollision); - if (multAll > 0) { - histos.fill(HIST("hMultiplicityMCgenAll"), multAll); + auto nTrk_multAll = 0; + for (const auto& GenParticle : GenParticles) { + if (!isGenTrackSelected(GenParticle)) { + continue; + } + nTrk_multAll++; + } + if (nTrk_multAll > 0) { + histos.fill(HIST("hMultiplicityMCgenAll"), nTrk_multAll); } bool atLeastOne = false; @@ -407,9 +413,15 @@ struct StudyPnch { if (atLeastOne) { histos.fill(HIST("MCEventHist"), 2); - auto multSel = countGenTracks(GenParticles, mcCollision); - if (multSel > 0) { - histos.fill(HIST("hMultiplicityMCgenSel"), multSel); + auto nTrk_multSel = 0; + for (const auto& GenParticle : GenParticles) { + if (!isGenTrackSelected(GenParticle)) { + continue; + } + nTrk_multSel++; + } + if (nTrk_multSel > 0) { + histos.fill(HIST("hMultiplicityMCgenSel"), nTrk_multSel); } } } @@ -417,7 +429,6 @@ struct StudyPnch { PROCESS_SWITCH(StudyPnch, processData, "process data CentFT0C", false); PROCESS_SWITCH(StudyPnch, processCorrelation, "do correlation study in data", false); PROCESS_SWITCH(StudyPnch, processMonteCarlo, "process MC CentFT0C", false); - PROCESS_SWITCH(StudyPnch, processModifiedMonteCarlo, "process MC CentFT0C", false); PROCESS_SWITCH(StudyPnch, processEvtLossSigLossMC, "process Signal Loss, Event Loss", false); };