From cf525b8965913d17bde3df229515f7ba896793ca Mon Sep 17 00:00:00 2001 From: nzardosh Date: Mon, 26 Oct 2020 17:37:45 +0100 Subject: [PATCH] HF PreSelection Changes further chnages for PreSel adding final touches adding histograms correct vertex reconstruction adding Jpsitoee adding debug option mend --- Analysis/Core/CMakeLists.txt | 2 + .../Core/include/Analysis/HFConfigurables.h | 61 ++ Analysis/Core/src/AnalysisCoreLinkDef.h | 1 + Analysis/Core/src/HFConfigurables.cxx | 14 + .../include/Analysis/HFSecondaryVertex.h | 22 +- .../Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx | 765 +++++++++++++----- 6 files changed, 642 insertions(+), 223 deletions(-) create mode 100644 Analysis/Core/include/Analysis/HFConfigurables.h create mode 100644 Analysis/Core/src/HFConfigurables.cxx diff --git a/Analysis/Core/CMakeLists.txt b/Analysis/Core/CMakeLists.txt index 63d608698c184..f931540f2cb3e 100644 --- a/Analysis/Core/CMakeLists.txt +++ b/Analysis/Core/CMakeLists.txt @@ -16,6 +16,7 @@ o2_add_library(AnalysisCore src/AnalysisCut.cxx src/AnalysisCompositeCut.cxx src/TriggerAliases.cxx + src/HFConfigurables.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisTools) o2_target_root_dictionary(AnalysisCore @@ -28,6 +29,7 @@ o2_target_root_dictionary(AnalysisCore include/Analysis/AnalysisCompositeCut.h include/Analysis/TriggerAliases.h include/Analysis/MC.h + include/Analysis/HFConfigurables.h LINKDEF src/AnalysisCoreLinkDef.h) if(FastJet_FOUND) diff --git a/Analysis/Core/include/Analysis/HFConfigurables.h b/Analysis/Core/include/Analysis/HFConfigurables.h new file mode 100644 index 0000000000000..ead5dde738c3d --- /dev/null +++ b/Analysis/Core/include/Analysis/HFConfigurables.h @@ -0,0 +1,61 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// HF Configurable Classes +// +// Authors: Nima Zardoshti + +#ifndef O2_ANALYSIS_HFCONFIGURABLES_H +#define O2_ANALYSIS_HFCONFIGURABLES_H + +#include + +class HFTrackIndexSkimsCreatorConfigs +{ + public: + HFTrackIndexSkimsCreatorConfigs() = default; + ~HFTrackIndexSkimsCreatorConfigs() = default; + + // 2-prong cuts D0 + double mPtD0Min = 0.; + double mInvMassD0Min = 1.46; + double mInvMassD0Max = 2.26; + double mCPAD0Min = 0.75; + double mImpParProductD0Max = -0.00005; + // 2-prong cuts Jpsi + double mPtJpsiMin = 0.; + double mInvMassJpsiMin = 2.5; + double mInvMassJpsiMax = 4.1; + double mCPAJpsiMin = -2; + double mImpParProductJpsiMax = 1000.; + // 3-prong cuts - D+ + double mPtDPlusMin = 0.; + double mInvMassDPlusMin = 1.7; + double mInvMassDPlusMax = 2.05; + double mCPADPlusMin = 0.5; + double mDecLenDPlusMin = 0.; + // 3-prong cuts - Lc + double mPtLcMin = 0.; + double mInvMassLcMin = 2.1; + double mInvMassLcMax = 2.5; + double mCPALcMin = 0.5; + double mDecLenLcMin = 0.; + // 3-prong cuts - Ds + double mPtDsMin = 0.; + double mInvMassDsMin = 1.7; + double mInvMassDsMax = 2.2; + double mCPADsMin = 0.5; + double mDecLenDsMin = 0.; + + private: + ClassDef(HFTrackIndexSkimsCreatorConfigs, 1); +}; + +#endif diff --git a/Analysis/Core/src/AnalysisCoreLinkDef.h b/Analysis/Core/src/AnalysisCoreLinkDef.h index ed2833df8a9f0..c550c9a41d98a 100644 --- a/Analysis/Core/src/AnalysisCoreLinkDef.h +++ b/Analysis/Core/src/AnalysisCoreLinkDef.h @@ -20,3 +20,4 @@ #pragma link C++ class HistogramManager + ; #pragma link C++ class AnalysisCut + ; #pragma link C++ class AnalysisCompositeCut + ; +#pragma link C++ class HFTrackIndexSkimsCreatorConfigs + ; diff --git a/Analysis/Core/src/HFConfigurables.cxx b/Analysis/Core/src/HFConfigurables.cxx new file mode 100644 index 0000000000000..fb283c8145718 --- /dev/null +++ b/Analysis/Core/src/HFConfigurables.cxx @@ -0,0 +1,14 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// HF Configurable Classes +// +// Authors: Nima Zardoshti +#include "Analysis/HFConfigurables.h" diff --git a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h index c3520064b443f..26d4b5726fddb 100644 --- a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h @@ -25,15 +25,13 @@ namespace o2::aod { namespace hf_seltrack { -DECLARE_SOA_COLUMN(IsSel2Prong, isSel2Prong, int); -DECLARE_SOA_COLUMN(IsSel3Prong, isSel3Prong, int); +DECLARE_SOA_COLUMN(IsSelProng, isSelProng, int); DECLARE_SOA_COLUMN(DCAPrim0, dcaPrim0, float); DECLARE_SOA_COLUMN(DCAPrim1, dcaPrim1, float); } // namespace hf_seltrack DECLARE_SOA_TABLE(HFSelTrack, "AOD", "HFSELTRACK", - hf_seltrack::IsSel2Prong, - hf_seltrack::IsSel3Prong, + hf_seltrack::IsSelProng, hf_seltrack::DCAPrim0, hf_seltrack::DCAPrim1); @@ -50,6 +48,13 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Index1, index1, int, BigTracks, "fIndex1"); DECLARE_SOA_INDEX_COLUMN_FULL(Index2, index2, int, BigTracks, "fIndex2"); DECLARE_SOA_INDEX_COLUMN_FULL(Index3, index3, int, BigTracks, "fIndex3"); DECLARE_SOA_COLUMN(HFflag, hfflag, uint8_t); + +DECLARE_SOA_COLUMN(D0ToKPiFlag, d0ToKPiFlag, uint8_t); +DECLARE_SOA_COLUMN(JpsiToEEFlag, jpsiToEEFlag, uint8_t); + +DECLARE_SOA_COLUMN(DPlusPiKPiFlag, dPlusPiKPiFlag, uint8_t); +DECLARE_SOA_COLUMN(LcPKPiFlag, lcPKPiFlag, uint8_t); +DECLARE_SOA_COLUMN(DsKKPiFlag, dsKKPiFlag, uint8_t); } // namespace hf_track_index DECLARE_SOA_TABLE(HfTrackIndexProng2, "AOD", "HFTRACKIDXP2", @@ -57,12 +62,21 @@ DECLARE_SOA_TABLE(HfTrackIndexProng2, "AOD", "HFTRACKIDXP2", hf_track_index::Index1Id, hf_track_index::HFflag); +DECLARE_SOA_TABLE(HfCutStatusProng2, "AOD", "HFCUTSTATUSP2", + hf_track_index::D0ToKPiFlag, + hf_track_index::JpsiToEEFlag); + DECLARE_SOA_TABLE(HfTrackIndexProng3, "AOD", "HFTRACKIDXP3", hf_track_index::Index0Id, hf_track_index::Index1Id, hf_track_index::Index2Id, hf_track_index::HFflag); +DECLARE_SOA_TABLE(HfCutStatusProng3, "AOD", "HFCUTSTATUSP3", + hf_track_index::DPlusPiKPiFlag, + hf_track_index::LcPKPiFlag, + hf_track_index::DsKKPiFlag); + // general decay properties namespace hf_cand { diff --git a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx index f203f9c6e9d17..8fc8e9261a6a6 100644 --- a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx +++ b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx @@ -13,12 +13,17 @@ /// /// \author Gian Michele Innocenti , CERN /// \author Vít Kučera , CERN +/// \author Nima Zardoshti , CERN #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "DetectorsVertexing/DCAFitterN.h" #include "Analysis/HFSecondaryVertex.h" #include "Analysis/trackUtilities.h" +#include "Analysis/HFConfigurables.h" +//#include "Analysis/Centrality.h" +#include "Framework/HistogramRegistry.h" +#include using namespace o2; using namespace o2::framework; @@ -32,99 +37,105 @@ struct SelectTracks { Configurable d_bz{"d_bz", 5., "bz field"}; // quality cut Configurable doCutQuality{"doCutQuality", true, "apply quality cuts"}; - Configurable d_tpcnclsfound{"d_tpcnclsfound", 70, "min. number of TPC clusters"}; + Configurable d_tpcnclsfound{"d_tpcnclsfound", 70, ">= min. number of TPC clusters needed"}; // 2-prong cuts - Configurable ptmintrack_2prong{"ptmintrack_2prong", -1., "min. track pT"}; - Configurable dcatoprimxymin_2prong{"dcatoprimxymin_2prong", 0., "min. DCAXY to prim. vtx."}; - Configurable etamax_2prong{"etamax_2prong", 999., "max. pseudorapidity"}; + Configurable ptmintrack_2prong{"ptmintrack_2prong", -1., "min. track pT for 2 prong candidate"}; + Configurable dcatoprimxy_2prong_maxpt{"dcatoprimxymin_2prong_minpt", 2., "max pt cut for min. DCAXY to prim. vtx. for 2 prong candidate"}; + Configurable dcatoprimxymin_2prong{"dcatoprimxymin_2prong", 0., "min. DCAXY to prim. vtx. for 2 prong candidate"}; + Configurable dcatoprimxymax_2prong{"dcatoprimxymax_2prong", 1.0, "max. DCAXY to prim. vtx. for 2 prong candidate"}; + Configurable etamax_2prong{"etamax_2prong", 4., "max. pseudorapidity for 2 prong candidate"}; // 3-prong cuts - Configurable ptmintrack_3prong{"ptmintrack_3prong", -1., "min. track pT"}; - Configurable dcatoprimxymin_3prong{"dcatoprimxymin_3prong", 0., "min. DCAXY to prim. vtx."}; - Configurable etamax_3prong{"etamax_3prong", 999., "max. pseudorapidity"}; - - OutputObj hpt_nocuts{TH1F("hpt_nocuts", "all tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", 100, 0., 10.)}; - // 2-prong histograms - OutputObj hpt_cuts_2prong{TH1F("hpt_cuts_2prong", "tracks selected for 2-prong vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", 100, 0., 10.)}; - OutputObj hdcatoprimxy_cuts_2prong{TH1F("hdcatoprimxy_cuts_2prong", "tracks selected for 2-prong vertexing;DCAxy to prim. vtx. (cm);entries", 100, -1., 1.)}; - OutputObj heta_cuts_2prong{TH1F("heta_cuts_2prong", "tracks selected for 2-prong vertexing;#it{#eta};entries", 100, -1., 1.)}; - // 3-prong histograms - OutputObj hpt_cuts_3prong{TH1F("hpt_cuts_3prong", "tracks selected for 3-prong vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", 100, 0., 10.)}; - OutputObj hdcatoprimxy_cuts_3prong{TH1F("hdcatoprimxy_cuts_3prong", "tracks selected for 3-prong vertexing;DCAxy to prim. vtx. (cm);entries", 100, -1., 1.)}; - OutputObj heta_cuts_3prong{TH1F("heta_cuts_3prong", "tracks selected for 3-prong vertexing;#it{#eta};entries", 100, -1., 1.)}; + Configurable ptmintrack_3prong{"ptmintrack_3prong", -1., "min. track pT for 3 prong candidate"}; + Configurable dcatoprimxy_3prong_maxpt{"dcatoprimxy_3prong_minpt", 2., "max pt cut for min. DCAXY to prim. vtx. for 3 prong candidate"}; + Configurable dcatoprimxymin_3prong{"dcatoprimxymin_3prong", 0., "min. DCAXY to prim. vtx. for 3 prong candidate"}; + Configurable dcatoprimxymax_3prong{"dcatoprimxymax_3prong", 1.0, "max. DCAXY to prim. vtx. for 3 prong candidate"}; + Configurable etamax_3prong{"etamax_3prong", 4., "max. pseudorapidity for 3 prong candidate"}; + + HistogramRegistry registry{ + "registry", + {{"hpt_nocuts", "all tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 100.}}}}, + // 2-prong histograms + {"hpt_cuts_2prong", "tracks selected for 2-prong vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 100.}}}}, + {"hdcatoprimxy_cuts_2prong", "tracks selected for 2-prong vertexing;DCAxy to prim. vtx. (cm);entries", {HistType::kTH1F, {{static_cast(1.2 * dcatoprimxymax_2prong * 100), -1.2 * dcatoprimxymax_2prong, 1.2 * dcatoprimxymax_2prong}}}}, + {"heta_cuts_2prong", "tracks selected for 2-prong vertexing;#it{#eta};entries", {HistType::kTH1F, {{static_cast(1.2 * etamax_2prong * 100), -1.2 * etamax_2prong, 1.2 * etamax_2prong}}}}, + // 3-prong histograms + {"hpt_cuts_3prong", "tracks selected for 3-prong vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 100.}}}}, + {"hdcatoprimxy_cuts_3prong", "tracks selected for 3-prong vertexing;DCAxy to prim. vtx. (cm);entries", {HistType::kTH1F, {{static_cast(1.2 * dcatoprimxymax_3prong) * 100, -1.2 * dcatoprimxymax_3prong, 1.2 * dcatoprimxymax_3prong}}}}, + {"heta_cuts_3prong", "tracks selected for 3-prong vertexing;#it{#eta};entries", {HistType::kTH1F, {{static_cast(1.2 * etamax_3prong * 100), -1.2 * etamax_3prong, 1.2 * etamax_3prong}}}}}}; void process(aod::Collision const& collision, soa::Join const& tracks) { math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); for (auto& track : tracks) { - int status_2prong = 1; // selection flag - int status_3prong = 1; // selection flag + int status_prong = 3; // selection flag , 2 bits on + + auto trackPt = track.pt(); if (b_dovalplots) { - hpt_nocuts->Fill(track.pt()); + registry.get("hpt_nocuts")->Fill(trackPt); } // pT cut - if (track.pt() < ptmintrack_2prong) { - status_2prong = 0; + if (trackPt < ptmintrack_2prong) { + status_prong = status_prong & ~(1 << 0); // the bitwise operation & ~(1 << n) will set the nth bit to 0 } - if (track.pt() < ptmintrack_3prong) { - status_3prong = 0; + if (trackPt < ptmintrack_3prong) { + status_prong = status_prong & ~(1 << 1); } + auto trackEta = track.eta(); // eta cut - if (status_2prong && abs(track.eta()) > etamax_2prong) { - status_2prong = 0; + if ((status_prong & (1 << 0)) && abs(trackEta) > etamax_2prong) { + status_prong = status_prong & ~(1 << 0); } - if (status_3prong && abs(track.eta()) > etamax_3prong) { - status_3prong = 0; + if ((status_prong & (1 << 1)) && abs(trackEta) > etamax_3prong) { + status_prong = status_prong & ~(1 << 1); } // quality cut - if (doCutQuality && (status_2prong || status_3prong)) { + if (doCutQuality && status_prong > 0) { // FIXME to make a more complete selection e.g track.flags() & o2::aod::track::TPCrefit && track.flags() & o2::aod::track::GoldenChi2 && UChar_t clustermap = track.itsClusterMap(); - bool isselected = track.tpcNClsFound() >= d_tpcnclsfound && - track.flags() & o2::aod::track::ITSrefit && - (TESTBIT(clustermap, 0) || TESTBIT(clustermap, 1)); - if (!isselected) { - status_2prong = 0; - status_3prong = 0; + if (!(track.tpcNClsFound() >= d_tpcnclsfound && + track.flags() & o2::aod::track::ITSrefit && + (TESTBIT(clustermap, 0) || TESTBIT(clustermap, 1)))) { + status_prong = 0; } } // DCA cut array dca; - if (status_2prong || status_3prong) { + if (status_prong > 0) { + double dcatoprimxymin_2prong_ptdep = dcatoprimxymin_2prong * TMath::Max(0., (1 - TMath::Floor(trackPt / dcatoprimxy_2prong_maxpt))); + double dcatoprimxymin_3prong_ptdep = dcatoprimxymin_3prong * TMath::Max(0., (1 - TMath::Floor(trackPt / dcatoprimxy_3prong_maxpt))); auto trackparvar0 = getTrackParCov(track); - bool isprop = trackparvar0.propagateParamToDCA(vtxXYZ, d_bz, &dca, 100.); // get impact parameters - if (!isprop) { - status_2prong = 0; - status_3prong = 0; + if (!trackparvar0.propagateParamToDCA(vtxXYZ, d_bz, &dca, 100.)) { // get impact parameters + status_prong = 0; } - if (status_2prong && abs(dca[0]) < dcatoprimxymin_2prong) { - status_2prong = 0; + if ((status_prong & (1 << 0)) && (abs(dca[0]) < dcatoprimxymin_2prong_ptdep || abs(dca[0]) > dcatoprimxymax_2prong)) { + status_prong = status_prong & ~(1 << 0); } - if (status_3prong && abs(dca[0]) < dcatoprimxymin_3prong) { - status_3prong = 0; + if ((status_prong & (1 << 1)) && (abs(dca[0]) < dcatoprimxymin_3prong_ptdep || abs(dca[0]) > dcatoprimxymax_3prong)) { + status_prong = status_prong & ~(1 << 1); } } // fill histograms if (b_dovalplots) { - if (status_2prong == 1) { - hpt_cuts_2prong->Fill(track.pt()); - hdcatoprimxy_cuts_2prong->Fill(dca[0]); - heta_cuts_2prong->Fill(track.eta()); + if (status_prong & (1 << 0)) { + registry.get("hpt_cuts_2prong")->Fill(trackPt); + registry.get("hdcatoprimxy_cuts_2prong")->Fill(dca[0]); + registry.get("heta_cuts_2prong")->Fill(trackEta); } - if (status_3prong == 1) { - hpt_cuts_3prong->Fill(track.pt()); - hdcatoprimxy_cuts_3prong->Fill(dca[0]); - heta_cuts_3prong->Fill(track.eta()); + if (status_prong & (1 << 1)) { + registry.get("hpt_cuts_3prong")->Fill(trackPt); + registry.get("hdcatoprimxy_cuts_3prong")->Fill(dca[0]); + registry.get("heta_cuts_3prong")->Fill(trackEta); } } // fill table row - rowSelectedTrack(status_2prong, status_3prong, dca[0], dca[1]); + rowSelectedTrack(status_prong, dca[0], dca[1]); } } }; @@ -132,68 +143,79 @@ struct SelectTracks { /// Pre-selection of 2-prong and 3-prong secondary vertices struct HFTrackIndexSkimsCreator { Produces rowTrackIndexProng2; + Produces rowProng2CutStatus; Produces rowTrackIndexProng3; + Produces rowProng3CutStatus; + //Configurable nCollsMax{"nCollsMax", -1, "Max collisions per file"}; //can be added to run over limited collisions per file - for tesing purposes Configurable b_dovalplots{"b_dovalplots", true, "fill histograms"}; Configurable do3prong{"do3prong", 0, "do 3 prong"}; // event selection Configurable triggerindex{"triggerindex", -1, "trigger index"}; // vertexing parameters - Configurable d_bz{"d_bz", 5., "magnetic field"}; + Configurable d_bz{"d_bz", 5., "magnetic field kG"}; 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 if chi2/chi2old > this"}; - // 2-prong cuts - Configurable cut2ProngPtCandMin{"cut2ProngPtCandMin", -1., "min. pT of the 2-prong candidate"}; - Configurable cut2ProngInvMassD0Min{"cut2ProngInvMassD0Min", -1., "min. D0 candidate invariant mass"}; - Configurable cut2ProngInvMassD0Max{"cut2ProngInvMassD0Max", -1., "max. D0 candidate invariant mass"}; - Configurable cut2ProngCPAMin{"cut2ProngCPAMin", -2., "min. cosine of pointing angle"}; - Configurable cut2ProngImpParProductMax{"cut2ProngImpParProductMax", 100., "max. product of imp. par. of D0 candidate prongs"}; - // 3-prong cuts - Configurable cut3ProngPtCandMin{"cut3ProngPtCandMin", -1., "min. pT of the 3-prong candidate"}; - Configurable cut3ProngInvMassDPlusMin{"cut3ProngInvMassDPlusMin", 1.5, "min. D+ candidate invariant mass"}; - Configurable cut3ProngInvMassDPlusMax{"cut3ProngInvMassDPlusMax", 2.1, "min. D+ candidate invariant mass"}; - Configurable cut3ProngCPAMin{"cut3ProngCPAMin", -2., "min. cosine of pointing angle"}; - Configurable cut3ProngDecLenMin{"cut3ProngDecLenMin", -1., "min. decay length"}; - - // 2-prong histograms - OutputObj hvtx2_x{TH1F("hvtx2_x", "2-prong candidates;#it{x}_{sec. vtx.} (cm);entries", 1000, -2., 2.)}; - OutputObj hvtx2_y{TH1F("hvtx2_y", "2-prong candidates;#it{y}_{sec. vtx.} (cm);entries", 1000, -2., 2.)}; - OutputObj hvtx2_z{TH1F("hvtx2_z", "2-prong candidates;#it{z}_{sec. vtx.} (cm);entries", 1000, -20., 20.)}; - OutputObj hmass2{TH1F("hmass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; - // 3-prong histograms - OutputObj hvtx3_x{TH1F("hvtx3_x", "3-prong candidates;#it{x}_{sec. vtx.} (cm);entries", 1000, -2., 2.)}; - OutputObj hvtx3_y{TH1F("hvtx3_y", "3-prong candidates;#it{y}_{sec. vtx.} (cm);entries", 1000, -2., 2.)}; - OutputObj hvtx3_z{TH1F("hvtx3_z", "3-prong candidates;#it{z}_{sec. vtx.} (cm);entries", 1000, -20., 20.)}; - OutputObj hmass3{TH1F("hmass3", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", 500, 1.6, 2.1)}; - - /* - // Counter histograms - OutputObj hNCand2ProngVsNTracks{TH2F("hNCand2ProngVsNTracks", "2-prong candidates;# of tracks;# of candidates;entries", 1000, 0., 1000., 1000, 0., 10000.)}; - OutputObj hNCand3ProngVsNTracks{TH2F("hNCand3ProngVsNTracks", "3-prong candidates;# of tracks;# of candidates;entries", 1000, 0., 1000., 1000, 0., 10000.)}; - OutputObj hNCand2Prong{TH1F("hNCand2Prong", "2-prong candidates;# of candidates;entries", 1000, 0., 10000.)}; - OutputObj hNCand3Prong{TH1F("hNCand3Prong", "3-prong candidates;# of candidates;entries", 1000, 0., 10000.)}; - OutputObj hNTracks{TH1F("hNTracks", ";# of tracks;entries", 1000, 0., 1000.)}; - */ - - Filter filterSelectTracks = (aod::hf_seltrack::isSel2Prong == 1); + Configurable configs{"configs", {}, "configurables"}; + Configurable b_debug{"b_debug", false, "debug mode"}; + + HistogramRegistry registry{ + "registry", + {{"hNTracks", ";# of tracks;entries", {HistType::kTH1F, {{2500, 0., 25000.}}}}, + //2prong histograms + {"hvtx2_x", "2-prong candidates;#it{x}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx2_y", "2-prong candidates;#it{y}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx2_z", "2-prong candidates;#it{z}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -20., 20.}}}}, + {"hNCand2Prong", "2-prong candidates preselected;# of candidates;entries", {HistType::kTH1F, {{2000, 0., 200000.}}}}, + {"hNCand2ProngVsNTracks", "2-prong candidates preselected;# of selected tracks;# of candidates;entries", {HistType::kTH2F, {{2500, 0., 25000.}, {2000, 0., 200000.}}}}, + {"hmassD0", "D0 candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + {"hmassJpsi", "Jpsi candidates;inv. mass (e+ e-) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + //3prong histograms + {"hvtx3_x", "3-prong candidates;#it{x}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx3_y", "3-prong candidates;#it{y}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx3_z", "3-prong candidates;#it{z}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -20., 20.}}}}, + {"hNCand3Prong", "3-prong candidates preselected;# of candidates;entries", {HistType::kTH1F, {{5000, 0., 500000.}}}}, + {"hNCand3ProngVsNTracks", "3-prong candidates preselected;# of selected tracks;# of candidates;entries", {HistType::kTH2F, {{2500, 0., 25000.}, {5000, 0., 500000.}}}}, + {"hmassDPlus", "D+ candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + {"hmassLc", "Lc candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + {"hmassDs", "Ds candidates;inv. mass (K K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}}}; + + Filter filterSelectTracks = (aod::hf_seltrack::isSelProng > 0); using SelectedTracks = soa::Filtered>; + // FIXME //Partition tracksPos = aod::track::signed1Pt > 0.f; //Partition tracksNeg = aod::track::signed1Pt < 0.f; double massPi = RecoDecay::getMassPDG(kPiPlus); double massK = RecoDecay::getMassPDG(kKPlus); - double mass2PiK{0.}; - double mass2KPi{0.}; - double mass3PiKPi{0.}; + double massProton = RecoDecay::getMassPDG(kProton); + double massElectron = RecoDecay::getMassPDG(kElectron); - void process(aod::Collision const& collision, - aod::BCs const& bcs, - SelectedTracks const& tracks) + // int nColls{0}; //can be added to run over limited collisions per file - for tesing purposes + + void process( //soa::Join::iterator const& collision, //FIXME add centrality when option for variations to the process function appears + aod::Collision const& collision, + aod::BCs const& bcs, + SelectedTracks const& tracks) { + + //can be added to run over limited collisions per file - for tesing purposes + /* + if (nCollsMax > -1){ + if (nColls == nCollMax){ + return; + //can be added to run over limited collisions per file - for tesing purposes + } + nColls++; + } + */ + + //auto centrality = collision.centV0M(); //FIXME add centrality when option for variations to the process function appears + int trigindex = int{triggerindex}; if (trigindex != -1) { uint64_t triggerMask = collision.bc().triggerMask(); @@ -203,6 +225,60 @@ struct HFTrackIndexSkimsCreator { } } + //FIXME move above process function + const int n2ProngDecays = 2; //number of two prong hadron types + const int n3ProngDecays = 3; //number of three prong hadron types + int n2ProngBit = TMath::Power(2, n2ProngDecays) - 1; //bit value for 2 prongs candidates where each candidiate is one bit and they are all set it 1 + int n3ProngBit = TMath::Power(2, n3ProngDecays) - 1; //bit value for 3 prongs candidates where each candidiate is one bit and they are all set it 1 + + //retrieve cuts from json - to be made pT dependent when option appears in json + const int nCuts2Prong = 4; //how mant different selections are made on 2 prongs + double cut2ProngPtCandMin[n2ProngDecays] = {configs->mPtD0Min, configs->mPtJpsiMin}; + double cut2ProngInvMassCandMin[n2ProngDecays] = {configs->mInvMassD0Min, configs->mInvMassJpsiMin}; + double cut2ProngInvMassCandMax[n2ProngDecays] = {configs->mInvMassD0Max, configs->mInvMassJpsiMax}; + double cut2ProngCPACandMin[n2ProngDecays] = {configs->mCPAD0Min, configs->mCPAJpsiMin}; + double cut2ProngImpParProductCandMax[n2ProngDecays] = {configs->mImpParProductD0Max, configs->mImpParProductJpsiMax}; + + const int nCuts3Prong = 4; //how mant different selections are made on 3 prongs + double cut3ProngPtCandMin[n3ProngDecays] = {configs->mPtDPlusMin, configs->mPtLcMin, configs->mPtDsMin}; + double cut3ProngInvMassCandMin[n3ProngDecays] = {configs->mInvMassDPlusMin, configs->mInvMassLcMin, configs->mInvMassDsMin}; + double cut3ProngInvMassCandMax[n3ProngDecays] = {configs->mInvMassDPlusMax, configs->mInvMassLcMax, configs->mInvMassDsMax}; + double cut3ProngCPACandMin[n3ProngDecays] = {configs->mCPADPlusMin, configs->mCPALcMin, configs->mCPADsMin}; + double cut3ProngDecLenCandMin[n3ProngDecays] = {configs->mDecLenDPlusMin, configs->mDecLenLcMin, configs->mDecLenDsMin}; + + bool cutStatus2Prong[n2ProngDecays][nCuts2Prong]; + bool cutStatus3Prong[n3ProngDecays][nCuts3Prong]; + int nCutStatus2ProngBit = TMath::Power(2, nCuts2Prong) - 1; //bit value for selection status for each 2 prongs candidate where each selection is one bit and they are all set it 1 + int nCutStatus3ProngBit = TMath::Power(2, nCuts3Prong) - 1; //bit value for selection status for each 2 prongs candidate where each selection is one bit and they are all set it 1 + + //set mass hypothesis 1 for 2 prong + auto arr2Mass1 = array{ + array{massPi, massK}, + array{massElectron, massElectron}}; + + //set mass hypothesis 2 for 2 prong + auto arr2Mass2 = array{ + array{massK, massPi}, + array{massElectron, massElectron}}; + + //set mass hypothesis 1 for 3 prong + auto arr3Mass1 = array{ + array{massPi, massK, massPi}, + array{massProton, massK, massPi}, + array{massK, massK, massPi}}; + + //set mass hypothesis 2 for 3 prong + auto arr3Mass2 = array{ + array{massPi, massK, massPi}, + array{massPi, massK, massProton}, + array{massPi, massK, massK}}; + + double mass2ProngHypo1[n2ProngDecays]; + double mass2ProngHypo2[n2ProngDecays]; + + double mass3ProngHypo1[n3ProngDecays]; + double mass3ProngHypo2[n3ProngDecays]; + // 2-prong vertex fitter o2::vertexing::DCAFitterN<2> df2; df2.setBz(d_bz); @@ -223,8 +299,9 @@ struct HFTrackIndexSkimsCreator { df3.setMinRelChi2Change(d_minrelchi2change); df3.setUseAbsDCA(true); - //auto nCand2 = rowTrackIndexProng2.lastIndex(); - //auto nCand3 = rowTrackIndexProng3.lastIndex(); + //used to calculate number of candidiates per event + auto nCand2 = rowTrackIndexProng2.lastIndex(); + auto nCand3 = rowTrackIndexProng3.lastIndex(); // first loop over positive tracks //for (auto trackPos1 = tracksPos.begin(); trackPos1 != tracksPos.end(); ++trackPos1) { @@ -232,6 +309,14 @@ struct HFTrackIndexSkimsCreator { if (trackPos1.signed1Pt() < 0) { continue; } + bool sel2ProngStatus = true; + bool sel3ProngStatusPos1 = trackPos1.isSelProng() & (1 << 1); + if (!(trackPos1.isSelProng() & (1 << 0))) { + sel2ProngStatus = false; + } + if (!sel2ProngStatus && !sel3ProngStatusPos1) { + continue; + } auto trackParVarPos1 = getTrackParCov(trackPos1); @@ -241,85 +326,152 @@ struct HFTrackIndexSkimsCreator { if (trackNeg1.signed1Pt() > 0) { continue; } + bool sel3ProngStatusNeg1 = trackNeg1.isSelProng() & (1 << 1); + if (sel2ProngStatus && !(trackNeg1.isSelProng() & (1 << 0))) { + sel2ProngStatus = false; + } + if (!sel2ProngStatus && !sel3ProngStatusNeg1) { + continue; + } auto trackParVarNeg1 = getTrackParCov(trackNeg1); - auto pVecCandProng2 = array{trackPos1.px() + trackNeg1.px(), - trackPos1.py() + trackNeg1.py(), - trackPos1.pz() + trackNeg1.pz()}; - bool isSelectedCand2Prong = true; - - // candidate pT cut - if (RecoDecay::Pt(pVecCandProng2) < cut2ProngPtCandMin) { - isSelectedCand2Prong = false; - } - if (isSelectedCand2Prong) { - // reconstruct the 2-prong secondary vertex - if (df2.process(trackParVarPos1, trackParVarNeg1) == 0) { - continue; - } + int isSelected2ProngCand = n2ProngBit; //bitmap for checking status of two-prong candidates (1 is true, 0 is rejected) - // imp. par. product cut - if (isSelectedCand2Prong && cut2ProngImpParProductMax < 100.) { - if (trackPos1.dcaPrim0() * trackNeg1.dcaPrim0() > cut2ProngImpParProductMax) { - isSelectedCand2Prong = false; + if (b_debug) { + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + for (int n2cut = 0; n2cut < nCuts2Prong; n2cut++) { + cutStatus2Prong[n2][n2cut] = true; } } - - // get secondary vertex - const auto& secondaryVertex2 = df2.getPCACandidate(); - // get track momenta - array pvec0; - array pvec1; - df2.getTrack(0).getPxPyPzGlo(pvec0); - df2.getTrack(1).getPxPyPzGlo(pvec1); - - // CPA cut - if (isSelectedCand2Prong && cut2ProngCPAMin > -2.) { - pVecCandProng2 = RecoDecay::PVec(pvec0, pvec1); - auto cpa = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex2, pVecCandProng2); - if (cpa < cut2ProngCPAMin) { - isSelectedCand2Prong = false; + } + int iDebugCut = 0; + + // 2prong invariant-mass cut + if (sel2ProngStatus > 0) { + auto arrMom = array{ + array{trackPos1.px(), trackPos1.py(), trackPos1.pz()}, + array{trackNeg1.px(), trackNeg1.py(), trackNeg1.pz()}}; + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + mass2ProngHypo1[n2] = RecoDecay::M(arrMom, arr2Mass1[n2]); + mass2ProngHypo2[n2] = RecoDecay::M(arrMom, arr2Mass2[n2]); + if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cut2ProngInvMassCandMin[n2] >= 0. && cut2ProngInvMassCandMax[n2] > 0.) { //no need to check isSelected2Prong but to avoid mistakes + if ((mass2ProngHypo1[n2] < cut2ProngInvMassCandMin[n2] || mass2ProngHypo1[n2] >= cut2ProngInvMassCandMax[n2]) && + (mass2ProngHypo2[n2] < cut2ProngInvMassCandMin[n2] || mass2ProngHypo2[n2] >= cut2ProngInvMassCandMax[n2])) { + isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); + if (b_debug) { + cutStatus2Prong[n2][iDebugCut] = false; + } + } } } + iDebugCut++; - if (isSelectedCand2Prong) { - // calculate invariant masses - auto arrMom = array{pvec0, pvec1}; - mass2PiK = RecoDecay::M(arrMom, array{massPi, massK}); - mass2KPi = RecoDecay::M(arrMom, array{massK, massPi}); - } + //secondary vertex reconstruction and further 2 prong selections + if (isSelected2ProngCand > 0 && df2.process(trackParVarPos1, trackParVarNeg1) > 0) { //should it be this or > 0 or are they equivalent + // get secondary vertex + const auto& secondaryVertex2 = df2.getPCACandidate(); + // get track momenta + array pvec0; + array pvec1; + df2.getTrack(0).getPxPyPzGlo(pvec0); + df2.getTrack(1).getPxPyPzGlo(pvec1); + + auto pVecCandProng2 = RecoDecay::PVec(pvec0, pvec1); - // invariant-mass cut - if (isSelectedCand2Prong && cut2ProngInvMassD0Min >= 0. && cut2ProngInvMassD0Max > 0.) { - if ((mass2PiK < cut2ProngInvMassD0Min || mass2PiK > cut2ProngInvMassD0Max) && - (mass2KPi < cut2ProngInvMassD0Min || mass2KPi > cut2ProngInvMassD0Max)) { - isSelectedCand2Prong = false; + // candidate pT cut + if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngPtCandMin), std::end(cut2ProngPtCandMin), [](double d) { return d >= 0.; }) > 0)) { + double cand2ProngPt = RecoDecay::Pt(pVecCandProng2); + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cand2ProngPt < cut2ProngPtCandMin[n2]) { + isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); + if (b_debug) { + cutStatus2Prong[n2][iDebugCut] = false; + } + } + } } - } + iDebugCut++; + + // imp. par. product cut + if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngImpParProductCandMax), std::end(cut2ProngImpParProductCandMax), [](double d) { return d < 100.; }) > 0)) { + auto impParProduct = trackPos1.dcaPrim0() * trackNeg1.dcaPrim0(); + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + if ((b_debug || (isSelected2ProngCand & 1 << n2)) && impParProduct > cut2ProngImpParProductCandMax[n2]) { + isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); + if (b_debug) { + cutStatus2Prong[n2][iDebugCut] = false; + } + } + } + } + iDebugCut++; - if (isSelectedCand2Prong) { - // fill table row - rowTrackIndexProng2(trackPos1.globalIndex(), - trackNeg1.globalIndex(), 1); + // CPA cut + if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngCPACandMin), std::end(cut2ProngCPACandMin), [](double d) { return d > -2.; }) > 0)) { + auto cpa = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex2, pVecCandProng2); + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cpa < cut2ProngCPACandMin[n2]) { + isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); + if (b_debug) { + cutStatus2Prong[n2][iDebugCut] = false; + } + } + } + } + iDebugCut++; + + if (isSelected2ProngCand > 0) { + // fill table row + rowTrackIndexProng2(trackPos1.globalIndex(), + trackNeg1.globalIndex(), isSelected2ProngCand); + if (b_debug) { + int Prong2CutStatus[n2ProngDecays]; + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + Prong2CutStatus[n2] = nCutStatus2ProngBit; + for (int n2cut = 0; n2cut < nCuts2Prong; n2cut++) { + if (!cutStatus2Prong[n2][n2cut]) { + Prong2CutStatus[n2] = Prong2CutStatus[n2] & ~(1 << n2cut); + } + } + } + rowProng2CutStatus(Prong2CutStatus[0], Prong2CutStatus[1]); //FIXME when we can do this by looping over n2ProngDecays + } - // fill histograms - if (b_dovalplots) { - hvtx2_x->Fill(secondaryVertex2[0]); - hvtx2_y->Fill(secondaryVertex2[1]); - hvtx2_z->Fill(secondaryVertex2[2]); - hmass2->Fill(mass2PiK); - hmass2->Fill(mass2KPi); + // fill histograms + if (b_dovalplots) { + + registry.get("hvtx2_x")->Fill(secondaryVertex2[0]); + registry.get("hvtx2_y")->Fill(secondaryVertex2[1]); + registry.get("hvtx2_z")->Fill(secondaryVertex2[2]); + arrMom = array{pvec0, pvec1}; + for (int n2 = 0; n2 < n2ProngDecays; n2++) { + if (isSelected2ProngCand & 1 << n2) { + if ((cut2ProngInvMassCandMin[n2] < 0. && cut2ProngInvMassCandMax[n2] <= 0.) || (mass2ProngHypo1[n2] >= cut2ProngInvMassCandMin[n2] && mass2ProngHypo1[n2] < cut2ProngInvMassCandMax[n2])) { + mass2ProngHypo1[n2] = RecoDecay::M(arrMom, arr2Mass1[n2]); + if (n2 == 0) { + registry.get("hmassD0")->Fill(mass2ProngHypo1[n2]); + } + if (n2 == 1) { + registry.get("hmassJpsi")->Fill(mass2ProngHypo1[n2]); + } + } + if ((cut2ProngInvMassCandMin[n2] < 0. && cut2ProngInvMassCandMax[n2] <= 0.) || (mass2ProngHypo2[n2] >= cut2ProngInvMassCandMin[n2] && mass2ProngHypo2[n2] < cut2ProngInvMassCandMax[n2])) { + mass2ProngHypo2[n2] = RecoDecay::M(arrMom, arr2Mass2[n2]); + if (n2 == 0) { + registry.get("hmassD0")->Fill(mass2ProngHypo1[n2]); + } + } + } + } + } } } } // 3-prong vertex reconstruction if (do3prong == 1) { - if (trackPos1.isSel3Prong() == 0) { - continue; - } - if (trackNeg1.isSel3Prong() == 0) { + if (!sel3ProngStatusPos1 || !sel3ProngStatusNeg1) { continue; } @@ -329,36 +481,50 @@ struct HFTrackIndexSkimsCreator { if (trackPos2.signed1Pt() < 0) { continue; } - if (trackPos2.isSel3Prong() == 0) { + if (!(trackPos2.isSelProng() & (1 << 1))) { continue; } - auto pVecCandProng3Pos = RecoDecay::PVec(pVecCandProng2, array{trackPos2.px(), trackPos2.py(), trackPos2.pz()}); + int isSelected3ProngCand = n3ProngBit; - // candidate pT cut - if (RecoDecay::Pt(pVecCandProng3Pos) < cut3ProngPtCandMin) { - continue; + if (b_debug) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { + cutStatus3Prong[n3][n3cut] = true; + } + } } - - // invariant-mass cut - if (cut3ProngInvMassDPlusMin >= 0. && cut3ProngInvMassDPlusMax > 0.) { - // calculate invariant mass - auto arr3Mom = array{ - array{trackPos1.px(), trackPos1.py(), trackPos1.pz()}, - array{trackNeg1.px(), trackNeg1.py(), trackNeg1.pz()}, - array{trackPos2.px(), trackPos2.py(), trackPos2.pz()}}; - mass3PiKPi = RecoDecay::M(std::move(arr3Mom), array{massPi, massK, massPi}); - if (mass3PiKPi < cut3ProngInvMassDPlusMin || mass3PiKPi > cut3ProngInvMassDPlusMax) { - continue; + int iDebugCut = 0; + + // 3prong invariant-mass cut + auto arr3Mom = array{ + array{trackPos1.px(), trackPos1.py(), trackPos1.pz()}, + array{trackNeg1.px(), trackNeg1.py(), trackNeg1.pz()}, + array{trackPos2.px(), trackPos2.py(), trackPos2.pz()}}; + + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + mass3ProngHypo1[n3] = RecoDecay::M(arr3Mom, arr3Mass1[n3]); + mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); + if ((isSelected3ProngCand & 1 << n3) && cut3ProngInvMassCandMin[n3] >= 0. && cut3ProngInvMassCandMax[n3] > 0.) { + if ((mass3ProngHypo1[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo1[n3] >= cut3ProngInvMassCandMax[n3]) && + (mass3ProngHypo2[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo2[n3] >= cut3ProngInvMassCandMax[n3])) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } } } + if (!b_debug && isSelected3ProngCand == 0) { + continue; + } + iDebugCut++; // reconstruct the 3-prong secondary vertex auto trackParVarPos2 = getTrackParCov(trackPos2); if (df3.process(trackParVarPos1, trackParVarNeg1, trackParVarPos2) == 0) { continue; } - // get secondary vertex const auto& secondaryVertex3 = df3.getPCACandidate(); // get track momenta @@ -369,36 +535,109 @@ struct HFTrackIndexSkimsCreator { df3.getTrack(1).getPxPyPzGlo(pvec1); df3.getTrack(2).getPxPyPzGlo(pvec2); + auto pVecCandProng3Pos = RecoDecay::PVec(pvec0, pvec1, pvec2); + + // candidate pT cut + if (std::count_if(std::begin(cut3ProngPtCandMin), std::end(cut3ProngPtCandMin), [](double d) { return d >= 0.; }) > 0) { + double cand3ProngPt = RecoDecay::Pt(pVecCandProng3Pos); + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if (cand3ProngPt < cut3ProngPtCandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + } + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + if (!b_debug && isSelected3ProngCand == 0) { + continue; //this and all further instances should be changed if 4 track loop is added + } + } + iDebugCut++; + // CPA cut - if (cut3ProngCPAMin > -2.) { - pVecCandProng3Pos = RecoDecay::PVec(pvec0, pvec1, pvec2); + if (std::count_if(std::begin(cut3ProngCPACandMin), std::end(cut3ProngCPACandMin), [](double d) { return d > -2.; }) > 0) { auto cpa = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex3, pVecCandProng3Pos); - if (cpa < cut3ProngCPAMin) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if ((isSelected3ProngCand & 1 << n3) && cpa < cut3ProngCPACandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + } + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + if (!b_debug && isSelected3ProngCand == 0) { continue; } } + iDebugCut++; // decay length cut - if (cut3ProngDecLenMin > 0.) { + if (std::count_if(std::begin(cut3ProngDecLenCandMin), std::end(cut3ProngDecLenCandMin), [](double d) { return d > 0.; }) > 0) { auto decayLength = RecoDecay::distance(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex3); - if (decayLength < cut3ProngDecLenMin) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if ((isSelected3ProngCand & 1 << n3) && decayLength < cut3ProngDecLenCandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + } + if (!b_debug && isSelected3ProngCand == 0) { continue; } } + iDebugCut++; // fill table row rowTrackIndexProng3(trackPos1.globalIndex(), trackNeg1.globalIndex(), - trackPos2.globalIndex(), 2); + trackPos2.globalIndex(), isSelected3ProngCand); + + if (b_debug) { + int Prong3CutStatus[n3ProngDecays]; + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + Prong3CutStatus[n3] = nCutStatus3ProngBit; + for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { + if (!cutStatus3Prong[n3][n3cut]) { + Prong3CutStatus[n3] = Prong3CutStatus[n3] & ~(1 << n3cut); + } + } + } + rowProng3CutStatus(Prong3CutStatus[0], Prong3CutStatus[1], Prong3CutStatus[2]); //FIXME when we can do this by looping over n3ProngDecays + } // fill histograms if (b_dovalplots) { - hvtx3_x->Fill(secondaryVertex3[0]); - hvtx3_y->Fill(secondaryVertex3[1]); - hvtx3_z->Fill(secondaryVertex3[2]); - // calculate invariant mass - hmass3->Fill(RecoDecay::M(array{pvec0, pvec1, pvec2}, array{massPi, massK, massPi})); + registry.get("hvtx3_x")->Fill(secondaryVertex3[0]); + registry.get("hvtx3_y")->Fill(secondaryVertex3[1]); + registry.get("hvtx3_z")->Fill(secondaryVertex3[2]); + arr3Mom = array{pvec0, pvec1, pvec2}; + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if (isSelected3ProngCand & 1 << n3) { + if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo1[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo1[n3] < cut3ProngInvMassCandMax[n3])) { + mass3ProngHypo1[n3] = RecoDecay::M(arr3Mom, arr3Mass1[n3]); + if (n3 == 0) { + registry.get("hmassDPlus")->Fill(mass3ProngHypo1[n3]); + } + if (n3 == 1) { + registry.get("hmassLc")->Fill(mass3ProngHypo1[n3]); + } + if (n3 == 2) { + registry.get("hmassDs")->Fill(mass3ProngHypo1[n3]); + } + } + if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo2[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo2[n3] < cut3ProngInvMassCandMax[n3])) { + mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); + if (n3 == 1) { + registry.get("hmassLc")->Fill(mass3ProngHypo2[n3]); + } + if (n3 == 2) { + registry.get("hmassDs")->Fill(mass3ProngHypo2[n3]); + } + } + } + } } } @@ -408,29 +647,44 @@ struct HFTrackIndexSkimsCreator { if (trackNeg2.signed1Pt() > 0) { continue; } - if (trackNeg2.isSel3Prong() == 0) { + if (!(trackNeg2.isSelProng() & (1 << 1))) { continue; } - auto pVecCandProng3Neg = RecoDecay::PVec(pVecCandProng2, array{trackNeg2.px(), trackNeg2.py(), trackNeg2.pz()}); + int isSelected3ProngCand = n3ProngBit; - // candidate pT cut - if (RecoDecay::Pt(pVecCandProng3Neg) < cut3ProngPtCandMin) { - continue; + if (b_debug) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { + cutStatus3Prong[n3][n3cut] = true; + } + } } - - // invariant-mass cut - if (cut3ProngInvMassDPlusMin >= 0. && cut3ProngInvMassDPlusMax > 0.) { - // calculate invariant mass - auto arr3Mom = array{ - array{trackNeg1.px(), trackNeg1.py(), trackNeg1.pz()}, - array{trackPos1.px(), trackPos1.py(), trackPos1.pz()}, - array{trackNeg2.px(), trackNeg2.py(), trackNeg2.pz()}}; - mass3PiKPi = RecoDecay::M(std::move(arr3Mom), array{massPi, massK, massPi}); - if (mass3PiKPi < cut3ProngInvMassDPlusMin || mass3PiKPi > cut3ProngInvMassDPlusMax) { - continue; + int iDebugCut = 0; + + // 3prong invariant-mass cut + auto arr3Mom = array{ + array{trackNeg1.px(), trackNeg1.py(), trackNeg1.pz()}, + array{trackPos1.px(), trackPos1.py(), trackPos1.pz()}, + array{trackNeg2.px(), trackNeg2.py(), trackNeg2.pz()}}; + + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + mass3ProngHypo1[n3] = RecoDecay::M(arr3Mom, arr3Mass1[n3]); + mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); + if ((isSelected3ProngCand & 1 << n3) && cut3ProngInvMassCandMin[n3] >= 0. && cut3ProngInvMassCandMax[n3] > 0.) { + if ((mass3ProngHypo1[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo1[n3] >= cut3ProngInvMassCandMax[n3]) && + (mass3ProngHypo2[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo2[n3] >= cut3ProngInvMassCandMax[n3])) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } } } + if (!b_debug && isSelected3ProngCand == 0) { + continue; + } + iDebugCut++; // reconstruct the 3-prong secondary vertex auto trackParVarNeg2 = getTrackParCov(trackNeg2); @@ -448,51 +702,124 @@ struct HFTrackIndexSkimsCreator { df3.getTrack(1).getPxPyPzGlo(pvec1); df3.getTrack(2).getPxPyPzGlo(pvec2); + auto pVecCandProng3Neg = RecoDecay::PVec(pvec0, pvec1, pvec2); + + // candidate pT cut + if (std::count_if(std::begin(cut3ProngPtCandMin), std::end(cut3ProngPtCandMin), [](double d) { return d >= 0.; }) > 0) { + double cand3ProngPt = RecoDecay::Pt(pVecCandProng3Neg); + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if (cand3ProngPt < cut3ProngPtCandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + } + if (!b_debug && isSelected3ProngCand == 0) { + continue; + } + } + iDebugCut++; + // CPA cut - if (cut3ProngCPAMin > -2.) { - pVecCandProng3Neg = RecoDecay::PVec(pvec0, pvec1, pvec2); + if (std::count_if(std::begin(cut3ProngCPACandMin), std::end(cut3ProngCPACandMin), [](double d) { return d > -2.; }) > 0) { auto cpa = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex3, pVecCandProng3Neg); - if (cpa < cut3ProngCPAMin) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if ((isSelected3ProngCand & 1 << n3) && cpa < cut3ProngCPACandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + } + if (!b_debug && isSelected3ProngCand == 0) { continue; } } + iDebugCut++; // decay length cut - if (cut3ProngDecLenMin > 0.) { + if (std::count_if(std::begin(cut3ProngDecLenCandMin), std::end(cut3ProngDecLenCandMin), [](double d) { return d > 0.; }) > 0) { auto decayLength = RecoDecay::distance(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex3); - if (decayLength < cut3ProngDecLenMin) { + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if ((isSelected3ProngCand & 1 << n3) && decayLength < cut3ProngDecLenCandMin[n3]) { + isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); + if (b_debug) { + cutStatus3Prong[n3][iDebugCut] = false; + } + } + } + if (!b_debug && isSelected3ProngCand == 0) { continue; } } + iDebugCut++; // fill table row rowTrackIndexProng3(trackNeg1.globalIndex(), trackPos1.globalIndex(), - trackNeg2.globalIndex(), 2); + trackNeg2.globalIndex(), isSelected3ProngCand); + + if (b_debug) { + int Prong3CutStatus[n3ProngDecays]; + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + Prong3CutStatus[n3] = nCutStatus3ProngBit; + for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { + if (!cutStatus3Prong[n3][n3cut]) { + Prong3CutStatus[n3] = Prong3CutStatus[n3] & ~(1 << n3cut); + } + } + } + rowProng3CutStatus(Prong3CutStatus[0], Prong3CutStatus[1], Prong3CutStatus[2]); //FIXME when we can do this by looping over n3ProngDecays + } // fill histograms if (b_dovalplots) { - hvtx3_x->Fill(secondaryVertex3[0]); - hvtx3_y->Fill(secondaryVertex3[1]); - hvtx3_z->Fill(secondaryVertex3[2]); - // calculate invariant mass - hmass3->Fill(RecoDecay::M(array{pvec0, pvec1, pvec2}, array{massPi, massK, massPi})); + registry.get("hvtx3_x")->Fill(secondaryVertex3[0]); + registry.get("hvtx3_y")->Fill(secondaryVertex3[1]); + registry.get("hvtx3_z")->Fill(secondaryVertex3[2]); + arr3Mom = array{pvec0, pvec1, pvec2}; + for (int n3 = 0; n3 < n3ProngDecays; n3++) { + if (isSelected3ProngCand & 1 << n3) { + if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo1[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo1[n3] < cut3ProngInvMassCandMax[n3])) { + mass3ProngHypo1[n3] = RecoDecay::M(arr3Mom, arr3Mass1[n3]); + if (n3 == 0) { + registry.get("hmassDPlus")->Fill(mass3ProngHypo1[n3]); + } + if (n3 == 1) { + registry.get("hmassLc")->Fill(mass3ProngHypo1[n3]); + } + if (n3 == 2) { + registry.get("hmassDs")->Fill(mass3ProngHypo1[n3]); + } + } + if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo2[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo2[n3] < cut3ProngInvMassCandMax[n3])) { + mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); + if (n3 == 1) { + registry.get("hmassLc")->Fill(mass3ProngHypo2[n3]); + } + if (n3 == 2) { + registry.get("hmassDs")->Fill(mass3ProngHypo2[n3]); + } + } + } + } } } } } } - /* - auto nTracks = tracks.size(); // number of tracks in this collision + + auto nTracks = tracks.size(); // number of tracks passing 2 and 3 prong selection in this collision nCand2 = rowTrackIndexProng2.lastIndex() - nCand2; // number of 2-prong candidates in this collision nCand3 = rowTrackIndexProng3.lastIndex() - nCand3; // number of 3-prong candidates in this collision - hNTracks->Fill(nTracks); - hNCand2Prong->Fill(nCand2); - hNCand3Prong->Fill(nCand3); - hNCand2ProngVsNTracks->Fill(nTracks, nCand2); - hNCand3ProngVsNTracks->Fill(nTracks, nCand3); - */ + + registry.get("hNTracks")->Fill(nTracks); + registry.get("hNCand2Prong")->Fill(nCand2); + registry.get("hNCand3Prong")->Fill(nCand3); + registry.get("hNCand2ProngVsNTracks")->Fill(nTracks, nCand2); + registry.get("hNCand3ProngVsNTracks")->Fill(nTracks, nCand3); } };