Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcDirectLaserReconstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcDirectLaserReconstruction.cc
1 
7 #include "TpcDirectLaserReconstruction.h"
8 
9 #include "TpcSpaceChargeMatrixContainerv1.h"
10 
11 #include <fun4all/Fun4AllReturnCodes.h>
12 #include <phool/getClass.h>
13 
14 #include <g4detectors/PHG4TpcCylinderGeom.h>
15 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
16 
17 #include <trackbase/ActsTrackingGeometry.h>
18 #include <trackbase_historic/SvtxTrack.h>
19 #include <trackbase_historic/SvtxTrackMap.h>
20 #include <trackbase_historic/SvtxTrackState_v1.h>
21 #include <trackbase/TrkrHitSet.h>
22 #include <trackbase/TrkrHitSetContainer.h>
23 #include <trackbase/TrkrHitv2.h> // for TrkrHit
24 #include <trackbase/TpcDefs.h>
25 
26 #include <TFile.h>
27 #include <TH1.h>
28 #include <TH2.h>
29 #include <TH3.h>
30 #include <TVector3.h>
31 #include <TNtuple.h>
32 
33 #include <cassert>
34 
35 namespace
36 {
37 
39  template<class T> class range_adaptor
40  {
41  public:
42  range_adaptor( const T& range ):m_range(range){}
43  inline const typename T::first_type& begin() {return m_range.first;}
44  inline const typename T::second_type& end() {return m_range.second;}
45  private:
46  T m_range;
47  };
48 
50  template<class T>
51  inline constexpr T square( const T& x ) { return x*x; }
52 
54  template<class T>
55  inline constexpr T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
56 
57  // calculate intersection between line and circle
58  std::pair<double, double> line_circle_intersection( const TVector3& p, const TVector3& d, double radius )
59  {
60  std::pair<double,double> ret = std::make_pair(-1, -1);
61 
62  const double A = square(d.x()) + square(d.y());
63  const double B = 2*p.x()*d.x() + 2*p.y()*d.y();
64  const double C = square(p.x()) + square(p.y()) - square(radius);
65  const double delta = square(B)-4*A*C;
66  if( delta < 0 ) return ret;
67 
68  // we want 1 intersection
69  const double tup = (-B + std::sqrt(delta))/(2*A);
70  const double tdn = (-B-sqrt(delta))/(2*A);
71  // std::cout << " tup " << tup << " tdn " << tdn << std::endl;
72  ret = std::make_pair(tup, tdn);
73 
74  return ret;
75 
76  /*
77  if( tup > 0 && tup < tdn)
78  return tup;
79  if(tdn > 0 && tdn < tup)
80  return tdn;
81 
82  // no valid extrapolation
83  return -1;
84  */
85  }
86 
88  inline std::ostream& operator << (std::ostream& out, const TVector3& vector )
89  {
90  out << "( " << vector.x() << ", " << vector.y() << ", " << vector.z() << ")";
91  return out;
92  }
93 
95  template< class T>
96  inline constexpr T delta_phi( const T& phi )
97  {
98  if( phi >= M_PI ) return phi - 2*M_PI;
99  else if( phi < -M_PI ) return phi + 2*M_PI;
100  else return phi;
101  }
102 
103  // phi range
104  static constexpr float m_phimin = 0;
105  static constexpr float m_phimax = 2.*M_PI;
106 
107  // TODO: could try to get the r and z range from TPC geometry
108  // r range
109  static constexpr float m_rmin = 20;
110  static constexpr float m_rmax = 78;
111 
112  // z range
113  static constexpr float m_zmin = -105.5;
114  static constexpr float m_zmax = 105.5;
115 
116 }
117 
118 //_____________________________________________________________________
119 TpcDirectLaserReconstruction::TpcDirectLaserReconstruction( const std::string& name ):
120  SubsysReco( name)
121  , PHParameterInterface(name)
122  , m_matrix_container( new TpcSpaceChargeMatrixContainerv1 )
123 {
124  InitializeParameters();
125 }
126 
127 //_____________________________________________________________________
128 int TpcDirectLaserReconstruction::Init(PHCompositeNode*)
129 {
130  m_total_hits = 0;
131  m_matched_hits = 0;
132  m_accepted_clusters = 0;
133 
134  if( m_savehistograms ) create_histograms();
135  return Fun4AllReturnCodes::EVENT_OK;
136 }
137 
138 //_____________________________________________________________________
139 int TpcDirectLaserReconstruction::InitRun(PHCompositeNode* )
140 {
141  UpdateParametersWithMacro();
142  m_max_dca = get_double_param( "directlaser_max_dca" );
143  m_max_drphi = get_double_param( "directlaser_max_drphi" );
144  m_max_dz = get_double_param( "directlaser_max_dz" );
145 
146  // print
147  if( Verbosity() )
148  {
149  std::cout
150  << "TpcDirectLaserReconstruction::InitRun\n"
151  << " m_outputfile: " << m_outputfile << "\n"
152  << " m_max_dca: " << m_max_dca << "\n"
153  << " m_max_drphi: " << m_max_drphi << "\n"
154  << " m_max_dz: " << m_max_dz << "\n"
155  << std::endl;
156 
157  // also identify the matrix container
158  m_matrix_container->identify();
159  }
160 
161  return Fun4AllReturnCodes::EVENT_OK;
162 }
163 
164 //_____________________________________________________________________
165 int TpcDirectLaserReconstruction::process_event(PHCompositeNode* topNode)
166 {
167  // load nodes
168  const auto res = load_nodes(topNode);
169  if( res != Fun4AllReturnCodes::EVENT_OK ) return res;
170 
171  process_tracks();
172  return Fun4AllReturnCodes::EVENT_OK;
173 }
174 
175 //_____________________________________________________________________
176 int TpcDirectLaserReconstruction::End(PHCompositeNode* )
177 {
178  // save matrix container in output file
179  if( m_matrix_container )
180  {
181  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
182  outputfile->cd();
183  m_matrix_container->Write( "TpcSpaceChargeMatrixContainer" );
184  }
185 
186  // write evaluation histograms to output
187  if( m_savehistograms && m_histogramfile )
188  {
189  m_histogramfile->cd();
190  //for(const auto& o:std::initializer_list<TObject*>({ h_dca_layer, h_deltarphi_layer_south,h_deltarphi_layer_north, h_deltaz_layer, h_deltar_r,h_deltheta_delphi,h_deltheta_delphi_1,h_deltheta_delphi_2,h_deltheta_delphi_3,h_deltheta_delphi_4,h_deltheta_delphi_5,h_deltheta_delphi_6,h_deltheta_delphi_7,h_deltheta_delphi_8, h_entries,h_relangle_lasrangle,h_relangle_theta_lasrangle,h_relangle_phi_lasrangle,h_GEMs_hit,h_layers_hit,h_xy,h_xz,h_xy_pca,h_xz_pca,h_dca_path,h_zr,h_zr_pca, h_dz_z }))
191 
192  for(const auto& o:std::initializer_list<TObject*>({ h_dca_layer, h_deltarphi_layer_south,h_deltarphi_layer_north, h_deltaz_layer, h_deltar_r,h_deltheta_delphi,h_deltheta_delphi_1,h_deltheta_delphi_2,h_deltheta_delphi_3,h_deltheta_delphi_4,h_deltheta_delphi_5,h_deltheta_delphi_6,h_deltheta_delphi_7,h_deltheta_delphi_8, h_entries,h_hits,h_hits_reco,h_adc,h_adc_reco,h_adc_sum,h_adc_sum_reco,h_adc_sum_ratio,h_adc_sum_ratio_true,h_adc_vs_DCA_true,h_adc_sum_ratio_lasrangle,h_num_sum,h_num_sum_reco,h_num_sum_ratio,
193 h_num_sum_ratio_true,h_adc_sum_ratio_true,h_adc_sum_ratio,h_num_sum_ratio_lasrangle,h_bright_hits_laser1,h_bright_hits_laser2,h_bright_hits_laser3,h_bright_hits_laser4,h_relangle_lasrangle,
194 h_relangle_theta_lasrangle,h_relangle_phi_lasrangle,h_GEMs_hit,h_layers_hit,h_xy,h_xz,h_xy_pca,h_xz_pca,h_dca_path,h_zr,h_zr_pca, h_dz_z }))
195 
196  { if( o ) o->Write(); }
197  m_histogramfile->Close();
198  }
199 
200 // add these back in later !!!
201 //,h_assoc_hits
202 //h_hits
203 //h_bright_hits_laser1
204 //h_bright_hits_laser2
205 //h_bright_hits_laser3
206 //h_bright_hits_laser4
207 //h_origins
208 //h_clusters
209 
210  // print counters
211  std::cout << "TpcDirectLaserReconstruction::End - m_total_hits: " << m_total_hits << std::endl;
212  std::cout << "TpcDirectLaserReconstruction::End - m_matched_hits: " << m_matched_hits << std::endl;
213  std::cout << "TpcDirectLaserReconstruction::End - m_accepted_clusters: " << m_accepted_clusters << std::endl;
214  std::cout << "TpcDirectLaserReconstruction::End - fraction cluster/hits: " << m_accepted_clusters/m_total_hits << std::endl;
215 
216  return Fun4AllReturnCodes::EVENT_OK;
217 }
218 
219 //___________________________________________________________________________
220 void TpcDirectLaserReconstruction::SetDefaultParameters()
221 {
222 
223  // DCA cut, to decide whether a cluster should be associated to a given laser track or not (set to 5 - Charles 10/24/23)
224  set_default_double_param( "directlaser_max_dca", 1.0 );
225 
226 
227 // // residual cuts, used to decide if a given cluster is used to fill SC reconstruction matrices
228 // set_default_double_param( "directlaser_max_drphi", 0.5 );
229 // set_default_double_param( "directlaser_max_dz", 0.5 );
230 
231  set_default_double_param( "directlaser_max_drphi", 2. );
232  set_default_double_param( "directlaser_max_dz", 2. );
233 }
234 
235 //_____________________________________________________________________
236 void TpcDirectLaserReconstruction::set_grid_dimensions( int phibins, int rbins, int zbins )
237 { m_matrix_container->set_grid_dimensions( phibins, rbins, zbins ); }
238 
239 //_____________________________________________________________________
240 int TpcDirectLaserReconstruction::load_nodes( PHCompositeNode* topNode )
241 {
242  m_geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
243  assert( m_geom_container );
244 
245  // acts geometry
246  m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
247  assert( m_tGeometry );
248 
249  // tracks
250  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
251  assert(m_track_map);
252 
253  // get node containing the digitized hits
254  m_hit_map = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
255  assert(m_hit_map);
256 
257  return Fun4AllReturnCodes::EVENT_OK;
258 }
259 
260 //_____________________________________________________________________
261 void TpcDirectLaserReconstruction::create_histograms()
262 {
263  std::cout << "TpcDirectLaserReconstruction::makeHistograms - writing evaluation histograms to: " << m_histogramfilename << std::endl;
264  m_histogramfile.reset( new TFile(m_histogramfilename.c_str(), "RECREATE") );
265  m_histogramfile->cd();
266 
267  // residuals vs layers
268  h_dca_layer = new TH2F( "dca_layer", ";radius; DCA (cm)", 78, 0, 78, 500, 0, 20 );
269  h_deltarphi_layer_north = new TH2F( "deltarphi_layer_north", ";radius; r.#Delta#phi_{track-cluster} (cm)", 78, 0, 78, 2000, -5, 5 );
270  h_deltarphi_layer_south = new TH2F( "deltarphi_layer_south", ";radius; r.#Delta#phi_{track-cluster} (cm)", 78, 0, 78, 2000, -5, 5 );
271  h_deltaz_layer = new TH2F( "deltaz_layer", ";radius; #Deltaz_{track-cluster} (cm)", 78, 0, 78, 2000, -20, 20 );
272  h_deltar_r = new TH2F("deltar_r",";radius;#Deltar_{track-cluster} (cm)", 78,0,78,2000,-3,3);
273 
274  h_xy = new TH2F("h_xy"," x vs y", 320,-80,80,320,-80,80);
275  h_xz = new TH2F("h_xz"," x vs z", 320,-80,80,440,-110,110);
276  h_xy_pca = new TH2F("h_xy_pca"," x vs y pca", 320,-80,80,320,-80,80);
277  h_xz_pca = new TH2F("h_xz_pca"," x vs z pca", 320,-80,80,440,-110,110);
278  h_dca_path = new TH2F("h_dca_path"," dca vs pathlength", 440,0,110,100,0,20);
279  h_zr = new TH2F("h_zr"," z vs r", 440,-110,110,1000,28,80);
280  h_zr->GetXaxis()->SetTitle("z");
281  h_zr->GetYaxis()->SetTitle("rad");
282  h_zr_pca = new TH2F("h_zr_pca"," z vs r pca", 440,-110,110,1000,28,80);
283  h_dz_z = new TH2F("h_dz_z"," dz vs z", 440,-110,110, 1000, -20, 20);
284 
285  h_hits = new TNtuple("hits","raw hits","x:y:z:adc");
286  h_hits_reco = new TNtuple("hits_reco","raw hits","x:y:z:adc");
287 
288  h_adc = new TH1F("adc","pedestal subtracted ADC spectra of ALL lasers (MCTRUTH direction)",1063,-19.5,1042.5);
289  h_adc->SetXTitle("ADC - PEDESTAL [ADU]");
290 
291  h_adc_reco = new TH1F("adc_reco","pedestal subtracted ADC spectra of ALL lasers (RECO DIRECTION)",1063,-19.5,1042.5);
292  h_adc_reco->SetXTitle("ADC - PEDESTAL [ADU]");
293 
294  h_adc_sum = new TH1F("adc_sum","non-pedestal subtracted sum of ADC for each laser (MCTRUTH DIRECTION)",1600,-0.5,159999.5);
295  h_adc_sum->SetXTitle("#Sigma ADC [ADU]");
296 
297  h_adc_sum_reco = new TH1F("adc_sum_reco","non-pedestal subtracted sum of ADC for each laser (RECO DIRECTION)",1600,-0.5,159999.5);
298  h_adc_sum_reco->SetXTitle("#Sigma ADC [ADU]");
299 
300  h_num_sum = new TH1F("num_sum","sum of hits for each laser (MCTRUTH DIRECTION)",1600,-0.5,7999.5);
301  h_num_sum->SetXTitle("#Sigma N_{hits}^{laser X} [arb. units]");
302 
303  h_num_sum_reco = new TH1F("num_sum_reco","sum of hits for each laser (RECO DIRECTION)",1600,-0.5,7999.5);
304  h_num_sum_reco->SetXTitle("#Sigma N_{hits}^{laser X} [arb. units]");
305 
306  h_adc_sum_ratio = new TH1F("adc_sum_ratio","RATIO of non-pedestal subtracted sum of ADC for each laser (RECO DIRECTION/MCTRUTH DIRECTION)",120,-0.5,5.5);
307  h_adc_sum_ratio->SetXTitle("#frac{#Sigma^{RECO} ADC}{#Sigma^{MCTRUTH} ADC} [arb. units]");
308 
309  h_adc_sum_ratio_true = new TH1F("adc_sum_ratio_true","RATIO of non-pedestal subtracted sum of ADC for each laser (MCTRUTH DIRECTION/ALL )",120,-0.5,5.5);
310  h_adc_sum_ratio_true->SetXTitle("#frac{#Sigma^{MCTRUTH} ADC}{#Sigma^{ALL} ADC} [arb. units]");
311 
312  h_num_sum_ratio = new TH1F("num_sum_ratio","RATIO of number of hits associated to each laser (RECO DIRECTION/MCTRUTH DIRECTION)",120,-0.5,5.5);
313  h_num_sum_ratio->SetXTitle("#frac{#Sigma^{RECO} N_{hits}^{laser X}}{#Sigma^{MCTRUTH} N_{hits}^{laser X}} [arb. units]");
314 
315  h_num_sum_ratio_true = new TH1F("num_sum_ratio_true","RATIO of number of hits associated to each laser (MCTRUTH DIRECTION/ALL )",120,-0.5,5.5);
316  h_num_sum_ratio_true->SetXTitle("#frac{#Sigma^{MCTRUTH} N_{hits}^{laser X}}{#Sigma^{ALL} N_{hits}^{laser X}} [arb. units]");
317 
318  h_adc_vs_DCA_true = new TH2F("adc_vs_DCA_true","PEDESTAL SUBTRACTED ADC vs DCA for ALL HITS from ALL LASERS",100,0,100,1024,-0.5,1023.5);
319  h_adc_vs_DCA_true->SetXTitle("DCA [mm]");
320  h_adc_vs_DCA_true->SetYTitle("ADC - PEDESTAL [ADU]");
321 
322  h_adc_sum_ratio_lasrangle = new TH2F("adc_sum_ratio_lasrangle","RATIO of non-pedestal subtracted sum of ADC for LASER 1 ONLY (RECO DIRECTION/MCTRUTH DIRECTION)",91,-0.5,90.5,361,-0.5,360.5);
323  h_adc_sum_ratio_lasrangle->SetXTitle("#theta_{laser} (deg.)");
324  h_adc_sum_ratio_lasrangle->SetYTitle("#phi_{laser} (deg.)");
325 
326  h_num_sum_ratio_lasrangle = new TH2F("num_sum_ratio_lasrangle","RATIO of sum of hits for LASER 1 ONLY (RECO DIRECTION/MCTRUTH DIRECTION)",91,-0.5,90.5,361,-0.5,360.5);
327  h_num_sum_ratio_lasrangle->SetXTitle("#theta_{laser} (deg.)");
328  h_num_sum_ratio_lasrangle->SetYTitle("#phi_{laser} (deg.)");
329 
330  h_bright_hits_laser1 = new TNtuple("bright_hits_laser1","bright hits relative to laser 1","x:y:z:deltheta:delphi");
331  h_bright_hits_laser2 = new TNtuple("bright_hits_laser2","bright hits relative to laser 2","x:y:z:deltheta:delphi");
332  h_bright_hits_laser3 = new TNtuple("bright_hits_laser3","bright hits relative to laser 3","x:y:z:deltheta:delphi");
333  h_bright_hits_laser4 = new TNtuple("bright_hits_laser4","bright hits relative to laser 4","x:y:z:deltheta:delphi");
334 
335  //h_assoc_hits = new TNtuple("assoc_hits","hits associated with tracks (dca cut)","x:y:z");
336  //h_clusters = new TNtuple("clusters","associated clusters","x:y:z");
337  //h_origins = new TNtuple("origins","track origins","x:y:z");
338 
339  //h_relangle_lasrangle = new TH2F("relangle_lasrangle","Difference in pointing angle from laser and relative angle histogram for ALL lasers",51,-25.5,25.5,51,-25.5,25.5);
340  h_relangle_lasrangle = new TH2F("relangle_lasrangle","Difference in pointing angle from laser and relative angle histogram for ALL lasers",102,-25.5,25.5,102,-25.5,25.5);
341  h_relangle_lasrangle->SetXTitle("#delta#theta (deg.)");
342  h_relangle_lasrangle->SetYTitle("#delta#phi (deg.)");
343 
344  h_relangle_theta_lasrangle = new TH2F("relangle_theta_lasrangle","#delta#theta from laser and relative angle histogram for LASER 1 ONLY, ",91,-0.5,90.5,361,-0.5,360.5);
345  h_relangle_theta_lasrangle->SetXTitle("#theta_{laser} (deg.)");
346  h_relangle_theta_lasrangle->SetYTitle("#phi_{laser} (deg.)");
347 
348  h_relangle_phi_lasrangle = new TH2F("relangle_phi_lasrangle","#delta#phi from laser and relative angle histogram for LASER 1 ONLY",91,-0.5,90.5,361,-0.5,360.5);
349  h_relangle_phi_lasrangle->SetXTitle("#theta_{laser} (deg.)");
350  h_relangle_phi_lasrangle->SetYTitle("#phi_{laser} (deg.)");
351 
352  h_deltheta_delphi = new TH2F("deltheta_delphi","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and ALL laser start points", 181, -0.5,180.5,361,-0.5,360.5);
353  h_deltheta_delphi->SetXTitle("#Delta#theta");
354  h_deltheta_delphi->SetYTitle("#Delta#phi");
355 
356  h_deltheta_delphi_1 = new TH2F("deltheta_delphi_1","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 1 only", 181, -0.5,180.5,361,-0.5,360.5);
357  h_deltheta_delphi_1->SetXTitle("#Delta#theta");
358  h_deltheta_delphi_1->SetYTitle("#Delta#phi");
359 
360  h_deltheta_delphi_2 = new TH2F("deltheta_delphi_2","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 2 only", 181, -0.5,180.5,361,-0.5,360.5);
361  h_deltheta_delphi_2->SetXTitle("#Delta#theta");
362  h_deltheta_delphi_2->SetYTitle("#Delta#phi");
363 
364  h_deltheta_delphi_3 = new TH2F("deltheta_delphi_3","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 3 only", 181, -0.5,180.5,361,-0.5,360.5);
365  h_deltheta_delphi_3->SetXTitle("#Delta#theta");
366  h_deltheta_delphi_3->SetYTitle("#Delta#phi");
367 
368  h_deltheta_delphi_4 = new TH2F("deltheta_delphi_4","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 4 only", 181, -0.5,180.5,361,-0.5,360.5);
369  h_deltheta_delphi_4->SetXTitle("#Delta#theta");
370  h_deltheta_delphi_4->SetYTitle("#Delta#phi");
371 
372  h_deltheta_delphi_5 = new TH2F("deltheta_delphi_5","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 5 only", 181, -0.5,180.5,361,-0.5,360.5);
373  h_deltheta_delphi_5->SetXTitle("#Delta#theta");
374  h_deltheta_delphi_5->SetYTitle("#Delta#phi");
375 
376  h_deltheta_delphi_6 = new TH2F("deltheta_delphi_6","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 6 only", 181, -0.5,180.5,361,-0.5,360.5);
377  h_deltheta_delphi_6->SetXTitle("#Delta#theta");
378  h_deltheta_delphi_6->SetYTitle("#Delta#phi");
379 
380  h_deltheta_delphi_7 = new TH2F("deltheta_delphi_7","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 7 only", 181, -0.5,180.5,361,-0.5,360.5);
381  h_deltheta_delphi_7->SetXTitle("#Delta#theta");
382  h_deltheta_delphi_7->SetYTitle("#Delta#phi");
383 
384  h_deltheta_delphi_8 = new TH2F("deltheta_delphi_8","#Delta#theta, #Delta#phi for separation b/w TPC volume hits and LASER 8 only", 181, -0.5,180.5,361,-0.5,360.5);
385  h_deltheta_delphi_8->SetXTitle("#Delta#theta");
386  h_deltheta_delphi_8->SetYTitle("#Delta#phi");
387 
388 
389  h_GEMs_hit = new TH1F("GEMS_hit","Number of Unique GEM Modules hit for each laser",8,0,8);
390  h_layers_hit = new TH1F("layers_hit","Number of Unique Layers hit for each laser",8,8,8);
391 
392  char GEM_bin_label[128];
393  for(int GEMhistiter = 0; GEMhistiter < 8; GEMhistiter++)
394  { // (pos z) laser 1 {0,60}, laser 2 {60,0}, laser 3 {0,-60}, laser 4 {-60,0}, (neg z) laser 5 {0,60}, laser 2 {60,0}, laser 3 {0,-60}, laser 4 {-60,0}
395  sprintf(GEM_bin_label,"laser %i",GEMhistiter+1);
396  h_GEMs_hit->GetXaxis()->SetBinLabel(GEMhistiter+1,GEM_bin_label);
397  h_layers_hit->GetXaxis()->SetBinLabel(GEMhistiter+1,GEM_bin_label);
398  }
399  h_GEMs_hit->SetYTitle("Number of Unique GEM Modules Hit");
400  h_layers_hit->SetYTitle("Number of Unique Layers Hit");
401 
402  // entries vs cell grid
403  /* histogram dimension and axis limits must match that of TpcSpaceChargeMatrixContainer */
404  if( m_matrix_container )
405  {
406  int phibins = 0;
407  int rbins = 0;
408  int zbins = 0;
409  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
410  h_entries = new TH3F( "entries", ";#phi;r (cm);z (cm)", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax );
411  }
412 }
413 
414 //_____________________________________________________________________
415 void TpcDirectLaserReconstruction::process_tracks()
416 {
417  if( !( m_track_map && m_hit_map ) ) return;
418 
419  // loop over tracks and process
420  for( auto iter = m_track_map->begin(); iter != m_track_map->end(); ++iter )
421  { process_track( iter->second ); }
422 }
423 
424 //_____________________________________________________________________
425 void TpcDirectLaserReconstruction::process_track( SvtxTrack* track )
426 {
427  std::multimap<unsigned int, std::pair<float, TVector3>> cluspos_map;
428  std::set<unsigned int> layer_bin_set;
429 
430  // get track parameters
431  const TVector3 origin( track->get_x(), track->get_y(), track->get_z() );
432  const TVector3 direction( track->get_px(), track->get_py(), track->get_pz() );
433  TVector3 direction2( track->get_px(), track->get_py(), track->get_pz() );
434 
435 /*
436  if(h_origins)
437  {
438  h_origins->Fill(track->get_x(), track->get_y(), track->get_z());
439  }
440 */
441  const unsigned int trkid = track->get_id();
442 
443  const float theta_orig = origin.Theta();
444  const float phi_orig = origin.Phi();
445 
446  float theta_trk = direction.Theta();
447  direction2.RotateZ(-phi_orig);
448  float phi_trk = direction2.Phi();
449 
450  if(theta_trk > M_PI/2.) theta_trk = M_PI - theta_trk;
451 
452  //if( phi_trk < 0) phi_trk*=-1;
453 
454  while( phi_trk < m_phimin ) phi_trk += 2.*M_PI;
455  while( phi_trk >= m_phimax ) phi_trk -= 2.*M_PI;
456 
457  //const double xo = track->get_x();
458  //const double yo = track->get_y();
459  //const double zo = track->get_z();
460 
461  //if( phi_orig < 0 )phi_orig += 2.*M_PI;
462 
463  if( track->get_id() > 3 ){ phi_trk = (2 * M_PI) - phi_trk;}
464 
465  //if( Verbosity() )
466  //{
467  std::cout << "---------- processing track " << track->get_id() << std::endl;
468  std::cout << "TpcDirectLaserReconstruction::process_track - position: " << origin << " direction: " << direction << std::endl;
469  std::cout << "Position Orientation - Theta: " << theta_orig*(180./M_PI) << ", Phi: " << phi_orig*(180./M_PI) << std::endl;
470  std::cout << "Track Angle Direction - Theta: " << theta_trk*(180./M_PI) << ", Phi: " << phi_trk*(180./M_PI) << std::endl;
471  //}
472 
473  int GEM_Mod_Arr[72] = {0};
474 
475  // loop over hits
476  TrkrHitSetContainer::ConstRange hitsetrange = m_hit_map->getHitSets(TrkrDefs::TrkrId::tpcId);
477 
478  float max_adc = 0.;
479  float sum_adc_truth=0;
480  float sum_adc_truth_all=0;
481 
482  float sum_n_hits_truth=0;
483  float sum_n_hits_truth_all=0;
484 
485  for(auto hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; ++hitsetitr)
486  {
487  const TrkrDefs::hitsetkey& hitsetkey = hitsetitr->first;
488  const int side = TpcDefs::getSide(hitsetkey);
489 
490  auto hitset = hitsetitr->second;
491  const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
492  const auto layergeom = m_geom_container->GetLayerCellGeom(layer);
493  const auto layer_center_radius = layergeom->get_radius();
494 
495  // maximum drift time.
496  /* it is needed to calculate a given hit position from its drift time */
497  static constexpr double AdcClockPeriod = 53.0; // ns
498  const unsigned short NTBins = (unsigned short)layergeom->get_zbins();
499  const float tdriftmax = AdcClockPeriod * NTBins / 2.0;
500 
501  // get corresponding hits
502  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
503 
504  for (auto hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
505  {
506  ++m_total_hits;
507  sum_n_hits_truth_all++;
508 
509  const unsigned short phibin = TpcDefs::getPad(hitr->first);
510  const unsigned short zbin = TpcDefs::getTBin(hitr->first);
511 
512  const double phi = layergeom->get_phicenter(phibin);
513  const double x = layer_center_radius * cos(phi);
514  const double y = layer_center_radius * sin(phi);
515 
516  const double zdriftlength = layergeom->get_zcenter(zbin)*m_tGeometry->get_drift_velocity();
517  double z = tdriftmax*m_tGeometry->get_drift_velocity() - zdriftlength;
518  if(side == 0) z *= -1;
519 
520  const TVector3 global(x,y,z);
521 
522  float adc = (hitr->second->getAdc()) - m_pedestal;
523  float adc_unsub = (hitr->second->getAdc());
524  sum_adc_truth_all += adc_unsub;
525 /*
526  if(h_hits && trkid == 0)
527  {
528  h_hits->Fill(x,y,z,adc);
529  }
530 */
531  if( adc > max_adc ){ max_adc = adc; }
532 
533  // calculate dca
534  // origin is track origin, direction is track direction
535  bool sameside = sameSign(global.z(),origin.z());
536  const TVector3 oc( global.x()-origin.x(), global.y()-origin.y(), global.z()-origin.z() ); // vector from track origin to cluster
537  auto t = direction.Dot( oc )/square( direction.Mag() );
538  auto om = direction*t; // vector from track origin to PCA
539  const auto dca = (oc-om).Mag();
540 
541  if (sameside){ h_adc_vs_DCA_true->Fill(dca,adc);}
542 
543  //rotate displacement vector by the rotation of the coordinates:
544  //oc2.RotateY(-theta_orig * trkzdir); //undoing rotation in PHG4TpcDirectLaser
545 
546  //global2.RotateZ(-phi_orig);
547  //origin2.RotateZ(-phi_orig);
548 
549  TVector3 oc2( global.x()-origin.x(), global.y()-origin.y(), global.z()-origin.z() ); // vector from track origin to cluster
550  //const double xt = x-xo;
551  //const double yt = y-yo;
552  //const double zt = z-zo;
553  oc2.RotateZ(-phi_orig);
554 
555 
556  //auto theta1 = origin.Theta()*(180./M_PI);
557  //auto theta2 = global.Theta()*(180./M_PI);
558 
559  //auto phi1 = origin.Phi()*(180./M_PI);
560  //auto phi2 = global.Phi()*(180./M_PI);
561 
562  //float deltheta = GetRelTheta( origin.x() , origin.y(), origin.z(), global.x() , global.y() , global.z() , origin.Theta(), origin.Phi()) *(180./M_PI);
563  //float delphi = GetRelPhi(origin.x() , origin.y() , global.x() , global.y() , origin.Phi()) *(180./M_PI);
564 
565  float deltheta = oc2.Theta();
566  float delphi = oc2.Phi();
567 
568  if (deltheta > M_PI/2.) deltheta = M_PI - deltheta;
569 
570  while( delphi < m_phimin ) delphi += 2.*M_PI;
571  while( delphi >= m_phimax ) delphi -= 2.*M_PI;
572 
573  if( track->get_id() > 3 ){ delphi = (2 * M_PI) - delphi;}
574 
575  //if (delphi < 0) delphi*=-1;
576 
577  deltheta *= (180./M_PI);
578  delphi *= (180./M_PI);
579 
580 /*
581  if( trkid < 4)
582  {
583  deltheta = 180 - theta2;
584  if ( trkid == 1) delphi = fabs(phi2 - phi);
585  else delphi = phi2 - phi1;
586  }
587  else
588  {
589  deltheta = theta2;
590  if (trkid == 7 ) delphi = 360. - fabs(phi2 - phi1);
591  else delphi = fabs(phi2 - phi1);
592  }
593 */
594 
595  // relative angle histogram - only fill for hits in the same side as origin
596 
597  if(sameside){
598 
599  //std::cout<<"Trk ID = "<<trkid<<std::endl;
600 
601  h_deltheta_delphi->Fill(deltheta, delphi);
602 
603  if(track->get_id() == 0) //only fill for 1 laser !!!
604  {
605  //std::cout<<"Trk ID = "<<trkid<<std::endl;
606  h_deltheta_delphi_1->Fill(deltheta, delphi) ;
607 /*
608  //if( h_bright_hits && deltheta > 80 && deltheta < 95 && delphi > 130 && delphi < 150)
609  //filling bright_hits for theta = 75, phi = 80 (simple debugging only - to be removed)
610  if( h_bright_hits_laser1 )
611  {
612  h_bright_hits_laser1->Fill(x,y,z,deltheta,delphi);
613  }
614 */
615  }
616  if(trkid == 1) //only fill for 1 laser !!!
617  {
618  h_deltheta_delphi_2->Fill(deltheta, delphi) ;
619 /*
620  if( h_bright_hits_laser2 )
621  {
622  h_bright_hits_laser2->Fill(x,y,z,deltheta,delphi);
623  }
624 */
625  }
626  if(trkid == 2) //only fill for 1 laser !!!
627  {
628  h_deltheta_delphi_3->Fill(deltheta, delphi) ;
629 /*
630  if( h_bright_hits_laser3 )
631  {
632  h_bright_hits_laser3->Fill(x,y,z,deltheta,delphi);
633  }
634 */
635  }
636  if(trkid == 3) //only fill for 1 laser !!!
637  {
638  h_deltheta_delphi_4->Fill(deltheta, delphi) ;
639 /*
640  if( h_bright_hits_laser4 )
641  {
642  h_bright_hits_laser4->Fill(x,y,z,deltheta,delphi);
643  }
644 */
645  }
646  if(trkid == 4) //only fill for 1 laser !!!
647  {
648  h_deltheta_delphi_5->Fill(deltheta, delphi) ;
649  }
650  if(trkid == 5) //only fill for 1 laser !!!
651  {
652  h_deltheta_delphi_6->Fill(deltheta, delphi) ;
653  }
654  if(trkid == 6) //only fill for 1 laser !!!
655  {
656  h_deltheta_delphi_7->Fill(deltheta, delphi) ;
657  }
658  if(trkid == 7) //only fill for 1 laser !!!
659  {
660  h_deltheta_delphi_8->Fill(deltheta, delphi) ;
661  }
662 
663  }
664 
665  // do not associate if dca is too large
666  if( dca > m_max_dca ) continue;
667 
668  if(sameside){sum_n_hits_truth++;h_adc->Fill(adc);sum_adc_truth += adc_unsub;} //increment truth BUT ONLY for hits on the same side as origin CHARLES 10.31.23
669 
670  ++m_matched_hits;
671 /*
672  if(h_assoc_hits){
673  h_assoc_hits->Fill( x,y,z);
674  }
675 */
676  //for locating the associated hits
677  float r2 = std::sqrt((x*x) + (y*y));
678  float phi3 = phi;
679  float z2 = z;
680  while( phi3 < m_phimin ) phi3 += 2.*M_PI;
681  while( phi3 >= m_phimax ) phi3 -= 2.*M_PI;
682 
683  const int locateid = Locate(r2 , phi3 , z2); //find where the cluster is
684 
685  if( z2 > m_zmin || z2 < m_zmax ) GEM_Mod_Arr[locateid-1]++; // the array ath the cluster location - counts the number of clusters in each array, (IFF its in the volume !!!)
686 
687  // bin hits by layer
688  const auto cluspos_pair = std::make_pair(adc, global);
689  cluspos_map.insert(std::make_pair(layer, cluspos_pair));
690  layer_bin_set.insert(layer);
691  } //end looping over hits
692  } //end looping over hitset
693 
694  h_adc_sum->Fill(sum_adc_truth);
695  h_num_sum->Fill(sum_n_hits_truth);
696 
697  if(sum_adc_truth_all > 0) // you MUST have hits in this to take ratio
698  {
699  h_adc_sum_ratio_true->Fill(sum_adc_truth/sum_adc_truth_all); //for a hits associated with the same laser fire take the ratio
700  }
701 
702  if(sum_n_hits_truth_all > 0) // you MUST have hits in this to take ratio
703  {
704  h_num_sum_ratio_true->Fill(sum_n_hits_truth/sum_n_hits_truth_all); //for a hits associated with the same laser fire take the ratio
705  }
706 
707  int maxbin;
708  int deltheta_max, delphi_max, dummy_z;
709 
710  float theta_reco =0;
711  float phi_reco =0;
712 
713  float direction_reco;
714  if( trkid < 4 ){direction_reco = -1;} //first four lasers are on positive z readout plane, and shoot towards negative z
715  else{direction_reco = 1;} // next four lasers are on negative z readout plane and shoot towards positive z
716 
717  TVector3 dir(0, 0, direction_reco); //unit vector starts out pointing along z axis
718 
719  if( (trkid == 0 && h_deltheta_delphi_1->GetEntries() > 0) && max_adc > 10)
720  {
721  h_deltheta_delphi_1->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
722  h_deltheta_delphi_1->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
723 
724  maxbin= h_deltheta_delphi_1->GetMaximumBin();
725 
726  h_deltheta_delphi_1->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
727 
728  std::cout << "Laser # "<<(trkid+1)<<" naieve deltheta max: "<< h_deltheta_delphi_1->GetXaxis()->GetBinCenter(deltheta_max) <<", naieve delphi max: "<< h_deltheta_delphi_1->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
729 
730  h_relangle_lasrangle->Fill(h_deltheta_delphi_1->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_1->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
731 
732  h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_1->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
733  h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_1->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
734 
735  theta_reco = h_deltheta_delphi_1->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
736  phi_reco = h_deltheta_delphi_1->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
737  }
738 
739  else if( (trkid == 1 && h_deltheta_delphi_2->GetEntries() > 0) && max_adc > 10 )
740  {
741  h_deltheta_delphi_2->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
742  h_deltheta_delphi_2->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
743 
744  maxbin= h_deltheta_delphi_2->GetMaximumBin();
745 
746  h_deltheta_delphi_2->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
747 
748  std::cout << "Laser # "<<(trkid+1)<<" naieve deltheta max: "<< h_deltheta_delphi_2->GetXaxis()->GetBinCenter(deltheta_max) <<", naive delphi max: "<< h_deltheta_delphi_2->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
749 
750  h_relangle_lasrangle->Fill(h_deltheta_delphi_2->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_2->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
751 
752  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_2->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
753  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_2->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
754 
755  theta_reco = h_deltheta_delphi_2->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
756  phi_reco = h_deltheta_delphi_2->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
757  }
758 
759  else if( (trkid == 2 && h_deltheta_delphi_3->GetEntries() > 0) && max_adc > 10 )
760  {
761  h_deltheta_delphi_3->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
762  h_deltheta_delphi_3->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
763 
764  maxbin= h_deltheta_delphi_3->GetMaximumBin();
765 
766  h_deltheta_delphi_3->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
767 
768  std::cout << "Laser # "<<(trkid+1)<<" naieve deltheta max: "<< h_deltheta_delphi_3->GetXaxis()->GetBinCenter(deltheta_max) <<", naive delphi max: "<< h_deltheta_delphi_3->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
769 
770  h_relangle_lasrangle->Fill(h_deltheta_delphi_3->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_3->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
771 
772  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_3->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
773  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_3->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
774 
775  theta_reco = h_deltheta_delphi_3->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
776  phi_reco = h_deltheta_delphi_3->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
777  }
778 
779  else if( (trkid == 3 && h_deltheta_delphi_4->GetEntries() > 0) && max_adc > 10 )
780  {
781  h_deltheta_delphi_4->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
782  h_deltheta_delphi_4->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
783 
784  maxbin= h_deltheta_delphi_4->GetMaximumBin();
785 
786  h_deltheta_delphi_4->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
787 
788  std::cout << "Laser # "<<(trkid+1)<<" naieve deltheta max: "<< h_deltheta_delphi_4->GetXaxis()->GetBinCenter(deltheta_max) <<", naieve delphi max: "<< h_deltheta_delphi_4->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
789 
790  h_relangle_lasrangle->Fill(h_deltheta_delphi_4->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_4->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
791 
792  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_4->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
793  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_4->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
794 
795  theta_reco = h_deltheta_delphi_4->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
796  phi_reco = h_deltheta_delphi_4->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
797  }
798 
799  else if( (trkid == 4 && h_deltheta_delphi_5->GetEntries() > 0) && max_adc > 10 )
800  {
801  h_deltheta_delphi_5->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
802  h_deltheta_delphi_5->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
803 
804  maxbin= h_deltheta_delphi_5->GetMaximumBin();
805 
806  h_deltheta_delphi_5->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
807 
808  std::cout << "Laser # "<<(trkid+1)<<" deltheta max: "<< h_deltheta_delphi_5->GetXaxis()->GetBinCenter(deltheta_max) <<", delphi max: "<< h_deltheta_delphi_5->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
809 
810  h_relangle_lasrangle->Fill(h_deltheta_delphi_5->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_5->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
811 
812  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_5->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
813  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_5->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
814 
815  theta_reco = h_deltheta_delphi_5->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
816  phi_reco = h_deltheta_delphi_5->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
817  }
818 
819  else if( (trkid == 5 && h_deltheta_delphi_6->GetEntries() > 0) && max_adc > 10 )
820  {
821  h_deltheta_delphi_6->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
822  h_deltheta_delphi_6->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
823 
824  maxbin= h_deltheta_delphi_6->GetMaximumBin();
825 
826  h_deltheta_delphi_6->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
827 
828  std::cout << "Laser # "<<(trkid+1)<<" deltheta max: "<< h_deltheta_delphi_6->GetXaxis()->GetBinCenter(deltheta_max) <<", delphi max: "<< h_deltheta_delphi_6->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
829 
830  h_relangle_lasrangle->Fill(h_deltheta_delphi_6->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_6->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
831 
832  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_6->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
833  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_6->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
834 
835  theta_reco = h_deltheta_delphi_6->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
836  phi_reco = h_deltheta_delphi_6->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
837  }
838 
839  else if( (trkid == 6 && h_deltheta_delphi_7->GetEntries() > 0) && max_adc > 10 )
840  {
841  h_deltheta_delphi_7->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
842  h_deltheta_delphi_7->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
843 
844  maxbin= h_deltheta_delphi_7->GetMaximumBin();
845 
846  h_deltheta_delphi_7->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
847 
848  std::cout << "Laser # "<<(trkid+1)<<" deltheta max: "<< h_deltheta_delphi_7->GetXaxis()->GetBinCenter(deltheta_max) <<", delphi max: "<< h_deltheta_delphi_7->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
849 
850  h_relangle_lasrangle->Fill(h_deltheta_delphi_7->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_7->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
851 
852  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_7->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
853  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_7->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
854 
855  theta_reco = h_deltheta_delphi_7->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
856  phi_reco = h_deltheta_delphi_7->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
857  }
858 
859  else if( (trkid == 7 && h_deltheta_delphi_8->GetEntries() > 0) && max_adc > 10 )
860  {
861  h_deltheta_delphi_8->GetXaxis()->SetRangeUser(-5+(theta_trk*(180./M_PI)),5+(theta_trk*(180./M_PI)));
862  h_deltheta_delphi_8->GetYaxis()->SetRangeUser(-5+(phi_trk*(180./M_PI)),5+(phi_trk*(180./M_PI)));
863 
864  maxbin= h_deltheta_delphi_8->GetMaximumBin();
865 
866  h_deltheta_delphi_8->GetBinXYZ(maxbin, deltheta_max, delphi_max, dummy_z);
867 
868  std::cout << "Laser # "<<(trkid+1)<<" deltheta max: "<< h_deltheta_delphi_8->GetXaxis()->GetBinCenter(deltheta_max) <<", delphi max: "<< h_deltheta_delphi_8->GetYaxis()->GetBinCenter(delphi_max) << std::endl;
869 
870  h_relangle_lasrangle->Fill(h_deltheta_delphi_8->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)),h_deltheta_delphi_8->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)));
871 
872  //h_relangle_theta_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_8->GetXaxis()->GetBinCenter(deltheta_max)-(theta_trk*(180./M_PI)) );
873  //h_relangle_phi_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),h_deltheta_delphi_8->GetYaxis()->GetBinCenter(delphi_max)-(phi_trk*(180./M_PI)) );
874 
875  theta_reco = h_deltheta_delphi_8->GetXaxis()->GetBinCenter(deltheta_max) * (M_PI/180); //convert to radians
876  phi_reco = h_deltheta_delphi_8->GetYaxis()->GetBinCenter(delphi_max) * (M_PI/180); //convert to radians
877  }
878 
879  // Determining Vector for reconstructed laser direction
880  dir.RotateY(theta_reco* direction_reco);
881  if(direction_reco == -1) dir.RotateZ(phi_reco); //if +z facing -z
882  else dir.RotateZ(-phi_reco); //if -z facting +z
883  dir.RotateZ(phi_orig); //rotate by the laser coordinate
884 
886 
887  TVector3 direction2_reco( dir.Px(), dir.Py(), dir.Pz() );
888 
889  float theta_trk_reco = dir.Theta();
890  direction2_reco.RotateZ(-phi_orig);
891  float phi_trk_reco = direction2_reco.Phi();
892 
893  if(theta_trk_reco > M_PI/2.) theta_trk_reco = M_PI - theta_trk_reco;
894 
895  //if( phi_trk < 0) phi_trk*=-1;
896 
897  while( phi_trk_reco < m_phimin ) phi_trk_reco += 2.*M_PI;
898  while( phi_trk_reco >= m_phimax ) phi_trk_reco -= 2.*M_PI;
899 
900  //const double xo = track->get_x();
901  //const double yo = track->get_y();
902  //const double zo = track->get_z();
903 
904  //if( phi_orig < 0 )phi_orig += 2.*M_PI;
905 
906  if( track->get_id() > 3 ){ phi_trk_reco = (2 * M_PI) - phi_trk_reco;}
907 
908  std::cout << "RECONSTRUCTED Track Angle Direction - Theta: " << theta_trk_reco*(180./M_PI) << ", Phi: " << phi_trk_reco*(180./M_PI) << std::endl;
909 
911 
912  //now loop over hits AGAIN and get ADC spectrum
913 
914  float sum_adc_reco=0;
915  float sum_n_hits_reco=0;
916 
917  for(auto hitsetitr = hitsetrange.first; hitsetitr != hitsetrange.second; ++hitsetitr)
918  {
919  const TrkrDefs::hitsetkey& hitsetkey_2 = hitsetitr->first;
920  const int side_2 = TpcDefs::getSide(hitsetkey_2);
921 
922  auto hitset_2 = hitsetitr->second;
923  const unsigned int layer_2 = TrkrDefs::getLayer(hitsetkey_2);
924  const auto layergeom_2 = m_geom_container->GetLayerCellGeom(layer_2);
925  const auto layer_center_radius_2 = layergeom_2->get_radius();
926 
927  // maximum drift time.
928  /* it is needed to calculate a given hit position from its drift time */
929  static constexpr double AdcClockPeriod_2 = 53.0; // ns
930  const unsigned short NTBins_2 = (unsigned short)layergeom_2->get_zbins();
931  const float tdriftmax_2 = AdcClockPeriod_2 * NTBins_2 / 2.0;
932 
933  // get corresponding hits
934  TrkrHitSet::ConstRange hitrangei_2 = hitset_2->getHits();
935 
936  for (auto hitr = hitrangei_2.first; hitr != hitrangei_2.second; ++hitr)
937  {
938 
939  const unsigned short phibin_2 = TpcDefs::getPad(hitr->first);
940  const unsigned short zbin_2 = TpcDefs::getTBin(hitr->first);
941 
942  const double phi_2 = layergeom_2->get_phicenter(phibin_2);
943  const double x_2 = layer_center_radius_2 * cos(phi_2);
944  const double y_2 = layer_center_radius_2 * sin(phi_2);
945 
946  const double zdriftlength_2 = layergeom_2->get_zcenter(zbin_2)*m_tGeometry->get_drift_velocity();
947  double z_2 = tdriftmax_2*m_tGeometry->get_drift_velocity() - zdriftlength_2;
948  if(side_2 == 0) z_2 *= -1;
949 
950  const TVector3 global_2(x_2,y_2,z_2);
951 
952  bool sameside_reco = sameSign(global_2.z(),origin.z());
953 
954  float adc_2 = (hitr->second->getAdc()) - m_pedestal;
955  float adc_2_unsub = (hitr->second->getAdc());
956 
957  // calculate dca
958  // origin is track origin, direction is track direction
959  const TVector3 oc_2( global_2.x()-origin.x(), global_2.y()-origin.y(), global_2.z()-origin.z() ); // vector from track origin to cluster
960  auto t_2 = dir.Dot( oc_2 )/square( dir.Mag() );
961  auto om_2 = dir*t_2; // vector from track origin to PCA
962  const auto dca_2 = (oc_2-om_2).Mag();
963 
964  if( dca_2 > m_max_dca ) continue; // do not record if hit is outside DCA, proceed to next hit
965  if(sameside_reco){sum_n_hits_reco++;h_adc_reco->Fill(adc_2);sum_adc_reco += adc_2_unsub;} //increment reco but only if hits are on the same side
966 /*
967  if(h_hits_reco && trkid == 0)
968  {
969  h_hits_reco->Fill(x_2,y_2,z_2,adc_2);
970  }
971 */
972  } //end loop over hits again
973  } //end loop over hitset again
974 
975  h_adc_sum_reco->Fill(sum_adc_reco);
976  h_num_sum_reco->Fill(sum_n_hits_reco);
977 
978  if(sum_adc_truth > 0) // you MUST have hits in this to take ratio
979  {
980  h_adc_sum_ratio->Fill(sum_adc_reco/sum_adc_truth); //for a hits associated with the same laser fire take the ratio
981  if( trkid == 0 ) //for laser 1 only, show me WHERE we are messing up reconstruction
982  {
983  h_adc_sum_ratio_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),(sum_adc_reco/sum_adc_truth));
984  }
985  }
986 
987  if(sum_n_hits_truth > 0) // you MUST have hits in this to take ratio
988  {
989  h_num_sum_ratio->Fill(sum_n_hits_reco/sum_n_hits_truth); //for a hits associated with the same laser fire take the ratio
990  if( trkid == 0 ) //for laser 1 only, show me WHERE we are messing up reconstruction
991  {
992  h_num_sum_ratio_lasrangle->Fill(theta_trk*(180./M_PI),phi_trk*(180./M_PI),(sum_n_hits_reco/sum_n_hits_truth));
993  }
994  }
995 
996  for(int GEMS_iter = 0; GEMS_iter < 72; GEMS_iter++)
997  {
998  if(GEM_Mod_Arr[GEMS_iter] > 0){ //laser 0
999  h_GEMs_hit->Fill(trkid + 0.5);
1000  }
1001  }
1002 
1003  h_deltheta_delphi_1->Reset("ICESM");
1004  h_deltheta_delphi_2->Reset("ICESM");
1005  h_deltheta_delphi_3->Reset("ICESM");
1006  h_deltheta_delphi_4->Reset("ICESM");
1007  h_deltheta_delphi_5->Reset("ICESM");
1008  h_deltheta_delphi_6->Reset("ICESM");
1009  h_deltheta_delphi_7->Reset("ICESM");
1010  h_deltheta_delphi_8->Reset("ICESM");
1011 
1012 
1013 
1014  // We will bring this back here when we fix the clustering
1015  //int GEM_Mod_Arr[72] = {0};
1016 
1017  // we have all of the residuals for this track added to the map, binned by layer
1018  // now we calculate the centroid of the clusters in each bin
1019 
1020  for(auto layer : layer_bin_set)
1021  {
1022  PHG4TpcCylinderGeom *layergeom = m_geom_container->GetLayerCellGeom(layer);
1023  const auto layer_center_radius = layergeom->get_radius();
1024  const auto layer_inner_radius = layer_center_radius - layergeom->get_thickness() / 2.0;
1025  const auto layer_outer_radius = layer_center_radius + layergeom->get_thickness() / 2.0;
1026 
1027  h_layers_hit->Fill(trkid + 0.5);
1028 
1029  // does track pass completely through this layer? If not do not use the hits
1030  auto tupdn = line_circle_intersection(origin, direction, layer_outer_radius);
1031  if( tupdn.first <= 0 && tupdn.second <= 0)
1032  {
1033  std::cout << " punt: layer " << layer << " layer outer radius " << layer_outer_radius << " tup " << tupdn.first << " tdn " << tupdn.second << std::endl;
1034  continue;
1035  }
1036  double layer_entry = tupdn.first;
1037  if(tupdn.second >= 0 && tupdn.second < tupdn.first) layer_entry = tupdn.second;
1038 
1039  tupdn = line_circle_intersection(origin, direction, layer_inner_radius);
1040  if( tupdn.first <= 0 && tupdn.second <= 0)
1041  {
1042  std::cout << " punt: layer " << layer << " layer inner radius " << layer_inner_radius << " tup " << tupdn.first << " tdn " << tupdn.second << std::endl;
1043  continue;
1044  }
1045  double layer_exit = tupdn.first;
1046  if(tupdn.second > 0 && tupdn.second < tupdn.first) layer_exit = tupdn.second;
1047  if(Verbosity() > 2) std::cout << " layer " << layer << " layer entry " << layer_entry << " layer exit " << layer_exit << std::endl;
1048 
1049  // calculate track intersection with layer center
1050  tupdn = line_circle_intersection(origin, direction, layer_center_radius);
1051  double tup = tupdn.first;
1052  double tdn = tupdn.second;
1053  // take 1 one if there are two intersections
1054  double t = tup;
1055  if(tdn > 0 && tdn < tup) t = tdn;
1056  if( t < 0 )
1057  {
1058  std::cout << " punt: layer " << layer << " layer center radius " << layer_center_radius << " t " << t << " tup " << tupdn.first << " tdn " << tupdn.second << std::endl;
1059  continue;
1060  }
1061 
1062  // get position of layer center on track
1063  auto om = direction*t;
1064 
1065  // get vector pointing to the track intersection with layer center
1066  const auto projection = origin + om;
1067 
1068  double zmax = -999.;
1069  double zmin = 999.;
1070 
1071  auto bin_residual = cluspos_map.equal_range(layer);
1072 
1073  TVector3 clus_centroid(0,0,0);
1074  float wt = 0;
1075  for( auto mit = bin_residual.first; mit != bin_residual.second; ++mit)
1076  {
1077  float adc = (float) mit->second.first;
1078  TVector3 cluspos =mit->second.second;
1079 
1080  if(fabs(cluspos.z() - projection.z()) > m_max_zrange) continue; // to reject clusters from possible second traverse of layer
1081 
1082  if(Verbosity() > 2)
1083  {
1084  std::cout << " layer " << layer << " adc " << adc << std::endl;
1085  std::cout << " cluspos " << cluspos.x() << " " << cluspos.y() << " " << cluspos.z()
1086  << " clus radius " << sqrt(square(cluspos.x()) + square(cluspos.y())) << std::endl;
1087  }
1088 
1089  clus_centroid += cluspos * adc;
1090  wt += adc;
1091 
1092  if(cluspos.z() < zmin) zmin = cluspos.z();
1093  if(cluspos.z() > zmax) zmax = cluspos.z();
1094  }
1095 
1096  clus_centroid.SetX(clus_centroid.x()/ wt);
1097  clus_centroid.SetY(clus_centroid.y()/ wt);
1098  clus_centroid.SetZ(clus_centroid.z()/ wt);
1099 
1100  double zrange = zmax-zmin;
1101  if(zrange > m_max_zrange)
1102  {
1103  std::cout << " exceeded max zrange: zrange " << zrange << " max zrange " << m_max_zrange << std::endl;
1104  continue;
1105  }
1106 
1107  // get the distance of the cluster centroid to the track-layer intersection point
1108 
1109  const TVector3 oc( clus_centroid.x()-origin.x(), clus_centroid.y()-origin.y(), clus_centroid.z()-origin.z() ); // vector from track origin to cluster
1110 
1111  const auto dca = (oc-om).Mag();
1112 
1113  // path length
1114  const auto pathlength = om.Mag();
1115 
1116  // Correct cluster z for the track transit time using the pathlength
1117  double ns_per_cm = 1e9 / 3e10;
1118  double dt = pathlength * ns_per_cm;
1119  double transit_dz = dt * m_tGeometry->get_drift_velocity();
1120  if(origin.z() > 0)
1121  clus_centroid.SetZ(clus_centroid.z() + transit_dz);
1122  else
1123  clus_centroid.SetZ(clus_centroid.z() - transit_dz);
1124 
1125  if(Verbosity() > 0)
1126  {
1127  std::cout << " layer " << layer << " radius " << layer_center_radius << " wt " << wt << " t " << t << " dca " << dca
1128  << " pathlength " << pathlength << " transit_dz " << transit_dz << std::endl;
1129  std::cout << " clus_centroid " << clus_centroid.x() << " " << clus_centroid.y() << " " << clus_centroid.z() << " zrange " << zrange << std::endl;
1130  std::cout << " projection " << projection.x() << " " << projection.y() << " " << projection.z() << " dz " << clus_centroid.z() - projection.z() << std::endl;
1131  }
1132 
1133  // create relevant state vector and assign to track
1134  SvtxTrackState_v1 state( pathlength );
1135  state.set_x( projection.x() );
1136  state.set_y( projection.y() );
1137  state.set_z( projection.z() );
1138 
1139  state.set_px( direction.x());
1140  state.set_py( direction.y());
1141  state.set_pz( direction.z());
1142  track->insert_state( &state );
1143 
1144  // also associate cluster to track
1145  //track->insert_cluster_key( key );
1146 
1147  // cluster r, phi and z
1148  const auto cluster_r = get_r(clus_centroid.x(), clus_centroid.y());
1149  const auto cluster_phi = std::atan2(clus_centroid.y(),clus_centroid.x());
1150  const auto cluster_z = clus_centroid.z();
1151 
1152  // We will bring this back when we fix clustering !!
1153  /*
1154  float r2 = cluster_r;
1155  float phi2 = cluster_phi;
1156  float z2 = cluster_z;
1157  while( phi2 < m_phimin ) phi2 += 2.*M_PI;
1158  while( phi2 >= m_phimax ) phi2 -= 2.*M_PI;
1159 
1160  const int locateid = Locate(r2 , phi2 , z2); //find where the cluster is
1161 
1162  if( z2 > m_zmin || z2 < m_zmax ) GEM_Mod_Arr[locateid-1]++; // the array ath the cluster location - counts the number of clusters in each array, (IFF its in the volume !!!)
1163  */
1164 
1165  // cluster errors
1166  const auto cluster_rphi_error = 0.015;
1167  const auto cluster_z_error = 0.075;
1168 
1169  // track position
1170  const auto track_phi = std::atan2( projection.y(), projection.x() );
1171  const auto track_z = projection.z();
1172 
1173  // track angles
1174  const auto cosphi( std::cos( track_phi ) );
1175  const auto sinphi( std::sin( track_phi ) );
1176  const auto track_pphi = -state.get_px()*sinphi + state.get_py()*cosphi;
1177  const auto track_pr = state.get_px()*cosphi + state.get_py()*sinphi;
1178  const auto track_pz = state.get_pz();
1179  const auto talpha = -track_pphi/track_pr;
1180  const auto tbeta = -track_pz/track_pr;
1181 
1182  // sanity check
1183  if( std::isnan(talpha) )
1184  {
1185  std::cout << "TpcDirectLaserReconstruction::process_track - talpha is nan" << std::endl;
1186  continue;
1187  }
1188 
1189  if( std::isnan(tbeta) )
1190  {
1191  std::cout << "TpcDirectLaserReconstruction::process_track - tbeta is nan" << std::endl;
1192  continue;
1193  }
1194 
1195  // residuals
1196  const auto drp = cluster_r*delta_phi( cluster_phi - track_phi );
1197  const auto dz = cluster_z - track_z;
1198 
1199  // sanity checks
1200  if( std::isnan(drp) )
1201  {
1202  std::cout << "TpcDirectLaserReconstruction::process_track - drp is nan" << std::endl;
1203  continue;
1204  }
1205 
1206  if( std::isnan(dz) )
1207  {
1208  std::cout << "TpcDirectLaserReconstruction::process_track - dz is nan" << std::endl;
1209  continue;
1210  }
1211 
1212  if(m_savehistograms)
1213  {
1214  const float r = get_r( projection.x(), projection.y() );
1215  const float dr = cluster_r - r;
1216  if(h_dca_layer) h_dca_layer->Fill(r, dca);
1217  if(h_deltarphi_layer_south && clus_centroid.z() < 0) h_deltarphi_layer_south->Fill(r, drp);
1218  if(h_deltarphi_layer_north && clus_centroid.z() > 0) h_deltarphi_layer_north->Fill(r, drp);
1219  if(h_deltaz_layer) h_deltaz_layer->Fill(r, dz);
1220  if(h_deltar_r) h_deltar_r->Fill(r, dr);
1221  if(h_entries)
1222  {
1223  auto phi = cluster_phi;
1224  while( phi < m_phimin ) phi += 2.*M_PI;
1225  while( phi >= m_phimax ) phi -= 2.*M_PI;
1226  h_entries->Fill( phi, cluster_r, cluster_z );
1227  }
1228  if(h_xy) h_xy->Fill(clus_centroid.x(), clus_centroid.y());
1229  if(h_xz) h_xz->Fill(clus_centroid.x(), clus_centroid.z());
1230  if(h_xy_pca) h_xy_pca->Fill(projection.x(), projection.y());
1231  if(h_xz_pca) h_xz_pca->Fill(projection.x(), projection.z());
1232  if(h_dca_path) h_dca_path->Fill(pathlength,dca);
1233  if(h_zr) h_zr->Fill(clus_centroid.z(), cluster_r);
1234  if(h_zr_pca) h_zr_pca->Fill(projection.z(), r);
1235  if(h_dz_z)h_dz_z->Fill(projection.z(), clus_centroid.z()- projection.z());
1236  //if(h_clusters)h_clusters->Fill(clus_centroid.x(),clus_centroid.y(),clus_centroid.z());
1237  }
1238 
1239  // // check against limits
1240  // if( std::abs( drp ) > m_max_drphi ) continue;
1241  // if( std::abs( dz ) > m_max_dz ) continue;
1242 
1243  // residual errors squared
1244  const auto erp = square(cluster_rphi_error);
1245  const auto ez = square(cluster_z_error);
1246 
1247  // sanity check
1248  if( std::isnan( erp ) )
1249  {
1250  std::cout << "TpcDirectLaserReconstruction::process_track - erp is nan" << std::endl;
1251  continue;
1252  }
1253 
1254  if( std::isnan( ez ) )
1255  {
1256  std::cout << "TpcDirectLaserReconstruction::process_track - ez is nan" << std::endl;
1257  continue;
1258  }
1259 
1260  // get cell
1261  const auto i = get_cell_index( clus_centroid );
1262  if( i < 0 )
1263  {
1264  if( Verbosity() )
1265  {
1266  std::cout << "TpcDirectLaserReconstruction::process_track - invalid cell index"
1267  << " r: " << cluster_r
1268  << " phi: " << cluster_phi
1269  << " z: " << cluster_z
1270  << std::endl;
1271  }
1272  continue;
1273  }
1274 
1275  // update matrices
1276  // see https://indico.bnl.gov/event/7440/contributions/43328/attachments/31334/49446/talk.pdf for details
1277  m_matrix_container->add_to_lhs(i, 0, 0, 1./erp );
1278  m_matrix_container->add_to_lhs(i, 0, 1, 0 );
1279  m_matrix_container->add_to_lhs(i, 0, 2, talpha/erp );
1280 
1281  m_matrix_container->add_to_lhs(i, 1, 0, 0 );
1282  m_matrix_container->add_to_lhs(i, 1, 1, 1./ez );
1283  m_matrix_container->add_to_lhs(i, 1, 2, tbeta/ez );
1284 
1285  m_matrix_container->add_to_lhs(i, 2, 0, talpha/erp );
1286  m_matrix_container->add_to_lhs(i, 2, 1, tbeta/ez );
1287  m_matrix_container->add_to_lhs(i, 2, 2, square(talpha)/erp + square(tbeta)/ez );
1288 
1289  m_matrix_container->add_to_rhs(i, 0, drp/erp );
1290  m_matrix_container->add_to_rhs(i, 1, dz/ez );
1291  m_matrix_container->add_to_rhs(i, 2, talpha*drp/erp + tbeta*dz/ez );
1292 
1293  // update entries in cell
1294  m_matrix_container->add_to_entries(i);
1295 
1296  // increment number of accepted clusters
1297  ++m_accepted_clusters;
1298  }
1299 
1300 // We will eventually bring this back when we fix clustering
1301 /*
1302  for(int GEMS_iter = 0; GEMS_iter < 72; GEMS_iter++)
1303  {
1304  if(GEM_Mod_Arr[GEMS_iter] > 0){ //laser 0
1305  h_GEMs_hit->Fill(trkid + 0.5);
1306  }
1307  }
1308 */
1309 
1310 }
1311 
1312 
1313 //_____________________________________________________________________
1314 int TpcDirectLaserReconstruction::get_cell_index( const TVector3& global ) const
1315 {
1316  // get grid dimensions from matrix container
1317  int phibins = 0;
1318  int rbins = 0;
1319  int zbins = 0;
1320  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
1321 
1322  // phi
1323  // bound check
1324  float phi = std::atan2( global.y(), global.x() );
1325  while( phi < m_phimin ) phi += 2.*M_PI;
1326  while( phi >= m_phimax ) phi -= 2.*M_PI;
1327  int iphi = phibins*(phi-m_phimin)/(m_phimax-m_phimin);
1328 
1329  // radius
1330  const float r = get_r( global.x(), global.y() );
1331  if( r < m_rmin || r >= m_rmax ) return -1;
1332  int ir = rbins*(r-m_rmin)/(m_rmax-m_rmin);
1333 
1334  // z
1335  const float z = global.z();
1336  if( z < m_zmin || z >= m_zmax ) return -1;
1337  int iz = zbins*(z-m_zmin)/(m_zmax-m_zmin);
1338 
1339  return m_matrix_container->get_cell_index( iphi, ir, iz );
1340 }
1341 //_____________________________________________________________________
1342 int TpcDirectLaserReconstruction::Locate(float r , float phi , float z)
1343 {
1344 
1346  // North Label Conventions //
1348  // 1 - 00_R1 16 - 05_R1 31 - 10_R1
1349  // 2 - 00_R2 17 - 05_R2 32 - 10_R2
1350  // 3 - 00_R3 18 - 05_R3 33 - 10_R3
1351  // 4 - 01_R1 19 - 06_R1 34 - 11_R1
1352  // 5 - 01_R2 20 - 06_R2 35 - 11_R2
1353  // 6 - 01_R3 21 - 06_R3 36 - 11_R3
1354  // 7 - 02_R1 22 - 07_R1
1355  // 8 - 02_R2 23 - 07_R2
1356  // 9 - 02_R3 24 - 07_R3
1357  // 10 - 03_R1 25 - 08_R1
1358  // 11 - 03_R2 26 - 08_R2
1359  // 12 - 03_R3 27 - 08_R3
1360  // 13 - 04_R1 28 - 09_R1
1361  // 14 - 04_R2 29 - 09_R2
1362  // 15 - 04_R3 30 - 09_R3
1363 
1365  // South Label Conventions //
1367  // 37 - 12_R1 52 - 17_R1 67 - 22_R1
1368  // 38 - 12_R2 53 - 17_R2 68 - 22_R2
1369  // 39 - 12_R3 54 - 17_R3 69 - 22_R3
1370  // 40 - 13_R1 55 - 18_R1 70 - 23_R1
1371  // 41 - 13_R2 56 - 18_R2 71 - 23_R2
1372  // 42 - 13_R3 57 - 18_R3 72 - 23_R3
1373  // 43 - 14_R1 58 - 19_R1
1374  // 44 - 14_R2 59 - 19_R2
1375  // 45 - 14_R3 60 - 19_R3
1376  // 46 - 15_R1 61 - 20_R1
1377  // 47 - 15_R2 62 - 20_R2
1378  // 48 - 15_R3 63 - 20_R3
1379  // 49 - 16_R1 64 - 21_R1
1380  // 50 - 16_R2 65 - 21_R2
1381  // 51 - 16_R3 66 - 21_R3
1382 
1383 
1384  int GEM_Mod_ID; //integer from 1 to 72
1385 
1386  const float Angle_Bins[13] = {23*M_PI/12,1*M_PI/12,3*M_PI/12,5*M_PI/12,7*M_PI/12,9*M_PI/12,11*M_PI/12,13*M_PI/12,15*M_PI/12,17*M_PI/12,19*M_PI/12,21*M_PI/12,23*M_PI/12}; //12 angle bins on each side
1387 
1388  const float r_bins[4] = {30,46,62,78}; //4 r bins
1389 
1390  int angle_id = 0; // default to 0
1391  int r_id = 0; //default to 0
1392  int side_id = 0; //default to 0 (north side)
1393 
1394  for(int r_iter = 0; r_iter < 3; r_iter++ ){
1395  if( r > r_bins[r_iter] && r < r_bins[r_iter+1] )
1396  {
1397  r_id = r_iter+1;
1398  break; //break out of the for loop if you found where the hit is in r
1399  }
1400  }
1401 
1402  for(int ang_iter = 0; ang_iter < 12; ang_iter++ ){
1403  if( phi > Angle_Bins[ang_iter] && phi < Angle_Bins[ang_iter+1] )
1404  {
1405  angle_id = ang_iter;
1406  break; //break out the for loop if you found where the hit is in phi
1407  }
1408  }
1409  if(z < 0)side_id=1; //(south side)
1410 
1411  return GEM_Mod_ID = (36*side_id) + (3*angle_id) + r_id;
1412 }
1413 
1414 float TpcDirectLaserReconstruction::GetRelPhi( float xorig, float yorig, float x, float y, float phiorig )
1415 {
1416 // getting the relative phi in the frame rotated to the ith laser origin
1417 //
1418 // inputs:
1419 // xorig = x coordinate of ith laser origin
1420 // yorig = y coordinate of ith laser origin
1421 // x = x coordinate of arbitrary point in TPC volume
1422 // y = y coordinate of arbitrary point in TPC volume
1423 // phiorig = phi (from + x axis) of ith laser origin
1424 //
1425 // output:
1426 // relphi = azimuthal displacement of arbitrary point in TPC volume from ith laser origin
1427 
1428  if (phiorig < 0) phiorig += 2. * M_PI;
1429  else if (phiorig > 2.*M_PI) phiorig -= 2. * M_PI;
1430 
1431  float dx = x - xorig;
1432  float dy = y - yorig;
1433  float relphi = atan2(dy, dx) - phiorig;
1434  if (relphi < 0) relphi += 2. * M_PI;
1435  else if (relphi > 2.*M_PI) relphi -= 2. * M_PI;
1436 
1437  return relphi;
1438 }
1439 
1440 float TpcDirectLaserReconstruction::GetRelTheta( float xorig, float yorig, float zorig, float x, float y, float z, float thetaorig, float phiorig )
1441 {
1442 // getting the relative theta in the frame rotated to the ith laser origin
1443 //
1444 // inputs:
1445 // xorig = x coordinate of ith laser origin
1446 // yorig = y coordinate of ith laser origin
1447 // zorig = z cooridante of ith laser origin
1448 // x = x coordinate of arbitrary point in TPC volume
1449 // y = y coordinate of arbitrary point in TPC volume
1450 // z = z coordinate of arbitrary point in TPC volume
1451 // thetaorig = theta (from + y axis?) of ithe laser origin
1452 // phiorig = phi (from + x axis) of ith laser origin
1453 //
1454 // output:
1455 // reltheta = polar displacement of arbitrary point in TPC volume from ith laser origin
1456 
1457  if (phiorig < 0) phiorig += 2. * M_PI;
1458  else if (phiorig > 2.*M_PI) phiorig -= 2. * M_PI;
1459 
1460  if (thetaorig < 0) thetaorig += 2. * M_PI;
1461  else if (thetaorig > 2.*M_PI) thetaorig -= 2. * M_PI;
1462 
1463 
1464  float dx = x - xorig;
1465  float dy = y - yorig;
1466  float dz = z - zorig;
1467  float r = sqrt(dx*dx + dy*dy + dz*dz);
1468 
1469  float cos_beta = (dx*cos(phiorig) + dy*sin(phiorig))/r;
1470  float sin_beta = dz/r;
1471 
1472  float reltheta = acos(cos_beta*cos(thetaorig) + sin_beta*sin(thetaorig)) - M_PI/2.;
1473  if (reltheta < 0) reltheta += 2. * M_PI;
1474  else if (reltheta > 2.*M_PI) reltheta -= 2. * M_PI;
1475 
1476  return reltheta;
1477 }
1478 
1479 bool TpcDirectLaserReconstruction::sameSign(float num1, float num2)
1480 {
1481  return (num1 >= 0 && num2 >= 0) || (num1 < 0 && num2 < 0);
1482 }
1483