diff --git a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc index 09de461fe9..4797d957ec 100644 --- a/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc +++ b/offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc @@ -390,15 +390,39 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay) m_reco_track_eta[index] = m_dst_track->get_eta(); m_reco_track_phi[index] = m_dst_track->get_phi(); m_reco_track_chi2nDoF[index] = m_dst_track->get_chisq() / m_dst_track->get_ndf(); - if (m_dst_track->get_silicon_seed()) - { - m_reco_track_silicon_seeds[index] = static_cast(m_dst_track->get_silicon_seed()->size_cluster_keys()); - } - else + m_reco_track_silicon_seeds[index] = 0; + m_reco_track_tpc_seeds[index] = 0; + + for (auto state_iter = m_dst_track->begin_states(); + state_iter != m_dst_track->end_states(); + ++state_iter) { - m_reco_track_silicon_seeds[index] = 0; + SvtxTrackState *tstate = state_iter->second; + if (tstate->get_pathlength() != 0) // The first track state is an extrapolation so has no cluster + { + auto stateckey = tstate->get_cluskey(); + if (stateckey == TrkrDefs::CLUSKEYMAX) + { + continue; + } + uint8_t id = TrkrDefs::getTrkrId(stateckey); + + switch (id) + { + case TrkrDefs::mvtxId: + [[fallthrough]]; + case TrkrDefs::inttId: + ++m_reco_track_silicon_seeds[index]; + break; + case TrkrDefs::tpcId: + ++m_reco_track_tpc_seeds[index]; + break; + default: + break; + } + } } - m_reco_track_tpc_seeds[index] = static_cast(m_dst_track->get_tpc_seed()->size_cluster_keys()); + m_min_reco_track_pT = std::min(m_reco_track_pT[index], m_min_reco_track_pT); m_max_reco_track_pT = std::max(m_reco_track_pT[index], m_max_reco_track_pT); @@ -513,8 +537,8 @@ void HFTrackEfficiency::resetBranches() m_reco_track_phi[iTrack] = std::numeric_limits::quiet_NaN(); m_true_track_PID[iTrack] = std::numeric_limits::quiet_NaN(); m_reco_track_chi2nDoF[iTrack] = std::numeric_limits::quiet_NaN(); - m_reco_track_silicon_seeds[iTrack] = 0; - m_reco_track_tpc_seeds[iTrack] = 0; + m_reco_track_silicon_seeds[iTrack] = -1; + m_reco_track_tpc_seeds[iTrack] = -1; } m_primary_vtx_x = std::numeric_limits::quiet_NaN(); diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.cc index 9bb76608ee..8877c5be00 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.cc @@ -327,6 +327,14 @@ std::vector KFParticle_Tools::makeAllDaughterParticles(PHCompositeNo } } + if (m_verbosity > 100) + { + printSelectionCheck("MVTX states", m_nMVTXStates, MVTX_states, 5); + printSelectionCheck("INTT states", m_nINTTStates, INTT_states, 5); + printSelectionCheck("TPC states", m_nTPCStates, TPC_states, 100); + printSelectionCheck("TPOT states", m_nTPOTStates, TPOT_states, 5); + } + if (MVTX_states < m_nMVTXStates) { continue; @@ -460,6 +468,23 @@ int KFParticle_Tools::getTracksFromVertex(PHCompositeNode *topNode, const KFPart { goodTrack = true; } + + + if (m_verbosity >= 10) + { + printSelectionCheck("This track", "passed", "failed", "the selection", goodTrack); + if (m_verbosity >= 11) + { + printSelectionCheck("Track pT", m_track_min_pt, pt, m_track_max_pt); + printSelectionCheck("Track pT chi^2", 0, ptchi2, m_track_ptchi2); + printSelectionCheck("IP", m_track_ip, min_ip, FLT_MAX); + printSelectionCheck("IP chi^2", m_track_ipchi2, min_ipchi2, FLT_MAX); + printSelectionCheck("IP xy", m_track_ip_xy, min_ip_xy, FLT_MAX); + printSelectionCheck("IP xy chi^2", m_track_ipchi2_xy, min_ipchi2_xy, FLT_MAX); + printSelectionCheck("Track chi^2/nDoF", 0, trackchi2ndof, m_track_chi2ndof); + } + } + return goodTrack; } @@ -513,7 +538,7 @@ std::vector KFParticle_Tools::findAllGoodTracks(const std::vector> KFParticle_Tools::findTwoProngs(std::vector daughterParticles, std::vector goodTrackIndex, int nTracks) const +std::vector> KFParticle_Tools::findTwoProngs(std::vector daughterParticles, std::vector goodTrackIndex, int nTracks) { std::vector> goodTracksThatMeet; @@ -526,6 +551,16 @@ std::vector> KFParticle_Tools::findTwoProngs(std::vector= 10) + { + printSelectionCheck("This track pair", "passed", "failed", "the DCA selection", (dca <= m_comb_DCA) && (dca_xy <= m_comb_DCA_xy)); + if (m_verbosity >= 11) + { + printSelectionCheck("Pair DCA", 0., dca, m_comb_DCA); + printSelectionCheck("Pair DCA xy", 0., dca_xy, m_comb_DCA_xy); + } + } + if (dca <= m_comb_DCA && dca_xy <= m_comb_DCA_xy) { KFVertex twoParticleVertex; @@ -535,6 +570,16 @@ std::vector> KFParticle_Tools::findTwoProngs(std::vector combination = {*i_it, *j_it}; + if (nTracks == 2 && m_verbosity >= 10) + { + printSelectionCheck("This track pair", "passed", "failed", "the quality and radius selection", (vertexchi2ndof <= m_vertex_chi2ndof) && (sv_radial_position >= m_min_radial_SV)); + if (m_verbosity >= 11) + { + printSelectionCheck("SV chi^2/nDoFA", 0., vertexchi2ndof, m_vertex_chi2ndof); + printSelectionCheck("SV radius", m_min_radial_SV, sv_radial_position, FLT_MAX); + } + } + if (nTracks == 2 && vertexchi2ndof > m_vertex_chi2ndof) { continue; @@ -581,6 +626,16 @@ std::vector> KFParticle_Tools::findNProngs(std::vector= 10) + { + printSelectionCheck("This track", "combined", "did not combine", "with a SV set", (dca <= m_comb_DCA) && (dca_xy <= m_comb_DCA_xy)); + if (m_verbosity >= 11) + { + printSelectionCheck("Pair DCA", 0., dca, m_comb_DCA); + printSelectionCheck("Pair DCA xy", 0., dca_xy, m_comb_DCA_xy); + } + } + if (dca > m_comb_DCA || dca_xy > m_comb_DCA_xy) { dcaMet = false; @@ -601,6 +656,16 @@ std::vector> KFParticle_Tools::findNProngs(std::vector= 10) + { + printSelectionCheck("This SV combination", "passed", "failed", "the quality and radius selection", (vertexchi2ndof <= m_vertex_chi2ndof) && (sv_radial_position >= m_min_radial_SV)); + if (m_verbosity >= 11) + { + printSelectionCheck("SV chi^2/nDoFA", 0., vertexchi2ndof, m_vertex_chi2ndof); + printSelectionCheck("SV radius", m_min_radial_SV, sv_radial_position, FLT_MAX); + } + } + if ((unsigned int) nRequiredTracks == nProngs && vertexchi2ndof > m_vertex_chi2ndof) { continue; @@ -825,6 +890,12 @@ std::tuple KFParticle_Tools::buildMother(KFParticle vDaughters float calculated_dEdx_value = get_dEdx(topNode, vDaughters[i]); double expected_dEdx_value = get_dEdx_fitValue((Int_t) vDaughters[i].GetQ() * vDaughters[i].GetP(), track_PDG_ID); bool accept_dEdx = isInRange((1 - m_dEdx_band_width) * expected_dEdx_value, calculated_dEdx_value, (1 + m_dEdx_band_width) * expected_dEdx_value); + + if (m_verbosity >= 11) + { + printSelectionCheck("dE/dx check", (1 - m_dEdx_band_width) * expected_dEdx_value, calculated_dEdx_value, (1 + m_dEdx_band_width) * expected_dEdx_value); + } + if (!accept_dEdx) { delete[] inputTracks; @@ -875,7 +946,9 @@ std::tuple KFParticle_Tools::buildMother(KFParticle vDaughters float calculated_mass; float calculated_mass_err; mother.GetMass(calculated_mass, calculated_mass_err); - float calculated_pt = mother.GetPt(); + float calculated_pt; + float calculated_pt_err; + mother.GetPt(calculated_pt, calculated_pt_err); float min_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].first : m_min_mass; float max_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].second : m_max_mass; @@ -909,6 +982,12 @@ std::tuple KFParticle_Tools::buildMother(KFParticle vDaughters { goodCandidate = false; } + + if (m_verbosity >= 11) + { + bool accept = crossings.size() == 1; + printSelectionCheck("", "All tracks are from the same BC", "Tracks are from different BC", "", accept); + } } // Check the requirements of an intermediate states against this mother and re-do goodCandidate @@ -923,6 +1002,28 @@ std::tuple KFParticle_Tools::buildMother(KFParticle vDaughters { goodCandidate = false; } + + if (m_verbosity >= 10) + { + printSelectionCheck("", "Accepted", "Rejected", "the intermediate selection", goodCandidate); + if (m_verbosity >= 11) + { + printSelectionCheck("Intermediate DIRA", m_intermediate_min_dira[k], intermediate_DIRA, FLT_MAX); + printSelectionCheck("Intermediate FD chi^2", m_intermediate_min_fdchi2[k], intermediate_FDchi2, FLT_MAX); + } + } + } + } + + if (m_verbosity >= 10) + { + printSelectionCheck("", "Accepted", "Rejected", "the mother selection", goodCandidate); + if (m_verbosity >= 11) + { + printSelectionCheck("Vertex charge is", "right", "wrong", "", chargeCheck); + printSelectionCheck("Invariant Mass", min_mass, calculated_mass, max_mass); + printSelectionCheck("Mother pT", min_pt, calculated_pt, FLT_MAX); + printSelectionCheck("Mother SV volume", 0., calculateEllipsoidVolume(mother), max_vertex_volume); } } delete[] inputTracks; @@ -970,6 +1071,28 @@ void KFParticle_Tools::constrainToVertex(KFParticle &particle, bool &goodCandida { goodCandidate = true; } + + if (m_verbosity >= 10) + { + printSelectionCheck("", "Passed", "Failed", "the PV constraint", goodCandidate); + if (m_verbosity >= 11) + { + printSelectionCheck("Mother DIRA", m_dira_min, calculated_dira, m_dira_max); + printSelectionCheck("Mother DIRA xy", m_dira_xy_min, calculated_dira_xy, m_dira_xy_max); + printSelectionCheck("Mother FD chi^2", m_fdchi2, calculated_fdchi2, FLT_MAX); + printSelectionCheck("Mother IP", 0, calculated_ip, m_mother_ip); + printSelectionCheck("Mother IP chi^2", 0., calculated_ipchi2, m_mother_ipchi2); + printSelectionCheck("Mother IP xy", 0., calculated_ip_xy, m_mother_ip_xy); + printSelectionCheck("Mother IP xy chi^2", 0., calculated_ipchi2_xy, m_mother_ipchi2_xy); + printSelectionCheck("Mother Decay Time", m_min_decayTime, calculated_decayTime, m_max_decayTime); + printSelectionCheck("Mother Decay Time Significance", m_mother_min_decay_time_significance, calculated_decay_time_significance, FLT_MAX); + printSelectionCheck("Mother Decay Time xy", m_min_decayTime_xy, calculated_decayTime_xy, m_max_decayTime_xy); + printSelectionCheck("Mother Decay Length", m_min_decayLength, calculated_decayLength, m_max_decayLength); + printSelectionCheck("Mother Decay Length Significance", m_mother_min_decay_length_significance, calculated_decay_length_significance, FLT_MAX); + printSelectionCheck("Mother Decay Length xy", m_min_decayLength_xy, calculated_decayLength_xy, m_max_decayLength_xy); + printSelectionCheck("Mother Decay Length xy Significance", m_mother_min_decay_length_xy_significance, calculated_decay_length_xy_significance, FLT_MAX); + } + } } std::tuple KFParticle_Tools::getCombination(KFParticle vDaughters[], int daughterOrder[], KFParticle vertex, bool constrain_to_vertex, bool isIntermediate, int intermediateNumber, int nTracks, bool constrainMass, float required_vertexID, PHCompositeNode *topNode) @@ -1210,7 +1333,10 @@ void KFParticle_Tools::init_dEdx_fits() if (m_use_local_PID_file) { - std::cout << PHWHERE << " using local file " << m_local_PID_filename << std::endl; + if (m_verbosity > 4) + { + std::cout << PHWHERE << " using local file " << m_local_PID_filename << std::endl; + } // new method is independent of charge filefit->GetObject("pi_band",f_pion_plus); filefit->GetObject("K_band",f_kaon_plus); @@ -1299,3 +1425,25 @@ bool KFParticle_Tools::checkTrackAndVertexMatch(KFParticle vDaughters[], int nTr return vertexAndTrackMatch; } + +void KFParticle_Tools::printSelectionCheck(const std::string ¶meter, float min, float val, float max) +{ + std::string trailer = "the " + parameter + " requirement\033[0m"; + std::string passOrFail = isInRange(min, val, max) ? "\033[1;" + accept_colour + "mPassed " + trailer + : "\033[1;" + reject_colour + "mFailed " + trailer; + std::cout << passOrFail << ". Lower bound = " << min << ", measured value = " << val << ", upper bound = " << max << std::endl; +} + +void KFParticle_Tools::printSelectionCheck(const std::string &start, const std::string &accept, const std::string &reject, const std::string &end, bool equality) +{ + std::string decision = equality ? accept : reject; + std::string colour = equality ? accept_colour : reject_colour; + std::string spacing = start.empty() ? "" : " "; + std::cout << "\033[1;" << colour << "m" << start << spacing << decision << " " << end << "\033[0m" << std::endl; +} + +void KFParticle_Tools::printSelectionCheck(const std::string &info, unsigned int value) +{ + std::string colour = value > 0 ? accept_colour : reject_colour; + std::cout << info << " = \033[1;" << colour << "m" << value << "\033[0m" << std::endl; +} diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.h b/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.h index feabd4aafc..434a64cf4d 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.h +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_Tools.h @@ -75,7 +75,7 @@ class KFParticle_Tools : protected KFParticle_MVA std::vector findAllGoodTracks(const std::vector &daughterParticles, const std::vector &primaryVertices); - std::vector> findTwoProngs(std::vector daughterParticles, std::vector goodTrackIndex, int nTracks) const; + std::vector> findTwoProngs(std::vector daughterParticles, std::vector goodTrackIndex, int nTracks); std::vector> findNProngs(std::vector daughterParticles, const std::vector &goodTrackIndex, @@ -126,6 +126,8 @@ class KFParticle_Tools : protected KFParticle_MVA void set_dont_use_global_vertex(bool set_variable) { m_dont_use_global_vertex = set_variable; } protected: + int m_verbosity = 0; + std::string m_mother_name_Tools; int m_num_intermediate_states{-1}; std::vector m_num_tracks_from_intermediate; @@ -281,6 +283,12 @@ class KFParticle_Tools : protected KFParticle_MVA void removeDuplicates(std::vector &v); void removeDuplicates(std::vector> &v); void removeDuplicates(std::vector> &v); + + void printSelectionCheck(const std::string ¶meter, float min, float val, float max); + void printSelectionCheck(const std::string &start, const std::string &accept, const std::string &reject, const std::string &end, bool equality); + void printSelectionCheck(const std::string &info, unsigned int value); + std::string accept_colour = "32"; + std::string reject_colour = "31"; }; #endif // KFPARTICLESPHENIX_KFPARTICLETOOLS_H diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_eventReconstruction.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_eventReconstruction.cc index 64d6b932fa..3c40fd4e58 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_eventReconstruction.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_eventReconstruction.cc @@ -76,6 +76,13 @@ void KFParticle_eventReconstruction::createDecay(PHCompositeNode* topNode, std:: std::vector goodTrackIndex = findAllGoodTracks(daughterParticles, primaryVertices); + if (m_verbosity >= 10) + { + printSelectionCheck("Number of daughters passing state selection", daughterParticles.size()); + printSelectionCheck("Number of daughters passing track selection", goodTrackIndex.size()); + printSelectionCheck("Number of PVs passing selection", primaryVertices.size()); + } + if (!m_has_intermediates) { buildBasicChain(selectedMother, selectedVertex, selectedDaughters, daughterParticles, goodTrackIndex, primaryVertices, topNode); @@ -102,6 +109,11 @@ void KFParticle_eventReconstruction::buildBasicChain(std::vector& se goodTracksThatMeet = findNProngs(daughterParticlesBasic, goodTrackIndexBasic, goodTracksThatMeet, m_num_tracks, p); } + if (m_verbosity >= 10) + { + printSelectionCheck("Number of SVs passing selection", goodTracksThatMeet.size()); + } + getCandidateDecay(selectedMotherBasic, selectedVertexBasic, selectedDaughtersBasic, daughterParticlesBasic, goodTracksThatMeet, primaryVerticesBasic, 0, m_num_tracks, false, 0, true, topNode); } @@ -136,6 +148,12 @@ void KFParticle_eventReconstruction::buildChain(std::vector& selecte goodTracksThatMeet, m_num_tracks_from_intermediate[i], p); } + + if (m_verbosity >= 10) + { + printSelectionCheck("Number of SVs passing selection", goodTracksThatMeet.size()); + } + getCandidateDecay(potentialIntermediates[i], vertices, potentialDaughters[i], daughterParticlesAdv, goodTracksThatMeet, primaryVerticesAdv, track_start, track_stop, true, i, m_constrain_int_mass, topNode); track_start += track_stop; diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc index b3b119e4f6..bc9f0b4c3a 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc @@ -303,9 +303,9 @@ void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode) m_tree->Branch("runNumber", &m_runNumber, "runNumber/I"); m_tree->Branch("eventNumber", &m_evtNumber, "eventNumber/I"); - m_tree->Branch("event_bco", &m_event_bco, "event_bco/L"); //adding for the current event BCO, not shifted - m_tree->Branch("BCO", &m_bco, "BCO/L"); //already there, this is shifted BCO - m_tree->Branch("last_event_bco", &m_last_event_bco, "last_event_bco/L"); //BCO for the last event + m_tree->Branch("Collision_BCO", &m_bco, "Collision_BCO/L"); //already there, this is shifted BCO + m_tree->Branch("GL1_BCO", &m_event_bco, "GL1_BCO/L"); //adding for the current event BCO, not shifted + m_tree->Branch("last_GL1_BCO", &m_last_event_bco, "last_GL1_BCO/L"); //BCO for the last event if (m_get_trigger_info) { diff --git a/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc b/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc index db46bffafe..1b597afcef 100644 --- a/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc +++ b/offline/packages/KFParticle_sPHENIX/KFParticle_sPHENIX.cc @@ -88,6 +88,7 @@ KFParticle_sPHENIX::KFParticle_sPHENIX(const std::string &name) int KFParticle_sPHENIX::Init(PHCompositeNode *topNode) { + m_verbosity = Verbosity(); if (m_save_output && Verbosity() >= VERBOSITY_SOME) { @@ -135,26 +136,13 @@ int KFParticle_sPHENIX::InitRun(PHCompositeNode *topNode) } int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) -{ - +{ std::vector mother; std::vector vertex_kfparticle; std::vector> daughters; std::vector> intermediates; int nPVs; int multiplicity; - - SvtxTrackMap *check_trackmap = findNode::getClass(topNode, m_trk_map_node_name); - multiplicity = check_trackmap->size(); - - if (check_trackmap->empty()) - { - if (Verbosity() >= VERBOSITY_SOME) - { - std::cout << "KFParticle: Event skipped as there are no tracks" << std::endl; - } - return Fun4AllReturnCodes::EVENT_OK; - } // Adding BCO Matching auto* evtHeader = findNode::getClass(topNode, "EventHeader"); // event header node @@ -162,57 +150,70 @@ int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) if (!gl1packet) { - gl1packet = findNode::getClass(topNode, "GL1Packet"); + gl1packet = findNode::getClass(topNode, "GL1Packet"); } if (evtHeader && gl1packet) { - const int64_t run = evtHeader->get_RunNumber(); - const int64_t evn = evtHeader->get_EvtSequence(); - m_this_event_bco = static_cast(gl1packet->lValue(0, "BCO")); - - if (Verbosity() >= VERBOSITY_SOME) - { - std::cout << "Event start | run: " << run << " event: " << evn << " this_event_bco: " << m_this_event_bco << std::endl; - } + const int64_t run = evtHeader->get_RunNumber(); + const int64_t evn = evtHeader->get_EvtSequence(); + m_this_event_bco = static_cast(gl1packet->lValue(0, "BCO")); - if (run != m_prev_runNumber || evn != m_prev_eventNumber) - { + if (Verbosity() >= VERBOSITY_A_LOT) + { + std::cout << "Event start | run: " << run << " event: " << evn << " this_event_bco: " << m_this_event_bco << std::endl; + } - if (Verbosity() >= VERBOSITY_SOME) + if (run != m_prev_runNumber || evn != m_prev_eventNumber) { - std::cout << "New event detected" << std::endl; - std::cout << "Previous event BCO: " << m_prev_event_bco << std::endl; - } - //m_last_event_bco = m_prev_event_bco; - m_last_event_bco = (run == m_prev_runNumber) ? m_prev_event_bco : -1; - m_prev_event_bco = m_this_event_bco; - m_prev_runNumber = run; - m_prev_eventNumber = evn; + if (Verbosity() >= VERBOSITY_A_LOT) + { + std::cout << "New event detected" << std::endl; + std::cout << "Previous event BCO: " << m_prev_event_bco << std::endl; + } + + m_last_event_bco = (run == m_prev_runNumber) ? m_prev_event_bco : -1; + m_prev_event_bco = m_this_event_bco; - if (Verbosity() >= VERBOSITY_SOME) - { - std::cout << "Updated values | last_event_bco: " << m_last_event_bco - << " stored_prev_event_bco: " << m_prev_event_bco - << std::endl; - } + m_prev_runNumber = run; + m_prev_eventNumber = evn; + + if (Verbosity() >= VERBOSITY_A_LOT) + { + std::cout << "Updated values | last_event_bco: " << m_last_event_bco + << " stored_prev_event_bco: " << m_prev_event_bco + << std::endl; + } + } } - } else { - if (Verbosity() >= VERBOSITY_SOME) + if (Verbosity() >= VERBOSITY_MORE) + { + std::cout << "KFParticle: EventHeader or GL1 packet not found" << std::endl; + } + m_this_event_bco = 0; + m_last_event_bco = 0; + m_prev_event_bco = 0; + m_prev_runNumber = 0; + m_prev_eventNumber = 0; + } + // End BCO matching here + + SvtxTrackMap *check_trackmap = findNode::getClass(topNode, m_trk_map_node_name); + + if (!check_trackmap || check_trackmap->empty()) { - std::cout << "EventHeader or GL1 packet not found" << std::endl; - } - m_this_event_bco = -1; - m_last_event_bco = -1; - m_prev_event_bco = -1; - m_prev_runNumber = -1; - m_prev_eventNumber = -1; + if (Verbosity() >= VERBOSITY_SOME) + { + std::cout << "KFParticle: Event skipped as there are no tracks" << std::endl; + } + return Fun4AllReturnCodes::EVENT_OK; } -// End BCO matching here + multiplicity = check_trackmap->size(); + if (!m_use_fake_pv) { if (m_use_mbd_vertex) @@ -241,6 +242,7 @@ int KFParticle_sPHENIX::process_event(PHCompositeNode *topNode) } } + createDecay(topNode, mother, vertex_kfparticle, daughters, intermediates, nPVs); if (!m_has_intermediates_sPHENIX) { diff --git a/offline/packages/decayfinder/DecayFinder.cc b/offline/packages/decayfinder/DecayFinder.cc index a3521df780..596cce9c63 100644 --- a/offline/packages/decayfinder/DecayFinder.cc +++ b/offline/packages/decayfinder/DecayFinder.cc @@ -93,14 +93,14 @@ int DecayFinder::Init(PHCompositeNode* topNode) int DecayFinder::process_event(PHCompositeNode* topNode) { - bool decayFound = findDecay(topNode); + int decayFound = findDecay(topNode); - if (decayFound && m_save_dst && Verbosity() >= VERBOSITY_MORE) + if (decayFound > 0 && m_save_dst && Verbosity() >= VERBOSITY_MORE) { printNode(topNode); } - if (m_triggerOnDecay && !decayFound) + if (m_triggerOnDecay && decayFound < 1) { if (Verbosity() >= VERBOSITY_MORE) { @@ -317,10 +317,10 @@ int DecayFinder::parseDecayDescriptor() * as decays wont enter the HepMC record * need a switch to go to Geant4 record */ -bool DecayFinder::findDecay(PHCompositeNode* topNode) +int DecayFinder::findDecay(PHCompositeNode* topNode) { bool decayWasFound = false; - bool reconstructableDecayWasFound = false; + int reconstructableDecayWasFound = 0; bool aTrackFailedPT = false; bool aTrackFailedETA = false; bool aMotherHasPhoton = false; @@ -384,6 +384,11 @@ bool DecayFinder::findDecay(PHCompositeNode* topNode) std::cout << "parent->pdg_id(): " << g4particle->get_pid() << std::endl; } + aTrackFailedPT = false; + aTrackFailedETA = false; + aMotherHasPhoton = false; + aMotherHasPi0 = false; + bool breakOut = false; correctMotherProducts.clear(); decayChain.clear(); @@ -416,7 +421,7 @@ bool DecayFinder::findDecay(PHCompositeNode* topNode) else { m_nCandReconstructable += 1; - reconstructableDecayWasFound = true; + ++reconstructableDecayWasFound; if (m_save_dst) { fillDecayNode(topNode, decayChain); @@ -450,7 +455,7 @@ bool DecayFinder::findDecay(PHCompositeNode* topNode) if (!m_genevt) { std::cout << "DecayFinder: Missing node PHHepMCGenEvent" << std::endl; - return false; + continue; } HepMC::GenEvent* theEvent = m_genevt->getEvent(); @@ -465,6 +470,11 @@ bool DecayFinder::findDecay(PHCompositeNode* topNode) std::cout << "parent->pdg_id(): " << (*p)->pdg_id() << std::endl; } + aTrackFailedPT = false; + aTrackFailedETA = false; + aMotherHasPhoton = false; + aMotherHasPi0 = false; + bool breakOut = false; correctMotherProducts.clear(); decayChain.clear(); @@ -505,7 +515,7 @@ bool DecayFinder::findDecay(PHCompositeNode* topNode) else { m_nCandReconstructable += 1; - reconstructableDecayWasFound = true; + ++reconstructableDecayWasFound; if (m_save_dst) { fillDecayNode(topNode, decayChain); diff --git a/offline/packages/decayfinder/DecayFinder.h b/offline/packages/decayfinder/DecayFinder.h index bd893546a1..3605733e9f 100644 --- a/offline/packages/decayfinder/DecayFinder.h +++ b/offline/packages/decayfinder/DecayFinder.h @@ -42,7 +42,7 @@ class DecayFinder : public SubsysReco int parseDecayDescriptor(); - bool findDecay(PHCompositeNode *topNode); + int findDecay(PHCompositeNode *topNode); bool findParticle(const std::string &particle); diff --git a/offline/packages/trackreco/PHSimpleVertexFinder.cc b/offline/packages/trackreco/PHSimpleVertexFinder.cc index 5f9340d8f6..954e790ba6 100644 --- a/offline/packages/trackreco/PHSimpleVertexFinder.cc +++ b/offline/packages/trackreco/PHSimpleVertexFinder.cc @@ -515,34 +515,13 @@ void PHSimpleVertexFinder::checkDCAs(SvtxTrackMap *track_map) { continue; } - if (_require_mvtx) + if (_require_mvtx && !passClusterRequirement(tr1, "MVTX")) { - unsigned int nmvtx = 0; - TrackSeed *siliconseed = tr1->get_silicon_seed(); - if (!siliconseed) - { - continue; - } - - for (auto clusit = siliconseed->begin_cluster_keys(); clusit != siliconseed->end_cluster_keys(); ++clusit) - { - if (TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId) - { - nmvtx++; - } - if (nmvtx >= _nmvtx_required) - { - break; - } - } - if (nmvtx < _nmvtx_required) - { - continue; - } - if (Verbosity() > 3) - { - std::cout << " tr1 id " << id1 << " has nmvtx at least " << nmvtx << std::endl; - } + continue; + } + if (_require_intt && !passClusterRequirement(tr1, "INTT")) + { + continue; } // look for close DCA matches with all other such tracks @@ -554,34 +533,13 @@ void PHSimpleVertexFinder::checkDCAs(SvtxTrackMap *track_map) { continue; } - if (_require_mvtx) + if (_require_mvtx && !passClusterRequirement(tr2, "MVTX")) { - unsigned int nmvtx = 0; - TrackSeed *siliconseed = tr2->get_silicon_seed(); - if (!siliconseed) - { - continue; - } - - for (auto clusit = siliconseed->begin_cluster_keys(); clusit != siliconseed->end_cluster_keys(); ++clusit) - { - if (TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId) - { - nmvtx++; - } - if (nmvtx >= _nmvtx_required) - { - break; - } - } - if (nmvtx < _nmvtx_required) - { - continue; - } - if (Verbosity() > 3) - { - std::cout << " tr2 id " << id2 << " has nmvtx at least " << nmvtx << std::endl; - } + continue; + } + if (_require_intt && !passClusterRequirement(tr2, "INTT")) + { + continue; } // find DCA of these two tracks @@ -616,13 +574,20 @@ void PHSimpleVertexFinder::checkDCAsZF(SvtxTrackMap *track_map) // tr1->identify(); TrackSeed *siliconseed = tr1->get_silicon_seed(); - if (_require_mvtx) - { - if (!siliconseed) - { - continue; - } - } + const bool needs_mvtx_seed = _require_mvtx && _nmvtx_required > 0; + const bool needs_intt_seed = _require_intt && _nintt_required > 0; + if ((needs_mvtx_seed || needs_intt_seed) && !siliconseed) + { + continue; + } + if (_require_mvtx && !passClusterRequirement(tr1, "MVTX")) + { + continue; + } + if (_require_intt && !passClusterRequirement(tr1, "INTT")) + { + continue; + } TrackSeed *tpcseed = tr1->get_tpc_seed(); std::vector global_vec; @@ -797,34 +762,13 @@ void PHSimpleVertexFinder::checkDCAs() { continue; } - if (_require_mvtx) + if (_require_mvtx && !passClusterRequirement(tr1, "MVTX")) { - unsigned int nmvtx = 0; - TrackSeed *siliconseed = tr1->get_silicon_seed(); - if (!siliconseed) - { - continue; - } - - for (auto clusit = siliconseed->begin_cluster_keys(); clusit != siliconseed->end_cluster_keys(); ++clusit) - { - if (TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId) - { - nmvtx++; - } - if (nmvtx >= _nmvtx_required) - { - break; - } - } - if (nmvtx < _nmvtx_required) - { - continue; - } - if (Verbosity() > 3) - { - std::cout << " tr1 id " << id1 << " has nmvtx at least " << nmvtx << std::endl; - } + continue; + } + if (_require_intt && !passClusterRequirement(tr1, "INTT")) + { + continue; } // look for close DCA matches with all other such tracks @@ -836,36 +780,14 @@ void PHSimpleVertexFinder::checkDCAs() { continue; } - if (_require_mvtx) + if (_require_mvtx && !passClusterRequirement(tr2, "MVTX")) { - unsigned int nmvtx = 0; - TrackSeed *siliconseed = tr2->get_silicon_seed(); - if (!siliconseed) - { - continue; - } - - for (auto clusit = siliconseed->begin_cluster_keys(); clusit != siliconseed->end_cluster_keys(); ++clusit) - { - if (TrkrDefs::getTrkrId(*clusit) == TrkrDefs::mvtxId) - { - nmvtx++; - } - if (nmvtx >= _nmvtx_required) - { - break; - } - } - if (nmvtx < _nmvtx_required) - { - continue; - } - if (Verbosity() > 3) - { - std::cout << " tr2 id " << id2 << " has nmvtx at least " << nmvtx << std::endl; - } + continue; + } + if (_require_intt && !passClusterRequirement(tr2, "INTT")) + { + continue; } - // find DCA of these two tracks if (Verbosity() > 3) { @@ -1324,3 +1246,53 @@ double PHSimpleVertexFinder::getAverage(std::vector &v) return avge; } + +bool PHSimpleVertexFinder::passClusterRequirement(SvtxTrack *track, const std::string &type) +{ + bool pass = false; + + std::vector acceptable_types = {"MVTX", "INTT"}; + bool accept_this_type = std::find(acceptable_types.begin(), acceptable_types.end(), type) != acceptable_types.end(); + + if (!accept_this_type) + { + if (Verbosity() > 3) + { + std::cout << "type " << type << " was not recognised" << std::endl; + } + return false; + } + + unsigned int _nclus_required = type == "MVTX" ? _nmvtx_required : _nintt_required; + if (_nclus_required == 0) + { + return true; + } + + TrackSeed *siliconseed = track->get_silicon_seed(); + if (!siliconseed) + { + return false; + } + + unsigned int nclus = 0; + for (auto clusit = siliconseed->begin_cluster_keys(); clusit != siliconseed->end_cluster_keys(); ++clusit) + { + uint8_t trkrId = type == "MVTX" ? TrkrDefs::mvtxId : TrkrDefs::inttId; + if (TrkrDefs::getTrkrId(*clusit) == trkrId) + { + nclus++; + } + if (nclus >= _nclus_required) + { + pass = true; + } + } + + if (Verbosity() > 3) + { + std::cout << " track id " << track->get_id() << " has " << nclus << " clusters for " << type << std::endl; + } + + return pass; +} diff --git a/offline/packages/trackreco/PHSimpleVertexFinder.h b/offline/packages/trackreco/PHSimpleVertexFinder.h index de5cff940b..6520174937 100644 --- a/offline/packages/trackreco/PHSimpleVertexFinder.h +++ b/offline/packages/trackreco/PHSimpleVertexFinder.h @@ -46,16 +46,18 @@ class PHSimpleVertexFinder : public SubsysReco void setBeamSpotCutY(const double cutlo, const double cuthi) { _beamline_y_cut_lo = cutlo; _beamline_y_cut_hi = cuthi; } void setDcaCut(const double cut) { _base_dcacut = cut; } void setTrackQualityCut(double cut) { _qual_cut = cut; } - void setRequireMVTX(bool set) { _require_mvtx = set; } + void setRequireMVTX(bool set = true) { _require_mvtx = set; } void setNmvtxRequired(unsigned int n) { _nmvtx_required = n; } + void setRequireINTT(bool set = true) { _require_intt = set; } + void setNinttRequired(unsigned int n) { _nintt_required = n; } void setTrackPtCut(const double cut) { _track_pt_cut = cut; } // void setUseTrackCovariance(bool set) {_use_track_covariance = set;} void setOutlierPairCut(const double cut) { _outlier_cut = cut; } void setTrackMapName(const std::string &name) { _track_map_name = name; } void setVertexMapName(const std::string &name) { _vertex_map_name = name; } - void zeroField(const bool flag) { _zero_field = flag; } + void zeroField(const bool flag = true) { _zero_field = flag; } void setTrkrClusterContainerName(const std::string &name){ m_clusterContainerName = name; } - void set_pp_mode(bool mode) { _pp_mode = mode; } + void set_pp_mode(bool mode = true) { _pp_mode = mode; } private: int GetNodes(PHCompositeNode *topNode); @@ -75,6 +77,7 @@ class PHSimpleVertexFinder : public SubsysReco void removeOutlierTrackPairs(); double getMedian(std::vector &v); double getAverage(std::vector &v); + bool passClusterRequirement(SvtxTrack *track, const std::string &type = "MVTX"); SvtxTrackMap *_track_map{nullptr}; TrkrClusterContainer* _cluster_map{nullptr}; @@ -91,7 +94,9 @@ class PHSimpleVertexFinder : public SubsysReco double _beamline_y_cut_hi = 0.2; double _qual_cut = 10.0; bool _require_mvtx = true; + bool _require_intt = false; unsigned int _nmvtx_required = 2; + unsigned int _nintt_required = 1; double _track_pt_cut = 0.0; double _outlier_cut = 0.015;