packages/tpccalib/TpcCentralMembraneMatching.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef TPCCALIB_TPCCENTRALMEMBRANEMATCHING_H
0004 #define TPCCALIB_TPCCENTRALMEMBRANEMATCHING_H
0005
0006 /**
0007  * \file TpcCentralMembraneMatching.h
0008  * \brief match reconstructed CM clusters to CM pads, calculate differences, store on the node tree and compute distortion reconstruction maps
0009  * \author Tony Frawley <frawley@fsunuc.physics.fsu.edu>, Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0010  */
0011
0012 #include <tpc/TpcDistortionCorrection.h>
0013 #include <tpc/TpcDistortionCorrectionContainer.h>
0014
0015 #include <trackbase/TrkrClusterContainer.h>
0016 #include <trackbase/TrkrDefs.h>
0017
0018 #include <fun4all/SubsysReco.h>
0019
0020 #include <TGraph.h>
0021 #include <TGraph2D.h>
0022
0023 #include <memory>
0024 #include <string>
0025
0026 class PHCompositeNode;
0027 class CMFlashDifferenceContainer;
0028 class LaserClusterContainer;
0029 class EventHeader;
0030
0031 class TF1;
0032 class TFile;
0033 class TH1;
0034 class TH2;
0035 class TGraph;
0036 class TGraph2D;
0037 class TNtuple;
0038 class TTree;
0039 class TVector3;
0040
0041 class TpcCentralMembraneMatching : public SubsysReco
0042 {
0043  public:
0044   TpcCentralMembraneMatching(const std::string &name = "TpcCentralMembraneMatching");
0045
0046   ~TpcCentralMembraneMatching() override = default;
0047
0048   /// set to true to store evaluation histograms and ttrees
0049   void setSavehistograms(bool value)
0050   {
0051     m_savehistograms = value;
0052   }
0053
0054   /// output file name for evaluation histograms
0055   void setHistogramOutputfile(const std::string &outputfile)
0056   {
0057     m_histogramfilename = outputfile;
0058   }
0059
0060   /// output file name for storing the space charge reconstruction matrices
0061   void setOutputfile(const std::string &outputfile)
0062   {
0063     m_outputfile = outputfile;
0064   }
0065
0066   void setDebugOutputFile(const std::string &debugfile)
0067   {
0068     m_debugfilename = debugfile;
0069   }
0070
0071   void set_nHitsInCluster_minimum(const unsigned int minHits)
0072   {
0073     m_nHitsInCuster_minimum = minHits;
0074   }
0075
0076   void set_fixShifts(bool fixShifts)
0077   {
0078     m_fixShifts = fixShifts;
0079   }
0080
0081   void set_fieldOn(bool fieldOn)
0082   {
0083     m_fieldOn = fieldOn;
0084   }
0085
0086   void set_doFancy(bool fancy)
0087   {
0088     m_doFancy = fancy;
0089   }
0090
0091   void set_doHadd(bool hadd)
0092   {
0093     m_doHadd = hadd;
0094   }
0095
0096   void set_averageMode(bool averageMode)
0097   {
0098     m_averageMode = averageMode;
0099   }
0100
0101   void set_event_sequence(int seq)
0102   {
0103     m_event_sequence = seq;
0104     m_event_index = 100 * seq;
0105   }
0106
0107   // void set_laminationFile(const std::string& filename)
0108   //{
0109   // m_lamfilename = filename;
0110   //}
0111
0112   void set_grid_dimensions(int phibins, int rbins);
0113
0114   //! run initialization
0115   int InitRun(PHCompositeNode *topNode) override;
0116
0117   //! event processing
0118   int process_event(PHCompositeNode *topNode) override;
0119
0120   //! end of process
0121   int End(PHCompositeNode *topNode) override;
0122
0123  private:
0124   EventHeader *eventHeader{nullptr};
0125
0126   int GetNodes(PHCompositeNode *topNode);
0127
0128   double getPhiRotation_smoothed(TH1 *hitHist, TH1 *clustHist, bool side);
0129
0130   std::vector<int> doGlobalRMatching(TH2 *r_phi, bool side);
0131
0132   void getRegionPhiRotation(bool side);
0133
0134   int getClusterRMatch(double clusterR, int side);
0135
0136   //! tpc distortion correction utility class
0137   TpcDistortionCorrection m_distortionCorrection;
0138
0139   //! CMFlashClusterContainer *m_corrected_CMcluster_map{nullptr};
0140   LaserClusterContainer *m_corrected_CMcluster_map{nullptr};
0141   CMFlashDifferenceContainer *m_cm_flash_diffs{nullptr};
0142
0143   //!@name distortion correction containers
0144   //@{
0145   /** used in input to correct CM clusters before calculating residuals */
0146   TpcDistortionCorrectionContainer *m_dcc_in_module_edge{nullptr};
0147   TpcDistortionCorrectionContainer *m_dcc_in_static{nullptr};
0148   TpcDistortionCorrectionContainer *m_dcc_in_average{nullptr};
0149   //@}
0150
0151   //! fluctuation distortion container
0152   /** used in output to write fluctuation distortions */
0153   TpcDistortionCorrectionContainer *m_dcc_out{nullptr};
0154
0155   /// local fluctuation distortion container
0156   /*
0157    * this one is used to aggregate multple CM events in a single map
0158    * it is not stored on the node tree but saved to a dedicated file at the end of the run
0159    */
0160   std::unique_ptr<TpcDistortionCorrectionContainer> m_dcc_out_aggregated;
0161
0162   /// output file, to which aggregated central membrane distortion corrections are stored
0163   std::string m_outputfile{"CMDistortionCorrections.root"};
0164
0165   ///@name evaluation histograms
0166   //@{
0167   bool m_savehistograms{false};
0168   std::string m_histogramfilename = "TpcCentralMembraneMatching.root";
0169
0170   TH2 *hxy_reco{nullptr};
0171   TH2 *hxy_truth{nullptr};
0172   TH2 *hdrdphi{nullptr};
0173   TH2 *hrdr{nullptr};
0174   TH2 *hrdphi{nullptr};
0175   TH1 *hdrphi{nullptr};
0176   TH1 *hdphi{nullptr};
0177   TH1 *hdr1_single{nullptr};
0178   TH1 *hdr2_single{nullptr};
0179   TH1 *hdr3_single{nullptr};
0180   TH1 *hdr1_double{nullptr};
0181   TH1 *hdr2_double{nullptr};
0182   TH1 *hdr3_double{nullptr};
0183   TH1 *hnclus{nullptr};
0184
0185   TFile *fout{};
0186
0187   TFile *m_debugfile{};
0188   std::string m_debugfilename{"CMMatcher.root"};
0189
0190   TH2 *truth_r_phi[2]{nullptr};
0191
0192   TH2 *reco_r_phi[2]{nullptr};
0193
0194   TH2 *m_matchResiduals[2]{nullptr};
0195
0196   //  TNtuple *match_ntup {nullptr};
0197   TTree *match_tree{nullptr};
0198   TTree *event_tree{nullptr};
0199
0200   bool m_useHeader{true};
0201   bool m_averageMode{false};
0202
0203   std::vector<bool> e_matched;
0204   std::vector<int> e_truthIndex;
0205   std::vector<float> e_truthR;
0206   std::vector<float> e_truthPhi;
0207   std::vector<float> e_recoR;
0208   std::vector<float> e_recoPhi;
0209   std::vector<bool> e_side;
0210   std::vector<bool> e_isSaturated;
0211   std::vector<bool> e_fitMode;
0212   std::vector<float> e_NNR;
0213   std::vector<float> e_NNPhi;
0214   std::vector<float> e_NNDistance;
0215   std::vector<int> e_NNIndex;
0216   std::vector<unsigned int> e_adc;
0217
0218   int m_event_index{0};
0219   int m_event_sequence{0};
0220   bool m_matched{false};
0221   int m_truthIndex{0};
0222   float m_truthR{0.0};
0223   float m_truthPhi{0.0};
0224   float m_recoR{0.0};
0225   float m_recoPhi{0.0};
0226   float m_recoZ{0.0};
0227   float m_rawR{0.0};
0228   float m_rawPhi{0.0};
0229   float m_staticR{0.0};
0230   float m_staticPhi{0.0};
0231   bool m_side{false};
0232   bool m_fitMode{false};
0233   bool m_isSaturated{false};
0234   unsigned int m_adc{0};
0235   unsigned int m_nhits{0};
0236   unsigned int m_nLayers{0};
0237   unsigned int m_nIPhi{0};
0238   unsigned int m_nIT{0};
0239   float m_layersSD{0.0};
0240   float m_IPhiSD{0.0};
0241   float m_ITSD{0.0};
0242   float m_layersWeightedSD{0.0};
0243   float m_IPhiWeightedSD{0.0};
0244   float m_ITWeightedSD{0.0};
0245   int m_lowShift{0};
0246   int m_highShift{0};
0247   float m_phiRotation{0.0};
0248   float m_distanceToTruth{0.0};
0249   float m_NNDistance{0.0};
0250   float m_NNR{0.0};
0251   float m_NNPhi{0.0};
0252   int m_NNIndex{0};
0253
0254   //@}
0255
0256   // double meanA[2] = {0.0, 0.0};
0257   // double meanB[2] = {0.0, 0.0};
0258   // double meanC[2] = {0.0, 0.0};
0259   // std::vector<double> lamPhis[2];
0260
0261   // std::string m_lamfilename = "";
0262
0263   unsigned int m_nHitsInCuster_minimum = 5;
0264
0265   /// radius cut for matching clusters to pad, for size 2 clusters
0266   //  double m_rad_cut{0}.5;
0267
0268   std::map<int, float> deltaRByIndex;
0269   std::map<int, float> deltaPhiByIndex;
0270   std::map<int, float> meanRByIndex;
0271   std::map<int, float> meanPhiByIndex;
0272   std::map<int, int> nByIndex;
0273
0274   TGraph2D *gr_dR[2]{nullptr, nullptr};
0275   TGraph2D *gr_dPhi[2]{nullptr, nullptr};
0276   TGraph *gr_points[2]{nullptr, nullptr};
0277
0278   /// phi cut for matching clusters to pad
0279   /** TODO: this will need to be adjusted to match beam-induced time averaged distortions */
0280   double m_phi_cut{0.025};
0281
0282   ///@name distortion correction histograms
0283   //@{
0284
0285   /// distortion correction grid size along phi
0286   // int m_phibins{24};
0287   int m_phibins{500};
0288
0289   static constexpr float m_phiMin{0};
0290   static constexpr float m_phiMax{2. * M_PI};
0291
0292   /// distortion correction grid size along r
0293   // int m_rbins{12};
0294   int m_rbins{500};
0295
0296   static constexpr float m_rMin{20};  // cm
0297   static constexpr float m_rMax{80};  // cm
0298
0299   //@}
0300
0301   ///@name central membrane pads definitions
0302   //@{
0303   static constexpr double mm{1.0};
0304   static constexpr double cm{10.0};
0305
0306   static constexpr int nRadii{8};
0307   static constexpr int nStripes_R1{6};
0308   static constexpr int nStripes_R2{8};
0309   static constexpr int nStripes_R3{12};
0310
0311   static constexpr int nPads_R1{6 * 16};
0312   static constexpr int nPads_R2{8 * 16};
0313   static constexpr int nPads_R3{12 * 16};
0314
0315   /// stripe radii
0316   static constexpr std::array<double, nRadii> R1_e{{227.0902789 * mm, 238.4100043 * mm, 249.7297296 * mm, 261.049455 * mm, 272.3691804 * mm, 283.6889058 * mm, 295.0086312 * mm, 306.3283566 * mm}};
0317   static constexpr std::array<double, nRadii> R1{{317.648082 * mm, 328.9678074 * mm, 340.2875328 * mm, 351.6072582 * mm, 362.9269836 * mm, 374.246709 * mm, 385.5664344 * mm, 396.8861597 * mm}};
0318   static constexpr std::array<double, nRadii> R2{{421.705532 * mm, 442.119258 * mm, 462.532984 * mm, 482.9467608 * mm, 503.36069 * mm, 523.774416 * mm, 544.188015 * mm, 564.601868 * mm}};
0319   static constexpr std::array<double, nRadii> R3{{594.6048725 * mm, 616.545823 * mm, 638.4867738 * mm, 660.4277246 * mm, 682.3686754 * mm, 704.3096262 * mm, 726.250577 * mm, 748.1915277 * mm}};
0320
0321   double cx1_e[nStripes_R1][nRadii]{};
0322   double cx1[nStripes_R1][nRadii]{};
0323   double cx2[nStripes_R2][nRadii]{};
0324   double cx3[nStripes_R3][nRadii]{};
0325
0326   double cy1_e[nStripes_R1][nRadii]{};
0327   double cy1[nStripes_R1][nRadii]{};
0328   double cy2[nStripes_R2][nRadii]{};
0329   double cy3[nStripes_R3][nRadii]{};
0330
0331   // Check which stripes get removed
0332   std::array<int, nRadii> nGoodStripes_R1_e = {};
0333   std::array<int, nRadii> nGoodStripes_R1 = {};
0334   std::array<int, nRadii> nGoodStripes_R2 = {};
0335   std::array<int, nRadii> nGoodStripes_R3 = {};
0336
0337   /// min stripe index
0338   static constexpr std::array<int, nRadii> keepThisAndAfter{{1, 0, 1, 0, 1, 0, 1, 0}};
0339
0340   /// max stripe index
0341   static constexpr std::array<int, nRadii> keepUntil_R1_e{{4, 4, 5, 4, 5, 5, 5, 5}};
0342   static constexpr std::array<int, nRadii> keepUntil_R1{{5, 5, 6, 5, 6, 5, 6, 5}};
0343   static constexpr std::array<int, nRadii> keepUntil_R2{{7, 7, 8, 7, 8, 8, 8, 8}};
0344   static constexpr std::array<int, nRadii> keepUntil_R3{{11, 10, 11, 11, 11, 11, 12, 11}};
0345
0346   std::array<int, nRadii> nStripesIn_R1_e{};
0347   std::array<int, nRadii> nStripesIn_R1{};
0348   std::array<int, nRadii> nStripesIn_R2{};
0349   std::array<int, nRadii> nStripesIn_R3{};
0350   std::array<int, nRadii> nStripesBefore_R1_e{};
0351   std::array<int, nRadii> nStripesBefore_R1{};
0352   std::array<int, nRadii> nStripesBefore_R2{};
0353   std::array<int, nRadii> nStripesBefore_R3{};
0354
0355   static constexpr int nStripesPerPetal{213};
0356   static constexpr int nPetals{18};
0357   static constexpr int nTotStripes{nStripesPerPetal * nPetals};
0358
0359   void CalculateCenters(
0360       int nPads,
0361       const std::array<double, nRadii> &R,
0362       std::array<int, nRadii> &nGoodStripes,
0363       const std::array<int, nRadii> &keepUntil,
0364       std::array<int, nRadii> &nStripesIn,
0365       std::array<int, nRadii> &nStripesBefore,
0366       double cx[][nRadii], double cy[][nRadii]);
0367
0368   /// store centers of all central membrane pads
0369   std::vector<TVector3> m_truth_pos;
0370   std::vector<int> m_truth_index;
0371
0372   std::vector<double> m_truth_RPeaks{22.709, 23.841, 24.973, 26.1049, 27.2369, 28.3689, 29.5009, 30.6328, 31.7648, 32.8968, 34.0288, 35.1607, 36.2927, 37.4247, 38.5566, 39.6886, 42.1706, 44.2119, 46.2533, 48.2947, 50.3361, 52.3774, 54.4188, 56.4602, 59.4605, 61.6546, 63.8487, 66.0428, 68.2369, 70.431, 72.6251, 74.8192};
0373
0374   //@}
0375
0376   std::vector<double> phiSpacing = {0.0749989, 0.0729917, 0.0711665, 0.0694996, 0.0679712, 0.0665648, 0.0652663, 0.0640638, 0.062947, 0.0619071, 0.0609364, 0.0600281, 0.0591765, 0.0583765, 0.0576234, 0.0569132, 0.0473084, 0.0462573, 0.045299, 0.0444217, 0.0436155, 0.0428722, 0.0421847, 0.0415468, 0.0325076, 0.0319331, 0.031398, 0.0308985, 0.0304311, 0.0299928, 0.029581, 0.0291934};
0377
0378   bool m_fixShifts{false};
0379   bool m_fieldOn{true};
0380   bool m_doFancy{false};
0381   bool m_doHadd{false};
0382
0383   std::vector<double> m_reco_RPeaks[2];
0384   double m_m[2]{};
0385   double m_b[2]{};
0386   int m_matchLow[2]{};
0387   int m_matchHigh[2]{};
0388   std::vector<int> m_reco_RMatches[2];
0389
0390   double m_recoRotation[2][3]{{-999, -999, -999}, {-999, -999, -999}};
0391 };
0392
0393 #endif  // PHTPCCENTRALMEMBRANEMATCHER_H