diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.cc b/offline/packages/TrackingDiagnostics/TrackResiduals.cc index 1a5a69fae0..4b613c5dd5 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.cc +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.cc @@ -46,6 +46,8 @@ #include #include +#include + #include #include #include @@ -57,6 +59,7 @@ #include #include #include +#include #include @@ -126,6 +129,7 @@ int TrackResiduals::InitRun(PHCompositeNode* topNode) void TrackResiduals::clearClusterStateVectors() { m_cluskeys.clear(); + m_clussize.clear(); m_clusphisize.clear(); m_cluszsize.clear(); m_idealsurfcenterx.clear(); @@ -183,7 +187,23 @@ void TrackResiduals::clearClusterStateVectors() m_statelzlocderivqop.clear(); m_clusedge.clear(); + m_clussledge.clear(); + m_clussredge.clear(); + m_clustledge.clear(); + m_clustredge.clear(); + m_clusdledge.clear(); + m_clusdredge.clear(); + m_clushledge.clear(); + m_clushredge.clear(); + m_clusslmix.clear(); + m_clussrmix.clear(); + m_clustlmix.clear(); + m_clustrmix.clear(); m_clusoverlap.clear(); + m_clusPadCen.clear(); + m_clusTBinCen.clear(); + m_clusPadMax.clear(); + m_clusTBinMax.clear(); m_cluslx.clear(); m_cluslz.clear(); m_cluselx.clear(); @@ -197,7 +217,14 @@ void TrackResiduals::clearClusterStateVectors() m_clusgzunmoved.clear(); m_clusAdc.clear(); m_clusMaxAdc.clear(); + m_clusCenAdc.clear(); m_cluslayer.clear(); + m_clusphibinlo.clear(); + m_clusphibinhi.clear(); + m_clustbinlo.clear(); + m_clustbinhi.clear(); + m_cluspadphase.clear(); + m_clustbinphase.clear(); m_statelx.clear(); m_statelz.clear(); @@ -301,6 +328,21 @@ int TrackResiduals::process_event(PHCompositeNode* topNode) } } + EventHeader* eventheader = findNode::getClass(topNode, "EventHeader"); + + if(eventheader) + { + m_evt_id = eventheader->get_EvtSequence(); + } + else + { + m_evt_id = -1; + } + + auto *rcs = recoConsts::instance(); + m_runnumber = rcs->get_IntFlag("RUNNUMBER"); + m_segment = rcs->get_IntFlag("RUNSEGMENT"); + m_ntpcclus = 0; if (Verbosity() > 1) { @@ -662,11 +704,37 @@ void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters, m_scluseta = acos(glob.z() / std::sqrt(square(glob.x()) + square(glob.y()) + square(glob.z()))); m_adc = cluster->getAdc(); m_clusmaxadc = cluster->getMaxAdc(); + m_cluscenadc = cluster->getCenAdc(); + m_padcen = cluster->getPadCen(); + m_tbincen = cluster->getTBinCen(); + m_padmax = cluster->getPadMax(); + m_tbinmax = cluster->getTBinMax(); m_scluslx = cluster->getLocalX(); m_scluslz = cluster->getLocalY(); + m_phibinlo = cluster->getPhiBinLo(); + m_phibinhi = cluster->getPhiBinHi(); + m_tbinlo = cluster->getTBinLo(); + m_tbinhi = cluster->getTBinHi(); + m_padphase = cluster->getPadPhase(); + m_tbinphase = cluster->getTBinPhase(); auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, m_sclusgr, key); + m_size = cluster->getRSize(); m_phisize = cluster->getPhiSize(); m_zsize = cluster->getZSize(); + m_overlap = cluster->getOverlap(); + m_nedge = cluster->getEdge(); + m_sledge = cluster->getSLEdge(); + m_sredge = cluster->getSREdge(); + m_tledge = cluster->getTLEdge(); + m_tredge = cluster->getTREdge(); + m_dledge = cluster->getDLEdge(); + m_dredge = cluster->getDREdge(); + m_hledge = cluster->getHLEdge(); + m_hredge = cluster->getHREdge(); + m_slmix = cluster->getSLMix(); + m_srmix = cluster->getSRMix(); + m_tlmix = cluster->getTLMix(); + m_trmix = cluster->getTRMix(); m_scluselx = std::sqrt(para_errors.first); m_scluselz = std::sqrt(para_errors.second); @@ -1087,7 +1155,19 @@ void TrackResiduals::fillClusterBranchesKF(TrkrDefs::cluskey ckey, SvtxTrack* tr //! have cluster and state, fill vectors m_clusedge.push_back(cluster->getEdge()); + m_clussledge.push_back(cluster->getSLEdge()); + m_clussredge.push_back(cluster->getSREdge()); + m_clustledge.push_back(cluster->getTLEdge()); + m_clustredge.push_back(cluster->getTREdge()); + m_clusdledge.push_back(cluster->getDLEdge()); + m_clusdredge.push_back(cluster->getDREdge()); + m_clushledge.push_back(cluster->getHLEdge()); + m_clushredge.push_back(cluster->getHREdge()); m_clusoverlap.push_back(cluster->getOverlap()); + m_clusslmix.push_back(cluster->getSLMix()); + m_clussrmix.push_back(cluster->getSRMix()); + m_clustlmix.push_back(cluster->getTLMix()); + m_clustrmix.push_back(cluster->getTRMix()); // get new local coords from moved cluster Surface surf = geometry->maps().getSurface(ckey, cluster); @@ -1160,9 +1240,21 @@ void TrackResiduals::fillClusterBranchesKF(TrkrDefs::cluskey ckey, SvtxTrack* tr m_clusgzunmoved.push_back(clusglob.z()); m_clusAdc.push_back(cluster->getAdc()); m_clusMaxAdc.push_back(cluster->getMaxAdc()); + m_clusCenAdc.push_back(cluster->getCenAdc()); + m_clusPadCen.push_back(cluster->getPadCen()); + m_clusTBinCen.push_back(cluster->getTBinCen()); + m_clusPadMax.push_back(cluster->getPadMax()); + m_clusTBinMax.push_back(cluster->getTBinMax()); m_cluslayer.push_back(TrkrDefs::getLayer(ckey)); + m_clussize.push_back(cluster->getRSize()); m_clusphisize.push_back(cluster->getPhiSize()); m_cluszsize.push_back(cluster->getZSize()); + m_clusphibinlo.push_back(cluster->getPhiBinLo()); + m_clusphibinhi.push_back(cluster->getPhiBinHi()); + m_clustbinlo.push_back(cluster->getTBinLo()); + m_clustbinhi.push_back(cluster->getTBinHi()); + m_cluspadphase.push_back(cluster->getPadPhase()); + m_clustbinphase.push_back(cluster->getTBinPhase()); auto misaligncenter = surf->center(geometry->geometry().getGeoContext()); auto misalignnorm = -1 * surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1)); @@ -1416,6 +1508,18 @@ void TrackResiduals::fillClusterBranchesSeeds(TrkrDefs::cluskey ckey, // SvtxTr //! have cluster and state, fill vectors m_clusedge.push_back(cluster->getEdge()); + m_clussledge.push_back(cluster->getSLEdge()); + m_clussredge.push_back(cluster->getSREdge()); + m_clustledge.push_back(cluster->getTLEdge()); + m_clustredge.push_back(cluster->getTREdge()); + m_clusdledge.push_back(cluster->getDLEdge()); + m_clusdredge.push_back(cluster->getDREdge()); + m_clushledge.push_back(cluster->getHLEdge()); + m_clushredge.push_back(cluster->getHREdge()); + m_clusslmix.push_back(cluster->getSLMix()); + m_clussrmix.push_back(cluster->getSRMix()); + m_clustlmix.push_back(cluster->getTLMix()); + m_clustrmix.push_back(cluster->getTRMix()); m_clusoverlap.push_back(cluster->getOverlap()); // This is the nominal position of the cluster in local coords, completely uncorrected - is that what we want? @@ -1438,9 +1542,21 @@ void TrackResiduals::fillClusterBranchesSeeds(TrkrDefs::cluskey ckey, // SvtxTr m_clusgzunmoved.push_back(clusglob.z()); m_clusAdc.push_back(cluster->getAdc()); m_clusMaxAdc.push_back(cluster->getMaxAdc()); + m_clusCenAdc.push_back(cluster->getCenAdc()); + m_clusPadCen.push_back(cluster->getPadCen()); + m_clusTBinCen.push_back(cluster->getTBinCen()); + m_clusPadMax.push_back(cluster->getPadMax()); + m_clusTBinMax.push_back(cluster->getTBinMax()); m_cluslayer.push_back(TrkrDefs::getLayer(ckey)); + m_clussize.push_back(cluster->getRSize()); m_clusphisize.push_back(cluster->getPhiSize()); m_cluszsize.push_back(cluster->getZSize()); + m_clusphibinlo.push_back(cluster->getPhiBinLo()); + m_clusphibinhi.push_back(cluster->getPhiBinHi()); + m_clustbinlo.push_back(cluster->getTBinLo()); + m_clustbinhi.push_back(cluster->getTBinHi()); + m_cluspadphase.push_back(cluster->getPadPhase()); + m_clustbinphase.push_back(cluster->getTBinPhase()); if (Verbosity() > 1) { @@ -1607,6 +1723,7 @@ void TrackResiduals::createBranches() m_eventtree->Branch("run", &m_runnumber, "m_runnumber/I"); m_eventtree->Branch("segment", &m_segment, "m_segment/I"); m_eventtree->Branch("event", &m_event, "m_event/I"); + m_eventtree->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_eventtree->Branch("gl1bco", &m_bco, "m_bco/I"); m_eventtree->Branch("nmvtx", &m_nmvtx_all, "m_nmvtx_all/I"); m_eventtree->Branch("nintt", &m_nintt_all, "m_nintt_all/I"); @@ -1627,6 +1744,7 @@ void TrackResiduals::createBranches() m_failedfits->Branch("segment", &m_segment, "m_segment/I"); m_failedfits->Branch("trackid", &m_trackid, "m_trackid/I"); m_failedfits->Branch("event", &m_event, "m_event/I"); + m_failedfits->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_failedfits->Branch("silseedx", &m_silseedx, "m_silseedx/F"); m_failedfits->Branch("silseedy", &m_silseedy, "m_silseedy/F"); m_failedfits->Branch("silseedz", &m_silseedz, "m_silseedz/F"); @@ -1653,6 +1771,7 @@ void TrackResiduals::createBranches() m_vertextree->Branch("run", &m_runnumber, "m_runnumber/I"); m_vertextree->Branch("segment", &m_segment, "m_segment/I"); m_vertextree->Branch("event", &m_event, "m_event/I"); + m_vertextree->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_vertextree->Branch("firedTriggers", &m_firedTriggers); m_vertextree->Branch("gl1BunchCrossing", &m_gl1BunchCrossing, "m_gl1BunchCrossing/l"); m_vertextree->Branch("gl1bco", &m_bco, "m_bco/l"); @@ -1675,6 +1794,7 @@ void TrackResiduals::createBranches() m_hittree->Branch("run", &m_runnumber, "m_runnumber/I"); m_hittree->Branch("segment", &m_segment, "m_segment/I"); m_hittree->Branch("event", &m_event, "m_event/I"); + m_hittree->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_hittree->Branch("gl1bco", &m_bco, "m_bco/l"); m_hittree->Branch("hitsetkey", &m_hitsetkey, "m_hitsetkey/i"); m_hittree->Branch("gx", &m_hitgx, "m_hitgx/F"); @@ -1704,6 +1824,7 @@ void TrackResiduals::createBranches() m_clustree->Branch("run", &m_runnumber, "m_runnumber/I"); m_clustree->Branch("segment", &m_segment, "m_segment/I"); m_clustree->Branch("event", &m_event, "m_event/I"); + m_clustree->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_clustree->Branch("gl1bco", &m_bco, "m_bco/l"); m_clustree->Branch("lx", &m_scluslx, "m_scluslx/F"); m_clustree->Branch("lz", &m_scluslz, "m_scluslz/F"); @@ -1713,11 +1834,37 @@ void TrackResiduals::createBranches() m_clustree->Branch("phi", &m_sclusphi, "m_sclusphi/F"); m_clustree->Branch("eta", &m_scluseta, "m_scluseta/F"); m_clustree->Branch("adc", &m_adc, "m_adc/F"); + m_clustree->Branch("size", &m_size, "m_size/I"); m_clustree->Branch("phisize", &m_phisize, "m_phisize/I"); m_clustree->Branch("zsize", &m_zsize, "m_zsize/I"); + m_clustree->Branch("phibinlo", &m_phibinlo, "m_phibinlo/F"); + m_clustree->Branch("phibinhi", &m_phibinhi, "m_phibinhi/F"); + m_clustree->Branch("tbinlo", &m_tbinlo, "m_tbinlo/F"); + m_clustree->Branch("tbinhi", &m_tbinhi, "m_tbinhi/F"); + m_clustree->Branch("padphase", &m_padphase, "m_padphase/F"); + m_clustree->Branch("tbinphase", &m_tbinphase, "m_tbinphase/F"); + m_clustree->Branch("overlap", &m_overlap, "m_overlap/B"); + m_clustree->Branch("nedge", &m_nedge, "m_nedge/B"); + m_clustree->Branch("sledge", &m_sledge, "m_sledge/B"); + m_clustree->Branch("sredge", &m_sredge, "m_sredge/B"); + m_clustree->Branch("tledge", &m_tledge, "m_tledge/B"); + m_clustree->Branch("tredge", &m_tredge, "m_tredge/B"); + m_clustree->Branch("dledge", &m_dledge, "m_dledge/B"); + m_clustree->Branch("dredge", &m_dredge, "m_dredge/B"); + m_clustree->Branch("hledge", &m_hledge, "m_hledge/B"); + m_clustree->Branch("hredge", &m_hredge, "m_hredge/B"); + m_clustree->Branch("slmix", &m_slmix, "m_slmix/B"); + m_clustree->Branch("srmix", &m_srmix, "m_srmix/B"); + m_clustree->Branch("tlmix", &m_tlmix, "m_tlmix/B"); + m_clustree->Branch("trmix", &m_trmix, "m_trmix/B"); m_clustree->Branch("erphi", &m_scluselx, "m_scluselx/F"); m_clustree->Branch("ez", &m_scluselz, "m_scluselz/F"); m_clustree->Branch("maxadc", &m_clusmaxadc, "m_clusmaxadc/F"); + m_clustree->Branch("cenadc", &m_cluscenadc, "m_cluscenadc/F"); + m_clustree->Branch("padcen", &m_padcen, "m_padcen/F"); + m_clustree->Branch("tbincen", &m_tbincen, "m_tbincen/F"); + m_clustree->Branch("padmax", &m_padmax, "m_padmax/F"); + m_clustree->Branch("tbinmax", &m_tbinmax, "m_tbinmax/F"); m_clustree->Branch("sector", &m_clussector, "m_clussector/I"); m_clustree->Branch("side", &m_side, "m_side/I"); m_clustree->Branch("stave", &m_staveid, "m_staveid/I"); @@ -1734,6 +1881,7 @@ void TrackResiduals::createBranches() m_tree->Branch("run", &m_runnumber, "m_runnumber/I"); m_tree->Branch("segment", &m_segment, "m_segment/I"); m_tree->Branch("event", &m_event, "m_event/I"); + m_tree->Branch("evt_id", &m_evt_id, "m_evt_id/I"); m_tree->Branch("mbdcharge",&m_totalmbd, "m_totalmbd/F"); m_tree->Branch("mbdzvtx", &m_mbdvtxz, "m_mbdvtxz/F"); m_tree->Branch("firedTriggers", &m_firedTriggers); @@ -1815,7 +1963,27 @@ void TrackResiduals::createBranches() m_tree->Branch("clusside", &m_clside); m_tree->Branch("cluskeys", &m_cluskeys); m_tree->Branch("clusedge", &m_clusedge); + m_tree->Branch("clussledge", &m_clussledge); + m_tree->Branch("clussredge", &m_clussredge); + m_tree->Branch("clustledge", &m_clustledge); + m_tree->Branch("clustredge", &m_clustredge); + m_tree->Branch("clusdledge", &m_clusdledge); + m_tree->Branch("clusdredge", &m_clusdredge); + m_tree->Branch("clushledge", &m_clushledge); + m_tree->Branch("clushredge", &m_clushredge); + m_tree->Branch("clusslmix", &m_clusslmix); + m_tree->Branch("clussrmix", &m_clussrmix); + m_tree->Branch("clustlmix", &m_clustlmix); + m_tree->Branch("clustrmix", &m_clustrmix); m_tree->Branch("clusoverlap", &m_clusoverlap); + m_tree->Branch("clusphibinlo", &m_clusphibinlo); + m_tree->Branch("clusphibinhi", &m_clusphibinhi); + m_tree->Branch("clustbinlo", &m_clustbinlo); + m_tree->Branch("clustbinhi", &m_clustbinhi); + m_tree->Branch("clusPadCen", &m_clusPadCen); + m_tree->Branch("clusTBinCen", &m_clusTBinCen); + m_tree->Branch("clusPadMax", &m_clusPadMax); + m_tree->Branch("clusTBinMax", &m_clusTBinMax); m_tree->Branch("cluslx", &m_cluslx); m_tree->Branch("cluslz", &m_cluslz); m_tree->Branch("cluselx", &m_cluselx); @@ -1824,6 +1992,8 @@ void TrackResiduals::createBranches() m_tree->Branch("clusgy", &m_clusgy); m_tree->Branch("clusgz", &m_clusgz); m_tree->Branch("clusgr", &m_clusgr); + m_tree->Branch("cluspadphase", &m_cluspadphase); + m_tree->Branch("clustbinphase", &m_clustbinphase); if (m_doAlignment) { m_tree->Branch("clusgxunmoved", &m_clusgxunmoved); @@ -1832,6 +2002,8 @@ void TrackResiduals::createBranches() } m_tree->Branch("clusAdc", &m_clusAdc); m_tree->Branch("clusMaxAdc", &m_clusMaxAdc); + m_tree->Branch("clusCenAdc", &m_clusCenAdc); + m_tree->Branch("clussize", &m_clussize); m_tree->Branch("clusphisize", &m_clusphisize); m_tree->Branch("cluszsize", &m_cluszsize); @@ -2235,6 +2407,7 @@ void TrackResiduals::fillEventTree(PHCompositeNode* topNode) if (Verbosity() > 1) { std::cout << " m_event:" << m_event << std::endl; + std::cout << " m_evt_id:" << m_evt_id << std::endl; std::cout << " m_ntpc_clus0:" << m_ntpc_clus0 << std::endl; std::cout << " m_ntpc_clus1: " << m_ntpc_clus1 << std::endl; std::cout << " m_nmvtx_all:" << m_nmvtx_all << std::endl; diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.h b/offline/packages/TrackingDiagnostics/TrackResiduals.h index 7e789faba5..2fb8b6f648 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.h +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.h @@ -128,6 +128,7 @@ class TrackResiduals : public SubsysReco bool m_doMicromegasOnly = false; int m_event = 0; + int m_evt_id = -1; int m_segment = std::numeric_limits::quiet_NaN(); int m_runnumber = std::numeric_limits::quiet_NaN(); int m_ntpcclus = std::numeric_limits::quiet_NaN(); @@ -246,8 +247,34 @@ class TrackResiduals : public SubsysReco float m_scluseta = std::numeric_limits::quiet_NaN(); float m_adc = std::numeric_limits::quiet_NaN(); float m_clusmaxadc = std::numeric_limits::quiet_NaN(); + float m_cluscenadc = std::numeric_limits::quiet_NaN(); + int m_size = std::numeric_limits::quiet_NaN(); int m_phisize = std::numeric_limits::quiet_NaN(); int m_zsize = std::numeric_limits::quiet_NaN(); + char m_overlap = std::numeric_limits::max(); + char m_nedge = std::numeric_limits::max(); + char m_sledge = std::numeric_limits::max(); + char m_sredge = std::numeric_limits::max(); + char m_tledge = std::numeric_limits::max(); + char m_tredge = std::numeric_limits::max(); + char m_dledge = std::numeric_limits::max(); + char m_dredge = std::numeric_limits::max(); + char m_hledge = std::numeric_limits::max(); + char m_hredge = std::numeric_limits::max(); + char m_slmix = std::numeric_limits::max(); + char m_srmix = std::numeric_limits::max(); + char m_tlmix = std::numeric_limits::max(); + char m_trmix = std::numeric_limits::max(); + float m_phibinlo = std::numeric_limits::quiet_NaN(); + float m_phibinhi = std::numeric_limits::quiet_NaN(); + float m_tbinlo = std::numeric_limits::quiet_NaN(); + float m_tbinhi = std::numeric_limits::quiet_NaN(); + float m_padphase = std::numeric_limits::quiet_NaN(); + float m_tbinphase = std::numeric_limits::quiet_NaN(); + float m_padcen = std::numeric_limits::quiet_NaN(); + float m_tbincen = std::numeric_limits::quiet_NaN(); + float m_padmax = std::numeric_limits::quiet_NaN(); + float m_tbinmax = std::numeric_limits::quiet_NaN(); float m_scluslx = std::numeric_limits::quiet_NaN(); float m_scluslz = std::numeric_limits::quiet_NaN(); float m_sclusgx = std::numeric_limits::quiet_NaN(); @@ -270,6 +297,11 @@ class TrackResiduals : public SubsysReco //! clusters on track information std::vector m_clusAdc; std::vector m_clusMaxAdc; + std::vector m_clusCenAdc; + std::vector m_clusPadCen; + std::vector m_clusTBinCen; + std::vector m_clusPadMax; + std::vector m_clusTBinMax; std::vector m_cluslx; std::vector m_cluslz; std::vector m_cluselx; @@ -288,10 +320,29 @@ class TrackResiduals : public SubsysReco std::vector m_clsector; std::vector m_clside; std::vector m_cluslayer; + std::vector m_clussize; std::vector m_clusphisize; std::vector m_cluszsize; - std::vector m_clusedge; std::vector m_clusoverlap; + std::vector m_clusedge; + std::vector m_clussledge; + std::vector m_clussredge; + std::vector m_clustledge; + std::vector m_clustredge; + std::vector m_clusdledge; + std::vector m_clusdredge; + std::vector m_clushledge; + std::vector m_clushredge; + std::vector m_clusslmix; + std::vector m_clussrmix; + std::vector m_clustlmix; + std::vector m_clustrmix; + std::vector m_clusphibinlo; + std::vector m_clusphibinhi; + std::vector m_clustbinlo; + std::vector m_clustbinhi; + std::vector m_cluspadphase; + std::vector m_clustbinphase; std::vector m_cluskeys; std::vector m_idealsurfcenterx; std::vector m_idealsurfcentery; diff --git a/offline/packages/TrackingDiagnostics/TrkrNtuplizer.cc b/offline/packages/TrackingDiagnostics/TrkrNtuplizer.cc index d5fddf6538..1c90c4688c 100644 --- a/offline/packages/TrackingDiagnostics/TrkrNtuplizer.cc +++ b/offline/packages/TrackingDiagnostics/TrkrNtuplizer.cc @@ -289,6 +289,11 @@ enum n_cluster // NOLINT(readability-enum-initial-value, performance-enum-size) nclue, ncluadc, nclumaxadc, + nclucenadc, + nclupadcen, + nclutbincen, + nclupadmax, + nclutbinmax, ncluthick, ncluafac, nclubfac, @@ -301,7 +306,25 @@ enum n_cluster // NOLINT(readability-enum-initial-value, performance-enum-size) ncluzsize, nclupedge, ncluredge, + nclusledge, + nclusredge, + nclutledge, + nclutredge, + ncludledge, + ncludredge, + ncluhledge, + ncluhredge, + ncluslmix, + nclusrmix, + nclutlmix, + nclutrmix, ncluovlp, + ncluphibinlo, + ncluphibinhi, + nclutbinlo, + nclutbinhi, + nclupadphase, + nclutbinphase, nclutrackID, ncluniter, clusize = ncluniter + 1 @@ -337,7 +360,7 @@ int TrkrNtuplizer::Init(PHCompositeNode* /*unused*/) std::string str_vertex = {"vertexID:vx:vy:vz:ntracks:chi2:ndof"}; std::string str_event = {"event:seed:run:seg:job"}; std::string str_hit = {"hitID:e:adc:layer:phielem:zelem:cellID:ecell:phibin:zbin:tbin:phi:r:x:y:z"}; - std::string str_cluster = {"locx:locy:x:y:z:r:phi:eta:theta:phibin:tbin:fee:chan:sampa:ex:ey:ez:ephi:pez:pephi:e:adc:maxadc:thick:afac:bfac:dcal:layer:phielem:zelem:size:phisize:zsize:pedge:redge:ovlp:trackID:niter"}; + std::string str_cluster = {"locx:locy:x:y:z:r:phi:eta:theta:phibin:tbin:fee:chan:sampa:ex:ey:ez:ephi:pez:pephi:e:adc:maxadc:cenadc:padcen:tbincen:padmax:tbinmax:thick:afac:bfac:dcal:layer:phielem:zelem:size:phisize:zsize:pedge:redge:sledge:sredge:tledge:tredge:dledge:dredge:hledge:hredge:slmix:srmix:tlmix:trmix:ovlp:phibinlo:phibinhi:tbinlo:tbinhi:padphase:tbinphase:trackID:niter"}; std::string str_seed = {"seedID:siter:spt:sptot:seta:sphi:syxint:srzint:sxyslope:srzslope:sX0:sY0:sdZ0:sR0:scharge:sdedx:spidedx:skdedx:sprdedx:sn1pix:snsil:sntpc:snhits"}; std::string str_residual = {"alpha:beta:resphio:resphi:resz"}; std::string str_track = {"trackID:crossing:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:charge:quality:chisq:ndf:nhits:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:dedx:pidedx:kdedx:prdedx:vertexID:vx:vy:vz:dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:hlxpt:hlxeta:hlxphi:hlxX0:hlxY0:hlxZ0:hlxcharge"}; @@ -1386,7 +1409,7 @@ void TrkrNtuplizer::fillOutputNtuples(PHCompositeNode* topNode) } //----------------------- - // fill the Vertex NTuple and fixed NaN placeholders + // fill the Vertex NTuple and fixed NaN placeholders //----------------------- bool doit = true; if (_ntp_vertex && doit) @@ -1409,8 +1432,7 @@ void TrkrNtuplizer::fillOutputNtuples(PHCompositeNode* topNode) for (auto & iter : *vertexmap) { SvtxVertex* vertex = iter.second; - if (!vertex) { continue; -} + if (!vertex) { continue; } float fx_vertex[n_vertex::vtxsize]; for (float& i : fx_vertex) @@ -1453,8 +1475,6 @@ void TrkrNtuplizer::fillOutputNtuples(PHCompositeNode* topNode) _timer->stop(); std::cout << "vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl; } - - //-------------------- // fill the Hit NTuple //-------------------- @@ -2284,6 +2304,11 @@ void TrkrNtuplizer::FillCluster(float fXcluster[n_cluster::clusize], TrkrDefs::c fXcluster[n_cluster::nclue] = cluster->getAdc(); fXcluster[n_cluster::ncluadc] = cluster->getAdc(); fXcluster[n_cluster::nclumaxadc] = cluster->getMaxAdc(); + fXcluster[n_cluster::nclucenadc] = cluster->getCenAdc(); + fXcluster[n_cluster::nclupadcen] = cluster->getPadCen(); + fXcluster[n_cluster::nclutbincen] = cluster->getTBinCen(); + fXcluster[n_cluster::nclupadmax] = cluster->getPadMax(); + fXcluster[n_cluster::nclutbinmax] = cluster->getTBinMax(); fXcluster[n_cluster::nclulayer] = layer_local; if (layer_local < 3) @@ -2311,7 +2336,7 @@ void TrkrNtuplizer::FillCluster(float fXcluster[n_cluster::clusize], TrkrDefs::c } } */ - fXcluster[n_cluster::nclusize] = cluster->getSize(); + fXcluster[n_cluster::nclusize] = cluster->getRSize(); fXcluster[n_cluster::ncluphisize] = cluster->getPhiSize(); fXcluster[n_cluster::ncluzsize] = cluster->getZSize(); fXcluster[n_cluster::nclupedge] = cluster->getEdge(); @@ -2321,8 +2346,25 @@ void TrkrNtuplizer::FillCluster(float fXcluster[n_cluster::clusize], TrkrDefs::c { fXcluster[n_cluster::ncluredge] = 1; } - - fXcluster[n_cluster::ncluovlp] = 3; // cluster->getOvlp(); + fXcluster[n_cluster::nclusledge] = cluster->getSLEdge(); + fXcluster[n_cluster::nclusredge] = cluster->getSREdge(); + fXcluster[n_cluster::nclutledge] = cluster->getTLEdge(); + fXcluster[n_cluster::nclutredge] = cluster->getTREdge(); + fXcluster[n_cluster::ncludledge] = cluster->getDLEdge(); + fXcluster[n_cluster::ncludredge] = cluster->getDREdge(); + fXcluster[n_cluster::ncluhledge] = cluster->getHLEdge(); + fXcluster[n_cluster::ncluhredge] = cluster->getHREdge(); + fXcluster[n_cluster::ncluslmix] = cluster->getSLMix(); + fXcluster[n_cluster::nclusrmix] = cluster->getSRMix(); + fXcluster[n_cluster::nclutlmix] = cluster->getTLMix(); + fXcluster[n_cluster::nclutrmix] = cluster->getTRMix(); + fXcluster[n_cluster::ncluovlp] = cluster->getOverlap(); + fXcluster[n_cluster::ncluphibinlo] = cluster->getPhiBinLo(); + fXcluster[n_cluster::ncluphibinhi] = cluster->getPhiBinHi(); + fXcluster[n_cluster::nclutbinlo] = cluster->getTBinLo(); + fXcluster[n_cluster::nclutbinhi] = cluster->getTBinHi(); + fXcluster[n_cluster::nclupadphase] = cluster->getPadPhase(); + fXcluster[n_cluster::nclutbinphase] = cluster->getTBinPhase(); fXcluster[n_cluster::nclutrackID] = std::numeric_limits::quiet_NaN(); fXcluster[n_cluster::ncluniter] = 0; diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index 687d5b58f9..32479b6e5f 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include // for hitkey, getLayer #include #include @@ -84,6 +85,37 @@ namespace unsigned short edge = 0; }; + // NOLINTBEGIN(misc-non-private-member-variables-in-classes) + struct ClusterCounters + { + int overlap = 0; + + int nedge = 0; // Total No. of Edges + + int sledge = 0; // Touching Left Sector Edge + int sredge = 0; // Touching Right Sector Edge + + int tledge = 0; // Touching Left Time Edge + int tredge = 0; // Touching Right Time Edge + + int dledge = 0; // Touching Left Dead Edge + int dredge = 0; // Touching Right Dead Edge + + int hledge = 0; // Touching Left Hot Edge + int hredge = 0; // Touching Right Hot Edge + + int slmix = 0; // Touching Cluster at Left in Phibin + int srmix = 0; // Touching Cluster at Right in Phibin + + int tlmix = 0; // Touching Cluster at Left in Timebin + int trmix = 0; // Touching Cluster at Right in Timebin + + void clear() + { + *this = ClusterCounters{}; + } + }; + // NOLINTEND(misc-non-private-member-variables-in-classes) using vec_dVerbose = std::vector>>; // Neural network parameters and modules @@ -129,6 +161,7 @@ namespace hitMaskTpcSet *hotMap = nullptr; bool maskDead = false; bool maskHot = false; + bool debug = false; std::vector association_vector; std::vector cluster_vector; @@ -176,13 +209,14 @@ namespace } } - void find_t_range(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, int &tdown, int &tup, int &touch, int &edge) + void find_t_range(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, int &tdown, int &tup, ClusterCounters &counts, bool &ttop_edge, bool &tbottom_edge) { const int FitRangeT = (int) my_data.maxHalfSizeT; const int NTBinsMax = (int) my_data.tbins; const int FixedWindow = (int) my_data.FixedWindow; tup = 0; tdown = 0; + if (FixedWindow != 0) { tup = FixedWindow; @@ -190,15 +224,16 @@ namespace if (tbin + tup >= NTBinsMax) { tup = NTBinsMax - tbin - 1; - edge++; + counts.nedge++; } if ((tbin - tdown) <= 0) { tdown = tbin; - edge++; + counts.nedge++; } return; } + for (int it = 0; it < FitRangeT; it++) { int ct = tbin + it; @@ -206,7 +241,12 @@ namespace if (ct <= 0 || ct >= NTBinsMax) { // tup = it; - edge++; + if (!ttop_edge) + { + counts.nedge++; + counts.tredge = 1; + ttop_edge = true; + } break; // truncate edge } @@ -216,7 +256,7 @@ namespace } if (adcval[phibin][ct] == USHRT_MAX) { - touch++; + counts.overlap++; break; } if (my_data.do_split) @@ -228,7 +268,7 @@ namespace adcval[phibin][ct + 2] + adcval[phibin][ct + 3]) { // rising again tup = it + 1; - touch++; + counts.overlap++; break; } } @@ -241,7 +281,12 @@ namespace if (ct <= 0 || ct >= NTBinsMax) { // tdown = it; - edge++; + if (!tbottom_edge) + { + counts.nedge++; + counts.tledge = 1; + tbottom_edge = true; + } break; // truncate edge } if (adcval[phibin][ct] <= 0) @@ -250,7 +295,7 @@ namespace } if (adcval[phibin][ct] == USHRT_MAX) { - touch++; + counts.overlap++; break; } if (my_data.do_split) @@ -261,7 +306,7 @@ namespace adcval[phibin][ct - 2] + adcval[phibin][ct - 3]) { // rising again tdown = it + 1; - touch++; + counts.overlap++; break; } } @@ -271,13 +316,14 @@ namespace return; } - void find_phi_range(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, int &phidown, int &phiup, int &touch, int &edge) + void find_phi_range(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, int &phidown, int &phiup, ClusterCounters &counts, bool &phitop_edge, bool &phibottom_edge) { int FitRangePHI = (int) my_data.maxHalfSizePhi; int NPhiBinsMax = (int) my_data.phibins; const int FixedWindow = (int) my_data.FixedWindow; phidown = 0; phiup = 0; + if (FixedWindow != 0) { phiup = FixedWindow; @@ -285,22 +331,28 @@ namespace if (phibin + phiup >= NPhiBinsMax) { phiup = NPhiBinsMax - phibin - 1; - edge++; + counts.nedge++; } if (phibin - phidown <= 0) { phidown = phibin; - edge++; + counts.nedge++; } return; } + for (int iphi = 0; iphi < FitRangePHI; iphi++) { int cphi = phibin + iphi; if (cphi < 0 || cphi >= NPhiBinsMax) { // phiup = iphi; - edge++; + if (!phitop_edge) + { + counts.nedge++; + counts.sredge = 1; + phitop_edge = true; + } break; // truncate edge } @@ -312,7 +364,7 @@ namespace } if (adcval[cphi][tbin] == USHRT_MAX) { - touch++; + counts.overlap++; break; } if (my_data.do_split) @@ -323,7 +375,7 @@ namespace adcval[cphi + 2][tbin] + adcval[cphi + 3][tbin]) { // rising again phiup = iphi + 1; - touch++; + counts.overlap++; break; } } @@ -337,7 +389,12 @@ namespace if (cphi < 0 || cphi >= NPhiBinsMax) { // phidown = iphi; - edge++; + if (!phibottom_edge) + { + counts.nedge++; + counts.sledge = 1; + phibottom_edge = true; + } break; // truncate edge } @@ -348,7 +405,7 @@ namespace } if (adcval[cphi][tbin] == USHRT_MAX) { - touch++; + counts.overlap++; break; } if (my_data.do_split) @@ -359,7 +416,7 @@ namespace adcval[cphi - 2][tbin] + adcval[cphi - 3][tbin]) { // rising again phidown = iphi + 1; - touch++; + counts.overlap++; break; } } @@ -369,6 +426,61 @@ namespace return; } + void check_cluster_touching(const std::vector& ihit_list, const std::vector>& adcval, int phibins, int tbins, ClusterCounters &counts) + { + // Encode (iphi, it) into single integer for fast lookup + std::unordered_set cluster_hits; + cluster_hits.reserve(ihit_list.size()); + + auto encode = [tbins](int phi, int t) + { + return phi * tbins + t; + }; + + for (const auto &hit : ihit_list) + { + cluster_hits.insert(encode(hit.iphi, hit.it)); + } + + for (const auto &hit : ihit_list) + { + int iphi = hit.iphi; + int it = hit.it; + + for (int dphi = -1; dphi <= 1; ++dphi) + { + for (int dt = -1; dt <= 1; ++dt) + { + if (dphi == 0 && dt == 0) { continue; } + + int nphi = iphi + dphi; + int nt = it + dt; + + if (nphi < 0 || nphi >= phibins || + nt < 0 || nt >= tbins) { + continue; + } + + // skip same cluster + if (cluster_hits.contains(encode(nphi, nt))) { continue; } + + // neighbor has signal → touching + if (adcval[nphi][nt] > 0 && + adcval[nphi][nt] != USHRT_MAX) + { + // Check Phi + if (dphi == -1) { counts.slmix = 1; } + if (dphi == 1) { counts.srmix = 1; } + + // Check Time + if (dt == -1) { counts.tlmix = 1; } + if (dt == 1) { counts.trmix = 1; } + } + } + } + } + } + int is_hit_isolated(int iphi, int it, int NPhiBinsMax, int NTBinsMax, const std::vector> &adcval) { // check isolated hits @@ -426,20 +538,25 @@ namespace return isiso; } - void get_cluster(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, std::vector &ihit_list, int &touch, int &edge) + void get_cluster(int phibin, int tbin, const thread_data &my_data, const std::vector> &adcval, std::vector &ihit_list, ClusterCounters &counts) { + bool ttop_edge = false; + bool tbottom_edge = false; + bool phitop_edge = false; + bool phibottom_edge = false; + // search along phi at the peak in t // const int NPhiBinsMax = (int) my_data.phibins; // const int NTBinsMax = (int) my_data.tbins; int tup = 0; int tdown = 0; - find_t_range(phibin, tbin, my_data, adcval, tdown, tup, touch, edge); + find_t_range(phibin, tbin, my_data, adcval, tdown, tup, counts, ttop_edge, tbottom_edge); // now we have the t extent of the cluster, go find the phi edges for (int it = tbin - tdown; it <= (tbin + tup); it++) { int phiup = 0; int phidown = 0; - find_phi_range(phibin, it, my_data, adcval, phidown, phiup, touch, edge); + find_phi_range(phibin, it, my_data, adcval, phidown, phiup, counts, phitop_edge, phibottom_edge); for (int iphi = (phibin - phidown); iphi <= (phibin + phiup); iphi++) { if (adcval[iphi][it] > 0 && adcval[iphi][it] != USHRT_MAX) @@ -456,7 +573,7 @@ namespace hit.it = it; hit.adc = adcval[iphi][it]; - if (touch > 0) + if (counts.overlap > 0) { if ((iphi == (phibin - phidown)) || (iphi == (phibin + phiup))) @@ -472,7 +589,7 @@ namespace } void calc_cluster_parameter(const int iphi_center, const int it_center, - const std::vector &ihit_list, thread_data &my_data, int ntouch, int nedge) + const std::vector &ihit_list, thread_data &my_data, ClusterCounters counts) { // // get z range from layer geometry @@ -488,6 +605,8 @@ namespace double iphi_sum = 0.0; double iphi2_sum = 0.0; + double it_sum = 0.0; + double radius = my_data.layergeom->get_radius(); // returns center of layer int phibinhi = -1; @@ -497,6 +616,12 @@ namespace int clus_size = ihit_list.size(); int max_adc = 0; + int phibinmax = -1; + int tbinmax = -1; + double cen_adc = 0; + + int size = 0; + if (clus_size <= my_data.min_clus_size) { return; @@ -521,14 +646,16 @@ namespace training_hits->phistep = my_data.layergeom->get_phistep(); training_hits->zstep = my_data.layergeom->get_zstep() * my_data.tGeometry->get_drift_velocity(); training_hits->layer = my_data.layer; - training_hits->ntouch = ntouch; - training_hits->nedge = nedge; + training_hits->ntouch = counts.overlap; + training_hits->nedge = counts.nedge; training_hits->v_adc.fill(0); } // std::cout << "process list" << std::endl; std::vector hitkeyvec; + std::map, double> adc_map; + // keep track of the hit locations in a given cluster std::map m_phi{}; std::map m_z{}; @@ -544,7 +671,18 @@ namespace continue; } - max_adc = std::max(max_adc, static_cast(std::round(adc))); // preserves rounding (0.5 -> 1) + size++; + + int adc_int = static_cast(std::round(adc)); + + if (adc_int > max_adc) + { + max_adc = adc_int; + phibinmax = iphi; + tbinmax = it; + } + + // max_adc = std::max(max_adc, static_cast(std::round(adc))); // preserves rounding (0.5 -> 1) phibinhi = std::max(iphi, phibinhi); phibinlo = std::min(iphi, phibinlo); tbinhi = std::max(it, tbinhi); @@ -565,8 +703,12 @@ namespace t_sum += t * adc; t2_sum += square(t) * adc; + it_sum += it * adc; + adc_sum += adc; + adc_map[{iphi, it}] += adc; + if (my_data.fillClusHitsVerbose) { auto pnew = m_phi.try_emplace(iphi, adc); @@ -623,13 +765,15 @@ namespace left_pad >= my_data.phioffset && deadset.contains(TpcDefs::genHitKey(left_pad, 0))) { - nedge++; + counts.nedge++; + counts.dledge = 1; } if (right_pad < (my_data.phibins + my_data.phioffset) && deadset.contains(TpcDefs::genHitKey(right_pad, 0))) { - nedge++; + counts.nedge++; + counts.dredge = 1; } } } @@ -646,24 +790,62 @@ namespace left_pad >= my_data.phioffset && hotset.contains(TpcDefs::genHitKey(left_pad, 0))) { - nedge++; + counts.nedge++; + counts.hledge = 1; } if (right_pad < (my_data.phibins + my_data.phioffset) && hotset.contains(TpcDefs::genHitKey(right_pad, 0))) { - nedge++; + counts.nedge++; + counts.hredge = 1; } } } - // This is the global position + // This is local position double clusiphi = iphi_sum / adc_sum; + double clusit = it_sum / adc_sum; + + // This is the global position double clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side); + double clust = t_sum / adc_sum; + + // ADC of centroid bin + int iphi_centroid = static_cast(std::floor(clusiphi)); + int it_centroid = static_cast(std::floor(clusit)); + + auto it_cent = adc_map.find({iphi_centroid, it_centroid}); + if (it_cent != adc_map.end()) + { + cen_adc = it_cent->second; + } + else + { + cen_adc = 0.0; // centroid may not land on a real hit + } + + // Max ADC position in global coordinates + double maxphi = my_data.layergeom->get_phi(phibinmax, my_data.side); + double maxt = my_data.layergeom->get_zcenter(tbinmax); + + // Phase relative to max ADC position + double padphase = 0.0; + double tbinphase = 0.0; + + if (my_data.layergeom->get_phistep() > 0) + { + padphase = (clusphi - maxphi) / my_data.layergeom->get_phistep(); + } + + if (my_data.layergeom->get_zstep() > 0) + { + tbinphase = (clust - maxt) / my_data.layergeom->get_zstep(); + } double clusx = radius * cos(clusphi); double clusy = radius * sin(clusphi); - double clust = t_sum / adc_sum; + // needed for surface identification double zdriftlength = clust * my_data.tGeometry->get_drift_velocity(); // convert z drift length to z position in the TPC @@ -702,6 +884,7 @@ namespace char tsize = tbinhi - tbinlo + 1; char phisize = phibinhi - phibinlo + 1; + char rsize = size; // std::cout << "phisize: " << (int) phisize << " phibinhi " << phibinhi << " phibinlo " << phibinlo << std::endl; // phi_cov = (weighted mean of dphi^2) - (weighted mean of dphi)^2, which is essentially the weighted mean of dphi^2. The error is then: // e_phi = sigma_dphi/sqrt(N) = sqrt( sigma_dphi^2 / N ) -- where N is the number of samples of the distribution with standard deviation sigma_dphi @@ -726,20 +909,54 @@ namespace // std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl; if (sqrt(phi_err_square) > my_data.min_err_squared) { - auto *clus = new TrkrClusterv5; - // auto clus = std::make_unique(); + TrkrCluster* clus = nullptr; + + if (my_data.debug) + { + clus = new TrkrClusterv6; + } + else + { + clus = new TrkrClusterv5; + } + clus_base = clus; - clus->setAdc(adc_sum); - clus->setMaxAdc(max_adc); - clus->setEdge(nedge); - clus->setPhiSize(phisize); - clus->setZSize(tsize); - clus->setSubSurfKey(subsurfkey); - clus->setOverlap(ntouch); clus->setLocalX(local(0)); clus->setLocalY(clust); + clus->setSubSurfKey(subsurfkey); + clus->setAdc(adc_sum); + clus->setMaxAdc(max_adc); + clus->setCenAdc(cen_adc); + clus->setPadCen(clusiphi); + clus->setTBinCen(clusit); + clus->setPadMax(phibinmax); + clus->setTBinMax(tbinmax); clus->setPhiError(sqrt(phi_err_square)); clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(), 2))); + clus->setRSize(rsize); + clus->setPhiSize(phisize); + clus->setZSize(tsize); + clus->setOverlap(counts.overlap); + clus->setEdge(counts.nedge); + clus->setSLEdge(counts.sledge); + clus->setSREdge(counts.sredge); + clus->setTLEdge(counts.tledge); + clus->setTREdge(counts.tredge); + clus->setDLEdge(counts.dledge); + clus->setDREdge(counts.dredge); + clus->setHLEdge(counts.hledge); + clus->setHREdge(counts.hredge); + clus->setSLMix(counts.slmix); + clus->setSRMix(counts.srmix); + clus->setTLMix(counts.tlmix); + clus->setTRMix(counts.trmix); + clus->setPhiBinLo(phibinlo); + clus->setPhiBinHi(phibinhi); + clus->setTBinLo(tbinlo); + clus->setTBinHi(tbinhi); + clus->setPadPhase(padphase); + clus->setTBinPhase(tbinphase); + my_data.cluster_vector.push_back(clus); b_made_cluster = true; } @@ -1037,6 +1254,9 @@ namespace } } */ + + std::vector> adcval_orig = adcval; + // std::cout << "done filling " << std::endl; while (!all_hit_map.empty()) { @@ -1064,9 +1284,10 @@ namespace // start with highest adc hit // -> cluster around it and get vector of hits std::vector ihit_list; - int ntouch = 0; - int nedge = 0; - get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge); + // Setting all the counters + ClusterCounters counts; + + get_cluster(iphi, it, *my_data, adcval, ihit_list, counts); if (my_data->FixedWindow > 0) { @@ -1102,11 +1323,16 @@ namespace my_data->FixedWindow = 0; // reset hit list and try again without fixed window ihit_list.clear(); - get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge); + // resetting all the counters + counts.clear(); + get_cluster(iphi, it, *my_data, adcval, ihit_list, counts); // std::cout << " stepdown size after " << ihit_list.size() << std::endl; my_data->FixedWindow = window_cache; } } + + check_cluster_touching(ihit_list, adcval_orig, my_data->phibins, my_data->tbins, counts); + if (ihit_list.size() <= 1) { remove_hits(ihit_list, all_hit_map, adcval); @@ -1117,7 +1343,7 @@ namespace // -> add hits to truth association // remove hits from all_hit_map // repeat untill all_hit_map empty - calc_cluster_parameter(iphi, it, ihit_list, *my_data, ntouch, nedge); + calc_cluster_parameter(iphi, it, ihit_list, *my_data, counts); remove_hits(ihit_list, all_hit_map, adcval); ihit_list.clear(); } @@ -1536,7 +1762,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) thread_pair.data.phioffset = PhiOffset; thread_pair.data.tbins = NTBinsSide; thread_pair.data.toffset = TOffset; - + thread_pair.data.debug = m_debug; thread_pair.data.radius = layergeom->get_radius(); thread_pair.data.drift_velocity = m_tGeometry->get_drift_velocity(); thread_pair.data.pads_per_sector = 0; @@ -1625,6 +1851,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) thread_pair.data.pedestal = pedestal; thread_pair.data.sector = sector; thread_pair.data.side = side; + thread_pair.data.debug = m_debug; thread_pair.data.do_assoc = do_hit_assoc; thread_pair.data.do_wedge_emulation = do_wedge_emulation; thread_pair.data.tGeometry = m_tGeometry; diff --git a/offline/packages/tpc/TpcClusterizer.h b/offline/packages/tpc/TpcClusterizer.h index 801207654e..c566f5b4f8 100644 --- a/offline/packages/tpc/TpcClusterizer.h +++ b/offline/packages/tpc/TpcClusterizer.h @@ -91,6 +91,11 @@ class TpcClusterizer : public SubsysReco m_hotChannelMapName = hmap; } + void DetailedClusterAnalysis() + { + m_debug = true; + } + private: bool is_in_sector_boundary(int phibin, int sector, PHG4TpcGeom *layergeom) const; bool record_ClusHitsVerbose{false}; @@ -135,6 +140,7 @@ class TpcClusterizer : public SubsysReco bool m_maskDeadChannels {false}; bool m_maskHotChannels {false}; bool m_maskFromFile {false}; + bool m_debug{false}; std::string m_deadChannelMapName; std::string m_hotChannelMapName; }; diff --git a/offline/packages/trackbase/TrkrCluster.h b/offline/packages/trackbase/TrkrCluster.h index f37440646b..2035e19796 100644 --- a/offline/packages/trackbase/TrkrCluster.h +++ b/offline/packages/trackbase/TrkrCluster.h @@ -52,21 +52,21 @@ class TrkrCluster : public PHObject // cluster position // virtual float getLocalX() const { return std::numeric_limits::quiet_NaN(); } - virtual void setLocalX(float) {} + virtual void setLocalX(const float) {} virtual float getLocalY() const { return std::numeric_limits::quiet_NaN(); } - virtual void setLocalY(float) {} + virtual void setLocalY(const float) {} // // cluster info // - virtual void setAdc(unsigned int) {} + virtual void setAdc(const unsigned int) {} virtual unsigned int getAdc() const { return UINT_MAX; } - virtual void setMaxAdc(uint16_t) {} + virtual void setMaxAdc(const uint16_t) {} virtual unsigned int getMaxAdc() const { return UINT_MAX; } virtual char getOverlap() const { return std::numeric_limits::max(); } - virtual void setOverlap(char) {} + virtual void setOverlap(const char) {} virtual char getEdge() const { return std::numeric_limits::max(); } - virtual void setEdge(char) {} + virtual void setEdge(const char) {} virtual void setTime(const float) {} virtual float getTime() const { return std::numeric_limits::quiet_NaN(); } virtual char getSize() const { return std::numeric_limits::max(); } @@ -82,8 +82,8 @@ class TrkrCluster : public PHObject virtual unsigned int getCenAdc() const { return UINT_MAX; } virtual float getPadCen() const { return std::numeric_limits::quiet_NaN(); } virtual float getTBinCen() const { return std::numeric_limits::quiet_NaN(); } - virtual float getPadMax() const { return std::numeric_limits::quiet_NaN(); } - virtual float getTBinMax() const { return std::numeric_limits::quiet_NaN(); } + virtual int getPadMax() const { return std::numeric_limits::max(); } + virtual int getTBinMax() const { return std::numeric_limits::max(); } virtual char getSLEdge() const { return std::numeric_limits::max(); } virtual char getSREdge() const { return std::numeric_limits::max(); } virtual char getTLEdge() const { return std::numeric_limits::max(); } @@ -92,23 +92,52 @@ class TrkrCluster : public PHObject virtual char getDREdge() const { return std::numeric_limits::max(); } virtual char getHLEdge() const { return std::numeric_limits::max(); } virtual char getHREdge() const { return std::numeric_limits::max(); } - virtual int getSLMix() const { return std::numeric_limits::max(); } - virtual int getSRMix() const { return std::numeric_limits::max(); } - virtual int getTLMix() const { return std::numeric_limits::max(); } - virtual int getTRMix() const { return std::numeric_limits::max(); } - virtual float getPhiBinLo() const { return std::numeric_limits::quiet_NaN(); } - virtual float getPhiBinHi() const { return std::numeric_limits::quiet_NaN(); } - virtual float getTBinLo() const { return std::numeric_limits::quiet_NaN(); } - virtual float getTBinHi() const { return std::numeric_limits::quiet_NaN(); } + virtual char getSLMix() const { return std::numeric_limits::max(); } + virtual char getSRMix() const { return std::numeric_limits::max(); } + virtual char getTLMix() const { return std::numeric_limits::max(); } + virtual char getTRMix() const { return std::numeric_limits::max(); } + virtual unsigned short getPhiBinLo() const { return std::numeric_limits::max(); } + virtual unsigned short getPhiBinHi() const { return std::numeric_limits::max(); } + virtual unsigned short getTBinLo() const { return std::numeric_limits::max(); } + virtual unsigned short getTBinHi() const { return std::numeric_limits::max(); } virtual float getPadPhase() const { return std::numeric_limits::quiet_NaN(); } virtual float getTBinPhase() const { return std::numeric_limits::quiet_NaN(); } virtual float getRSize() const { return std::numeric_limits::quiet_NaN(); } + virtual void setSLEdge(const char) {}; + virtual void setSREdge(const char) {}; + virtual void setTLEdge(const char) {}; + virtual void setTREdge(const char) {}; + virtual void setDLEdge(const char) {}; + virtual void setDREdge(const char) {}; + virtual void setHLEdge(const char) {}; + virtual void setHREdge(const char) {}; + virtual void setSLMix(const char) {}; + virtual void setSRMix(const char) {}; + virtual void setTLMix(const char) {}; + virtual void setTRMix(const char) {}; + virtual void setPhiBinLo(const unsigned short) {}; + virtual void setPhiBinHi(const unsigned short) {}; + virtual void setTBinLo(const unsigned short) {}; + virtual void setTBinHi(const unsigned short) {}; + virtual void setPadPhase(const float) {}; + virtual void setTBinPhase(const float) {}; + virtual void setRSize(const char) {}; + virtual void setCenAdc(const uint16_t) {}; + virtual void setPadCen(const float) {}; + virtual void setTBinCen(const float) {}; + virtual void setPadMax(const int) {}; + virtual void setTBinMax(const int) {}; + virtual void setPhiError(const float) {}; + virtual void setZError(const float) {}; + virtual void setPhiSize(const char) {}; + virtual void setZSize(const char) {}; + /// Acts functions, for Acts modules use only virtual void setActsLocalError(unsigned int /*i*/, unsigned int /*j*/, float /*value*/) {} virtual float getActsLocalError(unsigned int /*i*/, unsigned int /*j*/) const { return std::numeric_limits::quiet_NaN(); } virtual TrkrDefs::subsurfkey getSubSurfKey() const { return TrkrDefs::SUBSURFKEYMAX; } - virtual void setSubSurfKey(TrkrDefs::subsurfkey /*id*/) {} + virtual void setSubSurfKey(const TrkrDefs::subsurfkey /*id*/) {} // Global coordinate functions are deprecated, use local // coordinate functions only diff --git a/offline/packages/trackbase/TrkrClusterv4.cc b/offline/packages/trackbase/TrkrClusterv4.cc index 542c1530cc..8ee11bf797 100644 --- a/offline/packages/trackbase/TrkrClusterv4.cc +++ b/offline/packages/trackbase/TrkrClusterv4.cc @@ -13,7 +13,7 @@ namespace { // square convenience function template - inline constexpr T square(const T& x) + constexpr T square(const T& x) { return x * x; } diff --git a/offline/packages/trackbase/TrkrClusterv4.h b/offline/packages/trackbase/TrkrClusterv4.h index ec0bc4703c..439e3d9aef 100644 --- a/offline/packages/trackbase/TrkrClusterv4.h +++ b/offline/packages/trackbase/TrkrClusterv4.h @@ -160,10 +160,10 @@ class TrkrClusterv4 : public TrkrCluster // void setSize(char size) { m_size = size; } float getPhiSize() const override { return (float) m_phisize; } - void setPhiSize(char phisize) { m_phisize = phisize; } + void setPhiSize(char phisize) override { m_phisize = phisize; } float getZSize() const override { return (float) m_zsize; } - void setZSize(char zsize) { m_zsize = zsize; } + void setZSize(char zsize) override { m_zsize = zsize; } char getOverlap() const override { return m_overlap; } void setOverlap(char overlap) override { m_overlap = overlap; } diff --git a/offline/packages/trackbase/TrkrClusterv5.cc b/offline/packages/trackbase/TrkrClusterv5.cc index 58e08745ad..a0cc7fbe52 100644 --- a/offline/packages/trackbase/TrkrClusterv5.cc +++ b/offline/packages/trackbase/TrkrClusterv5.cc @@ -13,7 +13,7 @@ namespace { // square convenience function template - inline constexpr T square(const T& x) + constexpr T square(const T& x) { return x * x; } diff --git a/offline/packages/trackbase/TrkrClusterv5.h b/offline/packages/trackbase/TrkrClusterv5.h index ebb0bae961..ee3ac20755 100644 --- a/offline/packages/trackbase/TrkrClusterv5.h +++ b/offline/packages/trackbase/TrkrClusterv5.h @@ -54,12 +54,12 @@ class TrkrClusterv5 : public TrkrCluster float getPosition(int coor) const override { return m_local[coor]; } void setPosition(int coor, float xi) override { m_local[coor] = xi; } float getLocalX() const override { return m_local[0]; } - void setLocalX(float loc0) override { m_local[0] = loc0; } + void setLocalX(const float loc0) override { m_local[0] = loc0; } float getLocalY() const override { return m_local[1]; } - void setLocalY(float loc1) override { m_local[1] = loc1; } + void setLocalY(const float loc1) override { m_local[1] = loc1; } TrkrDefs::subsurfkey getSubSurfKey() const override { return m_subsurfkey; } - void setSubSurfKey(TrkrDefs::subsurfkey id) override { m_subsurfkey = id; } + void setSubSurfKey(const TrkrDefs::subsurfkey id) override { m_subsurfkey = id; } // // cluster info @@ -69,7 +69,7 @@ class TrkrClusterv5 : public TrkrCluster return m_adc; } - void setAdc(unsigned int adc) override + void setAdc(const unsigned int adc) override { m_adc = adc; } @@ -79,7 +79,7 @@ class TrkrClusterv5 : public TrkrCluster return m_maxadc; } - void setMaxAdc(uint16_t maxadc) override + void setMaxAdc(const uint16_t maxadc) override { m_maxadc = maxadc; } @@ -96,11 +96,11 @@ class TrkrClusterv5 : public TrkrCluster return m_zerr; } - void setPhiError(float phierror) + void setPhiError(const float phierror) override { m_phierr = phierror; } - void setZError(float zerror) + void setZError(const float zerror) override { m_zerr = zerror; } @@ -156,16 +156,16 @@ class TrkrClusterv5 : public TrkrCluster // void setSize(char size) { m_size = size; } float getPhiSize() const override { return (float) m_phisize; } - void setPhiSize(char phisize) { m_phisize = phisize; } + void setPhiSize(const char phisize) override { m_phisize = phisize; } float getZSize() const override { return (float) m_zsize; } - void setZSize(char zsize) { m_zsize = zsize; } + void setZSize(const char zsize) override { m_zsize = zsize; } char getOverlap() const override { return m_overlap; } - void setOverlap(char overlap) override { m_overlap = overlap; } + void setOverlap(const char overlap) override { m_overlap = overlap; } char getEdge() const override { return m_edge; } - void setEdge(char edge) override { m_edge = edge; } + void setEdge(const char edge) override { m_edge = edge; } // float getPhiSize() const override //{ std::cout << "Deprecated size function"<< std::endl; return NAN;} diff --git a/offline/packages/trackbase/TrkrClusterv6.h b/offline/packages/trackbase/TrkrClusterv6.h index 37ffdffbbf..683aeded6e 100644 --- a/offline/packages/trackbase/TrkrClusterv6.h +++ b/offline/packages/trackbase/TrkrClusterv6.h @@ -60,7 +60,7 @@ class TrkrClusterv6 : public TrkrCluster { return (coor >= 0 && coor < 2) ? m_local[coor] : std::numeric_limits::quiet_NaN(); } - void setPosition(int coor, float xi) override + void setPosition(const int coor, const float xi) override { if (coor >= 0 && coor < 2) { @@ -68,36 +68,36 @@ class TrkrClusterv6 : public TrkrCluster } } float getLocalX() const override { return m_local[0]; } - void setLocalX(float loc0) override { m_local[0] = loc0; } + void setLocalX(const float loc0) override { m_local[0] = loc0; } float getLocalY() const override { return m_local[1]; } - void setLocalY(float loc1) override { m_local[1] = loc1; } + void setLocalY(const float loc1) override { m_local[1] = loc1; } TrkrDefs::subsurfkey getSubSurfKey() const override { return m_subsurfkey; } - void setSubSurfKey(TrkrDefs::subsurfkey id) override { m_subsurfkey = id; } + void setSubSurfKey(const TrkrDefs::subsurfkey id) override { m_subsurfkey = id; } // // cluster info // unsigned int getAdc() const override { return m_adc; } - void setAdc(unsigned int adc) override { m_adc = adc; } + void setAdc(const unsigned int adc) override { m_adc = adc; } unsigned int getMaxAdc() const override { return m_maxadc; } - void setMaxAdc(uint16_t maxadc) override { m_maxadc = maxadc; } + void setMaxAdc(const uint16_t maxadc) override { m_maxadc = maxadc; } unsigned int getCenAdc() const override { return m_cenadc; } - void setCenAdc(uint16_t cenadc) { m_cenadc = cenadc; } + void setCenAdc(const uint16_t cenadc) override { m_cenadc = cenadc; } float getPadCen() const override { return m_padcen; } - void setPadCen(float padcen) { m_padcen = padcen; } + void setPadCen(const float padcen) override { m_padcen = padcen; } float getTBinCen() const override { return m_tbincen; } - void setTBinCen(float tbincen) { m_tbincen = tbincen; } + void setTBinCen(const float tbincen) override { m_tbincen = tbincen; } - float getPadMax() const override { return m_padmax; } - void setPadMax(float padmax) { m_padmax = padmax; } + int getPadMax() const override { return m_padmax; } + void setPadMax(const int padmax) override { m_padmax = padmax; } - float getTBinMax() const override { return m_tbinmax; } - void setTBinMax(float tbinmax) { m_tbinmax = tbinmax; } + int getTBinMax() const override { return m_tbinmax; } + void setTBinMax(const int tbinmax) override { m_tbinmax = tbinmax; } // // convenience interface @@ -105,117 +105,117 @@ class TrkrClusterv6 : public TrkrCluster float getRPhiError() const override { return m_phierr; } float getZError() const override { return m_zerr; } - void setPhiError(float phierror) { m_phierr = phierror; } - void setZError(float zerror) { m_zerr = zerror; } + void setPhiError(const float phierror) override { m_phierr = phierror; } + void setZError(const float zerror) override { m_zerr = zerror; } char getSize() const override { return m_phisize * m_zsize; } - // void setSize(char size) { m_size = size; } + // void setSize(const char size) { m_size = size; } float getRSize() const override { return (float) m_rsize; } - void setRSize(unsigned char rsize) { m_rsize = rsize; } + void setRSize(const char rsize) override { m_rsize = rsize; } float getPhiSize() const override { return (float) m_phisize; } - void setPhiSize(char phisize) { m_phisize = phisize; } + void setPhiSize(const char phisize) override { m_phisize = phisize; } float getZSize() const override { return (float) m_zsize; } - void setZSize(char zsize) { m_zsize = zsize; } + void setZSize(const char zsize) override { m_zsize = zsize; } char getOverlap() const override { return m_overlap; } - void setOverlap(char overlap) override { m_overlap = overlap; } + void setOverlap(const char overlap) override { m_overlap = overlap; } char getEdge() const override { return m_edge; } - void setEdge(char edge) override { m_edge = edge; } + void setEdge(const char edge) override { m_edge = edge; } char getSLEdge() const override { return m_sledge; } - void setSLEdge(char sledge) { m_sledge = sledge; } + void setSLEdge(const char sledge) override { m_sledge = sledge; } char getSREdge() const override { return m_sredge; } - void setSREdge(char sredge) { m_sredge = sredge; } + void setSREdge(const char sredge) override { m_sredge = sredge; } char getTLEdge() const override { return m_tledge; } - void setTLEdge(char tledge) { m_tledge = tledge; } + void setTLEdge(const char tledge) override { m_tledge = tledge; } char getTREdge() const override { return m_tredge; } - void setTREdge(char tredge) { m_tredge = tredge; } + void setTREdge(const char tredge) override { m_tredge = tredge; } char getDLEdge() const override { return m_dledge; } - void setDLEdge(char dledge) { m_dledge = dledge; } + void setDLEdge(const char dledge) override { m_dledge = dledge; } char getDREdge() const override { return m_dredge; } - void setDREdge(char dredge) { m_dredge = dredge; } + void setDREdge(const char dredge) override { m_dredge = dredge; } char getHLEdge() const override { return m_hledge; } - void setHLEdge(char hledge) { m_hledge = hledge; } + void setHLEdge(const char hledge) override { m_hledge = hledge; } char getHREdge() const override { return m_hredge; } - void setHREdge(char hredge) { m_hredge = hredge; } + void setHREdge(const char hredge) override { m_hredge = hredge; } - int getSLMix() const override { return m_slmix; } - void setSLMix(char slmix) { m_slmix = slmix; } + char getSLMix() const override { return m_slmix; } + void setSLMix(const char slmix) override { m_slmix = slmix; } - int getSRMix() const override { return m_srmix; } - void setSRMix(char srmix) { m_srmix = srmix; } + char getSRMix() const override { return m_srmix; } + void setSRMix(const char srmix) override { m_srmix = srmix; } - int getTLMix() const override { return m_tlmix; } - void setTLMix(char tlmix) { m_tlmix = tlmix; } + char getTLMix() const override { return m_tlmix; } + void setTLMix(const char tlmix) override { m_tlmix = tlmix; } - int getTRMix() const override { return m_trmix; } - void setTRMix(char trmix) { m_trmix = trmix; } + char getTRMix() const override { return m_trmix; } + void setTRMix(const char trmix) override { m_trmix = trmix; } - float getPhiBinLo() const override { return m_phibinlo; } - void setPhiBinLo(float phibinlo) { m_phibinlo = phibinlo; } + unsigned short getPhiBinLo() const override { return m_phibinlo; } + void setPhiBinLo(const unsigned short phibinlo) override { m_phibinlo = phibinlo; } - float getPhiBinHi() const override { return m_phibinhi; } - void setPhiBinHi(float phibinhi) { m_phibinhi = phibinhi; } + unsigned short getPhiBinHi() const override { return m_phibinhi; } + void setPhiBinHi(const unsigned short phibinhi) override { m_phibinhi = phibinhi; } - float getTBinLo() const override { return m_tbinlo; } - void setTBinLo(float tbinlo) { m_tbinlo = tbinlo; } + unsigned short getTBinLo() const override { return m_tbinlo; } + void setTBinLo(const unsigned short tbinlo) override { m_tbinlo = tbinlo; } - float getTBinHi() const override { return m_tbinhi; } - void setTBinHi(float tbinhi) { m_tbinhi = tbinhi; } + unsigned short getTBinHi() const override { return m_tbinhi; } + void setTBinHi(const unsigned short tbinhi) override { m_tbinhi = tbinhi; } float getPadPhase() const override { return m_padphase; } - void setPadPhase(float padphase) { m_padphase = padphase; } + void setPadPhase(const float padphase) override { m_padphase = padphase; } float getTBinPhase() const override { return m_tbinphase; } - void setTBinPhase(float tbinphase){ m_tbinphase = tbinphase; } + void setTBinPhase(const float tbinphase) override { m_tbinphase = tbinphase; } private: float m_local[2]{std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; //< 2D local position [cm] 2 * 32 64bit - cumul 1*64 TrkrDefs::subsurfkey m_subsurfkey {TrkrDefs::SUBSURFKEYMAX}; //< unique identifier for hitsetkey-surface maps 16 bit - float m_phierr{0}; - float m_zerr{0}; - unsigned short int m_adc{0}; //< cluster sum adc 16 - unsigned short int m_maxadc{0}; //< cluster max adc 16 - unsigned short int m_cenadc{0}; //< cluster centroid adc 16 - float m_padcen{0}; - float m_tbincen{0}; - float m_padmax{0}; - float m_tbinmax{0}; - unsigned char m_rsize{0}; // 8bit - char m_phisize{0}; // 8bit - char m_zsize{0}; // 8bit - char m_overlap{0}; // 8bit - char m_edge{0}; // 8bit - cumul 2*64 - char m_sledge{0}; // 8bit - char m_sredge{0}; // 8bit - char m_tledge{0}; // 8bit - char m_tredge{0}; // 8bit - char m_dledge{0}; // 8bit - char m_dredge{0}; // 8bit - char m_hledge{0}; // 8bit - char m_hredge{0}; // 8bit - char m_slmix{0}; // 8bit - char m_srmix{0}; // 8bit - char m_tlmix{0}; // 8bit - char m_trmix{0}; // 8bit - float m_phibinlo{0}; - float m_phibinhi{0}; - float m_tbinlo{0}; - float m_tbinhi{0}; - float m_padphase{0}; - float m_tbinphase{0}; + float m_phierr{std::numeric_limits::quiet_NaN()}; + float m_zerr{std::numeric_limits::quiet_NaN()}; + unsigned short m_adc{std::numeric_limits::max()}; //< cluster sum adc 16 + unsigned short m_maxadc{std::numeric_limits::max()}; //< cluster max adc 16 + unsigned short m_cenadc{std::numeric_limits::max()}; //< cluster centroid adc 16 + float m_padcen{std::numeric_limits::quiet_NaN()}; + float m_tbincen{std::numeric_limits::quiet_NaN()}; + int m_padmax{std::numeric_limits::max()}; + int m_tbinmax{std::numeric_limits::max()}; + char m_rsize{std::numeric_limits::max()}; + char m_phisize{std::numeric_limits::max()}; + char m_zsize{std::numeric_limits::max()}; + char m_overlap{std::numeric_limits::max()}; + char m_edge{std::numeric_limits::max()}; + char m_sledge{std::numeric_limits::max()}; + char m_sredge{std::numeric_limits::max()}; + char m_tledge{std::numeric_limits::max()}; + char m_tredge{std::numeric_limits::max()}; + char m_dledge{std::numeric_limits::max()}; + char m_dredge{std::numeric_limits::max()}; + char m_hledge{std::numeric_limits::max()}; + char m_hredge{std::numeric_limits::max()}; + char m_slmix{std::numeric_limits::max()}; + char m_srmix{std::numeric_limits::max()}; + char m_tlmix{std::numeric_limits::max()}; + char m_trmix{std::numeric_limits::max()}; + unsigned short m_phibinlo{std::numeric_limits::max()}; + unsigned short m_phibinhi{std::numeric_limits::max()}; + unsigned short m_tbinlo{std::numeric_limits::max()}; + unsigned short m_tbinhi{std::numeric_limits::max()}; + float m_padphase{std::numeric_limits::quiet_NaN()}; + float m_tbinphase{std::numeric_limits::quiet_NaN()}; ClassDefOverride(TrkrClusterv6, 1) };