From 594f6915b260766c9f35d255d9b495c3df706300 Mon Sep 17 00:00:00 2001 From: bkimelman Date: Tue, 23 May 2023 11:24:37 -0400 Subject: [PATCH 1/5] Personal changes, testing basics of using runNode --- .../PHTpcCentralMembraneClusterizer.cc | 7 +++-- .../tpccalib/PHTpcCentralMembraneMatcher.cc | 27 ++++++++++++++++--- .../tpccalib/PHTpcCentralMembraneMatcher.h | 7 +++-- 3 files changed, 33 insertions(+), 8 deletions(-) diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc index 6ce2370f0d..633835bc77 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc @@ -177,6 +177,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) mean_z_content_minus = mean_z_content_minus/hz_neg->GetNbinsX(); + std::cout << "mean_z_content_plus " << mean_z_content_plus << " mean_z_content_minus " << mean_z_content_minus << std::endl; + //use peak in z distributions to identify if there is a CM Flash. Peak should be >20% of entries (needs to be tuned) if(hz_pos->GetMaximum() < 5*mean_z_content_plus || hz_neg->GetMaximum() < 5*mean_z_content_minus) return Fun4AllReturnCodes::EVENT_OK; @@ -512,10 +514,12 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) } } else { + if(isAcrossGap[i]) continue; + if(_histos) hClustE[2]->Fill(energy[i]); // These single cluster cases have good phi, but do not have a good radius centroid estimate - may want to skip them, record nclusters and identify if across gap // if(layer[i] == 7) isAcrossGap[i] = true; - + aveenergy.push_back(energy[i]); avepos.push_back(pos[i]); nclusters.push_back(1); @@ -535,7 +539,6 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) { auto cmfc = new CMFlashClusterv3(); - cmfc->setX(avepos[iv].X()); cmfc->setY(avepos[iv].Y()); cmfc->setZ(avepos[iv].Z()); diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc index 59c6e9fa61..ab3504a237 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc @@ -267,7 +267,7 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) hdrphi = new TH1F("hdrphi","r * dphi", 200, -0.05, 0.05); hnclus = new TH1F("hnclus", " nclusters ", 3, 0., 3.); - fout2.reset ( new TFile(m_histogramfilename2.c_str(),"RECREATE") ); + m_debugfile.reset ( new TFile(m_debugfilename.c_str(),"RECREATE") ); match_ntup = new TNtuple("match_ntup","Match NTuple","event:truthR:truthPhi:recoR:recoPhi:recoZ:nclus:r1:phi1:e1:layer1:r2:phi2:e2:layer2"); } @@ -817,14 +817,16 @@ int PHTpcCentralMembraneMatcher::End(PHCompositeNode * /*topNode*/ ) fout->Close(); } - if(m_savehistograms && fout2) + if(m_savehistograms && m_debugfile) { - fout2->cd(); + m_debugfile->cd(); match_ntup->Write(); hit_r_phi->Write(); clust_r_phi_pos->Write(); clust_r_phi_neg->Write(); + + m_debugfile->Close(); } return Fun4AllReturnCodes::EVENT_OK; @@ -862,7 +864,24 @@ int PHTpcCentralMembraneMatcher::GetNodes(PHCompositeNode* topNode) // create node for results of matching std::cout << "Creating node CM_FLASH_DIFFERENCES" << std::endl; PHNodeIterator iter(topNode); + + // Looking for the DST node + PHCompositeNode *runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); + if (!runNode) + { + std::cout << PHWHERE << "RUN Node missing, doing nothing." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + PHNodeIterator runiter(runNode); + m_cm_flash_diffs = new CMFlashDifferenceContainerv1; + PHIODataNode *CMFlashDifferenceNode = + new PHIODataNode(m_cm_flash_diffs, "CM_FLASH_DIFFERENCES", "PHObject"); + runNode->addNode(CMFlashDifferenceNode); + + + + /* // Looking for the DST node PHCompositeNode *dstNode = dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); if (!dstNode) @@ -883,7 +902,7 @@ int PHTpcCentralMembraneMatcher::GetNodes(PHCompositeNode* topNode) PHIODataNode *CMFlashDifferenceNode = new PHIODataNode(m_cm_flash_diffs, "CM_FLASH_DIFFERENCES", "PHObject"); DetNode->addNode(CMFlashDifferenceNode); - + */ //// output tpc fluctuation distortion container //// this one is filled on the fly on a per-CM-event basis, and applied in the tracking chain diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h index 6d2d4f1308..7f7e7f7ab5 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h @@ -51,6 +51,9 @@ class PHTpcCentralMembraneMatcher : public SubsysReco void setOutputfile(const std::string &outputfile) {m_outputfile = outputfile;} + void setDebugOutputFile(const std::string &debugfile) + {m_debugfilename = debugfile;} + void setNMatchIter( int val ){ m_nMatchIter = val; } void set_useOnly_nClus2( bool val ){ m_useOnly_nClus2 = val; } @@ -118,8 +121,8 @@ class PHTpcCentralMembraneMatcher : public SubsysReco std::unique_ptr fout; - std::unique_ptr fout2; - std::string m_histogramfilename2 = "CMMatcher.root"; + std::unique_ptr m_debugfile; + std::string m_debugfilename = "CMMatcher.root"; TH2F *hit_r_phi; From dc2f600850cacce7962377b4a28dadd8126df6b0 Mon Sep 17 00:00:00 2001 From: bkimelman Date: Fri, 2 Jun 2023 08:27:27 -0400 Subject: [PATCH 2/5] Changes for tuning --- .../PHTpcCentralMembraneClusterizer.cc | 7 ++++-- .../tpccalib/PHTpcCentralMembraneMatcher.cc | 22 +++++++++++++++---- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc index 633835bc77..f8eac792a6 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc @@ -242,10 +242,10 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) //Find which phiMod bin cluster is in (with given layer) and if counts below threshold don't use if(z >= 0){ phiBin = hphi_reco_pos[lay-7]->GetXaxis()->FindBin(phiMod); - if(hphi_reco_pos[lay-7]->GetBinContent(phiBin) < m_moduloThreshold) aboveThreshold = false; + if(hphi_reco_pos[lay-7]->GetBinContent(phiBin) < (m_moduloThreshold + (mean_z_content_plus/18.0))) aboveThreshold = false; }else{ phiBin = hphi_reco_neg[lay-7]->GetXaxis()->FindBin(phiMod); - if(hphi_reco_neg[lay-7]->GetBinContent(phiBin) < m_moduloThreshold) aboveThreshold = false; + if(hphi_reco_neg[lay-7]->GetBinContent(phiBin) < (m_moduloThreshold + (mean_z_content_minus/18.0))) aboveThreshold = false; } if( !aboveThreshold ) continue; @@ -629,6 +629,9 @@ int PHTpcCentralMembraneClusterizer::End(PHCompositeNode * /*topNode*/ ) hDist2Adj->Write(); for(int i=0; i<47; i++){ + hphi_reco_pos[i]->Write(); + hphi_reco_neg[i]->Write(); + hphi_reco_pair_pos[i]->Write(); hphi_reco_pair_neg[i]->Write(); } diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc index ce718f0ffb..5f24033fc6 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc @@ -506,6 +506,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) R23Gap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin()+R12Gap_neg+1,clust_RGaps_neg.end())); } + int min_match_pos = 100; + int min_match_neg = 100; std::vector hitMatches_pos; //match cluster peaks to hit peaks using gap positions for(int i=0; i<(int)clust_RPeaks_pos.size(); i++){ @@ -518,6 +520,7 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) if(clust_RGaps_pos[R12Gap_pos] > 3.6) hitMatches_pos[i] -= 1; if(clust_RGaps_pos[R12Gap_pos] > 4.6) hitMatches_pos[i] -= 1; if(clust_RGaps_pos[R12Gap_pos] > 5.8) hitMatches_pos[i] -= 1; + if(hitMatches_pos[i] < min_match_pos) min_match_pos = hitMatches_pos[i]; } //module 1-2 gap is between 15 & 16 else if(i < (R23Gap_pos+1) && i >= (R12Gap_pos+1)) hitMatches_pos.push_back(15+1 + i - (R12Gap_pos+1)); @@ -534,6 +537,7 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) if(clust_RGaps_neg[R12Gap_neg] > 3.6) hitMatches_neg[i] -= 1; if(clust_RGaps_neg[R12Gap_neg] > 4.6) hitMatches_neg[i] -= 1; if(clust_RGaps_neg[R12Gap_neg] > 5.8) hitMatches_neg[i] -= 1; + if(hitMatches_neg[i] < min_match_neg) min_match_neg = hitMatches_neg[i]; } else if(i < (R23Gap_neg+1) && i >= (R12Gap_neg+1)) hitMatches_neg.push_back(15+1 + i - (R12Gap_neg+1)); else if(i >= R23Gap_neg+1) hitMatches_neg.push_back(23+1 + i - (R23Gap_neg+1)); @@ -676,6 +680,16 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) const auto n_reco_size2 = std::count_if( reco_nclusters.begin(), reco_nclusters.end(), []( const unsigned int& value ) { return value==2; } ); std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size: " << m_truth_pos.size() << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size, r>30cm: " << n_valid_truth << std::endl; + int pos_stripes_add = 0; + int neg_stripes_add = 0; + int nStr[8] = {3,4,4,4,4,5,4,5}; + for(int str=min_match_pos; str <= 7; str++){ + pos_stripes_add += nStr[str]*18; + } + for(int str=min_match_neg; str <= 7; str++){ + neg_stripes_add += nStr[str]*18; + } + std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size, r>30cm + rows below with match: " << n_valid_truth + pos_stripes_add + neg_stripes_add<< std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size: " << reco_pos.size() << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==1): " << n_reco_size1 << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==2): " << n_reco_size2 << std::endl; @@ -867,7 +881,7 @@ int PHTpcCentralMembraneMatcher::GetNodes(PHCompositeNode* topNode) PHNodeIterator iter(topNode); - // Looking for the DST node + // Looking for the RUN node PHCompositeNode *runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); if (!runNode) { @@ -911,9 +925,9 @@ int PHTpcCentralMembraneMatcher::GetNodes(PHCompositeNode* topNode) m_dcc_out = findNode::getClass(topNode,dcc_out_node_name); if( !m_dcc_out ) { - - /// Get the RUN node and check - auto runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); + + /// Get the RUN node and check + //auto runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); if (!runNode) { std::cout << "PHTpcCentralMembraneMatcher::InitRun - RUN Node missing, quitting" << std::endl; From 4ad48f222503ca1abb65126766776c7762fcffa7 Mon Sep 17 00:00:00 2001 From: bkimelman Date: Wed, 14 Jun 2023 11:20:56 -0400 Subject: [PATCH 3/5] Changes for tracking hit truth --- offline/packages/trackbase/Makefile.am | 16 +- offline/packages/trackbase/TrkrCluster.h | 5 + offline/packages/trackbase/TrkrClusterv6.cc | 83 +++++++++ offline/packages/trackbase/TrkrClusterv6.h | 163 ++++++++++++++++++ .../packages/trackbase/TrkrClusterv6LinkDef.h | 5 + offline/packages/trackbase/TrkrHit.h | 3 + offline/packages/trackbase/TrkrHitv3.cc | 52 ++++++ offline/packages/trackbase/TrkrHitv3.h | 59 +++++++ offline/packages/trackbase/TrkrHitv3LinkDef.h | 5 + .../g4tpc/PHG4TpcCentralMembrane.cc | 12 +- .../g4tpc/PHG4TpcCentralMembrane.h | 5 +- .../g4tpc/PHG4TpcElectronDrift.cc | 12 +- 12 files changed, 403 insertions(+), 17 deletions(-) create mode 100644 offline/packages/trackbase/TrkrClusterv6.cc create mode 100644 offline/packages/trackbase/TrkrClusterv6.h create mode 100644 offline/packages/trackbase/TrkrClusterv6LinkDef.h create mode 100644 offline/packages/trackbase/TrkrHitv3.cc create mode 100644 offline/packages/trackbase/TrkrHitv3.h create mode 100644 offline/packages/trackbase/TrkrHitv3LinkDef.h diff --git a/offline/packages/trackbase/Makefile.am b/offline/packages/trackbase/Makefile.am index 41db1bb5ed..c5a5abf1f7 100644 --- a/offline/packages/trackbase/Makefile.am +++ b/offline/packages/trackbase/Makefile.am @@ -85,6 +85,7 @@ pkginclude_HEADERS = \ TrkrClusterv3.h \ TrkrClusterv4.h \ TrkrClusterv5.h \ + TrkrClusterv6.h \ TrkrDefs.h \ TrkrHit.h \ TrkrHitSet.h \ @@ -94,7 +95,8 @@ pkginclude_HEADERS = \ TrkrHitTruthAssoc.h \ TrkrHitTruthAssocv1.h \ TrkrHitv1.h \ - TrkrHitv2.h + TrkrHitv2.h \ + TrkrHitv3.h ROOTDICTS = \ @@ -137,6 +139,7 @@ ROOTDICTS = \ TrkrClusterv3_Dict.cc \ TrkrClusterv4_Dict.cc \ TrkrClusterv5_Dict.cc \ + TrkrClusterv6_Dict.cc \ TrkrHitSetContainer_Dict.cc \ TrkrHitSetContainerv1_Dict.cc \ TrkrHitSet_Dict.cc \ @@ -145,7 +148,8 @@ ROOTDICTS = \ TrkrHitTruthAssocv1_Dict.cc \ TrkrHit_Dict.cc \ TrkrHitv1_Dict.cc \ - TrkrHitv2_Dict.cc + TrkrHitv2_Dict.cc \ + TrkrHitv3_Dict.cc pcmdir = $(libdir) @@ -189,6 +193,7 @@ nobase_dist_pcm_DATA = \ TrkrClusterv3_Dict_rdict.pcm \ TrkrClusterv4_Dict_rdict.pcm \ TrkrClusterv5_Dict_rdict.pcm \ + TrkrClusterv6_Dict_rdict.pcm \ TrkrHitSetContainer_Dict_rdict.pcm \ TrkrHitSetContainerv1_Dict_rdict.pcm \ TrkrHitSet_Dict_rdict.pcm \ @@ -197,7 +202,8 @@ nobase_dist_pcm_DATA = \ TrkrHitTruthAssocv1_Dict_rdict.pcm \ TrkrHit_Dict_rdict.pcm \ TrkrHitv1_Dict_rdict.pcm \ - TrkrHitv2_Dict_rdict.pcm + TrkrHitv2_Dict_rdict.pcm \ + TrkrHitv3_Dict_rdict.pcm # sources for io library libtrack_io_la_SOURCES = \ @@ -246,6 +252,7 @@ libtrack_io_la_SOURCES = \ TrkrClusterv3.cc \ TrkrClusterv4.cc \ TrkrClusterv5.cc \ + TrkrClusterv6.cc \ TrkrDefs.cc \ TrkrHitSet.cc \ TrkrHitSetContainer.cc \ @@ -253,7 +260,8 @@ libtrack_io_la_SOURCES = \ TrkrHitSetv1.cc \ TrkrHitTruthAssocv1.cc \ TrkrHitv1.cc \ - TrkrHitv2.cc + TrkrHitv2.cc \ + TrkrHitv3.cc libtrack_io_la_LIBADD = \ -lphool \ diff --git a/offline/packages/trackbase/TrkrCluster.h b/offline/packages/trackbase/TrkrCluster.h index 46edfc90ce..704c08f9f5 100644 --- a/offline/packages/trackbase/TrkrCluster.h +++ b/offline/packages/trackbase/TrkrCluster.h @@ -72,6 +72,11 @@ class TrkrCluster : public PHObject virtual float getTime() const { return NAN;} virtual char getSize() const {return std::numeric_limits::max(); } + virtual unsigned long long getParentID() const { return ULONG_LONG_MAX; } + virtual void setParentID(unsigned long long) {} + virtual double getPctParentID() const { return NAN; } + virtual void setPctParentID(double) {} + // // convenience interface diff --git a/offline/packages/trackbase/TrkrClusterv6.cc b/offline/packages/trackbase/TrkrClusterv6.cc new file mode 100644 index 0000000000..26e52adf9f --- /dev/null +++ b/offline/packages/trackbase/TrkrClusterv6.cc @@ -0,0 +1,83 @@ +/** + * @file trackbase/TrkrClusterv6.cc + * @author J. Osborn + * @date October 2021 + * @brief Implementation of TrkrClusterv6 + */ +#include "TrkrClusterv6.h" + +#include +#include // for swap + +namespace +{ + // square convenience function + template inline constexpr T square( const T& x ) + { return x*x; } +} + +TrkrClusterv6::TrkrClusterv6() + : m_subsurfkey(TrkrDefs::SUBSURFKEYMAX) + , m_phierr(0) + , m_zerr(0) + , m_adc(0) + , m_maxadc(0) + , m_phisize(0) + , m_zsize(0) + , m_overlap(0) + , m_edge(0) +{ + for (int i = 0; i < 2; i++) + { + m_local[i] = NAN; + } +} + +void TrkrClusterv6::identify(std::ostream& os) const +{ + os << "---TrkrClusterv6--------------------" << std::endl; + + os << " (rphi,z) = (" << getLocalX(); + os << ", " << getLocalY() << ") cm "; + + os << " valid = " << isValid() << std::endl; + + os << std::endl; + os << "-----------------------------------------------" << std::endl; + + return; +} + +int TrkrClusterv6::isValid() const +{ + for (int i = 0; i < 2; ++i) + { + if (std::isnan(getPosition(i))) { return 0; } + } + if (m_adc == 0xFFFF) { return 0; } + + return 1; +} + +void TrkrClusterv6::CopyFrom( const TrkrCluster& source ) +{ + // do nothing if copying onto oneself + if( this == &source ) return; + + // parent class method + TrkrCluster::CopyFrom( source ); + + setLocalX( source.getLocalX() ); + setLocalY( source.getLocalY() ); + setSubSurfKey( source.getSubSurfKey() ); + setAdc( source.getAdc() ); + setMaxAdc( source.getMaxAdc()); + setPhiError(source.getPhiError()); + setZError(source.getZError()); + setPhiSize(source.getPhiSize()); + setZSize(source.getZSize()); + setOverlap(source.getOverlap()); + setEdge(source.getEdge()); +} + + diff --git a/offline/packages/trackbase/TrkrClusterv6.h b/offline/packages/trackbase/TrkrClusterv6.h new file mode 100644 index 0000000000..5978b60ee9 --- /dev/null +++ b/offline/packages/trackbase/TrkrClusterv6.h @@ -0,0 +1,163 @@ +/** + * @file trackbase/TrkrClusterv6.h + * @author C. Roland + * @date March 2023 + * @brief Version 5 of TrkrCluster + */ +#ifndef TRACKBASE_TRKRCLUSTERV6_H +#define TRACKBASE_TRKRCLUSTERV6_H + +#include "TrkrCluster.h" +#include "TrkrDefs.h" +#include + +class PHObject; + +/** + * @brief Version 5 of TrkrCluster + * + * This version of TrkrCluster is blown up to contain a maximum of information + */ + +class TrkrClusterv6 : public TrkrCluster +{ + public: + + //! ctor + TrkrClusterv6(); + + //!dtor + ~TrkrClusterv6() override = default; + + // PHObject virtual overloads + + void identify(std::ostream& os = std::cout) const override; + void Reset() override {} + int isValid() const override; + PHObject* CloneMe() const override { return new TrkrClusterv6(*this); } + + //! import PHObject CopyFrom, in order to avoid clang warning + using PHObject::CopyFrom; + + //! copy content from base class + void CopyFrom( const TrkrCluster& ) override; + + //! copy content from base class + void CopyFrom( TrkrCluster* source ) override + { CopyFrom( *source ); } + + // + // cluster position + // + 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; } + float getLocalY() const override { return m_local[1]; } + void setLocalY(float loc1) override { m_local[1] = loc1; } + + TrkrDefs::subsurfkey getSubSurfKey() const override { return m_subsurfkey; } + void setSubSurfKey(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; + } + + unsigned int getMaxAdc() const override { + return m_maxadc; + } + + void setMaxAdc(uint16_t maxadc) override { + m_maxadc = maxadc; + } + + // + // convenience interface + // + 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;} + + /// deprecated global funtions with a warning + float getX() const override + { std::cout << "Deprecated getx trkrcluster function!"< + +TrkrHitv3::TrkrHitv3() +{ +} + + // these set and get the energy before digitization +void TrkrHitv3::addEnergy(const double edep) +{ + double max_adc = (double) USHRT_MAX; + + // overflow occurs if (sum of the existing + new ADC values) > USHRT_MAX + double ein = edep* TrkrDefs::EdepScaleFactor; + if( (double) m_adc + ein > max_adc) + { + m_adc = USHRT_MAX; + } + else + { + m_adc += (unsigned short) (ein); + } +} + +double TrkrHitv3::getEnergy() +{ + return ((double) m_adc) / TrkrDefs::EdepScaleFactor; +} + +void TrkrHitv3::setAdc(const unsigned int adc) + { + if(adc > USHRT_MAX) + m_adc = USHRT_MAX; + else + m_adc = (unsigned short) adc; + } + +unsigned int TrkrHitv3::getAdc() { + return (unsigned int) m_adc; + } + + +void TrkrHitv3::setParentID(const unsigned long long id){ + if(id > ULONG_LONG_MAX) + m_PHG4HitID = ULONG_LONG_MAX; + else + m_PHG4HitID = (unsigned long long) id; +} + +unsigned long long TrkrHitv3::getParentID() { + return (unsigned long long) m_PHG4HitID; +} diff --git a/offline/packages/trackbase/TrkrHitv3.h b/offline/packages/trackbase/TrkrHitv3.h new file mode 100644 index 0000000000..ecd766b031 --- /dev/null +++ b/offline/packages/trackbase/TrkrHitv3.h @@ -0,0 +1,59 @@ +/** + * @file trackbase/TrkrHit.h + * @author D. McGlinchey + * @date 4 June 2018 + * @brief Base class for hit object + */ +#ifndef TRACKBASE_TRKRHITV3_H +#define TRACKBASE_TRKRHITV3_H + +#include "TrkrDefs.h" +#include "TrkrHit.h" + +#include + +#include + +/** + * @brief Base class for hit object + * + * This is the empyt virtual base class for a hit object. + * Each subsystem should implement an inherited version + * which contains the actual storage information. + */ +class TrkrHitv3 : public TrkrHit +{ + public: + //! ctor + TrkrHitv3(); + + //! dtor + ~TrkrHitv3() override {} + // PHObject virtual overloads + void identify(std::ostream& os = std::cout) const override + { + os << "TrkrHitv3 class with adc = " << m_adc << std::endl; + } + void Reset() override {} + int isValid() const override { return 0; } + + // these set and get the energy before digitization + void addEnergy(const double edep) override; + double getEnergy() override; + + // after digitization, these are the adc values + void setAdc(const unsigned int adc) override; + unsigned int getAdc() override ; + + void setParentID(const unsigned long long id) override; + unsigned long long getParentID() override ; + + protected: + + unsigned short m_adc = 0; + unsigned long long m_PHG4HitID = 0; + + ClassDefOverride(TrkrHitv3, 1); +}; + +#endif //TRACKBASE_TRKRHITV2_H diff --git a/offline/packages/trackbase/TrkrHitv3LinkDef.h b/offline/packages/trackbase/TrkrHitv3LinkDef.h new file mode 100644 index 0000000000..38d2cf314b --- /dev/null +++ b/offline/packages/trackbase/TrkrHitv3LinkDef.h @@ -0,0 +1,5 @@ +#ifdef __CINT__ + +#pragma link C++ class TrkrHitv3+; + +#endif diff --git a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc index 55e8491f0a..ca30cd4d7d 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc @@ -99,7 +99,7 @@ int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */) * - adjust hit time and z, * - insert in container */ - auto adjust_hits = [&](PHG4Hit* source) { + auto adjust_hits = [&](PHG4Hitv1* source) { // adjust time to account for central membrane delay source->set_t(0, m_centralMembraneDelay); source->set_t(1, m_centralMembraneDelay); @@ -107,6 +107,8 @@ int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */) // assign to positive side source->set_z(0, 1.); source->set_z(1, 1.); + std::cout << "source hit id " << source->get_hit_id() << std::endl; + PHG4Hits.push_back(source); // clone @@ -114,6 +116,9 @@ int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */) auto copy = new PHG4Hitv1(source); copy->set_z(0, -1.); copy->set_z(1, -1.); + copy->set_hit_id( (PHG4HitDefs::keytype) (copy->get_hit_id() + (18*10000))); + std::cout << "copy hit id " << copy->get_hit_id() << std::endl; + PHG4Hits.push_back(copy); }; @@ -443,10 +448,10 @@ int PHG4TpcCentralMembrane::getSearchResult(double xcheck, double ycheck) const return result; } -PHG4Hit* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const +PHG4Hitv1* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const { //this function generates a PHG4 hit using coordinates from a stripe const double phi_petal = M_PI / 9.0; // angle span of one petal - PHG4Hit* hit; + PHG4Hitv1* hit; TVector3 dummyPos0, dummyPos1; //could put in some sanity checks here but probably not necessary since this is only really used within the class @@ -456,6 +461,7 @@ PHG4Hit* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, //from phg4tpcsteppingaction.cc hit = new PHG4Hitv1(); hit->set_layer(-1); // dummy number + hit->set_hit_id( (PHG4HitDefs::keytype) ((10000*petalID) + (100*radiusID) + stripeID)); //here we set the entrance values in cm if (moduleID == 0) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h index f8532284cb..af50cd14b6 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h @@ -11,6 +11,7 @@ class PHCompositeNode; class PHG4Hit; +class PHG4Hitv1; // all distances in mm, all angles in rad // class that generates stripes and dummy hit coordinates @@ -58,7 +59,7 @@ class PHG4TpcCentralMembrane : public SubsysReco, public PHParameterInterface /// g4hitnode name std::string hitnodename = "G4HIT_TPC"; - std::vector PHG4Hits; + std::vector PHG4Hits; int m_eventModulo = 10; int m_eventNum = 0; @@ -215,7 +216,7 @@ class PHG4TpcCentralMembrane : public SubsysReco, public PHParameterInterface const double x3b[][nRadii], const double y3b[][nRadii], double x, double y, const std::array& nGoodStripes) const; - PHG4Hit* GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const; + PHG4Hitv1* GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const; }; #endif diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index 29ef025a89..85d187dd7b 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -193,13 +193,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) seggeonodename = "CYLINDERCELLGEOM_SVTX"; // + detector; seggeo = findNode::getClass(topNode, seggeonodename); - if (!seggeo) - { - seggeo = new PHG4TpcCylinderGeomContainer(); - auto newNode = new PHIODataNode(seggeo, seggeonodename, "PHObject"); - runNode->addNode(newNode); - } - + assert( seggeo ); UpdateParametersWithMacro(); PHNodeIterator runIter(runNode); @@ -301,7 +295,6 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } padplane->InitRun(topNode); - padplane->CreateReadoutGeometry(topNode, seggeo); // print all layers radii if (Verbosity()) @@ -386,6 +379,9 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // clustering loopers in the same HitSetKey surfaces in multiple passes for (auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter) { + + std::cout << PHWHERE << " hit num: " << count_g4hits << " z: " << hiter->second->get_z(0) << " hitID: " << hiter->second->get_hit_id() << " hiter first " << hiter->first << std::endl; + count_g4hits++; dump_counter++; From c31683265fa996de31b26682a173d8133157a38b Mon Sep 17 00:00:00 2001 From: bkimelman Date: Mon, 3 Jul 2023 09:36:47 -0400 Subject: [PATCH 4/5] Changes to follow track id of hit on CM stripes through all reco to fluctuation distortion determination --- offline/packages/tpc/TpcClusterizer.cc | 56 +++++++- offline/packages/tpc/TpcClusterizer.h | 2 +- .../PHTpcCentralMembraneClusterizer.cc | 44 +++++- .../tpccalib/PHTpcCentralMembraneMatcher.cc | 39 ++++-- .../tpccalib/PHTpcCentralMembraneMatcher.h | 1 + offline/packages/trackbase/CMFlashCluster.h | 4 + .../packages/trackbase/CMFlashClusterv4.cc | 85 ++++++++++++ offline/packages/trackbase/CMFlashClusterv4.h | 127 ++++++++++++++++++ .../trackbase/CMFlashClusterv4LinkDef.h | 5 + offline/packages/trackbase/Makefile.am | 4 + offline/packages/trackbase/TrkrClusterv6.cc | 2 + .../g4tpc/PHG4TpcCentralMembrane.cc | 32 +++-- .../g4tpc/PHG4TpcCentralMembrane.h | 4 +- .../g4tpc/PHG4TpcElectronDrift.cc | 14 +- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 13 +- .../g4simulation/g4tpc/PHG4TpcSubsystem.cc | 4 + 16 files changed, 395 insertions(+), 41 deletions(-) create mode 100644 offline/packages/trackbase/CMFlashClusterv4.cc create mode 100644 offline/packages/trackbase/CMFlashClusterv4.h create mode 100644 offline/packages/trackbase/CMFlashClusterv4LinkDef.h diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index c93de17320..2854437e12 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include // for hitkey, getLayer #include @@ -63,6 +64,7 @@ namespace unsigned short it = 0; unsigned short adc = 0; unsigned short edge = 0; + unsigned long long parentID = -1; }; struct thread_data @@ -91,7 +93,7 @@ namespace unsigned short maxHalfSizePhi = 0; double m_tdriftmax = 0; double sampa_tbias = 0; - int cluster_version = 4; + int cluster_version = 6; std::vector association_vector; std::vector cluster_vector; int verbosity = 0; @@ -252,7 +254,7 @@ namespace return; } - 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, unsigned long long parentID, std::vector &ihit_list, int &touch, int &edge) { // search along phi at the peak in t @@ -271,6 +273,7 @@ namespace hit.iphi = iphi; hit.it = it; hit.adc = adcval[iphi][it]; + hit.parentID = parentID; if(touch>0){ if((iphi == (phibin - phidown))|| (iphi == (phibin + phiup))){ @@ -311,6 +314,8 @@ namespace int clus_size = ihit_list.size(); int max_adc = 0; if(clus_size == 1) return; + + std::vector parentIDs; // std::cout << "process list" << std::endl; std::vector hitkeyvec; for(auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){ @@ -325,6 +330,8 @@ namespace if(iphi < phibinlo) phibinlo = iphi; if(it > tbinhi) tbinhi = it; if(it < tbinlo) tbinlo = it; + + parentIDs.push_back(iter->parentID); // update phi sums // double phi_center = my_data.layergeom->get_phicenter(iphi); @@ -350,8 +357,27 @@ namespace // std::cout << "done process list" << std::endl; if (adc_sum < 10){ hitkeyvec.clear(); + parentIDs.clear(); return; // skip obvious noise "clusters" } + + unsigned long long commonParentID = -1; + int maxNpIDs = 0; + for(int pID=0; pID<(int)parentIDs.size(); pID++){ + + int nPIDs = 0; + for(int pID2=0; pID2<(int)parentIDs.size(); pID2++){ + if(pID == pID2) continue; + if(parentIDs[pID] == parentIDs[pID2]) nPIDs++; + } + if(nPIDs > maxNpIDs){ + maxNpIDs = nPIDs; + commonParentID = parentIDs[pID]; + } + } + double pctParentID = 1.0*maxNpIDs/parentIDs.size(); + + // This is the global position double clusiphi = iphi_sum / adc_sum; double clusphi = my_data.layergeom->get_phi(clusiphi); @@ -469,6 +495,26 @@ namespace clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(),2))); my_data.cluster_vector.push_back(clus); } + }else if(my_data.cluster_version==6){ + //std::cout << "ver5" << std::endl; + // std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl; + if(sqrt(phi_err_square) > 0.01){ + auto clus = new TrkrClusterv6; + //auto clus = std::make_unique(); + clus->setAdc(adc_sum); + clus->setMaxAdc(max_adc); + clus->setEdge(nedge); + clus->setPhiSize(phisize); + clus->setZSize(tsize); + clus->setSubSurfKey(subsurfkey); + clus->setLocalX(local(0)); + clus->setLocalY(clust); + clus->setPhiError(sqrt(phi_err_square)); + clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(),2))); + clus->setParentID(commonParentID); + clus->setPctParentID(pctParentID); + my_data.cluster_vector.push_back(clus); + } } //std::cout << "end clus out" << std::endl; @@ -521,7 +567,7 @@ namespace for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr){ - + if( TpcDefs::getPad(hitr->first) - phioffset < 0 ){ //std::cout << "WARNING phibin out of range: " << TpcDefs::getPad(hitr->first) - phioffset << " | " << phibins << std::endl; continue; @@ -556,6 +602,7 @@ namespace thisHit.it = tbin; thisHit.adc = adc; thisHit.edge = 0; + thisHit.parentID = hitr->second->getParentID(); all_hit_map.insert(std::make_pair(adc, thisHit)); } if(adc>my_data->threshold){ @@ -639,13 +686,14 @@ namespace ihit hiHit = iter->second; int iphi = hiHit.iphi; int it = hiHit.it; + unsigned long long iparentID = hiHit.parentID; //put all hits in the all_hit_map (sorted by adc) //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 ); + get_cluster(iphi, it, *my_data, adcval, iparentID, ihit_list, ntouch, nedge ); // -> calculate cluster parameters // -> add hits to truth association diff --git a/offline/packages/tpc/TpcClusterizer.h b/offline/packages/tpc/TpcClusterizer.h index b277b67275..005049a983 100644 --- a/offline/packages/tpc/TpcClusterizer.h +++ b/offline/packages/tpc/TpcClusterizer.h @@ -64,7 +64,7 @@ class TpcClusterizer : public SubsysReco double SectorFiducialCut = 0.5; unsigned short MaxClusterHalfSizePhi = 3; unsigned short MaxClusterHalfSizeT = 5; - int cluster_version = 4; + int cluster_version = 6; double m_tdriftmax = 0; double AdcClockPeriod = 53.0; // ns diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc index f8eac792a6..a18f09cc32 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc @@ -21,7 +21,8 @@ #include #include #include -#include +#include +#include #include #include #include @@ -110,11 +111,15 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) std::vectori_pair; //vector for pair matching std::vectorenergy;//vector for energy values of clusters std::vectorisAcrossGap; + std::vectorparentIDs; + std::vectorpctParentIDs; int nTpcClust = 0; double mean_z_content_plus = 0.0; double mean_z_content_minus = 0.0; + int clusCount = 0; + //first loop over clusters to make mod phi histograms of each layer and each pair of layers for(const auto& hitsetkey:_cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId)) { @@ -124,6 +129,12 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) { const auto& [cluskey, cluster] = *clusiter; + + // TrkrClusterv6 *clusterv6 = new TrkrClusterv6(); + // clusterv6->CopyFrom(clusiter->second); + + //std::cout << "cluster " << clusCount << " parentID: " << cluster->getParentID() << " pctParentID: " << cluster->getPctParentID() << std::endl; + auto glob = tgeometry->getGlobalPosition(cluskey, cluster); TVector3 tmp_pos(glob(0),glob(1),glob(2)); @@ -162,9 +173,9 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) if(layer < 54) hphi_reco_pair_neg[layer-7]->Fill(phiMod); if(layer > 7) hphi_reco_pair_neg[layer-8]->Fill(phiMod); } - - } - } + clusCount++; + }//end loop over clusters + }//end loop over hitsets for(int i=1; iGetNbinsX(); i++){ mean_z_content_plus += hz_pos->GetBinContent(i); @@ -256,7 +267,12 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) //require that z be within 1cm of peak in z distributions (separate for each side) if( (z>=0 && std::abs(z - hz_pos->GetBinCenter(hz_pos->GetMaximumBin())) > 1.0) || (z<0 && std::abs(z - hz_neg->GetBinCenter(hz_neg->GetMaximumBin())) > 1.0) ) continue; + //std::cout << "accepted cluster " << m_accepted_clusters << " parentID: " << cluster->getParentID() << " pctParentID: " << cluster->getPctParentID() << std::endl; + + ++m_accepted_clusters; + + i_pair.push_back(-1); energy.push_back(cluster->getAdc()); @@ -265,6 +281,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) layer.push_back((int)(TrkrDefs::getLayer(cluskey))); side.push_back(TpcDefs::getSide(cluskey)); isAcrossGap.push_back(false); + parentIDs.push_back(cluster->getParentID()); + pctParentIDs.push_back(cluster->getPctParentID()); if(Verbosity() > 0) std::cout << ":\t" << x << "\t" << y << "\t" << z < nclusters; std::vector isREdge; std::vector > pairNums; + std::vector parentID; + std::vector pctParentID; // int nR2 = 0; //int nR3 = 0; @@ -498,7 +518,15 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) nclusters.push_back(2); if(isAcrossGap[i] && isAcrossGap[i_pair[i]]) isREdge.push_back(true); else isREdge.push_back(false); - + + if(energy[i] > energy[i_pair[i]]){ + parentID.push_back(parentIDs[i]); + pctParentID.push_back(energy[i]*pctParentIDs[i]/(energy[i]+energy[i_pair[i]])); + }else{ + parentID.push_back(parentIDs[i_pair[i]]); + pctParentID.push_back(energy[i_pair[i]]*pctParentIDs[i_pair[i]]/(energy[i]+energy[i_pair[i]])); + } + tmp_pair.first = i; tmp_pair.second = i_pair[i]; pairNums.push_back(tmp_pair); @@ -524,6 +552,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) avepos.push_back(pos[i]); nclusters.push_back(1); isREdge.push_back(isAcrossGap[i]); + parentID.push_back(parentIDs[i]); + pctParentID.push_back(pctParentIDs[i]); tmp_pair.first = i; tmp_pair.second = -1; pairNums.push_back(tmp_pair); @@ -537,7 +567,7 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) for(unsigned int iv = 0; iv setX(avepos[iv].X()); cmfc->setY(avepos[iv].Y()); @@ -545,6 +575,8 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) cmfc->setAdc(aveenergy[iv]); cmfc->setNclusters(nclusters[iv]); cmfc->setIsRGap(isREdge[iv]); + cmfc->setParentID(parentID[iv]); + cmfc->setPctParentID(pctParentID[iv]); int pair1 = pairNums[iv].first; int pair2 = pairNums[iv].second; diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc index 5f24033fc6..af897aa2e0 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc @@ -11,7 +11,7 @@ #include #include #include -#include +#include #include #include #include @@ -267,7 +267,7 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) hnclus = new TH1F("hnclus", " nclusters ", 3, 0., 3.); m_debugfile.reset ( new TFile(m_debugfilename.c_str(),"RECREATE") ); - match_ntup = new TNtuple("match_ntup","Match NTuple","event:truthR:truthPhi:recoR:recoPhi:recoZ:nclus:r1:phi1:e1:layer1:r2:phi2:e2:layer2"); + match_ntup = new TNtuple("match_ntup","Match NTuple","event:truthR:truthPhi:recoR:recoPhi:recoZ:correctMatch:nclus:r1:phi1:e1:r2:phi2:e2:"); } hit_r_phi = new TH2F("hit_r_phi","hit r vs #phi;#phi (rad); r (cm)",360,-M_PI,M_PI,500,0,100); @@ -288,15 +288,17 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) * - assign proper z, * - insert in container */ - auto save_truth_position = [&](TVector3 source) + auto save_truth_position = [&](TVector3 source, unsigned long long parentID) { source.SetZ( +1 ); m_truth_pos.push_back( source ); + m_truth_parentID.push_back(parentID); hit_r_phi->Fill(source.Phi(), source.Perp()); source.SetZ( -1 ); m_truth_pos.push_back( source ); + m_truth_parentID.push_back(parentID + 180000); hit_r_phi->Fill(source.Phi(), source.Perp()); }; @@ -308,7 +310,7 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) { TVector3 dummyPos(cx1_e[i][j], cy1_e[i][j], 0.0); dummyPos.RotateZ(k * phi_petal); - save_truth_position(dummyPos); + save_truth_position(dummyPos, (k*10000) + (j*100) + i ); if(Verbosity() > 2) std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y() @@ -324,7 +326,7 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) { TVector3 dummyPos(cx1[i][j], cy1[i][j], 0.0); dummyPos.RotateZ(k * phi_petal); - save_truth_position(dummyPos); + save_truth_position(dummyPos, (k*10000) + ((j+8)*100) + i ); if(Verbosity() > 2) std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y() @@ -339,7 +341,7 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) { TVector3 dummyPos(cx2[i][j], cy2[i][j], 0.0); dummyPos.RotateZ(k * phi_petal); - save_truth_position(dummyPos); + save_truth_position(dummyPos, (k*10000) + ((j+16)*100) + i ); if(Verbosity() > 2) std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y() @@ -354,7 +356,8 @@ int PHTpcCentralMembraneMatcher::InitRun(PHCompositeNode *topNode) { TVector3 dummyPos(cx3[i][j], cy3[i][j], 0.0); dummyPos.RotateZ(k * phi_petal); - save_truth_position(dummyPos); + save_truth_position(dummyPos, (k*10000) + ((j+24)*100) + i ); + if(Verbosity() > 2) std::cout << " i " << i << " j " << j << " k " << k << " x1 " << dummyPos.X() << " y1 " << dummyPos.Y() @@ -379,6 +382,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) std::vector adc2; std::vector layer1; std::vector layer2; + std::vector reco_parentID; + std::vector reco_pctParentID; // reset output distortion correction container histograms for( const auto& harray:{m_dcc_out->m_hDRint, m_dcc_out->m_hDPint, m_dcc_out->m_hDZint, m_dcc_out->m_hentries} ) @@ -444,6 +449,11 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) layer1.push_back(cmclus->getLayer1()); layer2.push_back(cmclus->getLayer2()); + // std::cout << "parent ID: " << cmclus->getParentID() << " pctParentID: " << cmclus->getPctParentID() << std::endl; + reco_parentID.push_back(cmclus->getParentID()); + reco_pctParentID.push_back(cmclus->getPctParentID()); + + if(tmp_pos.Z() < 0) clust_r_phi_neg->Fill(tmp_pos.Phi(),tmp_pos.Perp()); else clust_r_phi_pos->Fill(tmp_pos.Phi(),tmp_pos.Perp()); @@ -678,6 +688,10 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) const auto n_valid_truth = std::count_if( m_truth_pos.begin(), m_truth_pos.end(), []( const TVector3& pos ) { return get_r( pos.x(), pos.y() ) > 30; } ); const auto n_reco_size1 = std::count_if( reco_nclusters.begin(), reco_nclusters.end(), []( const unsigned int& value ) { return value==1; } ); const auto n_reco_size2 = std::count_if( reco_nclusters.begin(), reco_nclusters.end(), []( const unsigned int& value ) { return value==2; } ); + int n_correct_match = 0; + for(int pairs=0; pairs<(int)matched_pair.size(); pairs++){ + if( m_truth_parentID[matched_pair[pairs].first] == reco_parentID[matched_pair[pairs].second]) n_correct_match++; + } std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size: " << m_truth_pos.size() << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - m_truth_pos size, r>30cm: " << n_valid_truth << std::endl; int pos_stripes_add = 0; @@ -694,6 +708,10 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==1): " << n_reco_size1 << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - reco_pos size (nclus==2): " << n_reco_size2 << std::endl; std::cout << "PHTpcCentralMembraneMatcher::process_event - matched_pair size: " << matched_pair.size() << std::endl; + std::cout << "PHTpcCentralMembraneMatcher::process_event - number of correct matches: " << n_correct_match << std::endl; + std::cout << "PHTpcCentralMembraneMatcher::process_event - number of correct matches over number of matches: " << 1.0*n_correct_match/matched_pair.size() << std::endl; + std::cout << "PHTpcCentralMembraneMatcher::process_event - number of correct matches over number of truth: " << 1.0*n_correct_match/m_truth_pos.size() << std::endl; + std::cout << "PHTpcCentralMembraneMatcher::process_event - number of correct matches over number of truth r>30cm: " << 1.0*n_correct_match/n_valid_truth << std::endl; } for(unsigned int ip = 0; ip < matched_pair.size(); ++ip) @@ -716,7 +734,12 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) m_cm_flash_diffs->addDifferenceSpecifyKey(key, cmdiff); - if(m_savehistograms) match_ntup->Fill(m_event_index,m_truth_pos[p.first].Perp(),m_truth_pos[p.first].Phi(),reco_pos[p.second].Perp(),reco_pos[p.second].Phi(),reco_pos[p.second].Z(),nclus,pos1[p.second].Perp(),pos1[p.second].Phi(),adc1[p.second],layer1[p.second],pos2[p.second].Perp(),pos2[p.second].Phi(),adc2[p.second],layer2[p.second]); + + double isCorrectMatch = 0.0; + + if( m_truth_parentID[p.first] == reco_parentID[p.second]) isCorrectMatch = 1.0; + + if(m_savehistograms) match_ntup->Fill(m_event_index,m_truth_pos[p.first].Perp(),m_truth_pos[p.first].Phi(),reco_pos[p.second].Perp(),reco_pos[p.second].Phi(),reco_pos[p.second].Z(),isCorrectMatch,nclus,pos1[p.second].Perp(),pos1[p.second].Phi(),adc1[p.second],pos2[p.second].Perp(),pos2[p.second].Phi(),adc2[p.second]); // store cluster position const double clus_r = reco_pos[p.second].Perp(); diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h index 7f7e7f7ab5..9d205bab56 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.h @@ -228,6 +228,7 @@ class PHTpcCentralMembraneMatcher : public SubsysReco /// store centers of all central membrane pads std::vector m_truth_pos; + std::vector m_truth_parentID; //@} diff --git a/offline/packages/trackbase/CMFlashCluster.h b/offline/packages/trackbase/CMFlashCluster.h index 0de1ed2bc3..b0958a290c 100644 --- a/offline/packages/trackbase/CMFlashCluster.h +++ b/offline/packages/trackbase/CMFlashCluster.h @@ -94,6 +94,10 @@ class CMFlashCluster : public PHObject virtual void setIsPhiGap(bool) {} virtual bool getIsPhiGap() const { return false; } + virtual unsigned long long getParentID() const {return ULONG_LONG_MAX;} + virtual void setParentID(unsigned long long) {} + virtual double getPctParentID() const {return NAN;} + virtual void setPctParentID(double) {} protected: CMFlashCluster() = default; diff --git a/offline/packages/trackbase/CMFlashClusterv4.cc b/offline/packages/trackbase/CMFlashClusterv4.cc new file mode 100644 index 0000000000..4eb7f9fba3 --- /dev/null +++ b/offline/packages/trackbase/CMFlashClusterv4.cc @@ -0,0 +1,85 @@ +/** + * @file trackbase/CMFlashClusterv3.cc + * @author Ben Kimelman + * @date March 2023 + * @brief Implementation of CMFlashClusterv4 + */ +#include "CMFlashClusterv4.h" + +#include +#include // for swap + +void CMFlashClusterv4::identify(std::ostream& os) const +{ + os << "---CMFlashClusterv4--------------------" << std::endl; + + os << " (x,y,z) = (" << m_pos[0]; + os << ", " << m_pos[1] << ", "; + os << m_pos[2] << ") cm"; + + os << " adc = " << getAdc() << std::endl; + + os << std::endl; + os << "-----------------------------------------------" << std::endl; + + return; +} + +int CMFlashClusterv4::isValid() const +{ + if(std::isnan(getX())) return 0; + if(std::isnan(getY())) return 0; + if(std::isnan(getZ())) return 0; + + + if(std::isnan(getX1())) return 0; + if(std::isnan(getY1())) return 0; + if(std::isnan(getZ1())) return 0; + + + if(std::isnan(getX2())) return 0; + if(std::isnan(getY2())) return 0; + if(std::isnan(getZ2())) return 0; + + if (m_adc == 0xFFFFFFFF) return 0; + if (m_adc1 == 0xFFFFFFFF) return 0; + if (m_adc2 == 0xFFFFFFFF) return 0; + + return 1; +} + +void CMFlashClusterv4::CopyFrom( const CMFlashCluster& source ) +{ + // do nothing if copying onto oneself + if( this == &source ) return; + + // parent class method + CMFlashCluster::CopyFrom( source ); + + setX( source.getX() ); + setY( source.getY() ); + setZ( source.getZ() ); + + + setX1( source.getX1() ); + setY1( source.getY1() ); + setZ1( source.getZ1() ); + + setX2( source.getX2() ); + setY2( source.getY2() ); + setZ2( source.getZ2() ); + + setLayer1( source.getLayer1() ); + setLayer2( source.getLayer2() ); + + setAdc( source.getAdc() ); + setAdc1( source.getAdc1() ); + setAdc2( source.getAdc2() ); + setIsRGap( source.getIsRGap() ); + setIsPhiGap( source.getIsPhiGap() ); + + setParentID( source.getParentID() ); + setPctParentID( source.getPctParentID() ); + +} + diff --git a/offline/packages/trackbase/CMFlashClusterv4.h b/offline/packages/trackbase/CMFlashClusterv4.h new file mode 100644 index 0000000000..792010a64c --- /dev/null +++ b/offline/packages/trackbase/CMFlashClusterv4.h @@ -0,0 +1,127 @@ +/** + * @file trackbase/CMFlashClusterv3.h + * @author Ben Kimelman + * @date March 2023 + * @brief Version 3 of CMFlashCluster + */ +#ifndef TRACKBASE_CMFLASHCLUSTERV4_H +#define TRACKBASE_CMFLASHCLUSTERV4_H + +#include "CMFlashCluster.h" + +#include + +class PHObject; + +/** + * @brief Version 3 of CMFlashCluster + * + *Adding variable to keep track of clusters + *put into metaclusters + * + */ +class CMFlashClusterv4 : public CMFlashCluster +{ + public: + //! ctor + CMFlashClusterv4() = default; + + // PHObject virtual overloads + void identify(std::ostream& os = std::cout) const override; + void Reset() override {} + int isValid() const override; + PHObject* CloneMe() const override { return new CMFlashClusterv4(*this); } + + //! copy content from base class + void CopyFrom( const CMFlashCluster& ) override; + + //! copy content from base class + void CopyFrom( CMFlashCluster* source ) override + { CopyFrom( *source ); } + + // + // cluster position + // + float getX() const override { return m_pos[0]; } + void setX(float x) override { m_pos[0] = x; } + float getY() const override { return m_pos[1]; } + void setY(float y) override { m_pos[1] = y; } + float getZ() const override { return m_pos[2]; } + void setZ(float z) override { m_pos[2] = z; } + + + float getX1() const override { return m_pos1[0]; } + void setX1(float x) override { m_pos1[0] = x; } + float getY1() const override { return m_pos1[1]; } + void setY1(float y) override { m_pos1[1] = y; } + float getZ1() const override { return m_pos1[2]; } + void setZ1(float z) override { m_pos1[2] = z; } + + float getX2() const override { return m_pos2[0]; } + void setX2(float x) override { m_pos2[0] = x; } + float getY2() const override { return m_pos2[1]; } + void setY2(float y) override { m_pos2[1] = y; } + float getZ2() const override { return m_pos2[2]; } + void setZ2(float z) override { m_pos2[2] = z; } + + unsigned int getLayer1() const override {return m_layer1;} + void setLayer1(unsigned int layer) override { m_layer1 = layer;} + unsigned int getLayer2() const override {return m_layer2;} + void setLayer2(unsigned int layer) override { m_layer2 = layer;} + + unsigned int getNclusters() const override {return m_nclusters;} + void setNclusters(unsigned int n) override { m_nclusters = n;} + bool getIsRGap() const override { return m_isRGap; } + void setIsRGap(bool isRGap) override { m_isRGap = isRGap;} + bool getIsPhiGap() const override { return m_isPhiGap; } + void setIsPhiGap(bool isPhiGap) override { m_isPhiGap = isPhiGap;} + + // + // cluster info + // + unsigned int getAdc() const override { return m_adc; } + void setAdc(unsigned int adc) override { m_adc = adc; } + + unsigned int getAdc1() const override { return m_adc1; } + void setAdc1(unsigned int adc) override { m_adc1 = adc; } + + unsigned int getAdc2() const override { return m_adc2; } + void setAdc2(unsigned int adc) override { m_adc2 = adc; } + + unsigned long long getParentID() const override { return m_parentID; } + void setParentID(unsigned long long ID) override { m_parentID = ID; } + + double getPctParentID() const override { return m_pctParentID; } + void setPctParentID(double pct) override { m_pctParentID = pct; } + + protected: + + /// mean cluster position + float m_pos[3] = {NAN, NAN, NAN}; + + float m_pos1[3] = {NAN, NAN, NAN}; + float m_pos2[3] = {NAN, NAN, NAN}; + + /// cluster sum adc + unsigned int m_adc = 0xFFFFFFFF; + unsigned int m_adc1 = 0xFFFFFFFF; + unsigned int m_adc2 = 0xFFFFFFFF; + + unsigned int m_layer1 = UINT_MAX; + unsigned int m_layer2 = UINT_MAX; + + /// number of TPC clusters used to create this central mebrane cluster + unsigned int m_nclusters = UINT_MAX; + + /// bools to identify if meta-cluster is across sector/module gaps + bool m_isRGap = false; + bool m_isPhiGap = false; + + //parent id + unsigned long long m_parentID = ULONG_LONG_MAX; + double m_pctParentID = NAN; + + ClassDefOverride(CMFlashClusterv4, 1) +}; + +#endif //TRACKBASE_CMFLASHCLUSTERV4_H diff --git a/offline/packages/trackbase/CMFlashClusterv4LinkDef.h b/offline/packages/trackbase/CMFlashClusterv4LinkDef.h new file mode 100644 index 0000000000..2192a0eeea --- /dev/null +++ b/offline/packages/trackbase/CMFlashClusterv4LinkDef.h @@ -0,0 +1,5 @@ +#ifdef __CINT__ + +#pragma link C++ class CMFlashClusterv4+; + +#endif diff --git a/offline/packages/trackbase/Makefile.am b/offline/packages/trackbase/Makefile.am index c5a5abf1f7..08fdd62602 100644 --- a/offline/packages/trackbase/Makefile.am +++ b/offline/packages/trackbase/Makefile.am @@ -44,6 +44,7 @@ pkginclude_HEADERS = \ CMFlashClusterv1.h \ CMFlashClusterv2.h \ CMFlashClusterv3.h \ + CMFlashClusterv4.h \ CMFlashDifference.h \ CMFlashDifferenceContainer.h \ CMFlashDifferenceContainerv1.h \ @@ -107,6 +108,7 @@ ROOTDICTS = \ CMFlashClusterv1_Dict.cc \ CMFlashClusterv2_Dict.cc \ CMFlashClusterv3_Dict.cc \ + CMFlashClusterv4_Dict.cc \ CMFlashDifferenceContainer_Dict.cc \ CMFlashDifferenceContainerv1_Dict.cc \ CMFlashDifference_Dict.cc \ @@ -161,6 +163,7 @@ nobase_dist_pcm_DATA = \ CMFlashClusterv1_Dict_rdict.pcm \ CMFlashClusterv2_Dict_rdict.pcm \ CMFlashClusterv3_Dict_rdict.pcm \ + CMFlashClusterv4_Dict_rdict.pcm \ CMFlashDifferenceContainer_Dict_rdict.pcm \ CMFlashDifferenceContainerv1_Dict_rdict.pcm \ CMFlashDifference_Dict_rdict.pcm \ @@ -216,6 +219,7 @@ libtrack_io_la_SOURCES = \ CMFlashClusterv1.cc \ CMFlashClusterv2.cc \ CMFlashClusterv3.cc \ + CMFlashClusterv4.cc \ CMFlashDifferenceContainerv1.cc \ CMFlashDifferencev1.cc \ ClusterErrorPara.cc \ diff --git a/offline/packages/trackbase/TrkrClusterv6.cc b/offline/packages/trackbase/TrkrClusterv6.cc index 26e52adf9f..ab4d7561aa 100644 --- a/offline/packages/trackbase/TrkrClusterv6.cc +++ b/offline/packages/trackbase/TrkrClusterv6.cc @@ -78,6 +78,8 @@ void TrkrClusterv6::CopyFrom( const TrkrCluster& source ) setZSize(source.getZSize()); setOverlap(source.getOverlap()); setEdge(source.getEdge()); + setParentID(source.getParentID()); + setPctParentID(source.getPctParentID()); } diff --git a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc index ca30cd4d7d..a3cb9eca64 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.cc @@ -99,7 +99,7 @@ int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */) * - adjust hit time and z, * - insert in container */ - auto adjust_hits = [&](PHG4Hitv1* source) { + auto adjust_hits = [&](PHG4Hit* source) { // adjust time to account for central membrane delay source->set_t(0, m_centralMembraneDelay); source->set_t(1, m_centralMembraneDelay); @@ -107,17 +107,17 @@ int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */) // assign to positive side source->set_z(0, 1.); source->set_z(1, 1.); - std::cout << "source hit id " << source->get_hit_id() << std::endl; + //std::cout << "source hit id " << source->get_hit_id() << std::endl; PHG4Hits.push_back(source); // clone // assign to negative side and insert in list - auto copy = new PHG4Hitv1(source); + PHG4Hit *copy = new PHG4Hitv1(source); copy->set_z(0, -1.); copy->set_z(1, -1.); copy->set_hit_id( (PHG4HitDefs::keytype) (copy->get_hit_id() + (18*10000))); - std::cout << "copy hit id " << copy->get_hit_id() << std::endl; + //std::cout << "copy hit id " << copy->get_hit_id() << std::endl; PHG4Hits.push_back(copy); }; @@ -181,10 +181,24 @@ int PHG4TpcCentralMembrane::process_event(PHCompositeNode* topNode) // copy all hits from G4hits vector into container for (const auto& hit : PHG4Hits) { - auto copy = new PHG4Hitv1(hit); - g4hitcontainer->AddHit(detId, copy); + PHG4Hit *copy = new PHG4Hitv1(hit); + // g4hitcontainer->AddHit(detId, copy); + g4hitcontainer->AddHit(copy); } + + /* + PHG4HitContainer::ConstRange hit_begin_end = g4hitcontainer->getHits(); + int hitCount = 0; + for (auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter){ + + std::cout << PHWHERE << " hit num: " << hitCount << " z: " << hiter->second->get_z(0) << " hitID: " << hiter->second->get_hit_id() << " hiter first " << hiter->first << std::endl; + + hitCount++; + + } + */ + m_eventNum++; return Fun4AllReturnCodes::EVENT_OK; @@ -448,10 +462,10 @@ int PHG4TpcCentralMembrane::getSearchResult(double xcheck, double ycheck) const return result; } -PHG4Hitv1* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const +PHG4Hit* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const { //this function generates a PHG4 hit using coordinates from a stripe const double phi_petal = M_PI / 9.0; // angle span of one petal - PHG4Hitv1* hit; + PHG4Hit* hit; TVector3 dummyPos0, dummyPos1; //could put in some sanity checks here but probably not necessary since this is only really used within the class @@ -461,7 +475,7 @@ PHG4Hitv1* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleI //from phg4tpcsteppingaction.cc hit = new PHG4Hitv1(); hit->set_layer(-1); // dummy number - hit->set_hit_id( (PHG4HitDefs::keytype) ((10000*petalID) + (100*radiusID) + stripeID)); + hit->set_hit_id( (PHG4HitDefs::keytype) ((10000*petalID) + (100*(radiusID+(8*moduleID))) + stripeID)); //here we set the entrance values in cm if (moduleID == 0) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h index af50cd14b6..9cf7b4a7a0 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcCentralMembrane.h @@ -59,7 +59,7 @@ class PHG4TpcCentralMembrane : public SubsysReco, public PHParameterInterface /// g4hitnode name std::string hitnodename = "G4HIT_TPC"; - std::vector PHG4Hits; + std::vector PHG4Hits; int m_eventModulo = 10; int m_eventNum = 0; @@ -216,7 +216,7 @@ class PHG4TpcCentralMembrane : public SubsysReco, public PHParameterInterface const double x3b[][nRadii], const double y3b[][nRadii], double x, double y, const std::array& nGoodStripes) const; - PHG4Hitv1* GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const; + PHG4Hit* GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const; }; #endif diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index 85d187dd7b..e06d356827 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -6,6 +6,7 @@ #include "PHG4TpcPadPlane.h" // for PHG4TpcPadPlane #include +#include #include #include #include @@ -16,7 +17,7 @@ #include #include // for TrkrHitTruthA... #include -#include +#include #include #include @@ -380,7 +381,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) for (auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter) { - std::cout << PHWHERE << " hit num: " << count_g4hits << " z: " << hiter->second->get_z(0) << " hitID: " << hiter->second->get_hit_id() << " hiter first " << hiter->first << std::endl; + // std::cout << PHWHERE << " hit num: " << count_g4hits << " z: " << hiter->second->get_z(0) << " hitID: " << hiter->second->get_hit_id() << " hiter first " << hiter->first << std::endl; count_g4hits++; dump_counter++; @@ -591,7 +592,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) const int sector = TpcDefs::getSectorId(node_hitsetkey); const int side = TpcDefs::getSide(node_hitsetkey); - if (Verbosity() > 8) + if (Verbosity() > 100) std::cout << " hitsetkey " << node_hitsetkey << " layer " << layer << " sector " << sector << " side " << side << std::endl; // get all of the hits from the single hitset TrkrHitSet::ConstRange single_hit_range = single_hitset_iter->second->getHits(); @@ -658,12 +659,13 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) if (!node_hit) { // Otherwise, create a new one - node_hit = new TrkrHitv2(); + node_hit = new TrkrHitv3(); node_hitsetit->second->addHitSpecificKey(temp_hitkey, node_hit); } // Either way, add the energy to it node_hit->addEnergy(temp_tpchit->getEnergy()); + node_hit->setParentID(temp_tpchit->getParentID()); } // end loop over temp hits @@ -700,7 +702,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer TrkrDefs::hitsetkey hitsetkey = hitset_iter->first; const unsigned int layer = TrkrDefs::getLayer(hitsetkey); - if (layer != print_layer) continue; + //if (layer != print_layer) continue; const int sector = TpcDefs::getSectorId(hitsetkey); const int side = TpcDefs::getSide(hitsetkey); @@ -716,7 +718,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) TrkrDefs::hitkey hitkey = hit_iter->first; TrkrHit *tpchit = hit_iter->second; std::cout << " hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey) << " z bin " << TpcDefs::getTBin(hitkey) - << " energy " << tpchit->getEnergy() << std::endl; + << " energy " << tpchit->getEnergy() << " parentID " << tpchit->getParentID() << std::endl; } } } diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 6cd28bb40a..d41f9bdda2 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -17,7 +17,7 @@ #include // for TrkrHit #include #include -#include // for TrkrHit +#include // for TrkrHit #include #include @@ -178,9 +178,9 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( layernum = LayerGeom->get_layer(); /* pass_data.layerGeom = LayerGeom; */ /* pass_data.layer = layernum; */ - if (Verbosity() > 1000) + if (Verbosity() > 20) std::cout << " g4hit id " << hiter->first << " rad_gem " << rad_gem << " rad_low " << rad_low << " rad_high " << rad_high - << " layer " << hiter->second->get_layer() << " want to change to " << layernum << std::endl; + << " layer " << hiter->second->get_layer() << " want to change to " << layernum << std::endl; hiter->second->set_layer(layernum); // have to set here, since the stepping action knows nothing about layers } @@ -333,11 +333,13 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( if (!hit) { // create a new one - hit = new TrkrHitv2(); + hit = new TrkrHitv3(); hitsetit->second->addHitSpecificKey(hitkey, hit); } // Either way, add the energy to it -- adc values will be added at digitization hit->addEnergy(neffelectrons); + // std::cout << "parent id first " << hiter->first << " parent id get " << hiter->second->get_hit_id() << std::endl; + hit->setParentID(hiter->first); tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons); @@ -348,11 +350,12 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( if (!single_hit) { // create a new one - single_hit = new TrkrHitv2(); + single_hit = new TrkrHitv3(); single_hitsetit->second->addHitSpecificKey(hitkey, single_hit); } // Either way, add the energy to it -- adc values will be added at digitization single_hit->addEnergy(neffelectrons); + single_hit->setParentID(hiter->first); /* if (Verbosity() > 0) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc index 5cc669efce..4f8b0bb688 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc @@ -39,13 +39,17 @@ PHG4TpcSubsystem::~PHG4TpcSubsystem() //_______________________________________________________________________ int PHG4TpcSubsystem::InitRunSubsystem(PHCompositeNode *topNode) { + + std::cout << PHWHERE << " starting init" << std::endl; PHNodeIterator iter(topNode); PHCompositeNode *dstNode = dynamic_cast(iter.findFirst("PHCompositeNode", "DST")); // create display settings before detector (detector adds its volumes to it) m_DisplayAction = new PHG4TpcDisplayAction(Name()); // create detector + std::cout << "making detector" << std::endl; m_Detector = new PHG4TpcDetector(this, topNode, GetParams(), Name()); + std::cout << "made detector" << std::endl; m_Detector->SuperDetector(SuperDetector()); m_Detector->OverlapCheck(CheckOverlap()); std::set nodes; From f303f107edcc2efefe7b925ddafad6b24a27ea32 Mon Sep 17 00:00:00 2001 From: bkimelman Date: Mon, 3 Jul 2023 17:43:16 -0400 Subject: [PATCH 5/5] Extra cout statements for debugging when reading prdf's --- .../PHTpcCentralMembraneClusterizer.cc | 31 ++++++++++++----- .../tpccalib/PHTpcCentralMembraneMatcher.cc | 34 +++++++++++++++++-- 2 files changed, 55 insertions(+), 10 deletions(-) diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc index a18f09cc32..84e32ab7b6 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneClusterizer.cc @@ -60,8 +60,8 @@ int PHTpcCentralMembraneClusterizer::InitRun(PHCompositeNode *topNode) hxy = new TH2F("hxy","cluster x:y",800,-100,100,800,-80,80); hz = new TH1F("hz","cluster z", 220, -2,2); - hz_pos = new TH1F("hz_pos","cluster z>0", 400, 0,110); - hz_neg = new TH1F("hz_neg","cluster z<0", 400, -110,0); + hz_pos = new TH1F("hz_pos","cluster z>0", 750, 0,200); + hz_neg = new TH1F("hz_neg","cluster z<0", 750, -200,0); hClustE[0]= new TH1F("hRawClusterEnergy","Cluster Energy Before Merging;E[?]",200,0,2000); hClustE[1] = new TH1F("hMatchedClusterEnergy","Pair Cluster Energy After Merging;E[?]",200,0,2000); @@ -259,14 +259,25 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) if(hphi_reco_neg[lay-7]->GetBinContent(phiBin) < (m_moduloThreshold + (mean_z_content_minus/18.0))) aboveThreshold = false; } - if( !aboveThreshold ) continue; - - if(cluster->getAdc() < _min_adc_value) continue; - if( std::abs(z) < _min_z_value) continue; + if( !aboveThreshold ){ + std::cout << "below threshold" << std::endl; + continue; + } - //require that z be within 1cm of peak in z distributions (separate for each side) - if( (z>=0 && std::abs(z - hz_pos->GetBinCenter(hz_pos->GetMaximumBin())) > 1.0) || (z<0 && std::abs(z - hz_neg->GetBinCenter(hz_neg->GetMaximumBin())) > 1.0) ) continue; + if(cluster->getAdc() < _min_adc_value){ + std::cout << "adc too low" << std::endl; + continue; + } + if( std::abs(z) < _min_z_value){ + std::cout << "z too low" << std::endl; + continue; + } + //require that z be within 1cm of peak in z distributions (separate for each side) + if( (z>=0 && std::abs(z - hz_pos->GetBinCenter(hz_pos->GetMaximumBin())) > 1.0) || (z<0 && std::abs(z - hz_neg->GetBinCenter(hz_neg->GetMaximumBin())) > 1.0) ){ + std::cout << "z not close enough to peak" << std::endl; + continue; + } //std::cout << "accepted cluster " << m_accepted_clusters << " parentID: " << cluster->getParentID() << " pctParentID: " << cluster->getPctParentID() << std::endl; @@ -326,9 +337,13 @@ int PHTpcCentralMembraneClusterizer::process_event(PHCompositeNode *topNode) int layerMatch = -1; + /* int nRowsMatch = 2; if(layer[i] >= 39) nRowsMatch = 4; else if(layer[i] >= 23) nRowsMatch = 3; + */ + + int nRowsMatch = 1; if( pos[i].Z() >= 0 ){ if(layer[i] == 7){ diff --git a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc index af897aa2e0..c0110239a7 100644 --- a/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc +++ b/offline/packages/tpccalib/PHTpcCentralMembraneMatcher.cc @@ -385,6 +385,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) std::vector reco_parentID; std::vector reco_pctParentID; + std::cout << PHWHERE << " process_event" << std::endl; + // reset output distortion correction container histograms for( const auto& harray:{m_dcc_out->m_hDRint, m_dcc_out->m_hDPint, m_dcc_out->m_hDZint, m_dcc_out->m_hentries} ) { @@ -392,7 +394,7 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) clust_r_phi_pos->Reset(); clust_r_phi_neg->Reset(); - if(!m_corrected_CMcluster_map || m_corrected_CMcluster_map->size() < 100){ + if(!m_corrected_CMcluster_map || m_corrected_CMcluster_map->size() < 10){//was100 m_event_index++; return Fun4AllReturnCodes::EVENT_OK; } @@ -459,7 +461,7 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) - if(Verbosity()) + if(Verbosity() > 0) { double raw_rad = sqrt( cmclus->getX()*cmclus->getX() + cmclus->getY()*cmclus->getY() ); double corr_rad = sqrt( tmp_pos.X()*tmp_pos.X() + tmp_pos.Y()*tmp_pos.Y() ); @@ -475,20 +477,38 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) //get global phi rotation for each module + std::cout << "pos R1" << std::endl; m_clustRotation_pos[0] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR1",151,206),clust_r_phi_pos->ProjectionX("cR1_pos",151,206)); + std::cout << "pos R2" << std::endl; m_clustRotation_pos[1] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR2",206,290),clust_r_phi_pos->ProjectionX("cR2_pos",206,290)); + std::cout << "pos R3" << std::endl; m_clustRotation_pos[2] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR3",290,499),clust_r_phi_pos->ProjectionX("cR3_pos",290,499)); + /* m_clustRotation_neg[0] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR1",151,206),clust_r_phi_neg->ProjectionX("cR1_neg",151,206)); m_clustRotation_neg[1] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR2",206,290),clust_r_phi_neg->ProjectionX("cR2_neg",206,290)); m_clustRotation_neg[2] = getPhiRotation_smoothed(hit_r_phi->ProjectionX("hR3",290,499),clust_r_phi_neg->ProjectionX("cR3_neg",290,499)); + */ + std::cout << "neg R1" << std::endl; + m_clustRotation_neg[0] = 0.0; + std::cout << "neg R2" << std::endl; + m_clustRotation_neg[1] = 0.0; + std::cout << "neg R3" << std::endl; + m_clustRotation_neg[2] = 0.0; + + + std::cout << "getting R peaks" << std::endl; //get hit and cluster peaks std::vector hit_RPeaks = getRPeaks(hit_r_phi); + std::cout << "getting R peaks pos" << std::endl; std::vector clust_RPeaks_pos = getRPeaks(clust_r_phi_pos); + std::cout << "getting R peaks neg" << std::endl; std::vector clust_RPeaks_neg = getRPeaks(clust_r_phi_neg); + std::cout << "finding gaps pos" << std::endl; + //identify where gaps between module 1&2 and 2&3 are by finding largest distances between peaks std::vector clust_RGaps_pos; int R12Gap_pos = -1; @@ -503,9 +523,13 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) R23Gap_pos = std::distance(clust_RGaps_pos.begin(),std::max_element(clust_RGaps_pos.begin()+R12Gap_pos+1,clust_RGaps_pos.end())); } + std::cout << "finding gaps neg" << std::endl; + + std::vector clust_RGaps_neg; int R12Gap_neg = -1; int R23Gap_neg = -1; + /* for(int i=0; i<(int)clust_RPeaks_neg.size()-1; i++) clust_RGaps_neg.push_back(clust_RPeaks_neg[i+1] - clust_RPeaks_neg[i]); int tmpGap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin(),clust_RGaps_neg.end())); if(tmpGap_neg > (int)clust_RGaps_neg.size()/2){ @@ -515,6 +539,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) R12Gap_neg = tmpGap_neg; R23Gap_neg = std::distance(clust_RGaps_neg.begin(),std::max_element(clust_RGaps_neg.begin()+R12Gap_neg+1,clust_RGaps_neg.end())); } +*/ + std::cout << "pos peak matching" << std::endl; int min_match_pos = 100; int min_match_neg = 100; @@ -538,6 +564,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) else if(i >= R23Gap_pos+1) hitMatches_pos.push_back(23+1 + i - (R23Gap_pos+1)); } + std::cout << "neg peak matching" << std::endl; + std::vector hitMatches_neg; for(int i=0; i<(int)clust_RPeaks_neg.size(); i++){ @@ -554,6 +582,8 @@ int PHTpcCentralMembraneMatcher::process_event(PHCompositeNode * /*topNode*/) } + std::cout << "starting cluster matching" << std::endl; + // Match reco and truth positions //std::map matched_pair; std::vector> matched_pair;