Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcCentralMembraneMatcher.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcCentralMembraneMatcher.h
1 // Tell emacs that this is a C++ source
2 // -*- C++ -*-.
3 #ifndef PHTPCCENTRALMEMBRANEMATCHER_H
4 #define PHTPCCENTRALMEMBRANEMATCHER_H
5 
12 #include <fun4all/SubsysReco.h>
13 #include <tpc/TpcDistortionCorrectionContainer.h>
14 #include <tpc/TpcDistortionCorrection.h>
15 #include <trackbase/TrkrDefs.h>
16 #include <trackbase/TrkrClusterContainer.h>
17 
18 #include <memory>
19 #include <string>
20 
21 class PHCompositeNode;
22 class CMFlashClusterContainer;
23 class CMFlashDifferenceContainer;
24 
25 class TF1;
26 class TFile;
27 class TH1F;
28 class TH1D;
29 class TH2F;
30 class TGraph;
31 class TNtuple;
32 class TVector3;
33 
34 class PHTpcCentralMembraneMatcher : public SubsysReco
35 {
36  public:
37 
38  PHTpcCentralMembraneMatcher(const std::string &name = "PHTpcCentralMembraneMatcher");
39 
40  ~PHTpcCentralMembraneMatcher() override = default;
41 
43  void setSavehistograms( bool value )
44  {m_savehistograms = value;}
45 
47  void setHistogramOutputfile(const std::string &outputfile)
48  {m_histogramfilename = outputfile;}
49 
51  void setOutputfile(const std::string &outputfile)
52  {m_outputfile = outputfile;}
53 
54  void setNMatchIter( int val ){ m_nMatchIter = val; }
55 
56  void set_useOnly_nClus2( bool val ){ m_useOnly_nClus2 = val; }
57 
58  void set_grid_dimensions( int phibins, int rbins );
59 
61  int InitRun(PHCompositeNode *topNode) override;
62 
64  int process_event(PHCompositeNode *topNode) override;
65 
67  int End(PHCompositeNode * topNode) override;
68 
69  private:
70 
71  int GetNodes(PHCompositeNode* topNode);
72 
74  TpcDistortionCorrection m_distortionCorrection;
75 
76  CMFlashClusterContainer *m_corrected_CMcluster_map{nullptr};
77  CMFlashDifferenceContainer *m_cm_flash_diffs{nullptr};
78 
80 
81  TpcDistortionCorrectionContainer* m_dcc_in_static{nullptr};
82  TpcDistortionCorrectionContainer* m_dcc_in_average{nullptr};
83 
85 
86  TpcDistortionCorrectionContainer* m_dcc_out{nullptr};
87 
89  /*
90  * this one is used to aggregate multple CM events in a single map
91  * it is not stored on the node tree but saved to a dedicated file at the end of the run
92  */
93  std::unique_ptr<TpcDistortionCorrectionContainer> m_dcc_out_aggregated;
94 
96  std::string m_outputfile = "CMDistortionCorrections.root";
97 
99 
100  bool m_savehistograms = false;
101  std::string m_histogramfilename = "PHTpcCentralMembraneMatcher.root";
102 
103  TH2F *hxy_reco = nullptr;
104  TH2F *hxy_truth = nullptr;
105  TH2F *hdrdphi = nullptr;
106  TH2F *hrdr = nullptr;
107  TH2F *hrdphi = nullptr;
108  TH1F *hdrphi = nullptr;
109  TH1F *hdphi = nullptr;
110  TH1F *hdr1_single = nullptr;
111  TH1F *hdr2_single = nullptr;
112  TH1F *hdr3_single = nullptr;
113  TH1F *hdr1_double = nullptr;
114  TH1F *hdr2_double = nullptr;
115  TH1F *hdr3_double = nullptr;
116  TH1F *hnclus = nullptr;
117 
118  std::unique_ptr<TFile> fout;
119 
120 
121  std::unique_ptr<TFile> fout2;
122  std::string m_histogramfilename2 = "CMMatcher.root";
123 
124  TH2F *hit_r_phi;
125 
126  TH2F *clust_r_phi_pos;
127  TH2F *clust_r_phi_neg;
128 
129  TNtuple *match_ntup = nullptr;
130 
131  int m_event_index = 0;
132 
134 
136  // double m_rad_cut= 0.5;
137 
139 
140  double m_phi_cut= 0.02;
141 
143 
144 
146  int m_phibins = 24;
147 
148  static constexpr float m_phiMin = 0;
149  static constexpr float m_phiMax = 2.*M_PI;
150 
152  int m_rbins = 12;
153 
154  static constexpr float m_rMin = 20; // cm
155  static constexpr float m_rMax = 80; // cm
156 
158 
160 
161  static constexpr double mm = 1.0;
162  static constexpr double cm = 10.0;
163 
164  static constexpr int nRadii = 8;
165  static constexpr int nStripes_R1 = 6;
166  static constexpr int nStripes_R2 = 8;
167  static constexpr int nStripes_R3 = 12;
168 
169  static constexpr int nPads_R1 = 6 * 16;
170  static constexpr int nPads_R2 = 8 * 16;
171  static constexpr int nPads_R3 = 12 * 16;
172 
174  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}};
175  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}};
176  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}};
177  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}};
178 
179  double cx1_e[nStripes_R1][nRadii];
180  double cx1[nStripes_R1][nRadii];
181  double cx2[nStripes_R2][nRadii];
182  double cx3[nStripes_R3][nRadii];
183 
184  double cy1_e[nStripes_R1][nRadii];
185  double cy1[nStripes_R1][nRadii];
186  double cy2[nStripes_R2][nRadii];
187  double cy3[nStripes_R3][nRadii];
188 
189  //Check which stripes get removed
190  std::array<int, nRadii> nGoodStripes_R1_e = {};
191  std::array<int, nRadii> nGoodStripes_R1 = {};
192  std::array<int, nRadii> nGoodStripes_R2 = {};
193  std::array<int, nRadii> nGoodStripes_R3 = {};
194 
196  static constexpr std::array<int, nRadii> keepThisAndAfter = {{1, 0, 1, 0, 1, 0, 1, 0}};
197 
199  static constexpr std::array<int, nRadii> keepUntil_R1_e = {{4, 4, 5, 4, 5, 5, 5, 5}};
200  static constexpr std::array<int, nRadii> keepUntil_R1 = {{5, 5, 6, 5, 6, 5, 6, 5}};
201  static constexpr std::array<int, nRadii> keepUntil_R2 = {{7, 7, 8, 7, 8, 8, 8, 8}};
202  static constexpr std::array<int, nRadii> keepUntil_R3 = {{11, 10, 11, 11, 11, 11, 12, 11}};
203 
204  std::array<int, nRadii> nStripesIn_R1_e = {};
205  std::array<int, nRadii> nStripesIn_R1 = {};
206  std::array<int, nRadii> nStripesIn_R2 = {};
207  std::array<int, nRadii> nStripesIn_R3 = {};
208  std::array<int, nRadii> nStripesBefore_R1_e = {};
209  std::array<int, nRadii> nStripesBefore_R1 = {};
210  std::array<int, nRadii> nStripesBefore_R2 = {};
211  std::array<int, nRadii> nStripesBefore_R3 = {};
212 
213  static constexpr int nStripesPerPetal = 213;
214  static constexpr int nPetals = 18;
215  static constexpr int nTotStripes = nStripesPerPetal * nPetals;
216 
217  void CalculateCenters(
218  int nPads,
219  const std::array<double, nRadii>& R,
220  std::array<int, nRadii>& nGoodStripes,
221  const std::array<int, nRadii>& keepUntil,
222  std::array<int, nRadii>& nStripesIn,
223  std::array<int, nRadii>& nStripesBefore,
224  double cx[][nRadii], double cy[][nRadii] );
225 
227  std::vector<TVector3> m_truth_pos;
228 
230 
231  bool m_useOnly_nClus2 = false;
232 
233  int m_nMatchIter = 2;
234 
235  double m_clustRotation_pos[3];
236  double m_clustRotation_neg[3];
237 
238  double getPhiRotation_smoothed( TH1D *hitHist, TH1D *clustHist );
239 
240  std::vector<double> getRPeaks(TH2F *r_phi);
241 
242  int getClusterRMatch( std::vector<int> hitMatches, std::vector<double> clusterPeaks, double clusterR);
243 };
244 
245 #endif // PHTPCCENTRALMEMBRANEMATCHER_H