diff --git a/offline/packages/tpccalib/PHTpcResiduals.cc b/offline/packages/tpccalib/PHTpcResiduals.cc index 788e4ae06d..7d8ec8e095 100644 --- a/offline/packages/tpccalib/PHTpcResiduals.cc +++ b/offline/packages/tpccalib/PHTpcResiduals.cc @@ -128,6 +128,8 @@ namespace PHTpcResiduals::PHTpcResiduals(const std::string& name) : SubsysReco(name) , m_matrix_container(new TpcSpaceChargeMatrixContainerv2) + , m_matrix_container_pos(new TpcSpaceChargeMatrixContainerv2) + , m_matrix_container_neg(new TpcSpaceChargeMatrixContainerv2) { } @@ -139,6 +141,7 @@ int PHTpcResiduals::Init(PHCompositeNode* /*topNode*/) std::cout << "PHTpcResiduals::Init - m_maxTBeta: " << m_maxTBeta << std::endl; std::cout << "PHTpcResiduals::Init - m_maxResidualDrphi: " << m_maxResidualDrphi << " cm" << std::endl; std::cout << "PHTpcResiduals::Init - m_maxResidualDz: " << m_maxResidualDz << " cm" << std::endl; + std::cout << "PHTpcResiduals::Init - m_ignoreEdgeClusters: " << m_ignoreEdgeClusters << std::endl; std::cout << "PHTpcResiduals::Init - m_minRPhiErr: " << m_minRPhiErr << " cm" << std::endl; std::cout << "PHTpcResiduals::Init - m_minZErr: " << m_minZErr << " cm" << std::endl; std::cout << "PHTpcResiduals::Init - m_minPt: " << m_minPt << " GeV/c" << std::endl; @@ -186,11 +189,25 @@ int PHTpcResiduals::End(PHCompositeNode* /*topNode*/) std::cout << "PHTpcResiduals::End - writing matrices to " << m_outputfile << std::endl; // save matrix container in output file - if (m_matrix_container) + if (m_matrix_container || m_matrix_container_pos || m_matrix_container_neg) { std::unique_ptr outputfile(TFile::Open(m_outputfile.c_str(), "RECREATE")); outputfile->cd(); - m_matrix_container->Write("TpcSpaceChargeMatrixContainer"); + + if( m_matrix_container ) { + std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer object. entries: " << m_matrix_container->get_entries() << std::endl; + m_matrix_container->Write("TpcSpaceChargeMatrixContainer"); + } + + if( m_matrix_container_pos ) { + std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer_pos object. entries: " << m_matrix_container_pos->get_entries() << std::endl; + m_matrix_container_pos->Write("TpcSpaceChargeMatrixContainer_pos"); + } + + if( m_matrix_container_neg ) { + std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer_neg object. entries: " << m_matrix_container_neg->get_entries() << std::endl; + m_matrix_container_neg->Write("TpcSpaceChargeMatrixContainer_neg"); + } } // print counters @@ -456,8 +473,6 @@ void PHTpcResiduals::processTrack(SvtxTrack* track) for (const auto& cluskey : get_cluster_keys(track)) { - // increment counter - ++m_total_clusters; // make sure cluster is from TPC const auto detId = TrkrDefs::getTrkrId(cluskey); @@ -466,6 +481,10 @@ void PHTpcResiduals::processTrack(SvtxTrack* track) continue; } + // increment counter + /* only counts TPC clusters */ + ++m_total_clusters; + // find matching track state const auto stateiter = std::find_if( track->begin_states(), track->end_states(), [&cluskey]( const auto& state_pair ) { return state_pair.second->get_cluskey() == cluskey; } ); @@ -478,6 +497,14 @@ void PHTpcResiduals::processTrack(SvtxTrack* track) // calculate residuals with respect to cluster auto* const cluster = m_clusterContainer->findCluster(cluskey); + + // check cluster + if( !cluster ) { continue; } + + // check if cluster is an edge + if( m_ignoreEdgeClusters && cluster->getEdge() > 0 && cluster->getEdge() < std::numeric_limits::max() ) + { continue; } + const auto globClusPos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluskey, cluster, crossing); const double clusR = get_r(globClusPos(0), globClusPos(1)); const double clusPhi = std::atan2(globClusPos(1), globClusPos(0)); @@ -650,43 +677,58 @@ void PHTpcResiduals::processTrack(SvtxTrack* track) continue; } - // Fill distortion matrices - m_matrix_container->add_to_lhs(index, 0, 0, square(clusR) / erp); - m_matrix_container->add_to_lhs(index, 0, 1, 0); - m_matrix_container->add_to_lhs(index, 0, 2, clusR * trackAlpha / erp); + std::vector containers; + containers.emplace_back( m_matrix_container.get() ); + if( track->get_positive_charge() ) { + containers.emplace_back( m_matrix_container_pos.get() ); + } else { + containers.emplace_back( m_matrix_container_neg.get() ); + } + - m_matrix_container->add_to_lhs(index, 1, 0, 0); - m_matrix_container->add_to_lhs(index, 1, 1, 1. / ez); - m_matrix_container->add_to_lhs(index, 1, 2, trackBeta / ez); + for( auto& container:containers ) + { + + if( !container ) { continue; } + + // Fill distortion matrices + container->add_to_lhs(index, 0, 0, square(clusR) / erp); + container->add_to_lhs(index, 0, 1, 0); + container->add_to_lhs(index, 0, 2, clusR * trackAlpha / erp); - m_matrix_container->add_to_lhs(index, 2, 0, clusR * trackAlpha / erp); - m_matrix_container->add_to_lhs(index, 2, 1, trackBeta / ez); - m_matrix_container->add_to_lhs(index, 2, 2, square(trackAlpha) / erp + square(trackBeta) / ez); + container->add_to_lhs(index, 1, 0, 0); + container->add_to_lhs(index, 1, 1, 1. / ez); + container->add_to_lhs(index, 1, 2, trackBeta / ez); - m_matrix_container->add_to_rhs(index, 0, clusR * drphi / erp); - m_matrix_container->add_to_rhs(index, 1, dz / ez); - m_matrix_container->add_to_rhs(index, 2, trackAlpha * drphi / erp + trackBeta * dz / ez); + container->add_to_lhs(index, 2, 0, clusR * trackAlpha / erp); + container->add_to_lhs(index, 2, 1, trackBeta / ez); + container->add_to_lhs(index, 2, 2, square(trackAlpha) / erp + square(trackBeta) / ez); - // also update rphi reduced matrices - m_matrix_container->add_to_lhs_rphi(index, 0, 0, square(clusR) / erp); - m_matrix_container->add_to_lhs_rphi(index, 0, 1, clusR * trackAlpha / erp); - m_matrix_container->add_to_lhs_rphi(index, 1, 0, clusR * trackAlpha / erp); - m_matrix_container->add_to_lhs_rphi(index, 1, 1, square(trackAlpha) / erp); + container->add_to_rhs(index, 0, clusR * drphi / erp); + container->add_to_rhs(index, 1, dz / ez); + container->add_to_rhs(index, 2, trackAlpha * drphi / erp + trackBeta * dz / ez); - m_matrix_container->add_to_rhs_rphi(index, 0, clusR * drphi / erp); - m_matrix_container->add_to_rhs_rphi(index, 1, trackAlpha * drphi / erp); + // also update rphi reduced matrices + container->add_to_lhs_rphi(index, 0, 0, square(clusR) / erp); + container->add_to_lhs_rphi(index, 0, 1, clusR * trackAlpha / erp); + container->add_to_lhs_rphi(index, 1, 0, clusR * trackAlpha / erp); + container->add_to_lhs_rphi(index, 1, 1, square(trackAlpha) / erp); - // also update z reduced matrices - m_matrix_container->add_to_lhs_z(index, 0, 0, 1. / ez); - m_matrix_container->add_to_lhs_z(index, 0, 1, trackBeta / ez); - m_matrix_container->add_to_lhs_z(index, 1, 0, trackBeta / ez); - m_matrix_container->add_to_lhs_z(index, 1, 1, square(trackBeta) / ez); + container->add_to_rhs_rphi(index, 0, clusR * drphi / erp); + container->add_to_rhs_rphi(index, 1, trackAlpha * drphi / erp); - m_matrix_container->add_to_rhs_z(index, 0, dz / ez); - m_matrix_container->add_to_rhs_z(index, 1, trackBeta * dz / ez); + // also update z reduced matrices + container->add_to_lhs_z(index, 0, 0, 1. / ez); + container->add_to_lhs_z(index, 0, 1, trackBeta / ez); + container->add_to_lhs_z(index, 1, 0, trackBeta / ez); + container->add_to_lhs_z(index, 1, 1, square(trackBeta) / ez); - // update entries in cell - m_matrix_container->add_to_entries(index); + container->add_to_rhs_z(index, 0, dz / ez); + container->add_to_rhs_z(index, 1, trackBeta * dz / ez); + + // update entries in cell + container->add_to_entries(index); + } // increment number of accepted clusters ++m_accepted_clusters; @@ -776,5 +818,7 @@ int PHTpcResiduals::getNodes(PHCompositeNode* topNode) //____________________________________________________________________________ void PHTpcResiduals::setGridDimensions(const int phiBins, const int rBins, const int zBins) { - m_matrix_container->set_grid_dimensions(phiBins, rBins, zBins); + if( m_matrix_container ) { m_matrix_container->set_grid_dimensions(phiBins, rBins, zBins); } + if( m_matrix_container_pos ) { m_matrix_container_pos->set_grid_dimensions(phiBins, rBins, zBins); } + if( m_matrix_container_neg ) { m_matrix_container_neg->set_grid_dimensions(phiBins, rBins, zBins); } } diff --git a/offline/packages/tpccalib/PHTpcResiduals.h b/offline/packages/tpccalib/PHTpcResiduals.h index 1e750aad3a..3db64d5f4c 100644 --- a/offline/packages/tpccalib/PHTpcResiduals.h +++ b/offline/packages/tpccalib/PHTpcResiduals.h @@ -63,11 +63,20 @@ class PHTpcResiduals : public SubsysReco } //@} + /// true to remove "edge" clusters + /** these are clusters that touch the edge of a detector and are considered pathological */ + void setIgnoreEdgeClusters( bool value ) + { m_ignoreEdgeClusters = value; } + + /// minimum value for RPhi error. + /** a too small rphi error is usually a sign of pathological TPC cluster */ void setMinRPhiErr(float minRPhiErr) { m_minRPhiErr = minRPhiErr; } + /// minimum value for z error. + /** a too small z error is usually a sign of pathological TPC cluster */ void setMinZErr(float minZErr) { m_minZErr = minZErr; @@ -171,6 +180,9 @@ class PHTpcResiduals : public SubsysReco float m_maxTBeta = 1.5; float m_maxResidualDz = 0.5; // cm + /// ignore edge clusters + bool m_ignoreEdgeClusters = false; + float m_minRPhiErr = 0.005; // 0.005cm -- 50um float m_minZErr = 0.01; // 0.01cm -- 100um @@ -193,6 +205,12 @@ class PHTpcResiduals : public SubsysReco /// matrix container std::unique_ptr m_matrix_container; + /// matrix container positive charges only + std::unique_ptr m_matrix_container_pos; + + /// matrix container negative charges only + std::unique_ptr m_matrix_container_neg; + // TODO: check if needed int m_event = 0; diff --git a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainer.h b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainer.h index 701bc147d0..767e0f5135 100644 --- a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainer.h +++ b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainer.h @@ -45,6 +45,10 @@ class TpcSpaceChargeMatrixContainer : public PHObject virtual int get_cell_index( int /*iphibin*/, int /*irbin*/, int /*izbin*/ ) const { return -1; } + /// get all entries + virtual int get_entries() const + { return 0; } + /// get entries for a given cell virtual int get_entries( int /*cell_index*/ ) const { return 0; } diff --git a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.cc b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.cc index 3fa25bbbf3..e13c9baffd 100644 --- a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.cc +++ b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.cc @@ -8,6 +8,8 @@ #include "TpcSpaceChargeMatrixContainerv2.h" +#include + //___________________________________________________________ TpcSpaceChargeMatrixContainerv2::TpcSpaceChargeMatrixContainerv2() { @@ -56,6 +58,10 @@ int TpcSpaceChargeMatrixContainerv2::get_cell_index(int iphi, int ir, int iz) co return iz + m_zbins * (ir + m_rbins * iphi); } +//___________________________________________________________ +int TpcSpaceChargeMatrixContainerv2::get_entries() const +{ return std::accumulate( m_entries.begin(), m_entries.end(), 0); } + //___________________________________________________________ int TpcSpaceChargeMatrixContainerv2::get_entries(int cell_index) const { diff --git a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.h b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.h index 8005e4f1f3..5942a8c126 100644 --- a/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.h +++ b/offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.h @@ -39,6 +39,9 @@ class TpcSpaceChargeMatrixContainerv2 : public TpcSpaceChargeMatrixContainer /// get grid index for given sub-indexes int get_cell_index(int iphibin, int irbin, int izbin) const override; + /// get all entries + int get_entries() const override; + /// get entries for a given cell int get_entries(int cell_index) const override; diff --git a/offline/packages/tpccalib/TpcSpaceChargeMatrixInversion.cc b/offline/packages/tpccalib/TpcSpaceChargeMatrixInversion.cc index 9c787730aa..5fa957ab3e 100644 --- a/offline/packages/tpccalib/TpcSpaceChargeMatrixInversion.cc +++ b/offline/packages/tpccalib/TpcSpaceChargeMatrixInversion.cc @@ -165,6 +165,14 @@ bool TpcSpaceChargeMatrixInversion::add_from_file(const std::string& shortfilena return false; } + if( Verbosity() ) { + std::cout << "TpcSpaceChargeMatrixInversion::add_from_file -" + << " file: " << filename + << " objectname: " << objectname + << " entries: " << source->get_entries() + << std::endl; + } + // add object return add(*source); } @@ -200,6 +208,8 @@ void TpcSpaceChargeMatrixInversion::calculate_distortion_corrections(const Inver exit(1); } + std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortion_corrections - entries: " << m_matrix_container->get_entries() << std::endl; + // get grid dimensions from matrix container int phibins = 0; int rbins = 0; diff --git a/offline/packages/trackreco/PHMicromegasTpcTrackMatching.h b/offline/packages/trackreco/PHMicromegasTpcTrackMatching.h index 7f9f8e2611..23bd75854f 100644 --- a/offline/packages/trackreco/PHMicromegasTpcTrackMatching.h +++ b/offline/packages/trackreco/PHMicromegasTpcTrackMatching.h @@ -48,7 +48,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco int InitRun(PHCompositeNode* topNode) override; int process_event(PHCompositeNode*) override; int End(PHCompositeNode*) override; - + // deprecated calls inline void set_sc_calib_mode(const bool) {} inline void set_collision_rate(const double) {} @@ -77,7 +77,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco unsigned int _max_tpc_layer = 55; // pt cut for field-on data - float _pt_cut = 0.5; + float _pt_cut = 0.2; // delta_phi window between the last cluster in the tracklet and the projection float _dphi_cut = 0.9; @@ -90,7 +90,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco TrackSeedContainer* _si_track_map{nullptr}; std::string _clustermap_name = "TRKR_CLUSTER"; - + //! default rphi search window for each layer std::array _rphi_search_win{0.25, 13.0};