diff --git a/Analysis/ALICE3/include/ALICE3Analysis/MID.h b/Analysis/ALICE3/include/ALICE3Analysis/MID.h new file mode 100644 index 0000000000000..7d932b90802d3 --- /dev/null +++ b/Analysis/ALICE3/include/ALICE3Analysis/MID.h @@ -0,0 +1,44 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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 MID.h +/// \author Vít Kučera +/// \note Based on RICH.h +/// \brief Set of tables for the ALICE3 MID information +/// + +#ifndef O2_ANALYSIS_ALICE3_MID_H_ +#define O2_ANALYSIS_ALICE3_MID_H_ + +// O2 includes +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace alice3mid +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! +DECLARE_SOA_INDEX_COLUMN(Track, track); //! +DECLARE_SOA_COLUMN(MIDIsMuon, midIsMuon, uint8_t); //! FIXME: To be changed to bool once bool columns are groupable. +} // namespace alice3mid + +DECLARE_SOA_TABLE(MIDs, "AOD", "MID", //! + o2::soa::Index<>, + alice3mid::CollisionId, + alice3mid::TrackId, + alice3mid::MIDIsMuon); + +using MID = MIDs::iterator; + +} // namespace o2::aod + +#endif // O2_ANALYSIS_ALICE3_MID_H_ diff --git a/Analysis/Core/include/AnalysisCore/HFSelectorCuts.h b/Analysis/Core/include/AnalysisCore/HFSelectorCuts.h index f2264a1a94c97..66a6f09428765 100644 --- a/Analysis/Core/include/AnalysisCore/HFSelectorCuts.h +++ b/Analysis/Core/include/AnalysisCore/HFSelectorCuts.h @@ -29,7 +29,8 @@ enum Code { kDPlus = 411, kLambdaCPlus = 4122, kXiCPlus = 4232, - kJpsi = 443 + kJpsi = 443, + kX3872 = 9920443 }; } // namespace pdg @@ -99,7 +100,7 @@ auto pTBinsVec = std::vector{pTBins, pTBins + npTBins + 1}; // default values for the cuts constexpr double cuts[npTBins][nCutVars] = {{1.65, 2.15, 0.5, 100.}, /* 1 < pt < 5 */ - {1.65, 2.15, 0.5, 100.}}; /* 5 < pt > 1000 */ + {1.65, 2.15, 0.5, 100.}}; /* 5 < pt < 1000 */ // row labels static const std::vector pTBinLabels{}; @@ -383,7 +384,7 @@ static const std::vector cutVarLabels = {"m", "pT p", "pT K", "pT P namespace hf_cuts_jpsi_toee { static constexpr int npTBins = 9; -static constexpr int nCutVars = 4; +static constexpr int nCutVars = 5; // default values for the pT bin edges (can be used to configure histogram axis) // offset by 1 from the bin numbers in cuts array constexpr double pTBins[npTBins + 1] = { @@ -401,15 +402,15 @@ constexpr double pTBins[npTBins + 1] = { auto pTBins_v = std::vector{pTBins, pTBins + npTBins + 1}; // default values for the cuts -constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.2, 0.4, 1}, /* 0 < pT < 0.5 */ - {0.5, 0.2, 0.4, 1}, /* 0.5 < pT < 1 */ - {0.5, 0.2, 0.4, 1}, /* 1 < pT < 2 */ - {0.5, 0.2, 0.4, 1}, /* 2 < pT < 3 */ - {0.5, 0.2, 0.4, 1}, /* 3 < pT < 4 */ - {0.5, 0.2, 0.4, 1}, /* 4 < pT < 5 */ - {0.5, 0.2, 0.4, 1}, /* 5 < pT < 7 */ - {0.5, 0.2, 0.4, 1}, /* 7 < pT < 10 */ - {0.5, 0.2, 0.4, 1}}; /* 10 < pT < 15 */ +constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.2, 0.4, 1, 1.}, /* 0 < pT < 0.5 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 0.5 < pT < 1 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 1 < pT < 2 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 2 < pT < 3 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 3 < pT < 4 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 4 < pT < 5 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 5 < pT < 7 */ + {0.5, 0.2, 0.4, 1, 1.}, /* 7 < pT < 10 */ + {0.5, 0.2, 0.4, 1, 1.}}; /* 10 < pT < 15 */ // row labels static const std::vector pTBinLabels = { @@ -424,8 +425,54 @@ static const std::vector pTBinLabels = { "pT bin 8"}; // column labels -static const std::vector cutVarLabels = {"m", "DCA_xy", "DCA_z", "pT El"}; +static const std::vector cutVarLabels = {"m", "DCA_xy", "DCA_z", "pT El", "chi2PCA"}; } // namespace hf_cuts_jpsi_toee + +namespace hf_cuts_x_tojpsipipi +{ +static constexpr int npTBins = 9; +static constexpr int nCutVars = 7; +// default values for the pT bin edges (can be used to configure histogram axis) +// offset by 1 from the bin numbers in cuts array +constexpr double pTBins[npTBins + 1] = { + 0, + 0.5, + 1.0, + 2.0, + 3.0, + 4.0, + 5.0, + 7.0, + 10.0, + 15.0, +}; +auto pTBins_v = std::vector{pTBins, pTBins + npTBins + 1}; + +// default values for the cuts +// m CPA d0Jpsi d0Pi pTJpsi pTPi chi2PCA +constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.}, /* 0 pTBinLabels = { + "pT bin 0", + "pT bin 1", + "pT bin 2", + "pT bin 3", + "pT bin 4", + "pT bin 5", + "pT bin 6", + "pT bin 7", + "pT bin 8"}; +// column labels +static const std::vector cutVarLabels = {"m", "CPA", "d0 Jpsi", "d0 Pi", "pT Jpsi", "pT Pi", "chi2PCA"}; +} // namespace hf_cuts_x_tojpsipipi } // namespace o2::analysis #endif // HF_SELECTOR_CUTS_H_ diff --git a/Analysis/Core/include/AnalysisCore/TrackSelectorPID.h b/Analysis/Core/include/AnalysisCore/TrackSelectorPID.h index 3eda33cd74708..6c0267c3bab41 100644 --- a/Analysis/Core/include/AnalysisCore/TrackSelectorPID.h +++ b/Analysis/Core/include/AnalysisCore/TrackSelectorPID.h @@ -252,9 +252,159 @@ class TrackSelectorPID } } - /// Returns status of combined PID selection for a given track. + // RICH + + /// Set pT range where RICH PID is applicable. + void setRangePtRICH(float ptMin, float ptMax) + { + mPtRICHMin = ptMin; + mPtRICHMax = ptMax; + } + + /// Set RICH nσ range in which a track should be accepted. + void setRangeNSigmaRICH(float nsMin, float nsMax) + { + mNSigmaRICHMin = nsMin; + mNSigmaRICHMax = nsMax; + } + + /// Set RICH nσ range in which a track should be conditionally accepted if combined with TOF. + void setRangeNSigmaRICHCondTOF(float nsMin, float nsMax) + { + mNSigmaRICHMinCondTOF = nsMin; + mNSigmaRICHMaxCondTOF = nsMax; + } + + /// Checks if track is OK for RICH PID. + /// \param track track + /// \return true if track is OK for RICH PID + template + bool isValidTrackPIDRICH(const T& track) + { + if (track.richId() < 0) { + return false; + } + auto pt = track.pt(); + return mPtRICHMin <= pt && pt <= mPtRICHMax; + } + + /// Checks if track is compatible with given particle species hypothesis within given RICH nσ range. + /// \param track track + /// \param conditionalTOF variable to store the result of selection with looser cuts for conditional accepting of track if combined with TOF + /// \return true if track satisfies RICH PID hypothesis for given RICH nσ range + template + bool isSelectedTrackPIDRICH(const T& track, bool& conditionalTOF) + { + // Accept if selection is disabled via large values. + if (mNSigmaRICHMin < -999. && mNSigmaRICHMax > 999.) { + return true; + } + + // Get nσ for a given particle hypothesis. + double nSigma = 100.; + switch (mPdg) { + case kElectron: { + nSigma = track.rich().richNsigmaEl(); + break; + } + case kMuonMinus: { + nSigma = track.rich().richNsigmaMu(); + break; + } + case kPiPlus: { + nSigma = track.rich().richNsigmaPi(); + break; + } + case kKPlus: { + nSigma = track.rich().richNsigmaKa(); + break; + } + case kProton: { + nSigma = track.rich().richNsigmaPr(); + break; + } + default: { + LOGF(error, "ERROR: RICH PID not implemented for PDG %d", mPdg); + assert(false); + } + } + + if (mNSigmaRICHMinCondTOF < -999. && mNSigmaRICHMaxCondTOF > 999.) { + conditionalTOF = true; + } else { + conditionalTOF = mNSigmaRICHMinCondTOF <= nSigma && nSigma <= mNSigmaRICHMaxCondTOF; + } + return mNSigmaRICHMin <= nSigma && nSigma <= mNSigmaRICHMax; + } + + /// Returns status of RICH PID selection for a given track. + /// \param track track + /// \return RICH selection status (see TrackSelectorPID::Status) + template + int getStatusTrackPIDRICH(const T& track) + { + if (isValidTrackPIDRICH(track)) { + bool condTOF = false; + if (isSelectedTrackPIDRICH(track, condTOF)) { + return Status::PIDAccepted; // accepted + } else if (condTOF) { + return Status::PIDConditional; // potential to be accepted if combined with TOF + } else { + return Status::PIDRejected; // rejected + } + } else { + return Status::PIDNotApplicable; // PID not applicable + } + } + + // MID + + /// Checks if track is OK for MID PID. + /// \param track track + /// \return true if track is OK for MID PID + template + bool isValidTrackPIDMID(const T& track) + { + return track.midId() > -1; + } + + /// Checks if track is compatible with muon hypothesis in the MID detector. + /// \param track track + /// \return true if track has been identified as muon by the MID detector + template + bool isSelectedTrackPIDMID(const T& track) + { + if (mPdg != kMuonMinus) { + return false; + } + return track.mid().midIsMuon() == 1; // FIXME: change to return track.midIsMuon() once the column is bool. + } + + /// Returns status of MID PID selection for a given track. + /// \param track track + /// \return MID selection status (see TrackSelectorPID::Status) + template + int getStatusTrackPIDMID(const T& track) + { + if (mPdg != kMuonMinus) { + return Status::PIDRejected; + } + if (isValidTrackPIDMID(track)) { + if (isSelectedTrackPIDMID(track)) { + return Status::PIDAccepted; // accepted + } else { + return Status::PIDRejected; // rejected + } + } else { + return Status::PIDNotApplicable; // PID not applicable + } + } + + // Combined selection (TPC + TOF) + + /// Returns status of combined PID (TPC + TOF) selection for a given track. /// \param track track - /// \return combined-selection status (see TrackSelectorPID::Status) + /// \return status of combined PID (TPC + TOF) (see TrackSelectorPID::Status) template int getStatusTrackPIDAll(const T& track) { @@ -273,6 +423,60 @@ class TrackSelectorPID return Status::PIDNotApplicable; // (NotApplicable for one detector) and (NotApplicable or Conditional for the other) } + /// Checks whether a track is identified as electron and rejected as pion by TOF or RICH. + /// \param track track + /// \param useTOF switch to use TOF + /// \param useRICH switch to use RICH + /// \return true if track is selected by TOF or RICH + /// \note Ported from https://github.com/feisenhu/ALICE3-LoI-LMee/blob/main/efficiency/macros/anaEEstudy.cxx + template + bool isElectronAndNotPion(const T& track, bool useTOF = true, bool useRICH = true) + { + bool isSelTOF = false; + bool isSelRICH = false; + bool hasRICH = track.richId() > -1; + bool hasTOF = isValidTrackPIDTOF(track); + auto nSigmaTOFEl = track.tofNSigmaEl(); + auto nSigmaTOFPi = track.tofNSigmaPi(); + auto nSigmaRICHEl = hasRICH ? track.rich().richNsigmaEl() : -1000.; + auto nSigmaRICHPi = hasRICH ? track.rich().richNsigmaPi() : -1000.; + auto p = track.p(); + + // TOF + if (useTOF && hasTOF && (p < 0.6)) { + if (p > 0.4 && hasRICH) { + if ((std::abs(nSigmaTOFEl) < mNSigmaTOFMax) && (std::abs(nSigmaRICHEl) < mNSigmaRICHMax)) { + isSelTOF = true; // is selected as electron by TOF and RICH + } + } else if (p <= 0.4) { + if (std::abs(nSigmaTOFEl) < mNSigmaTOFMax) { + isSelTOF = true; // is selected as electron by TOF + } + } else { + isSelTOF = false; // This is rejecting all the heavier particles which do not have a RICH signal in the p area of 0.4-0.6 GeV/c + } + if (std::abs(nSigmaTOFPi) < mNSigmaTOFMax) { + isSelTOF = false; // is selected as pion by TOF + } + } else { + isSelTOF = false; + } + + // RICH + if (useRICH && hasRICH) { + if (std::abs(nSigmaRICHEl) < mNSigmaRICHMax) { + isSelRICH = true; // is selected as electron by RICH + } + if ((std::abs(nSigmaRICHPi) < mNSigmaRICHMax) && (p > 1.0) && (p < 2.0)) { + isSelRICH = false; // is selected as pion by RICH + } + } else { + isSelRICH = false; + } + + return isSelRICH || isSelTOF; + } + private: uint mPdg = kPiPlus; ///< PDG code of the expected particle @@ -291,6 +495,14 @@ class TrackSelectorPID float mNSigmaTOFMax = 3.; ///< maximum number of TOF σ float mNSigmaTOFMinCondTPC = -1000.; ///< minimum number of TOF σ if combined with TPC float mNSigmaTOFMaxCondTPC = 1000.; ///< maximum number of TOF σ if combined with TPC + + // RICH + float mPtRICHMin = 0.; ///< minimum pT for RICH PID [GeV/c] + float mPtRICHMax = 100.; ///< maximum pT for RICH PID [GeV/c] + float mNSigmaRICHMin = -3.; ///< minimum number of RICH σ + float mNSigmaRICHMax = 3.; ///< maximum number of RICH σ + float mNSigmaRICHMinCondTOF = -1000.; ///< minimum number of RICH σ if combined with TOF + float mNSigmaRICHMaxCondTOF = 1000.; ///< maximum number of RICH σ if combined with TOF }; #endif // O2_ANALYSIS_TRACKSELECTORPID_H_ diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h index 284646648de2a..57621f8695f02 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h @@ -40,9 +40,26 @@ DECLARE_SOA_TABLE(HFSelLcCandidate, "AOD", "HFSELLCCAND", //! namespace hf_selcandidate_jpsi { DECLARE_SOA_COLUMN(IsSelJpsiToEE, isSelJpsiToEE, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToMuMu, isSelJpsiToMuMu, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToEETopol, isSelJpsiToEETopol, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToEETpc, isSelJpsiToEETpc, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToEETof, isSelJpsiToEETof, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToEERich, isSelJpsiToEERich, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToEETofRich, isSelJpsiToEETofRich, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToMuMuTopol, isSelJpsiToMuMuTopol, int); //! +DECLARE_SOA_COLUMN(IsSelJpsiToMuMuMid, isSelJpsiToMuMuMid, int); //! } // namespace hf_selcandidate_jpsi -DECLARE_SOA_TABLE(HFSelJpsiToEECandidate, "AOD", "HFSELJPSICAND", //! - hf_selcandidate_jpsi::IsSelJpsiToEE); +DECLARE_SOA_TABLE(HFSelJpsiCandidate, "AOD", "HFSELJPSICAND", //! + hf_selcandidate_jpsi::IsSelJpsiToEE, + hf_selcandidate_jpsi::IsSelJpsiToMuMu, + hf_selcandidate_jpsi::IsSelJpsiToEETopol, + hf_selcandidate_jpsi::IsSelJpsiToEETpc, + hf_selcandidate_jpsi::IsSelJpsiToEETof, + hf_selcandidate_jpsi::IsSelJpsiToEERich, + hf_selcandidate_jpsi::IsSelJpsiToEETofRich, + hf_selcandidate_jpsi::IsSelJpsiToMuMuTopol, + hf_selcandidate_jpsi::IsSelJpsiToMuMuMid); + namespace hf_selcandidate_lc_k0sp { DECLARE_SOA_COLUMN(IsSelLcK0sP, isSelLcK0sP, int); @@ -50,6 +67,12 @@ DECLARE_SOA_COLUMN(IsSelLcK0sP, isSelLcK0sP, int); DECLARE_SOA_TABLE(HFSelLcK0sPCandidate, "AOD", "HFSELLCK0SPCAND", //! hf_selcandidate_lc_k0sp::IsSelLcK0sP); +namespace hf_selcandidate_x +{ +DECLARE_SOA_COLUMN(IsSelXToJpsiPiPi, isSelXToJpsiPiPi, int); //! +} // namespace hf_selcandidate_x +DECLARE_SOA_TABLE(HFSelXToJpsiPiPiCandidate, "AOD", "HFSELXCAND", //! + hf_selcandidate_x::IsSelXToJpsiPiPi); } // namespace o2::aod namespace o2::aod diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h index 6614b667f54dc..98ea974d9e5b6 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h @@ -65,8 +65,9 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Index3, index3, int, Tracks, "_3"); //! DECLARE_SOA_INDEX_COLUMN_FULL(IndexV0, indexV0, int, aod::V0Datas, "_V0"); //! DECLARE_SOA_COLUMN(HFflag, hfflag, uint8_t); //! -DECLARE_SOA_COLUMN(D0ToKPiFlag, d0ToKPiFlag, uint8_t); //! -DECLARE_SOA_COLUMN(JpsiToEEFlag, jpsiToEEFlag, uint8_t); //! +DECLARE_SOA_COLUMN(D0ToKPiFlag, d0ToKPiFlag, uint8_t); //! +DECLARE_SOA_COLUMN(JpsiToEEFlag, jpsiToEEFlag, uint8_t); //! +DECLARE_SOA_COLUMN(JpsiToMuMuFlag, jpsiToMuMuFlag, uint8_t); //! DECLARE_SOA_COLUMN(DPlusPiKPiFlag, dPlusPiKPiFlag, uint8_t); //! DECLARE_SOA_COLUMN(LcPKPiFlag, lcPKPiFlag, uint8_t); //! @@ -86,7 +87,8 @@ DECLARE_SOA_TABLE(HfTrackIndexCasc, "AOD", "HFTRACKIDXCASC", //! DECLARE_SOA_TABLE(HfCutStatusProng2, "AOD", "HFCUTSTATUSP2", //! hf_track_index::D0ToKPiFlag, - hf_track_index::JpsiToEEFlag); + hf_track_index::JpsiToEEFlag, + hf_track_index::JpsiToMuMuFlag); DECLARE_SOA_TABLE(HfTrackIndexProng3, "AOD", "HFTRACKIDXP3", //! hf_track_index::Index0Id, @@ -227,6 +229,7 @@ DECLARE_SOA_COLUMN(OriginMCGen, originMCGen, int8_t); //! particle origin, // mapping of decay types enum DecayType { D0ToPiK = 0, JpsiToEE, + JpsiToMuMu, N2ProngDecays }; //always keep N2ProngDecays at the end // functions for specific particles @@ -275,7 +278,7 @@ auto CosThetaStarD0bar(const T& candidate) return candidate.cosThetaStar(array{RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}, RecoDecay::getMassPDG(pdg::Code::kD0), 0); } -// J/ψ → e+ e− +// J/ψ template auto CtJpsi(const T& candidate) @@ -295,11 +298,21 @@ auto EJpsi(const T& candidate) return candidate.e(RecoDecay::getMassPDG(pdg::Code::kJpsi)); } +// J/ψ → e+ e− + template auto InvMassJpsiToEE(const T& candidate) { return candidate.m(array{RecoDecay::getMassPDG(kElectron), RecoDecay::getMassPDG(kElectron)}); } + +// J/ψ → μ+ μ− + +template +auto InvMassJpsiToMuMu(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kMuonPlus), RecoDecay::getMassPDG(kMuonMinus)}); +} } // namespace hf_cand_prong2 // general columns @@ -325,6 +338,7 @@ auto InvMassJpsiToEE(const T& candidate) // 2-prong decay candidate table DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE", //! + o2::soa::Index<>, // general columns HFCAND_COLUMNS, // 2-prong specific columns @@ -595,6 +609,33 @@ auto InvMassXicToPiKP(const T& candidate) { return InvMassLcpiKp(candidate); } + +// X → Jpsi π+ π- +// TODO: add pdg code for X (9920443), temporarily hardcode mass here: +float massX = 3.872; // replace this with: "RecoDecay::getMassPDG(9920443)" when pdg is added +template +auto CtX(const T& candidate) +{ + return candidate.ct(massX); +} + +template +auto YX(const T& candidate) +{ + return candidate.y(massX); +} + +template +auto EX(const T& candidate) +{ + return candidate.e(massX); +} +template +auto InvMassXToJpsiPiPi(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(443), RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kPiPlus)}); +} + } // namespace hf_cand_prong3 // 3-prong decay candidate table @@ -653,6 +694,75 @@ DECLARE_SOA_TABLE(HfCandProng3MCGen, "AOD", "HFCANDP3MCGEN", //! hf_cand_prong3::OriginMCGen, hf_cand_prong3::FlagMCDecayChanGen); +// specific X candidate properties +namespace hf_cand_x +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Index0, index0, int, HfCandProng2, "_0"); // Jpsi index +// MC matching result: +DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); // generator level +DECLARE_SOA_COLUMN(OriginMCRec, originMCRec, int8_t); // particle origin, reconstruction level +DECLARE_SOA_COLUMN(OriginMCGen, originMCGen, int8_t); // particle origin, generator level +DECLARE_SOA_COLUMN(FlagMCDecayChanRec, flagMCDecayChanRec, int8_t); // resonant decay channel flag, reconstruction level +DECLARE_SOA_COLUMN(FlagMCDecayChanGen, flagMCDecayChanGen, int8_t); // resonant decay channel flag, generator level +// mapping of decay types +enum DecayType { XToJpsiPiPi = 0 }; // move this to a dedicated cascade namespace in the future? +} // namespace hf_cand_x + +// declare dedicated X candidate table +DECLARE_SOA_TABLE(HfCandXBase, "AOD", "HFCANDXBASE", + // general columns + HFCAND_COLUMNS, + // 3-prong specific columns + hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, + hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, + hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2, + hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2, + hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, hf_cand::ErrorImpactParameter2, + hf_cand_x::Index0Id, hf_track_index::Index1Id, hf_track_index::Index2Id, // note the difference between Jpsi and pion indices + hf_track_index::HFflag, + /* dynamic columns */ + hf_cand_prong3::M, + hf_cand_prong3::M2, + /* prong 2 */ + hf_cand::PtProng2, + hf_cand::Pt2Prong2, + hf_cand::PVectorProng2, + /* dynamic columns that use candidate momentum components */ + hf_cand::Pt, + hf_cand::Pt2, + hf_cand::P, + hf_cand::P2, + hf_cand::PVector, + hf_cand::CPA, + hf_cand::CPAXY, + hf_cand::Ct, + hf_cand::ImpactParameterXY, + hf_cand_prong3::MaxNormalisedDeltaIP, + hf_cand::Eta, + hf_cand::Phi, + hf_cand::Y, + hf_cand::E, + hf_cand::E2); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(HfCandXExt, HfCandXBase, "HFCANDXEXT", + hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz); + +using HfCandX = HfCandXExt; + +// table with results of reconstruction level MC matching +DECLARE_SOA_TABLE(HfCandXMCRec, "AOD", "HFCANDXMCREC", //! + hf_cand_x::FlagMCMatchRec, + hf_cand_x::OriginMCRec, + hf_cand_x::FlagMCDecayChanRec); + +// table with results of generator level MC matching +DECLARE_SOA_TABLE(HfCandXMCGen, "AOD", "HFCANDXMCGEN", //! + hf_cand_x::FlagMCMatchGen, + hf_cand_x::OriginMCGen, + hf_cand_x::FlagMCDecayChanGen); + // definition of columns and tables for D-Dbar correlation pairs namespace hf_correlation_ddbar { diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index 3af51a5d77c50..f2f41075fac9d 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -9,6 +9,11 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +o2_add_dpl_workflow(qa-rejection + SOURCES qaPidRejection.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::ALICE3Analysis + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-track-index-skims-creator SOURCES HFTrackIndexSkimsCreator.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils ROOT::EG @@ -34,6 +39,11 @@ o2_add_dpl_workflow(hf-candidate-creator-3prong PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-candidate-creator-x + SOURCES HFCandidateCreatorX.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-tree-creator-lc-topkpi SOURCES HFTreeCreatorLcToPKPi.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG @@ -54,9 +64,9 @@ o2_add_dpl_workflow(hf-lc-candidate-selector PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) -o2_add_dpl_workflow(hf-jpsi-toee-candidate-selector - SOURCES HFJpsiToEECandidateSelector.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing +o2_add_dpl_workflow(hf-jpsi-candidate-selector + SOURCES HFJpsiCandidateSelector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::ALICE3Analysis COMPONENT_NAME Analysis) o2_add_dpl_workflow(hf-xic-topkpi-candidate-selector @@ -69,6 +79,11 @@ o2_add_dpl_workflow(hf-lc-tok0sp-candidate-selector PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-x-tojpsipipi-candidate-selector + SOURCES HFXToJpsiPiPiCandidateSelector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-task-d0 SOURCES taskD0.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing diff --git a/Analysis/Tasks/PWGHF/HFCandidateCreator2Prong.cxx b/Analysis/Tasks/PWGHF/HFCandidateCreator2Prong.cxx index 0715698e2733e..6508638e5715f 100644 --- a/Analysis/Tasks/PWGHF/HFCandidateCreator2Prong.cxx +++ b/Analysis/Tasks/PWGHF/HFCandidateCreator2Prong.cxx @@ -173,12 +173,21 @@ struct HFCandidateCreator2ProngMC { // J/ψ → e+ e− if (flag == 0) { //Printf("Checking J/ψ → e+ e−"); - indexRec = RecoDecay::getMatchedMCRec(particlesMC, std::move(arrayDaughters), pdg::Code::kJpsi, array{+kElectron, -kElectron}, true); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdg::Code::kJpsi, array{+kElectron, -kElectron}, true); if (indexRec > -1) { flag = 1 << DecayType::JpsiToEE; } } + // J/ψ → μ+ μ− + if (flag == 0) { + //Printf("Checking J/ψ → μ+ μ−"); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdg::Code::kJpsi, array{+kMuonPlus, -kMuonPlus}, true); + if (indexRec > -1) { + flag = 1 << DecayType::JpsiToMuMu; + } + } + // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { auto particle = particlesMC.iteratorAt(indexRec); @@ -208,6 +217,14 @@ struct HFCandidateCreator2ProngMC { } } + // J/ψ → μ+ μ− + if (flag == 0) { + //Printf("Checking J/ψ → μ+ μ−"); + if (RecoDecay::isMatchedMCGen(particlesMC, particle, pdg::Code::kJpsi, array{+kMuonPlus, -kMuonPlus}, true)) { + flag = 1 << DecayType::JpsiToMuMu; + } + } + // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); diff --git a/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx b/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx index 0a0a3ade42814..51c62cfc0433d 100644 --- a/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx +++ b/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx @@ -211,7 +211,7 @@ struct HFCandidateCreator3ProngMC { // Ξc± → p± K∓ π± if (flag == 0) { //Printf("Checking Ξc± → p± K∓ π±"); - indexRec = RecoDecay::getMatchedMCRec(particlesMC, std::move(arrayDaughters), pdg::Code::kXiCPlus, array{+kProton, -kKPlus, +kPiPlus}, true, &sign); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdg::Code::kXiCPlus, array{+kProton, -kKPlus, +kPiPlus}, true, &sign); if (indexRec > -1) { flag = sign * (1 << DecayType::XicToPKPi); } diff --git a/Analysis/Tasks/PWGHF/HFCandidateCreatorX.cxx b/Analysis/Tasks/PWGHF/HFCandidateCreatorX.cxx new file mode 100644 index 0000000000000..0d62fa63f4f9e --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFCandidateCreatorX.cxx @@ -0,0 +1,301 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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 HFCandidateCreatorX.cxx +/// \brief Reconstruction of X(3872) candidates +/// \note Adapted from HFCandidateCreator3Prong +/// +/// \author Rik Spijkers , Utrecht University + +#include "Framework/AnalysisTask.h" +#include "DetectorsVertexing/DCAFitterN.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisCore/trackUtilities.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/V0.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::aod::hf_cand_x; +using namespace o2::framework::expressions; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"doMC", VariantType::Bool, false, {"Perform MC matching."}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +/// Reconstruction of X candidates +struct HFCandidateCreatorX { + Produces rowCandidateBase; + + Configurable magneticField{"magneticField", 5., "magnetic field"}; + Configurable b_propdca{"b_propdca", true, "create tracks version propagated to PCA"}; + Configurable d_maxr{"d_maxr", 200., "reject PCA's above this radius"}; + Configurable d_maxdzini{"d_maxdzini", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable d_minparamchange{"d_minparamchange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable d_minrelchi2change{"d_minrelchi2change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable b_dovalplots{"b_dovalplots", true, "do validation plots"}; + + OutputObj hMassJpsi{TH1F("hMassJpsi", "J/#psi candidates;inv. mass (#e+ e-) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; + OutputObj hPtJpsi{TH1F("hPtJpsi", "J/#psi candidates;candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hCPAJpsi{TH1F("hCPAJpsi", "J/#psi candidates;cosine of pointing angle;entries", 110, -1.1, 1.1)}; + OutputObj hMassX{TH1F("hMassX", "3-prong candidates;inv. mass (J/#psi #pi+ #pi-) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; + OutputObj hCovPVXX{TH1F("hCovPVXX", "3-prong candidates;XX element of cov. matrix of prim. vtx position (cm^{2});entries", 100, 0., 1.e-4)}; + OutputObj hCovSVXX{TH1F("hCovSVXX", "3-prong candidates;XX element of cov. matrix of sec. vtx position (cm^{2});entries", 100, 0., 0.2)}; + + double massPi = RecoDecay::getMassPDG(kPiPlus); + double massJpsi = RecoDecay::getMassPDG(443); + double massJpsiPiPi; + + Configurable d_selectionFlagJpsi{"d_selectionFlagJpsi", 1, "Selection Flag for Jpsi"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= d_selectionFlagJpsi); + + void process(aod::Collision const& collision, + soa::Filtered> const& jpsiCands, + aod::BigTracks const& tracks) + { + // 2-prong vertex fitter (to rebuild Jpsi vertex) + o2::vertexing::DCAFitterN<2> df2; + df2.setBz(magneticField); + df2.setPropagateToPCA(b_propdca); + df2.setMaxR(d_maxr); + df2.setMaxDZIni(d_maxdzini); + df2.setMinParamChange(d_minparamchange); + df2.setMinRelChi2Change(d_minrelchi2change); + df2.setUseAbsDCA(true); + + // 3-prong vertex fitter + o2::vertexing::DCAFitterN<3> df3; + df3.setBz(magneticField); + df3.setPropagateToPCA(b_propdca); + df3.setMaxR(d_maxr); + df3.setMaxDZIni(d_maxdzini); + df3.setMinParamChange(d_minparamchange); + df3.setMinRelChi2Change(d_minrelchi2change); + df3.setUseAbsDCA(true); + + // loop over Jpsi candidates + for (auto& jpsiCand : jpsiCands) { + if (!(jpsiCand.hfflag() & 1 << JpsiToEE)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YJpsi(jpsiCand)) > cutYCandMax) { + continue; + } + hMassJpsi->Fill(InvMassJpsiToEE(jpsiCand)); + hPtJpsi->Fill(jpsiCand.pt()); + hCPAJpsi->Fill(jpsiCand.cpa()); + // create Jpsi track to pass to DCA fitter; use cand table + rebuild vertex + const std::array vertexJpsi = {jpsiCand.xSecondaryVertex(), jpsiCand.ySecondaryVertex(), jpsiCand.zSecondaryVertex()}; + array pvecJpsi = {jpsiCand.px(), jpsiCand.py(), jpsiCand.pz()}; + auto prong0 = jpsiCand.index0_as(); + auto prong1 = jpsiCand.index1_as(); + auto prong0TrackParCov = getTrackParCov(prong0); + auto prong1TrackParCov = getTrackParCov(prong1); + + if (df2.process(prong0TrackParCov, prong1TrackParCov) == 0) { + continue; + } + + // propogate prong tracks to Jpsi vertex + prong0TrackParCov.propagateTo(jpsiCand.xSecondaryVertex(), magneticField); + prong1TrackParCov.propagateTo(jpsiCand.xSecondaryVertex(), magneticField); + const std::array covJpsi = df2.calcPCACovMatrixFlat(); + // define the Jpsi track + auto trackJpsi = o2::dataformats::V0(vertexJpsi, pvecJpsi, covJpsi, prong0TrackParCov, prong1TrackParCov, {0, 0}, {0, 0}); //FIXME: also needs covxyz??? + + // used to check that prongs used for Jpsi and X reco are not the same prongs + int index0Jpsi = jpsiCand.index0Id(); + int index1Jpsi = jpsiCand.index1Id(); + + // loop over pi+ candidates + for (auto& trackPos : tracks) { + if (trackPos.sign() < 0) { // select only positive tracks - use partitions? + continue; + } + if (trackPos.globalIndex() == index0Jpsi) { + continue; + } + + // loop over pi- candidates + for (auto& trackNeg : tracks) { + if (trackNeg.sign() > 0) { // select only negative tracks - use partitions? + continue; + } + if (trackNeg.globalIndex() == index1Jpsi) { + continue; + } + + auto trackParVarPos = getTrackParCov(trackPos); + auto trackParVarNeg = getTrackParCov(trackNeg); + array pvecPos; + array pvecNeg; + + // reconstruct the 3-prong X vertex + if (df3.process(trackJpsi, trackParVarPos, trackParVarNeg) == 0) { + continue; + } + + // calculate relevant properties + const auto& XsecondaryVertex = df3.getPCACandidate(); + auto chi2PCA = df3.getChi2AtPCACandidate(); + auto covMatrixPCA = df3.calcPCACovMatrix().Array(); + hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. + + df3.propagateTracksToVertex(); // propagate the pions and Jpsi to the X vertex + df3.getTrack(0).getPxPyPzGlo(pvecJpsi); // update momentum of Jpsi at the X vertex + df3.getTrack(1).getPxPyPzGlo(pvecPos); // momentum of pi+ at the X vertex + df3.getTrack(2).getPxPyPzGlo(pvecNeg); // momentum of pi- at the X vertex + + // get track impact parameters + // This modifies track momenta! + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + hCovPVXX->Fill(covMatrixPV[0]); + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + o2::dataformats::DCA impactParameter2; + trackJpsi.propagateToDCA(primaryVertex, magneticField, &impactParameter0); + trackParVarPos.propagateToDCA(primaryVertex, magneticField, &impactParameter1); + trackParVarNeg.propagateToDCA(primaryVertex, magneticField, &impactParameter2); + + // get uncertainty of the decay length + double phi, theta; + getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, XsecondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + int hfFlag = 1 << XToJpsiPiPi; + + // fill the candidate table for the X here: + rowCandidateBase(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + XsecondaryVertex[0], XsecondaryVertex[1], XsecondaryVertex[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pvecJpsi[0], pvecJpsi[1], pvecJpsi[2], + pvecPos[0], pvecPos[1], pvecPos[2], + pvecNeg[0], pvecNeg[1], pvecNeg[2], + impactParameter0.getY(), impactParameter1.getY(), impactParameter2.getY(), + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), std::sqrt(impactParameter2.getSigmaY2()), + jpsiCand.globalIndex(), trackPos.globalIndex(), trackNeg.globalIndex(), + hfFlag); + + // calculate invariant mass + auto arrayMomenta = array{pvecJpsi, pvecPos, pvecNeg}; + massJpsiPiPi = RecoDecay::M(std::move(arrayMomenta), array{massJpsi, massPi, massPi}); + hMassX->Fill(massJpsiPiPi); + } // pi- loop + } // pi+ loop + } // Jpsi loop + } // process +}; // struct + +/// Extends the base table with expression columns. +struct HFCandidateCreatorXExpressions { + Spawns rowCandidateX; + void init(InitContext const&) {} +}; + +/// Performs MC matching. +struct HFCandidateCreatorXMC { + Produces rowMCMatchRec; + Produces rowMCMatchGen; + + void process(aod::HfCandX const& candidates, + aod::HfCandProng2, + aod::BigTracksMC const& tracks, + aod::McParticles const& particlesMC) + { + int indexRec = -1; + int8_t sign = 0; + int8_t flag = 0; + int8_t origin = 0; + int8_t channel = 0; + + // Match reconstructed candidates. + for (auto& candidate : candidates) { + //Printf("New rec. candidate"); + flag = 0; + origin = 0; + channel = 0; + auto jpsiTrack = candidate.index0(); + auto arrayDaughters = array{candidate.index1_as(), + candidate.index2_as(), + jpsiTrack.index0_as(), + jpsiTrack.index1_as()}; + + // X → J/ψ π+ π- + //Printf("Checking X → J/ψ π+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, 9920443, array{+kPiPlus, -kPiPlus, +kElectron, -kElectron}, true, &sign, 2); + if (indexRec > -1) { + flag = 1 << XToJpsiPiPi; + } + + // Check whether the particle is non-prompt (from a b quark). + if (flag != 0) { + auto particle = particlesMC.iteratorAt(indexRec); + origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt); + } + + rowMCMatchRec(flag, origin, channel); + } + + // Match generated particles. + for (auto& particle : particlesMC) { + //Printf("New gen. candidate"); + flag = 0; + origin = 0; + channel = 0; + + // X → J/ψ π+ π- + //Printf("Checking X → J/ψ π+ π-"); + if (RecoDecay::isMatchedMCGen(particlesMC, particle, 9920443, array{443, +kPiPlus, -kPiPlus}, true)) { + // Match J/psi --> e+e- + std::vector arrDaughter; + RecoDecay::getDaughters(particlesMC, particle, &arrDaughter, array{443}, 1); + auto jpsiCandMC = particlesMC.iteratorAt(arrDaughter[0]); + if (RecoDecay::isMatchedMCGen(particlesMC, jpsiCandMC, 443, array{+kElectron, -kElectron}, true)) { + flag = 1 << XToJpsiPiPi; + } + } + + // Check whether the particle is non-prompt (from a b quark). + if (flag != 0) { + origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt); + } + + rowMCMatchGen(flag, origin, channel); + } // candidate loop + } // process +}; // struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-x"}), + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-x-expressions"})}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-x-mc"})); + } + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/HFJpsiCandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFJpsiCandidateSelector.cxx new file mode 100644 index 0000000000000..2d4ae314256c2 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFJpsiCandidateSelector.cxx @@ -0,0 +1,334 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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 HFJpsiCandidateSelector.cxx +/// \brief J/ψ → e+ e−, μ+ μ− selection task +/// +/// \author Biao Zhang , CCNU +/// \author Nima Zardoshti , CERN +/// \author Vít Kučera , CERN + +#include "Framework/AnalysisTask.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" +#include "AnalysisCore/TrackSelectorPID.h" +#include "ALICE3Analysis/RICH.h" +#include "ALICE3Analysis/MID.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::analysis::hf_cuts_jpsi_toee; + +namespace o2::aod +{ +namespace hf_track_index_alice3_pid +{ +DECLARE_SOA_INDEX_COLUMN(Track, track); //! +DECLARE_SOA_INDEX_COLUMN(RICH, rich); //! +DECLARE_SOA_INDEX_COLUMN(MID, mid); //! +} // namespace hf_track_index_alice3_pid + +DECLARE_SOA_INDEX_TABLE_USER(HfTrackIndexALICE3PID, Tracks, "HFTRKIDXA3PID", //! + hf_track_index_alice3_pid::TrackId, + hf_track_index_alice3_pid::RICHId, + hf_track_index_alice3_pid::MIDId); +} // namespace o2::aod + +struct Alice3PidIndexBuilder { + Builds index; + void init(o2::framework::InitContext&) {} +}; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec isAlice3{"isAlice3", VariantType::Bool, false, {"Switch between ALICE 2 and ALICE 3 detector setup"}}; + workflowOptions.push_back(isAlice3); +} + +#include "Framework/runDataProcessing.h" + +/// Struct for applying J/ψ → e+ e−, μ+ μ− selection cuts +struct HfJpsiCandidateSelector { + Produces hfSelJpsiCandidate; + + Configurable selectENotPi{"selectENotPi", true, "Apply combined TOF + RICH e/π selection"}; + Configurable d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"}; + Configurable d_pTCandMax{"d_pTCandMax", 50., "Upper bound of candidate pT"}; + // TPC + Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 0.15, "Lower bound of track pT for TPC PID"}; + Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 10., "Upper bound of track pT for TPC PID"}; + Configurable d_nSigmaTPC{"d_nSigmaTPC", 3., "Nsigma cut on TPC only"}; + // TOF + Configurable d_pidTOFMinpT{"d_pidTOFMinpT", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 5., "Upper bound of track pT for TOF PID"}; + Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; + Configurable d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 5., "Nsigma cut on TOF combined with TPC"}; + // RICH + Configurable d_pidRICHMinpT{"d_pidRICHMinpT", 0.15, "Lower bound of track pT for RICH PID"}; + Configurable d_pidRICHMaxpT{"d_pidRICHMaxpT", 10., "Upper bound of track pT for RICH PID"}; + Configurable d_nSigmaRICH{"d_nSigmaRICH", 3., "Nsigma cut on RICH only"}; + Configurable d_nSigmaRICHCombinedTOF{"d_nSigmaRICHCombinedTOF", 5., "Nsigma cut on RICH combined with TOF"}; + // topological cuts + Configurable> pTBins{"pTBins", std::vector{hf_cuts_jpsi_toee::pTBins_v}, "pT bin limits"}; + Configurable> cuts{"Jpsi_to_ee_cuts", {hf_cuts_jpsi_toee::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Jpsi candidate selection per pT bin"}; + + /// Conjugate-independent topological cuts + /// \param candidate is candidate + /// \param trackPos is the positive track + /// \param trackNeg is the negative track + /// \return true if candidate passes all cuts + template + bool selectionTopol(const T1& candidate, const T2& trackPos, const T2& trackNeg, int& selEE, int& selMuMu) + { + auto candpT = candidate.pt(); + auto pTBin = findBin(pTBins, candpT); + if (pTBin == -1) { + return false; + } + + // check that the candidate pT is within the analysis range + if (candpT < d_pTCandMin || candpT >= d_pTCandMax) { + return false; + } + + // cut on e+ e− invariant mass + if (std::abs(InvMassJpsiToEE(candidate) - RecoDecay::getMassPDG(pdg::Code::kJpsi)) > cuts->get(pTBin, "m")) { + selEE = 0; + } + + // cut on μ+ μ− invariant mass + if (std::abs(InvMassJpsiToMuMu(candidate) - RecoDecay::getMassPDG(pdg::Code::kJpsi)) > cuts->get(pTBin, "m")) { + selMuMu = 0; + } + + if (selEE == 0 && selMuMu == 0) { + return false; + } + + // cut on daughter pT (same cut used for both channels) + if (trackNeg.pt() < cuts->get(pTBin, "pT El") || trackPos.pt() < cuts->get(pTBin, "pT El")) { + return false; + } + + // cut on daughter DCA - need to add secondary vertex constraint here + if (std::abs(trackNeg.dcaPrim0()) > cuts->get(pTBin, "DCA_xy") || std::abs(trackPos.dcaPrim0()) > cuts->get(pTBin, "DCA_xy")) { + return false; + } + + // cut on daughter DCA - need to add secondary vertex constraint here + if (std::abs(trackNeg.dcaPrim1()) > cuts->get(pTBin, "DCA_z") || std::abs(trackPos.dcaPrim1()) > cuts->get(pTBin, "DCA_z")) { + return false; + } + + // cut on chi2 point of closest approach + if (std::abs(candidate.chi2PCA()) > cuts->get(pTBin, "chi2PCA")) { + return false; + } + return true; + } + + using TracksPID = soa::Join; + + void processAlice2(aod::HfCandProng2 const& candidates, aod::BigTracksPID const&) + { + TrackSelectorPID selectorElectron(kElectron); + selectorElectron.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); + selectorElectron.setRangeNSigmaTPC(-d_nSigmaTPC, d_nSigmaTPC); + selectorElectron.setRangePtTOF(d_pidTOFMinpT, d_pidTOFMaxpT); + selectorElectron.setRangeNSigmaTOF(-d_nSigmaTOF, d_nSigmaTOF); + selectorElectron.setRangeNSigmaTOFCondTPC(-d_nSigmaTOFCombined, d_nSigmaTOFCombined); + selectorElectron.setRangePtRICH(d_pidRICHMinpT, d_pidRICHMaxpT); + selectorElectron.setRangeNSigmaRICH(-d_nSigmaRICH, d_nSigmaRICH); + selectorElectron.setRangeNSigmaRICHCondTOF(-d_nSigmaRICHCombinedTOF, d_nSigmaRICHCombinedTOF); + + // looping over 2-prong candidates + for (auto& candidate : candidates) { + + if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE) && !(candidate.hfflag() & 1 << DecayType::JpsiToMuMu)) { + hfSelJpsiCandidate(0, 0, 0, 0, 0, 0, 0, 0, 0); + //hfSelJpsiCandidate(0, 0); + continue; + } + + auto trackPos = candidate.index0_as(); // positive daughter + auto trackNeg = candidate.index1_as(); // negative daughter + + int selectedEETopol = 1; + int selectedEETpc = 1; + int selectedEETof = 1; + int selectedEERich = 1; + int selectedEETofRich = 1; + int selectedMuMuTopol = 1; + int selectedMuMuMid = 1; + int selectedEE = 1; + int selectedMuMu = 1; + + // track selection level need to add special cuts (additional cuts on decay length and d0 norm) + + if (!selectionTopol(candidate, trackPos, trackNeg, selectedEETopol, selectedMuMuTopol)) { + selectedEETopol = 0; + selectedMuMuTopol = 0; + selectedEE = 0; + selectedMuMu = 0; + //if (!selectionTopol(candidate, trackPos, trackNeg, selectedEE, selectedMuMu)) { + //hfSelJpsiCandidate(0, 0); + //continue; + } + + // track-level electron PID TOF selection + if (selectorElectron.getStatusTrackPIDTOF(trackPos) == TrackSelectorPID::Status::PIDRejected || + selectorElectron.getStatusTrackPIDTOF(trackNeg) == TrackSelectorPID::Status::PIDRejected) { + selectedEETof = 0; + selectedEE = 0; + //if (selectedMuMu == 0) { + // hfSelJpsiCandidate(0, 0); + // continue; + //} + } + + // track-level electron PID TPC selection + if (selectorElectron.getStatusTrackPIDTPC(trackPos) == TrackSelectorPID::Status::PIDRejected || + selectorElectron.getStatusTrackPIDTPC(trackNeg) == TrackSelectorPID::Status::PIDRejected) { + selectedEETpc = 0; + selectedEE = 0; + } + + hfSelJpsiCandidate(selectedEE, + selectedMuMu, + selectedEETopol, + selectedEETpc, + selectedEETof, + selectedEERich, + selectedEETofRich, + selectedMuMuTopol, + selectedMuMuMid); + //hfSelJpsiCandidate(selectedEE, selectedMuMu); + } + } + + PROCESS_SWITCH(HfJpsiCandidateSelector, processAlice2, "Use ALICE 2 detector setup", true); + + void processAlice3(aod::HfCandProng2 const& candidates, TracksPID const&, aod::RICHs const&, aod::MIDs const&) + { + TrackSelectorPID selectorElectron(kElectron); + selectorElectron.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); + selectorElectron.setRangeNSigmaTPC(-d_nSigmaTPC, d_nSigmaTPC); + selectorElectron.setRangePtTOF(d_pidTOFMinpT, d_pidTOFMaxpT); + selectorElectron.setRangeNSigmaTOF(-d_nSigmaTOF, d_nSigmaTOF); + selectorElectron.setRangeNSigmaTOFCondTPC(-d_nSigmaTOFCombined, d_nSigmaTOFCombined); + selectorElectron.setRangePtRICH(d_pidRICHMinpT, d_pidRICHMaxpT); + selectorElectron.setRangeNSigmaRICH(-d_nSigmaRICH, d_nSigmaRICH); + selectorElectron.setRangeNSigmaRICHCondTOF(-d_nSigmaRICHCombinedTOF, d_nSigmaRICHCombinedTOF); + + TrackSelectorPID selectorMuon(kMuonMinus); + + // looping over 2-prong candidates + for (auto& candidate : candidates) { + + if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE) && !(candidate.hfflag() & 1 << DecayType::JpsiToMuMu)) { + hfSelJpsiCandidate(0, 0, 0, 0, 0, 0, 0, 0, 0); + //hfSelJpsiCandidate(0, 0); + continue; + } + + auto trackPos = candidate.index0_as(); // positive daughter + auto trackNeg = candidate.index1_as(); // negative daughter + + int selectedEETopol = 1; + int selectedEETpc = 1; + int selectedEETof = 1; + int selectedEERich = 1; + int selectedEETofRich = 1; + int selectedMuMuTopol = 1; + int selectedMuMuMid = 1; + int selectedEE = 1; + int selectedMuMu = 1; + + // track selection level need to add special cuts (additional cuts on decay length and d0 norm) + + if (!selectionTopol(candidate, trackPos, trackNeg, selectedEETopol, selectedMuMuTopol)) { + selectedEETopol = 0; + selectedMuMuTopol = 0; + selectedEE = 0; + selectedMuMu = 0; + //if (!selectionTopol(candidate, trackPos, trackNeg, selectedEE, selectedMuMu)) { + //hfSelJpsiCandidate(0, 0); + //continue; + } + + //if (selectENotPi) { + // combined TOF + RICH e selection with π rejection + if (!selectorElectron.isElectronAndNotPion(trackPos) || + !selectorElectron.isElectronAndNotPion(trackNeg)) { + selectedEETofRich = 0; + selectedEE = 0; + } + //} else { + // track-level electron PID TOF selection + if (selectorElectron.getStatusTrackPIDTOF(trackPos) == TrackSelectorPID::Status::PIDRejected || + selectorElectron.getStatusTrackPIDTOF(trackNeg) == TrackSelectorPID::Status::PIDRejected) { + selectedEETof = 0; + selectedEE = 0; + //if (selectedMuMu == 0) { + // hfSelJpsiCandidate(0, 0); + // continue; + //} + } + + // track-level electron PID RICH selection + if (selectorElectron.getStatusTrackPIDRICH(trackPos) == TrackSelectorPID::Status::PIDRejected || + selectorElectron.getStatusTrackPIDRICH(trackNeg) == TrackSelectorPID::Status::PIDRejected) { + selectedEERich = 0; + selectedEE = 0; + } + //} + + //if (selectedEE == 0 && selectedMuMu == 0) { + // hfSelJpsiCandidate(0, 0); + // continue; + //} + + // track-level muon PID MID selection + if (selectorMuon.getStatusTrackPIDMID(trackPos) != TrackSelectorPID::Status::PIDAccepted || + selectorMuon.getStatusTrackPIDMID(trackNeg) != TrackSelectorPID::Status::PIDAccepted) { + selectedMuMuMid = 0; + selectedMuMu = 0; + } + + hfSelJpsiCandidate(selectedEE, + selectedMuMu, + selectedEETopol, + selectedEETpc, + selectedEETof, + selectedEERich, + selectedEETofRich, + selectedMuMuTopol, + selectedMuMuMid); + //hfSelJpsiCandidate(selectedEE, selectedMuMu); + } + } + + PROCESS_SWITCH(HfJpsiCandidateSelector, processAlice3, "Use ALICE 3 detector setup", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + const bool isAlice3 = cfgc.options().get("isAlice3"); + if (isAlice3) { + workflow.push_back(adaptAnalysisTask(cfgc)); + workflow.push_back(adaptAnalysisTask(cfgc, SetDefaultProcesses{{{"processAlice2", false}, {"processAlice3", true}}})); + } else { + workflow.push_back(adaptAnalysisTask(cfgc)); + } + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/HFJpsiToEECandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFJpsiToEECandidateSelector.cxx deleted file mode 100644 index ee6346c92507b..0000000000000 --- a/Analysis/Tasks/PWGHF/HFJpsiToEECandidateSelector.cxx +++ /dev/null @@ -1,151 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// 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 HFJpsiToEECandidateSelector.cxx -/// \brief J/ψ → e+ e− selection task -/// -/// \author Biao Zhang , CCNU -/// \author Nima Zardoshti , CERN -/// \author Vít Kučera , CERN - -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "AnalysisDataModel/HFSecondaryVertex.h" -#include "AnalysisDataModel/HFCandidateSelectionTables.h" -#include "AnalysisCore/TrackSelectorPID.h" - -using namespace o2; -using namespace o2::framework; -using namespace o2::aod::hf_cand_prong2; -using namespace o2::analysis::hf_cuts_jpsi_toee; - -/// Struct for applying J/ψ → e+ e− selection cuts -struct HFJpsiToEECandidateSelector { - Produces hfSelJpsiToEECandidate; - - Configurable d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"}; - Configurable d_pTCandMax{"d_pTCandMax", 50., "Upper bound of candidate pT"}; - // TPC - Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 0.15, "Lower bound of track pT for TPC PID"}; - Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 10., "Upper bound of track pT for TPC PID"}; - Configurable d_nSigmaTPC{"d_nSigmaTPC", 3., "Nsigma cut on TPC only"}; - //Configurable d_TPCNClsFindablePIDCut{"d_TPCNClsFindablePIDCut", 70., "Lower bound of TPC findable clusters for good PID"}; - // topological cuts - Configurable> pTBins{"pTBins", std::vector{hf_cuts_jpsi_toee::pTBins_v}, "pT bin limits"}; - Configurable> cuts{"Jpsi_to_ee_cuts", {hf_cuts_jpsi_toee::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Jpsi candidate selection per pT bin"}; - - /* - /// Selection on goodness of daughter tracks - /// \note should be applied at candidate selection - /// \param track is daughter track - /// \return true if track is good - template - bool daughterSelection(const T& track) - { - if (track.tpcNClsFound() == 0) { - return false; //is it clusters findable or found - need to check - } - return true; - } - */ - - /// Conjugate-independent topological cuts - /// \param candidate is candidate - /// \param trackPositron is the track with the positron hypothesis - /// \param trackElectron is the track with the electron hypothesis - /// \return true if candidate passes all cuts - template - bool selectionTopol(const T1& candidate, const T2& trackPositron, const T2& trackElectron) - { - auto candpT = candidate.pt(); - auto pTBin = findBin(pTBins, candpT); - if (pTBin == -1) { - return false; - } - - // check that the candidate pT is within the analysis range - if (candpT < d_pTCandMin || candpT >= d_pTCandMax) { - return false; - } - - // cut on invariant mass - if (std::abs(InvMassJpsiToEE(candidate) - RecoDecay::getMassPDG(pdg::Code::kJpsi)) > cuts->get(pTBin, "m")) { - return false; - } - - // cut on daughter pT - if (trackElectron.pt() < cuts->get(pTBin, "pT El") || trackPositron.pt() < cuts->get(pTBin, "pT El")) { - return false; - } - - // cut on daughter DCA - need to add secondary vertex constraint here - if (std::abs(trackElectron.dcaPrim0()) > cuts->get(pTBin, "DCA_xy") || std::abs(trackPositron.dcaPrim0()) > cuts->get(pTBin, "DCA_xy")) { - return false; - } - - // cut on daughter DCA - need to add secondary vertex constraint here - if (std::abs(trackElectron.dcaPrim1()) > cuts->get(pTBin, "DCA_z") || std::abs(trackPositron.dcaPrim1()) > cuts->get(pTBin, "DCA_z")) { - return false; - } - - return true; - } - - void process(aod::HfCandProng2 const& candidates, aod::BigTracksPID const&) - { - TrackSelectorPID selectorElectron(kElectron); - selectorElectron.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); - selectorElectron.setRangeNSigmaTPC(-d_nSigmaTPC, d_nSigmaTPC); - - // looping over 2-prong candidates - for (auto& candidate : candidates) { - - auto trackPos = candidate.index0_as(); // positive daughter - auto trackNeg = candidate.index1_as(); // negative daughter - - if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE)) { - hfSelJpsiToEECandidate(0); - continue; - } - - /* - // daughter track validity selection - if (!daughterSelection(trackPos) || !daughterSelection(trackNeg)) { - hfSelJpsiToEECandidate(0); - continue; - } - */ - - // implement filter bit 4 cut - should be done before this task at the track selection level - // need to add special cuts (additional cuts on decay length and d0 norm) - - if (!selectionTopol(candidate, trackPos, trackNeg)) { - hfSelJpsiToEECandidate(0); - continue; - } - - // track-level PID selection - if (selectorElectron.getStatusTrackPIDTPC(trackPos) == TrackSelectorPID::Status::PIDRejected || - selectorElectron.getStatusTrackPIDTPC(trackNeg) == TrackSelectorPID::Status::PIDRejected) { - hfSelJpsiToEECandidate(0); - continue; - } - - hfSelJpsiToEECandidate(1); - } - } -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"hf-jpsi-toee-candidate-selector"})}; -} diff --git a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx index fd2439189db41..46481073e7ea4 100644 --- a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx +++ b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx @@ -48,6 +48,7 @@ static const double massPi = RecoDecay::getMassPDG(kPiPlus); static const double massK = RecoDecay::getMassPDG(kKPlus); static const double massProton = RecoDecay::getMassPDG(kProton); static const double massElectron = RecoDecay::getMassPDG(kElectron); +static const double massMuon = RecoDecay::getMassPDG(kMuonPlus); void customize(std::vector& workflowOptions) { @@ -76,7 +77,6 @@ using MyTracks = soa::Join rowSelectedCollision; Configurable fillHistograms{"fillHistograms", true, "fill histograms"}; @@ -143,7 +143,6 @@ struct HfTagSelCollisions { /// Track selection struct HfTagSelTracks { - Produces rowSelectedTrack; Configurable fillHistograms{"fillHistograms", true, "fill histograms"}; @@ -452,9 +451,12 @@ struct HfTrackIndexSkimsCreator { // D0 cuts Configurable> pTBinsD0ToPiK{"pTBinsD0ToPiK", std::vector{hf_cuts_presel_2prong::pTBinsVec}, "pT bin limits for D0->piK pT-depentend cuts"}; Configurable> cutsD0ToPiK{"cutsD0ToPiK", {hf_cuts_presel_2prong::cuts[0], hf_cuts_presel_2prong::npTBins, hf_cuts_presel_2prong::nCutVars, hf_cuts_presel_2prong::pTBinLabels, hf_cuts_presel_2prong::cutVarLabels}, "D0->piK selections per pT bin"}; - // Jpsi cuts + // Jpsi -> ee cuts Configurable> pTBinsJpsiToEE{"pTBinsJpsiToEE", std::vector{hf_cuts_presel_2prong::pTBinsVec}, "pT bin limits for Jpsi->ee pT-depentend cuts"}; Configurable> cutsJpsiToEE{"cutsJpsiToEE", {hf_cuts_presel_2prong::cuts[0], hf_cuts_presel_2prong::npTBins, hf_cuts_presel_2prong::nCutVars, hf_cuts_presel_2prong::pTBinLabels, hf_cuts_presel_2prong::cutVarLabels}, "Jpsi->ee selections per pT bin"}; + // Jpsi -> mumu cuts + Configurable> pTBinsJpsiToMuMu{"pTBinsJpsiToMuMu", std::vector{hf_cuts_presel_2prong::pTBinsVec}, "pT bin limits for Jpsi->mumu pT-depentend cuts"}; + Configurable> cutsJpsiToMuMu{"cutsJpsiToMuMu", {hf_cuts_presel_2prong::cuts[0], hf_cuts_presel_2prong::npTBins, hf_cuts_presel_2prong::nCutVars, hf_cuts_presel_2prong::pTBinLabels, hf_cuts_presel_2prong::cutVarLabels}, "Jpsi->mumu selections per pT bin"}; // D+ cuts Configurable> pTBinsDPlusToPiKPi{"pTBinsDPlusToPiKPi", std::vector{hf_cuts_presel_3prong::pTBinsVec}, "pT bin limits for D+->piKpi pT-depentend cuts"}; Configurable> cutsDPlusToPiKPi{"cutsDPlusToPiKPi", {hf_cuts_presel_3prong::cuts[0], hf_cuts_presel_3prong::npTBins, hf_cuts_presel_3prong::nCutVars, hf_cuts_presel_3prong::pTBinLabels, hf_cuts_presel_3prong::cutVarLabels}, "D+->piKpi selections per pT bin"}; @@ -479,6 +481,7 @@ struct HfTrackIndexSkimsCreator { {"hNCand2ProngVsNTracks", "2-prong candidates preselected;# of selected tracks;# of candidates;entries", {HistType::kTH2F, {{2500, 0., 25000.}, {2000, 0., 200000.}}}}, {"hmassD0ToPiK", "D^{0} candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, {"hmassJpsiToEE", "J/#psi candidates;inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + {"hmassJpsiToMuMu", "J/#psi candidates;inv. mass (#mu^{#plus} #mu^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, // 3-prong histograms {"hVtx3ProngX", "3-prong candidates;#it{x}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, {"hVtx3ProngY", "3-prong candidates;#it{y}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, @@ -512,6 +515,9 @@ struct HfTrackIndexSkimsCreator { arrMass2Prong[hf_cand_prong2::DecayType::JpsiToEE] = array{array{massElectron, massElectron}, array{massElectron, massElectron}}; + arrMass2Prong[hf_cand_prong2::DecayType::JpsiToMuMu] = array{array{massMuon, massMuon}, + array{massMuon, massMuon}}; + arrMass3Prong[hf_cand_prong3::DecayType::DPlusToPiKPi] = array{array{massPi, massK, massPi}, array{massPi, massK, massPi}}; @@ -525,8 +531,8 @@ struct HfTrackIndexSkimsCreator { array{massPi, massK, massProton}}; // cuts for 2-prong decays retrieved by json. the order must be then one in hf_cand_prong2::DecayType - cut2Prong = {cutsD0ToPiK, cutsJpsiToEE}; - pTBins2Prong = {pTBinsD0ToPiK, pTBinsJpsiToEE}; + cut2Prong = {cutsD0ToPiK, cutsJpsiToEE, cutsJpsiToMuMu}; + pTBins2Prong = {pTBinsD0ToPiK, pTBinsJpsiToEE, pTBinsJpsiToMuMu}; // cuts for 3-prong decays retrieved by json. the order must be then one in hf_cand_prong3::DecayType cut3Prong = {cutsDPlusToPiKPi, cutsLcToPKPi, cutsDsToPiKK, cutsXicToPKPi}; pTBins3Prong = {pTBinsDPlusToPiKPi, pTBinsLcToPKPi, pTBinsDsToPiKK, pTBinsXicToPKPi}; @@ -864,7 +870,7 @@ struct HfTrackIndexSkimsCreator { } } } - rowProng2CutStatus(Prong2CutStatus[0], Prong2CutStatus[1]); //FIXME when we can do this by looping over n2ProngDecays + rowProng2CutStatus(Prong2CutStatus[0], Prong2CutStatus[1], Prong2CutStatus[2]); //FIXME when we can do this by looping over n2ProngDecays } // fill histograms @@ -877,10 +883,16 @@ struct HfTrackIndexSkimsCreator { if (TESTBIT(isSelected2ProngCand, iDecay2P)) { if (whichHypo2Prong[iDecay2P] == 1 || whichHypo2Prong[iDecay2P] == 3) { auto mass2Prong = RecoDecay::M(arrMom, arrMass2Prong[iDecay2P][0]); - if (iDecay2P == hf_cand_prong2::DecayType::D0ToPiK) { - registry.get(HIST("hmassD0ToPiK"))->Fill(mass2Prong); - } else if (iDecay2P == hf_cand_prong2::DecayType::JpsiToEE) { - registry.get(HIST("hmassJpsiToEE"))->Fill(mass2Prong); + switch (iDecay2P) { + case hf_cand_prong2::DecayType::D0ToPiK: + registry.get(HIST("hmassD0ToPiK"))->Fill(mass2Prong); + break; + case hf_cand_prong2::DecayType::JpsiToEE: + registry.get(HIST("hmassJpsiToEE"))->Fill(mass2Prong); + break; + case hf_cand_prong2::DecayType::JpsiToMuMu: + registry.get(HIST("hmassJpsiToMuMu"))->Fill(mass2Prong); + break; } } if (whichHypo2Prong[iDecay2P] >= 2) { diff --git a/Analysis/Tasks/PWGHF/HFXToJpsiPiPiCandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFXToJpsiPiPiCandidateSelector.cxx new file mode 100644 index 0000000000000..a2b2fb5f72ea1 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFXToJpsiPiPiCandidateSelector.cxx @@ -0,0 +1,233 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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 HFXToJpsiPiPiCandidateSelector.cxx +/// \brief X(3872) selection task. +/// \note Adapted from HFJpsiCandidateSelector.cxx +/// \author Rik Spijkers , Utrecht University + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "AnalysisCore/HFSelectorCuts.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_x; +using namespace o2::analysis; +using namespace o2::analysis::hf_cuts_x_tojpsipipi; + +/// Struct for applying Jpsi selection cuts +struct HFXToJpsiPiPiCandidateSelector { + + Produces hfSelXToJpsiPiPiCandidate; + + Configurable d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"}; + Configurable d_pTCandMax{"d_pTCandMax", 50., "Upper bound of candidate pT"}; + + Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 0.15, "Lower bound of track pT for TPC PID"}; + Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 10., "Upper bound of track pT for TPC PID"}; + Configurable d_pidTOFMinpT{"d_pidTOFMinpT", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 10., "Upper bound of track pT for TOF PID"}; + + Configurable d_TPCNClsFindablePIDCut{"d_TPCNClsFindablePIDCut", 70., "Lower bound of TPC findable clusters for good PID"}; + Configurable d_nSigmaTPC{"d_nSigmaTPC", 3., "Nsigma cut on TPC only"}; + Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; + + Configurable> pTBins{"pTBins", std::vector{hf_cuts_x_tojpsipipi::pTBins_v}, "pT bin limits"}; + Configurable> cuts{"X_to_jpsipipi_cuts", {hf_cuts_x_tojpsipipi::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Jpsi candidate selection per pT bin"}; + + /// Selection on goodness of daughter tracks + /// \note should be applied at candidate selection + /// \param track is daughter track + /// \return true if track is good + template + bool daughterSelection(const T& track) + { + return true; + } + + /// Conjugate independent toplogical cuts + /// \param hfCandX is candidate + /// \param trackNeutral is the track with the pi+ hypothesis + /// \param trackPos is the track with the pi+ hypothesis + /// \param trackNeg is the track with the pi- hypothesis + /// \return true if candidate passes all cuts + template + bool selectionTopol(const T1& hfCandX, const T2& hfCandJpsi, const T3& trackPos, const T3& trackNeg) + { + auto candpT = hfCandX.pt(); + int pTBin = findBin(pTBins, candpT); + if (pTBin == -1) { + // Printf("X topol selection failed at getpTBin"); + return false; + } + + if (candpT < d_pTCandMin || candpT >= d_pTCandMax) { + // Printf("X topol selection failed at cand pT check"); + return false; //check that the candidate pT is within the analysis range + } + + // TODO: replace hardcoded mass with "RecoDecay::getMassPDG(9920443)" + if (TMath::Abs(InvMassXToJpsiPiPi(hfCandX) - 3.87168) > cuts->get(pTBin, "m")) { + // Printf("X topol selection failed at mass diff check"); + return false; //check that mass difference is within bounds + } + + if ((hfCandJpsi.pt() < cuts->get(pTBin, "pT Jpsi")) || (trackNeg.pt() < cuts->get(pTBin, "pT Pi")) || (trackPos.pt() < cuts->get(pTBin, "pT Pi"))) { + // Printf("X topol selection failed at daughter pT check"); + return false; //cut on daughter pT + } + + if (hfCandX.cpa() < cuts->get(pTBin, "CPA")) { + return false; // CPA check + } + + if ((TMath::Abs(hfCandX.impactParameter0()) > cuts->get(pTBin, "d0 Jpsi")) || + (TMath::Abs(hfCandX.impactParameter1()) > cuts->get(pTBin, "d0 Pi")) || + (TMath::Abs(hfCandX.impactParameter2()) > cuts->get(pTBin, "d0 Pi"))) { + return false; // DCA check on daughters + } + + // add more cuts: d0 product? PCA? + + return true; + } + + /// Check if track is ok for TPC PID + /// \param track is the track + /// \note function to be expanded + /// \return true if track is ok for TPC PID + template + bool validTPCPID(const T& track) + { + if (TMath::Abs(track.pt()) < d_pidTPCMinpT || TMath::Abs(track.pt()) >= d_pidTPCMaxpT) { + return false; + } + //if (track.TPCNClsFindable() < d_TPCNClsFindablePIDCut) return false; + return true; + } + + /// Check if track is ok for TOF PID + /// \param track is the track + /// \note function to be expanded + /// \return true if track is ok for TOF PID + template + bool validTOFPID(const T& track) + { + if (TMath::Abs(track.pt()) < d_pidTOFMinpT || TMath::Abs(track.pt()) >= d_pidTOFMaxpT) { + return false; + } + return true; + } + + /// Check if track is compatible with given TPC Nsigma cut for the pion hypothesis + /// \param track is the track + /// \param nSigmaCut is the nsigma threshold to test against + /// \return true if track satisfies TPC pion hypothesis for given Nsigma cut + template + bool selectionPIDTPC(const T& track, int nSigmaCut) + { + if (nSigmaCut > 999.) { + return true; + } + return track.tpcNSigmaPi() < nSigmaCut; + } + + /// Check if track is compatible with given TOF NSigma cut for the pion hypothesis + /// \param track is the track + /// \param nSigmaCut is the nSigma threshold to test against + /// \return true if track satisfies TOF pion hypothesis for given NSigma cut + template + bool selectionPIDTOF(const T& track, double nSigmaCut) + { + if (nSigmaCut > 999.) { + return true; + } + return track.tofNSigmaPi() < nSigmaCut; + } + + /// PID selection on daughter track + /// \param track is the daughter track + /// \return 1 if successful PID match, 0 if successful PID rejection, -1 if no PID info + template + int selectionPID(const T& track) + { // use both TPC and TOF here; in run5 only TOF makes sense. add some flag for run3/run5 data later? + if (validTOFPID(track)) { + if (!selectionPIDTOF(track, d_nSigmaTOF)) { + return 0; //rejected by PID + } else { + return 1; //positive PID + } + } else { + return -1; //no PID info + } + + // no tpc in run5, so for now comment it out + // if (validTPCPID(track)) { + // if (!selectionPIDTPC(track, d_nSigmaTPC)) { + // return 0; //rejected by PID + // } else { + // return 1; //positive PID + // } + // } else { + // return -1; //no PID info + // } + } + + void process(aod::HfCandX const& hfCandXs, aod::HfCandProng2, aod::BigTracksPID const& tracks) + { + for (auto& hfCandX : hfCandXs) { //looping over X candidates + // note the difference between Jpsi (index0) and pions (index1,2) + auto candJpsi = hfCandX.index0(); + auto trackPos = hfCandX.index1_as(); //positive daughter + auto trackNeg = hfCandX.index2_as(); //negative daughter + + // check if flagged as X --> Jpsi Pi Pi + if (!(hfCandX.hfflag() & 1 << XToJpsiPiPi)) { + hfSelXToJpsiPiPiCandidate(0); + // Printf("X candidate selection failed at hfflag check"); + continue; + } + + // daughter track validity selection + if (!daughterSelection(trackPos) || !daughterSelection(trackNeg)) { + hfSelXToJpsiPiPiCandidate(0); + // Printf("X candidate selection failed at daughter selection"); + continue; + } + + //implement filter bit 4 cut - should be done before this task at the track selection level + //need to add special cuts (additional cuts on decay length and d0 norm) + + if (!selectionTopol(hfCandX, candJpsi, trackPos, trackNeg)) { + hfSelXToJpsiPiPiCandidate(0); + // Printf("X candidate selection failed at selection topology"); + continue; + } + + if (selectionPID(trackPos) == 0 || selectionPID(trackNeg) == 0) { + hfSelXToJpsiPiPiCandidate(0); + // Printf("X candidate selection failed at selection PID"); + continue; + } + + hfSelXToJpsiPiPiCandidate(1); + // Printf("X candidate selection successful, candidate should be selected"); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"hf-x-tojpsipipi-candidate-selector"})}; +} diff --git a/Analysis/Tasks/PWGHF/qaPidRejection.cxx b/Analysis/Tasks/PWGHF/qaPidRejection.cxx new file mode 100644 index 0000000000000..7a25b5f83ecdc --- /dev/null +++ b/Analysis/Tasks/PWGHF/qaPidRejection.cxx @@ -0,0 +1,523 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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. + +/// \author Peter Hristov , CERN +/// \author Gian Michele Innocenti , CERN +/// \author Henrique J C Zanoli , Utrecht University +/// \author Nicolo' Jacazio , CERN + +// O2 includes +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "ReconstructionDataFormats/DCA.h" +#include "AnalysisCore/trackUtilities.h" +#include "AnalysisCore/MC.h" +#include "AnalysisDataModel/TrackSelectionTables.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisCore/TrackSelectorPID.h" +#include "ALICE3Analysis/RICH.h" +#include "ALICE3Analysis/MID.h" + +using namespace o2; + +using namespace o2::framework; + +namespace o2::aod +{ + +namespace hf_track_index_alice3_pid +{ +DECLARE_SOA_INDEX_COLUMN(Track, track); //! +DECLARE_SOA_INDEX_COLUMN(RICH, rich); //! +DECLARE_SOA_INDEX_COLUMN(MID, mid); //! +} // namespace hf_track_index_alice3_pid + +DECLARE_SOA_INDEX_TABLE_USER(HfTrackIndexALICE3PID, Tracks, "HFTRKIDXA3PID", //! + hf_track_index_alice3_pid::TrackId, + hf_track_index_alice3_pid::RICHId, + hf_track_index_alice3_pid::MIDId); +} // namespace o2::aod + +struct Alice3PidIndexBuilder { + Builds index; + void init(o2::framework::InitContext&) {} +}; + +void customize(std::vector& workflowOptions) +{ + std::vector options{ + {"rej-el", VariantType::Int, 1, {"Efficiency for the Electron PDG code"}}, + {"rej-mu", VariantType::Int, 1, {"Efficiency for the Muon PDG code"}}, + {"rej-pi", VariantType::Int, 1, {"Efficiency for the Pion PDG code"}}, + {"rej-ka", VariantType::Int, 1, {"Efficiency for the Kaon PDG code"}}, + {"rej-pr", VariantType::Int, 0, {"Efficiency for the Proton PDG code"}}}; + std::swap(workflowOptions, options); +} + +#include "Framework/runDataProcessing.h" + +// ROOT includes +#include "TPDGCode.h" +#include "TEfficiency.h" +#include "TList.h" + +/// Task to QA the efficiency of a particular particle defined by particlePDG +template +struct QaTrackingRejection { + static constexpr PDG_t PDGs[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static_assert(particle < 5 && "Maximum of particles reached"); + static constexpr int particlePDG = PDGs[particle]; + // Particle selection + Configurable etaBins{"etaBins", 40, "Number of eta bins"}; + Configurable etaMin{"etaMin", -2.f, "Lower limit in eta"}; + Configurable etaMax{"etaMax", 2.f, "Upper limit in eta"}; + Configurable ptBins{"ptBins", 400, "Number of pT bins"}; + Configurable ptMin{"ptMin", 0.f, "Lower limit in pT"}; + Configurable ptMax{"ptMax", 20.f, "Upper limit in pT"}; + // TPC + Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 9999., "Lower bound of track pT for TPC PID"}; + Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 999999., "Upper bound of track pT for TPC PID"}; + Configurable d_nSigmaTPC{"d_nSigmaTPC", 99999, "Nsigma cut on TPC only"}; + // TOF + Configurable d_pidTOFMinpT{"d_pidTOFMinpT", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 5., "Upper bound of track pT for TOF PID"}; + Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; + Configurable d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 0., "Nsigma cut on TOF combined with TPC"}; + // RICH + Configurable d_pidRICHMinpT{"d_pidRICHMinpT", 0.15, "Lower bound of track pT for RICH PID"}; + Configurable d_pidRICHMaxpT{"d_pidRICHMaxpT", 10., "Upper bound of track pT for RICH PID"}; + Configurable d_nSigmaRICH{"d_nSigmaRICH", 3., "Nsigma cut on RICH only"}; + Configurable d_nSigmaRICHCombinedTOF{"d_nSigmaRICHCombinedTOF", 0., "Nsigma cut on RICH combined with TOF"}; + HistogramRegistry histos{"HistogramsRejection"}; + + void init(InitContext&) + { + AxisSpec ptAxis{ptBins, ptMin, ptMax}; + AxisSpec etaAxis{etaBins, etaMin, etaMax}; + + TString commonTitle = ""; + if (particlePDG != 0) { + commonTitle += Form("PDG %i", particlePDG); + } + + const TString pt = "#it{p}_{T} [GeV/#it{c}]"; + const TString p = "#it{p} [GeV/#it{c}]"; + const TString eta = "#it{#eta}"; + const TString phi = "#it{#varphi} [rad]"; + + histos.add("tracking/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselElectron/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselPion/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselKaon/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselProton/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselElectron/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselPion/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselKaon/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselProton/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingMIDselMuon/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("tracking/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselElectron/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselPion/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselKaon/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingTOFselProton/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselElectron/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselPion/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselKaon/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingRICHselProton/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("trackingMIDselMuon/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + } + + using TracksPID = soa::Join; + + void process(const o2::soa::Join& collisions, + const o2::soa::Join& tracks, + const o2::aod::McCollisions& mcCollisions, + const o2::aod::McParticles& mcParticles, aod::RICHs const&, aod::MIDs const&) + { + TrackSelectorPID selectorElectron(kElectron); + selectorElectron.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); + selectorElectron.setRangeNSigmaTPC(-d_nSigmaTPC, d_nSigmaTPC); + selectorElectron.setRangePtTOF(d_pidTOFMinpT, d_pidTOFMaxpT); + selectorElectron.setRangeNSigmaTOF(-d_nSigmaTOF, d_nSigmaTOF); + selectorElectron.setRangeNSigmaTOFCondTPC(-d_nSigmaTOFCombined, d_nSigmaTOFCombined); + selectorElectron.setRangePtRICH(d_pidRICHMinpT, d_pidRICHMaxpT); + selectorElectron.setRangeNSigmaRICH(-d_nSigmaRICH, d_nSigmaRICH); + selectorElectron.setRangeNSigmaRICHCondTOF(-d_nSigmaRICHCombinedTOF, d_nSigmaRICHCombinedTOF); + + auto selectorPion(selectorElectron); + selectorPion.setPDG(kPiPlus); + + auto selectorKaon(selectorElectron); + selectorKaon.setPDG(kKPlus); + + auto selectorProton(selectorElectron); + selectorProton.setPDG(kProton); + + TrackSelectorPID selectorMuon(kMuonPlus); + + std::vector recoEvt(collisions.size()); + std::vector recoTracks(tracks.size()); + LOGF(info, "%d", particlePDG); + for (const auto& track : tracks) { + const auto mcParticle = track.mcParticle(); + if (particlePDG != 0 && mcParticle.pdgCode() != particlePDG) { // Checking PDG code + continue; + } + bool isTOFhpElectron = !(selectorElectron.getStatusTrackPIDTOF(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpElectron = !(selectorElectron.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isTOFhpPion = !(selectorPion.getStatusTrackPIDTOF(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpPion = !(selectorPion.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isTOFhpKaon = !(selectorKaon.getStatusTrackPIDTOF(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpKaon = !(selectorKaon.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isTOFhpProton = !(selectorProton.getStatusTrackPIDTOF(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpProton = !(selectorProton.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isMIDhpMuon = (selectorMuon.getStatusTrackPIDMID(track) == TrackSelectorPID::Status::PIDAccepted); + + if (MC::isPhysicalPrimary(mcParticle)) { + histos.fill(HIST("tracking/pteta"), track.pt(), track.eta()); + if (isTOFhpElectron) { + histos.fill(HIST("trackingTOFselElectron/pteta"), track.pt(), track.eta()); + } + if (isTOFhpPion) { + histos.fill(HIST("trackingTOFselPion/pteta"), track.pt(), track.eta()); + } + if (isTOFhpKaon) { + histos.fill(HIST("trackingTOFselKaon/pteta"), track.pt(), track.eta()); + } + if (isTOFhpProton) { + histos.fill(HIST("trackingTOFselProton/pteta"), track.pt(), track.eta()); + } + if (isRICHhpElectron) { + histos.fill(HIST("trackingRICHselElectron/pteta"), track.pt(), track.eta()); + } + if (isRICHhpPion) { + histos.fill(HIST("trackingRICHselPion/pteta"), track.pt(), track.eta()); + } + if (isRICHhpKaon) { + histos.fill(HIST("trackingRICHselKaon/pteta"), track.pt(), track.eta()); + } + if (isRICHhpProton) { + histos.fill(HIST("trackingRICHselProton/pteta"), track.pt(), track.eta()); + } + if (isMIDhpMuon) { + histos.fill(HIST("trackingMIDselMuon/pteta"), track.pt(), track.eta()); + } + + histos.fill(HIST("tracking/peta"), track.p(), track.eta()); + if (isTOFhpElectron) { + histos.fill(HIST("trackingTOFselElectron/peta"), track.p(), track.eta()); + } + if (isTOFhpPion) { + histos.fill(HIST("trackingTOFselPion/peta"), track.p(), track.eta()); + } + if (isTOFhpKaon) { + histos.fill(HIST("trackingTOFselKaon/peta"), track.p(), track.eta()); + } + if (isTOFhpProton) { + histos.fill(HIST("trackingTOFselProton/peta"), track.p(), track.eta()); + } + if (isRICHhpElectron) { + histos.fill(HIST("trackingRICHselElectron/peta"), track.p(), track.eta()); + } + if (isRICHhpPion) { + histos.fill(HIST("trackingRICHselPion/peta"), track.p(), track.eta()); + } + if (isRICHhpKaon) { + histos.fill(HIST("trackingRICHselKaon/peta"), track.p(), track.eta()); + } + if (isRICHhpProton) { + histos.fill(HIST("trackingRICHselProton/peta"), track.p(), track.eta()); + } + if (isMIDhpMuon) { + histos.fill(HIST("trackingMIDselMuon/peta"), track.p(), track.eta()); + } + } + } + } +}; + +struct QaRejectionGeneral { + static constexpr PDG_t PDGs[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + // Cuts + Configurable etaMaxSel{"etaMaxSel", 1.44, "Max #eta single track"}; + Configurable ptMinSel{"ptMinSel", 0.6, "p_{T} min single track"}; + // Particle selection + Configurable etaBins{"etaBins", 40, "Number of eta bins"}; + Configurable etaMin{"etaMin", -2.f, "Lower limit in eta"}; + Configurable etaMax{"etaMax", 2.f, "Upper limit in eta"}; + Configurable ptBins{"ptBins", 400, "Number of pT bins"}; + Configurable ptMin{"ptMin", 0.f, "Lower limit in pT"}; + Configurable ptMax{"ptMax", 20.f, "Upper limit in pT"}; + // TPC + Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 9999., "Lower bound of track pT for TPC PID"}; + Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 999999., "Upper bound of track pT for TPC PID"}; + Configurable d_nSigmaTPC{"d_nSigmaTPC", 99999, "Nsigma cut on TPC only"}; + // TOF + Configurable d_pidTOFMinpT{"d_pidTOFMinpT", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 5., "Upper bound of track pT for TOF PID"}; + Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; + Configurable d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 0., "Nsigma cut on TOF combined with TPC"}; + // RICH + Configurable d_pidRICHMinpT{"d_pidRICHMinpT", 0.15, "Lower bound of track pT for RICH PID"}; + Configurable d_pidRICHMaxpT{"d_pidRICHMaxpT", 10., "Upper bound of track pT for RICH PID"}; + Configurable d_nSigmaRICH{"d_nSigmaRICH", 3., "Nsigma cut on RICH only"}; + Configurable d_nSigmaRICHCombinedTOF{"d_nSigmaRICHCombinedTOF", 0., "Nsigma cut on RICH combined with TOF"}; + HistogramRegistry histos{"HistogramsRejection"}; + + void init(InitContext&) + { + AxisSpec ptAxis{ptBins, ptMin, ptMax}; + AxisSpec etaAxis{etaBins, etaMin, etaMax}; + + TString commonTitle = ""; + + const TString pt = "#it{p}_{T} [GeV/#it{c}]"; + const TString p = "#it{p} [GeV/#it{c}]"; + const TString eta = "#it{#eta}"; + const TString phi = "#it{#varphi} [rad]"; + + histos.add("hAllNoSel/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionNoSel/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronNoSel/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonNoSel/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonNoSel/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllRICHSelHpElTight/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionRICHSelHpElTight/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronRICHSelHpElTight/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonRICHSelHpElTight/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonRICHSelHpElTight/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllRICHSelHpElTightAlt/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionRICHSelHpElTightAlt/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronRICHSelHpElTightAlt/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonRICHSelHpElTightAlt/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonRICHSelHpElTightAlt/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllRICHSelHpElTightAltDiff/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionRICHSelHpElTightAltDiff/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronRICHSelHpElTightAltDiff/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonRICHSelHpElTightAltDiff/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonRICHSelHpElTightAltDiff/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllRICHSelHpElLoose/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionRICHSelHpElLoose/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronRICHSelHpElLoose/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonRICHSelHpElLoose/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonRICHSelHpElLoose/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllMID/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronMID/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonMID/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionMID/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonMID/pteta", commonTitle + " Primary;" + pt, kTH2D, {ptAxis, etaAxis}); + + histos.add("hAllMID/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("hElectronMID/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("hMuonMID/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("hPionMID/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + histos.add("hKaonMID/peta", commonTitle + " Primary;" + p, kTH2D, {ptAxis, etaAxis}); + } + + using TracksPID = soa::Join; + + void process(const o2::soa::Join& collisions, + const o2::soa::Join& tracks, + const o2::aod::McCollisions& mcCollisions, + const o2::aod::McParticles& mcParticles, aod::RICHs const&, aod::MIDs const&) + { + TrackSelectorPID selectorElectron(kElectron); + selectorElectron.setRangePtTPC(d_pidTPCMinpT, d_pidTPCMaxpT); + selectorElectron.setRangeNSigmaTPC(-d_nSigmaTPC, d_nSigmaTPC); + selectorElectron.setRangePtTOF(d_pidTOFMinpT, d_pidTOFMaxpT); + selectorElectron.setRangeNSigmaTOF(-d_nSigmaTOF, d_nSigmaTOF); + selectorElectron.setRangeNSigmaTOFCondTPC(-d_nSigmaTOFCombined, d_nSigmaTOFCombined); + selectorElectron.setRangePtRICH(d_pidRICHMinpT, d_pidRICHMaxpT); + selectorElectron.setRangeNSigmaRICH(-d_nSigmaRICH, d_nSigmaRICH); + selectorElectron.setRangeNSigmaRICHCondTOF(-d_nSigmaRICHCombinedTOF, d_nSigmaRICHCombinedTOF); + + auto selectorPion(selectorElectron); + selectorPion.setPDG(kPiPlus); + + auto selectorKaon(selectorElectron); + selectorKaon.setPDG(kKPlus); + + auto selectorProton(selectorElectron); + selectorProton.setPDG(kProton); + + TrackSelectorPID selectorMuon(kMuonPlus); + + for (const auto& track : tracks) { + + if (std::abs(track.eta()) > etaMaxSel || track.pt() < ptMinSel) { + continue; + } + const auto mcParticle = track.mcParticle(); + histos.fill(HIST("hAllNoSel/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronNoSel/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonNoSel/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionNoSel/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonNoSel/pteta"), track.pt(), track.eta()); + } + + bool isRICHhpElectron = !(selectorElectron.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpPion = !(selectorPion.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpKaon = !(selectorKaon.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isRICHhpProton = !(selectorProton.getStatusTrackPIDRICH(track) == TrackSelectorPID::Status::PIDRejected); + bool isMIDhpMuon = (selectorMuon.getStatusTrackPIDMID(track) == TrackSelectorPID::Status::PIDAccepted); + + bool isRICHElLoose = isRICHhpElectron; + + bool isRICHElTight = true; + bool isrichel = true; + bool istofel = true; + bool istofpi = false; + bool isrichpi = false; + + if (track.p() < 0.6) { + istofel = std::abs(track.tofNSigmaEl()) < d_pidTOFMaxpT; + istofpi = std::abs(track.tofNSigmaPi()) < d_pidTOFMaxpT; + isrichel = std::abs(track.rich().richNsigmaEl()) < d_pidRICHMaxpT; + } else if (track.p() > 0.6 && track.p() < 1.0) { + isrichel = std::abs(track.rich().richNsigmaEl()) < d_pidRICHMaxpT; + } else if (track.p() > 1.0 && track.p() < 2.0) { + isrichel = std::abs(track.rich().richNsigmaEl()) < d_pidRICHMaxpT; + isrichpi = std::abs(track.rich().richNsigmaPi()) < d_pidRICHMaxpT; + } else { + isrichel = std::abs(track.rich().richNsigmaEl()) < d_pidRICHMaxpT; + } + isRICHElTight = isrichel && !isrichpi && istofel && !istofpi; + + auto isRICHElTightAlt = selectorElectron.isElectronAndNotPion(track); + + if (isRICHElLoose) { + histos.fill(HIST("hAllRICHSelHpElLoose/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronRICHSelHpElLoose/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonRICHSelHpElLoose/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionRICHSelHpElLoose/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonRICHSelHpElLoose/pteta"), track.pt(), track.eta()); + } + } + if (isRICHElTight) { + histos.fill(HIST("hAllRICHSelHpElTight/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronRICHSelHpElTight/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonRICHSelHpElTight/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionRICHSelHpElTight/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonRICHSelHpElTight/pteta"), track.pt(), track.eta()); + } + } + if (isRICHElTightAlt) { + histos.fill(HIST("hAllRICHSelHpElTightAlt/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronRICHSelHpElTightAlt/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonRICHSelHpElTightAlt/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionRICHSelHpElTightAlt/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonRICHSelHpElTightAlt/pteta"), track.pt(), track.eta()); + } + } + if (isRICHElTightAlt != isRICHElTight) { + histos.fill(HIST("hAllRICHSelHpElTightAltDiff/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronRICHSelHpElTightAltDiff/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonRICHSelHpElTightAltDiff/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionRICHSelHpElTightAltDiff/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonRICHSelHpElTightAltDiff/pteta"), track.pt(), track.eta()); + } + } + if (isMIDhpMuon) { + histos.fill(HIST("hAllMID/pteta"), track.pt(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronMID/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonMID/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionMID/pteta"), track.pt(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonMID/pteta"), track.pt(), track.eta()); + } + + histos.fill(HIST("hAllMID/peta"), track.p(), track.eta()); + if (mcParticle.pdgCode() == kElectron) { + histos.fill(HIST("hElectronMID/peta"), track.p(), track.eta()); + } + if (mcParticle.pdgCode() == kMuonPlus) { + histos.fill(HIST("hMuonMID/peta"), track.p(), track.eta()); + } + if (mcParticle.pdgCode() == kPiPlus) { + histos.fill(HIST("hPionMID/peta"), track.p(), track.eta()); + } + if (mcParticle.pdgCode() == kKPlus) { + histos.fill(HIST("hKaonMID/peta"), track.p(), track.eta()); + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec w; + w.push_back(adaptAnalysisTask(cfgc)); + if (cfgc.options().get("rej-el")) { + w.push_back(adaptAnalysisTask>(cfgc, TaskName{"qa-tracking-rejection-electron"})); + } + if (cfgc.options().get("rej-ka")) { + w.push_back(adaptAnalysisTask>(cfgc, TaskName{"qa-tracking-rejection-kaon"})); + } + if (cfgc.options().get("rej-pr")) { + w.push_back(adaptAnalysisTask>(cfgc, TaskName{"qa-tracking-rejection-proton"})); + } + if (cfgc.options().get("rej-mu")) { + w.push_back(adaptAnalysisTask>(cfgc, TaskName{"qa-tracking-rejection-mu"})); + } + if (cfgc.options().get("rej-pi")) { + w.push_back(adaptAnalysisTask>(cfgc, TaskName{"qa-tracking-rejection-pion"})); + } + w.push_back(adaptAnalysisTask(cfgc)); + return w; +} diff --git a/Analysis/Tasks/PWGHF/taskJpsi.cxx b/Analysis/Tasks/PWGHF/taskJpsi.cxx index c7434e554c0c1..7b1a8070f0592 100644 --- a/Analysis/Tasks/PWGHF/taskJpsi.cxx +++ b/Analysis/Tasks/PWGHF/taskJpsi.cxx @@ -39,42 +39,65 @@ void customize(std::vector& workflowOptions) struct TaskJpsi { HistogramRegistry registry{ "registry", - {{"hptcand", "2-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, - {"hptprong0", "2-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, - {"hptprong1", "2-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}}}; + {{"hptcand", "2-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, + {"hptprong0", "2-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, + {"hptprong1", "2-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}}}; - Configurable d_selectionFlagJpsi{"d_selectionFlagJpsi", 1, "Selection Flag for Jpsi"}; + Configurable d_selectionFlagJpsi{"d_selectionFlagJpsi", 0, "Selection Flag for Jpsi"}; + Configurable d_modeJpsiToMuMu{"d_modeJpsiToMuMu", false, "Perform Jpsi to mu+mu- analysis"}; Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; Configurable> bins{"pTBins", std::vector{hf_cuts_jpsi_toee::pTBins_v}, "pT bin limits"}; void init(o2::framework::InitContext&) { - registry.add("hmass", "2-prong candidates;inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0., 5.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hdeclength", "2-prong candidates;decay length (cm);entries", {HistType::kTH2F, {{200, 0., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hdeclengthxy", "2-prong candidates;decay length xy (cm);entries", {HistType::kTH2F, {{200, 0., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hd0Prong0", "2-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hd0Prong1", "2-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hd0d0", "2-prong candidates;product of DCAxy to prim. vertex (cm^{2});entries", {HistType::kTH2F, {{500, -1., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hCPA", "2-prong candidates;cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + if (d_modeJpsiToMuMu) { + registry.add("hmass", "2-prong candidates;inv. mass (#mu^{#plus} #mu^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } else { + registry.add("hmass", "2-prong candidates;inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + registry.add("hdeclength", "2-prong candidates;decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclengthxy", "2-prong candidates;decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0", "2-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1", "2-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0d0", "2-prong candidates;product of DCAxy to prim. vertex (cm^{2});entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPA", "2-prong candidates;cosine of pointing angle;entries", {HistType::kTH2F, {{1000, 0.5, 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); registry.add("hEta", "2-prong candidates;candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hImpParErr", "2-prong candidates;impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hDecLenErr", "2-prong candidates;decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hDecLenXYErr", "2-prong candidates;decay length xy error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hImpParErr", "2-prong candidates;impact parameter error (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenErr", "2-prong candidates;decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenXYErr", "2-prong candidates;decay length xy error (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); } - Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= d_selectionFlagJpsi); + Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEETopol >= d_selectionFlagJpsi || aod::hf_selcandidate_jpsi::isSelJpsiToMuMuTopol >= d_selectionFlagJpsi); - void process(soa::Filtered> const& candidates) + void process(soa::Filtered> const& candidates) { + int decayMode = d_modeJpsiToMuMu ? DecayType::JpsiToMuMu : DecayType::JpsiToEE; + for (auto& candidate : candidates) { - if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE)) { + + if (!(candidate.hfflag() & 1 << decayMode)) { continue; } + if (d_selectionFlagJpsi > 0) { + if (d_modeJpsiToMuMu) { + if (candidate.isSelJpsiToMuMuTopol() <= 0) { + continue; + } + } else { + if (candidate.isSelJpsiToEETopol() <= 0) { + continue; + } + } + } if (cutYCandMax >= 0. && std::abs(YJpsi(candidate)) > cutYCandMax) { continue; } - registry.fill(HIST("hmass"), InvMassJpsiToEE(candidate), candidate.pt()); + if (d_modeJpsiToMuMu) { + registry.fill(HIST("hmass"), InvMassJpsiToMuMu(candidate), candidate.pt()); + } else { + registry.fill(HIST("hmass"), InvMassJpsiToEE(candidate), candidate.pt()); + } registry.fill(HIST("hptcand"), candidate.pt()); registry.fill(HIST("hptprong0"), candidate.ptProng0()); registry.fill(HIST("hptprong1"), candidate.ptProng1()); @@ -97,10 +120,10 @@ struct TaskJpsi { struct TaskJpsiMC { HistogramRegistry registry{ "registry", - {{"hPtRecSig", "2-prong candidates (rec. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, - {"hPtRecBg", "2-prong candidates (rec. unmatched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, - {"hPtGen", "2-prong candidates (gen. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, - {"hPtGenSig", "2-prong candidates (rec. matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {{"hPtRecSig", "2-prong candidates (rec. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, + {"hPtRecBg", "2-prong candidates (rec. unmatched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, + {"hPtGen", "2-prong candidates (gen. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, + {"hPtGenSig", "2-prong candidates (rec. matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 20.}}}}, {"hCPARecSig", "2-prong candidates (rec. matched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, {"hCPARecBg", "2-prong candidates (rec. unmatched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, {"hEtaRecSig", "2-prong candidates (rec. matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, @@ -108,54 +131,126 @@ struct TaskJpsiMC { {"hEtaGen", "2-prong candidates (gen. matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}}}; Configurable d_selectionFlagJpsi{"d_selectionFlagJpsi", 1, "Selection Flag for Jpsi"}; + Configurable d_modeJpsiToMuMu{"d_modeJpsiToMuMu", false, "Perform Jpsi to mu+mu- analysis"}; Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; Configurable> bins{"pTBins", std::vector{hf_cuts_jpsi_toee::pTBins_v}, "pT bin limits"}; void init(o2::framework::InitContext&) { - registry.add("hmassSig", "2-prong candidates (rec matched);inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0., 5.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hmassBkg", "2-prong candidates (rec unmatched);inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0., 5.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + if (d_modeJpsiToMuMu) { + registry.add("hmassSig", "2-prong candidates (rec matched);inv. mass (#mu^{#plus} #mu^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hmassBg", "2-prong candidates (rec unmatched);inv. mass (#mu^{#plus} #mu^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } else { + registry.add("hmassSig", "2-prong candidates (rec matched);inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hmassBg", "2-prong candidates (rec unmatched);inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + registry.add("hdeclengthSig", "2-prong candidates (rec matched);decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclengthBg", "2-prong candidates (rec unmatched);decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclengthxySig", "2-prong candidates (rec matched);decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclengthxyBg", "2-prong candidates (rec unmatched);decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0Sig", "2-prong candidates (rec matched);prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0Bg", "2-prong candidates (rec unmatched);prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1Sig", "2-prong candidates (rec matched);prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1Bg", "2-prong candidates (rec unmatched);prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0d0Sig", "2-prong candidates (rec matched);product of DCAxy to prim. vertex (cm^{2});entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0d0Bg", "2-prong candidates (rec unmatched);product of DCAxy to prim. vertex (cm^{2});entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hChi2PCASig", "2-prong candidates (rec. matched);chi2 PCA (cm);entries", {HistType::kTH2F, {{1000, 0., 0.0001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hChi2PCABg", "2-prong candidates (rec. unmatched);chi2 PCA (cm);entries", {HistType::kTH2F, {{1000, 0., 0.0001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCtSig", "2-prong candidates (rec. matched);proper lifetime X(3872) * #it{c} (cm);entries", {HistType::kTH2F, {{400, 0., 0.001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCtBg", "2-prong candidates (rec. unmatched);proper lifetime X(3872) * #it{c} (cm);entries", {HistType::kTH2F, {{400, 0., 0.001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYGenSig", "2-prong candidates (rec. matched);candidate rapidity;entries", {HistType::kTH2F, {{10, -2., 2.}, {(std::vector)bins, "#it{p}^{gen}_{T} (GeV/#it{c})"}}}); + registry.add("hYSig", "2-prong candidates (rec. matched);candidate rapidity;entries", {HistType::kTH2F, {{10, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYBg", "2-prong candidates (rec. unmatched);candidate rapidity;entries", {HistType::kTH2F, {{10, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYGen", "2-prong MC particles (gen. matched);candidate rapidity;entries", {HistType::kTH2F, {{10, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtGenProng0", "2-prong candidates (gen. matched);prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtGenProng1", "2-prong candidates (gen. matched);prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); } - Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= d_selectionFlagJpsi); + Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEETopol >= d_selectionFlagJpsi || aod::hf_selcandidate_jpsi::isSelJpsiToMuMuTopol >= d_selectionFlagJpsi); - void process(soa::Filtered> const& candidates, - soa::Join const& particlesMC, aod::BigTracksMC const& tracks) + using McParticlesHf = soa::Join; + + void process(soa::Filtered> const& candidates, + McParticlesHf const& particlesMC, aod::BigTracksMC const& tracks) { // MC rec. //Printf("MC Candidates: %d", candidates.size()); + int decayMode = d_modeJpsiToMuMu ? DecayType::JpsiToMuMu : DecayType::JpsiToEE; + for (auto& candidate : candidates) { - if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE)) { + + if (!(candidate.hfflag() & 1 << decayMode)) { continue; } + if (d_selectionFlagJpsi > 0) { + if (d_modeJpsiToMuMu) { + if (candidate.isSelJpsiToMuMuTopol() <= 0) { + continue; + } + } else { + if (candidate.isSelJpsiToEETopol() <= 0) { + continue; + } + } + } + if (cutYCandMax >= 0. && std::abs(YJpsi(candidate)) > cutYCandMax) { continue; } - if (candidate.flagMCMatchRec() == 1 << DecayType::JpsiToEE) { + if (candidate.flagMCMatchRec() == 1 << decayMode) { //Get the corresponding MC particle. - auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as().mcParticle_as>(), 443, true); + auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as().mcParticle_as(), pdg::Code::kJpsi, true); auto particleMother = particlesMC.iteratorAt(indexMother); registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT registry.fill(HIST("hPtRecSig"), candidate.pt()); // rec. level pT registry.fill(HIST("hCPARecSig"), candidate.cpa()); registry.fill(HIST("hEtaRecSig"), candidate.eta()); - registry.fill(HIST("hmassSig"), InvMassJpsiToEE(candidate), candidate.pt()); + if (d_modeJpsiToMuMu) { + registry.fill(HIST("hmassSig"), InvMassJpsiToMuMu(candidate), candidate.pt()); + } else { + registry.fill(HIST("hmassSig"), InvMassJpsiToEE(candidate), candidate.pt()); + } + registry.fill(HIST("hdeclengthSig"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hdeclengthxySig"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hd0Prong0Sig"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1Sig"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hd0d0Sig"), candidate.impactParameterProduct(), candidate.pt()); + registry.fill(HIST("hChi2PCASig"), candidate.chi2PCA(), candidate.pt()); + registry.fill(HIST("hCtSig"), CtJpsi(candidate), candidate.pt()); + registry.fill(HIST("hYSig"), YJpsi(candidate), candidate.pt()); + registry.fill(HIST("hYGenSig"), RecoDecay::Y(array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::getMassPDG(particleMother.pdgCode())), particleMother.pt()); + } else { registry.fill(HIST("hPtRecBg"), candidate.pt()); registry.fill(HIST("hCPARecBg"), candidate.cpa()); registry.fill(HIST("hEtaRecBg"), candidate.eta()); - registry.fill(HIST("hmassBkg"), InvMassJpsiToEE(candidate), candidate.pt()); + if (d_modeJpsiToMuMu) { + registry.fill(HIST("hmassBg"), InvMassJpsiToMuMu(candidate), candidate.pt()); + } else { + registry.fill(HIST("hmassBg"), InvMassJpsiToEE(candidate), candidate.pt()); + } + registry.fill(HIST("hdeclengthBg"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hdeclengthxyBg"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hd0Prong0Bg"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1Bg"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hd0d0Bg"), candidate.impactParameterProduct(), candidate.pt()); + registry.fill(HIST("hChi2PCABg"), candidate.chi2PCA(), candidate.pt()); + registry.fill(HIST("hCtBg"), CtJpsi(candidate), candidate.pt()); + registry.fill(HIST("hYBg"), YJpsi(candidate), candidate.pt()); } } // MC gen. //Printf("MC Particles: %d", particlesMC.size()); for (auto& particle : particlesMC) { - if (particle.flagMCMatchGen() == 1 << DecayType::JpsiToEE) { + if (particle.flagMCMatchGen() == 1 << decayMode) { if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle.px(), particle.py(), particle.pz()}, RecoDecay::getMassPDG(particle.pdgCode()))) > cutYCandMax) { continue; } registry.fill(HIST("hPtGen"), particle.pt()); registry.fill(HIST("hEtaGen"), particle.eta()); + registry.fill(HIST("hYGen"), RecoDecay::Y(array{particle.px(), particle.py(), particle.pz()}, RecoDecay::getMassPDG(particle.pdgCode())), particle.pt()); + registry.fill(HIST("hPtGenProng0"), particle.daughter0_as().pt(), particle.pt()); + registry.fill(HIST("hPtGenProng1"), particle.daughter1_as().pt(), particle.pt()); } } } diff --git a/Analysis/Tasks/PWGHF/taskX.cxx b/Analysis/Tasks/PWGHF/taskX.cxx index 64de1cb9f7db5..7280ed78e0d58 100644 --- a/Analysis/Tasks/PWGHF/taskX.cxx +++ b/Analysis/Tasks/PWGHF/taskX.cxx @@ -17,11 +17,16 @@ #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" +#include "AnalysisCore/HFSelectorCuts.h" #include "AnalysisDataModel/HFSecondaryVertex.h" #include "AnalysisDataModel/HFCandidateSelectionTables.h" using namespace o2; +using namespace o2::analysis; +using namespace o2::analysis::hf_cuts_x_tojpsipipi; using namespace o2::framework; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::aod::hf_cand_x; using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_prong2; @@ -37,60 +42,196 @@ void customize(std::vector& workflowOptions) struct TaskX { HistogramRegistry registry{ "registry", - {{"hMassJpsi", "2-prong candidates;inv. mass (e+ e-) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, - {"hPtCand", "X candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}}}; + {{"hPtProng0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtProng1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtProng2", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtCand", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{150, 0., 15.}}}}}}; - Configurable selectionFlagJpsi{"selectionFlagJpsi", 1, "Selection Flag for Jpsi"}; + Configurable d_selectionFlagX{"d_selectionFlagX", 1, "Selection Flag for X"}; Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable> bins{"pTBins", std::vector{hf_cuts_x_tojpsipipi::pTBins_v}, "pT bin limits"}; - Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= selectionFlagJpsi); + void init(o2::framework::InitContext&) + { + registry.add("hMass", "3-prong candidates;inv. mass (J/#psi #pi+ #pi-) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 3., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclength", "3-prong candidates;decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hdeclengthxy", "3-prong candidates;decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong2", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPA", "3-prong candidates;cosine of pointing angle;entries", {HistType::kTH2F, {{500, 0.9, 1.0}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEta", "3-prong candidates;candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hImpParErr", "3-prong candidates;impact parameter error (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenErr", "3-prong candidates;decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenXYErr", "3-prong candidates;decay length xy error (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_x::isSelXToJpsiPiPi >= d_selectionFlagX); - /// aod::BigTracks is not soa::Filtered, should be added when filters are added - void process(aod::Collision const&, aod::BigTracks const& tracks, soa::Filtered> const& candidates) + void process(soa::Filtered> const& candidates) { for (auto& candidate : candidates) { - if (!(candidate.hfflag() & 1 << DecayType::JpsiToEE)) { + if (!(candidate.hfflag() & 1 << XToJpsiPiPi)) { continue; } - if (cutYCandMax >= 0. && std::abs(YJpsi(candidate)) > cutYCandMax) { + if (cutYCandMax >= 0. && std::abs(YX(candidate)) > cutYCandMax) { continue; } - registry.fill(HIST("hMassJpsi"), InvMassJpsiToEE(candidate)); - - int index0jpsi = candidate.index0Id(); - int index1jpsi = candidate.index1Id(); - for (auto& track1 : tracks) { - int signTrack1 = track1.sign(); - int indexTrack1 = track1.globalIndex(); - if (signTrack1 > 0) { - if (indexTrack1 == index0jpsi) { - continue; - } - } else if (indexTrack1 == index1jpsi) { + + registry.fill(HIST("hMass"), InvMassXToJpsiPiPi(candidate), candidate.pt()); + registry.fill(HIST("hPtCand"), candidate.pt()); + registry.fill(HIST("hPtProng0"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate.ptProng1()); + registry.fill(HIST("hPtProng2"), candidate.ptProng2()); + registry.fill(HIST("hdeclength"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hdeclengthxy"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hd0Prong0"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hd0Prong2"), candidate.impactParameter2(), candidate.pt()); + registry.fill(HIST("hCPA"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEta"), candidate.eta(), candidate.pt()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter0(), candidate.pt()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1(), candidate.pt()); + registry.fill(HIST("hDecLenErr"), candidate.errorDecayLength(), candidate.pt()); + registry.fill(HIST("hDecLenXYErr"), candidate.errorDecayLengthXY(), candidate.pt()); + } // candidate loop + } // process +}; // struct + +struct TaskXMC { + HistogramRegistry registry{ + "registry", + {{"hPtRecSig", "3-prong candidates (rec. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{150, 0., 15.}}}}, + {"hPtRecBg", "3-prong candidates (rec. unmatched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{150, 0., 15.}}}}, + {"hPtGen", "3-prong candidates (gen. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{150, 0., 15.}}}}, + {"hPtGenSig", "3-prong candidates (rec. matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{150, 0., 15.}}}}}}; + + Configurable d_selectionFlagX{"d_selectionFlagX", 1, "Selection Flag for X"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable> bins{"pTBins", std::vector{hf_cuts_x_tojpsipipi::pTBins_v}, "pT bin limits"}; + + void init(o2::framework::InitContext&) + { + registry.add("hCPARecSig", "3-prong candidates (rec. matched);cosine of pointing angle;entries", {HistType::kTH2F, {{500, 0.9, 1.0}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPARecBg", "3-prong candidates (rec. unmatched);cosine of pointing angle;entries", {HistType::kTH2F, {{500, 0.9, 1.0}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaRecSig", "3-prong candidates (rec. matched);#it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaRecBg", "3-prong candidates (rec. unmatched);#it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaGen", "3-prong candidates (gen. matched);#it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + + registry.add("hPtProng0RecSig", "3-prong candidates (rec. matched);prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng1RecSig", "3-prong candidates (rec. matched);prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng2RecSig", "3-prong candidates (rec. matched);prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng0RecBg", "3-prong candidates (rec. unmatched);prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng1RecBg", "3-prong candidates (rec. unmatched);prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng2RecBg", "3-prong candidates (rec. unmatched);prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtGenProng0", "3-prong candidates (gen. matched);prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtGenProng1", "3-prong candidates (gen. matched);prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtGenProng2", "3-prong candidates (gen. matched);prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + + registry.add("hMassRecSig", "3-prong candidates (rec. matched);inv. mass (J/#psi #pi+ #pi-) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 3., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassRecBg", "3-prong candidates (rec. unmatched);inv. mass (J/#psi #pi+ #pi-) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 3., 4.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0RecSig", "3-prong candidates (rec. matched);prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1RecSig", "3-prong candidates (rec. matched);prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong2RecSig", "3-prong candidates (rec. matched);prong 2 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0RecBg", "3-prong candidates (rec. unmatched);prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1RecBg", "3-prong candidates (rec. unmatched);prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong2RecBg", "3-prong candidates (rec. unmatched);prong 2 DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{400, -0.002, 0.002}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDeclengthRecSig", "3-prong candidates (rec. matched);decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDeclengthRecBg", "3-prong candidates (rec. unmatched);decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.01}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + + registry.add("hChi2PCASig", "3-prong candidates (rec. matched);chi2 PCA (cm);entries", {HistType::kTH2F, {{1000, 0., 0.0001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hChi2PCABg", "3-prong candidates (rec. unmatched);chi2 PCA (cm);entries", {HistType::kTH2F, {{1000, 0., 0.0001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCtSig", "3-prong candidates (rec. matched);proper lifetime X(3872) * #it{c} (cm);entries", {HistType::kTH2F, {{400, 0., 0.001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCtBg", "3-prong candidates (rec. unmatched);proper lifetime X(3872) * #it{c} (cm);entries", {HistType::kTH2F, {{400, 0., 0.001}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYSig", "3-prong candidates (rec. matched);candidate rapidity;entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYBg", "3-prong candidates (rec. unmatched);candidate rapidity;entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_x::isSelXToJpsiPiPi >= d_selectionFlagX); + + void process(soa::Filtered> const& candidates, + soa::Join const& particlesMC, aod::BigTracksMC const& tracks) + { + // MC rec. + //Printf("MC Candidates: %d", candidates.size()); + for (auto& candidate : candidates) { + if (!(candidate.hfflag() & 1 << XToJpsiPiPi)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YX(candidate)) > cutYCandMax) { + continue; + } + if (candidate.flagMCMatchRec() == 1 << XToJpsiPiPi) { + auto indexMother = RecoDecay::getMother(particlesMC, candidate.index1_as().mcParticle_as>(), 9920443, true); + auto particleMother = particlesMC.iteratorAt(indexMother); + registry.fill(HIST("hPtGenSig"), particleMother.pt()); + registry.fill(HIST("hPtRecSig"), candidate.pt()); + registry.fill(HIST("hCPARecSig"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEtaRecSig"), candidate.eta(), candidate.pt()); + + registry.fill(HIST("hDeclengthRecSig"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hMassRecSig"), InvMassXToJpsiPiPi(candidate), candidate.pt()); + registry.fill(HIST("hd0Prong0RecSig"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1RecSig"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hd0Prong2RecSig"), candidate.impactParameter2(), candidate.pt()); + registry.fill(HIST("hPtProng0RecSig"), candidate.ptProng0(), candidate.pt()); + registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1(), candidate.pt()); + registry.fill(HIST("hPtProng2RecSig"), candidate.ptProng2(), candidate.pt()); + registry.fill(HIST("hChi2PCASig"), candidate.chi2PCA(), candidate.pt()); + registry.fill(HIST("hCtSig"), CtX(candidate), candidate.pt()); + registry.fill(HIST("hYSig"), YX(candidate), candidate.pt()); + } else { + registry.fill(HIST("hPtRecBg"), candidate.pt()); + registry.fill(HIST("hCPARecBg"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEtaRecBg"), candidate.eta(), candidate.pt()); + + registry.fill(HIST("hDeclengthRecBg"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hMassRecBg"), InvMassXToJpsiPiPi(candidate), candidate.pt()); + registry.fill(HIST("hd0Prong0RecBg"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1RecBg"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hd0Prong2RecBg"), candidate.impactParameter2(), candidate.pt()); + registry.fill(HIST("hPtProng0RecBg"), candidate.ptProng0(), candidate.pt()); + registry.fill(HIST("hPtProng1RecBg"), candidate.ptProng1(), candidate.pt()); + registry.fill(HIST("hPtProng2RecBg"), candidate.ptProng2(), candidate.pt()); + registry.fill(HIST("hChi2PCABg"), candidate.chi2PCA(), candidate.pt()); + registry.fill(HIST("hCtBg"), CtX(candidate), candidate.pt()); + registry.fill(HIST("hYBg"), YX(candidate), candidate.pt()); + } + } // rec + // MC gen. + //Printf("MC Particles: %d", particlesMC.size()); + for (auto& particle : particlesMC) { + if (particle.flagMCMatchGen() == 1 << XToJpsiPiPi) { + // TODO: add X(3872) mass such that we can use the getMassPDG function instead of hardcoded mass + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle.px(), particle.py(), particle.pz()}, 3.87168)) > cutYCandMax) { + // Printf("MC Gen.: Y rejection: %g", RecoDecay::Y(array{particle.px(), particle.py(), particle.pz()}, 3.87168)); continue; } - for (auto track2 = track1 + 1; track2 != tracks.end(); ++track2) { - if (signTrack1 == track2.sign()) { - continue; - } - int indexTrack2 = track2.globalIndex(); - if (signTrack1 > 0) { - if (indexTrack2 == index1jpsi) { - continue; - } - } else if (indexTrack2 == index0jpsi) { - continue; - } - registry.fill(HIST("hPtCand"), candidate.pt() + track1.pt() + track2.pt()); - } // track2 loop (pion) - } // track1 loop (pion) - } // Jpsi loop - } // process -}; // struct + registry.fill(HIST("hPtGen"), particle.pt()); + registry.fill(HIST("hEtaGen"), particle.eta(), particle.pt()); + + // properties of gen matched X(3872), to get a first look at some cuts + float ptProngs[3]; + int counter = 0; + for (int iD = particle.daughter0Id(); iD <= particle.daughter1Id(); ++iD) { + ptProngs[counter] = particlesMC.iteratorAt(iD).pt(); + counter++; + } + registry.fill(HIST("hPtGenProng0"), ptProngs[0], particle.pt()); + registry.fill(HIST("hPtGenProng1"), ptProngs[1], particle.pt()); + registry.fill(HIST("hPtGenProng2"), ptProngs[2], particle.pt()); + } + } // gen + } // process +}; // struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{ adaptAnalysisTask(cfgc, TaskName{"hf-task-x"})}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-task-x-mc"})); + } return workflow; }