ECCE @ EIC 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 
8 
10 
12 #include <phool/getClass.h>
15 #include <trackbase/TrkrCluster.h>
21 
22 #include <TFile.h>
23 #include <TH1.h>
24 #include <TH2.h>
25 #include <TH3.h>
26 #include <TVector3.h>
27 
28 #include <cassert>
29 
30 namespace
31 {
32 
34  template<class T> class range_adaptor
35  {
36  public:
37  range_adaptor( const T& range ):m_range(range){}
38  inline const typename T::first_type& begin() {return m_range.first;}
39  inline const typename T::second_type& end() {return m_range.second;}
40  private:
41  T m_range;
42  };
43 
45  template<class T>
46  inline constexpr T square( const T& x ) { return x*x; }
47 
49  template<class T>
50  inline constexpr T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
51 
52  // calculate intersection between line and circle
53  double line_circle_intersection( const TVector3& p, const TVector3& d, double radius )
54  {
55  const double A = square(d.x()) + square(d.y());
56  const double B = 2*p.x()*d.x() + 2*p.y()*d.y();
57  const double C = square(p.x()) + square(p.y()) - square(radius);
58  const double delta = square(B)-4*A*C;
59  if( delta < 0 ) return -1;
60 
61  // check first intersection
62  const double tup = (-B + std::sqrt(delta))/(2*A);
63  if( tup >= 0 ) return tup;
64 
65  // check second intersection
66  const double tdn = (-B-sqrt(delta))/(2*A);
67  if( tdn >= 0 ) return tdn;
68 
69  // no valid extrapolation
70  return -1;
71  }
72 
74  inline std::ostream& operator << (std::ostream& out, const TVector3& vector )
75  {
76  out << "( " << vector.x() << ", " << vector.y() << ", " << vector.z() << ")";
77  return out;
78  }
79 
81  template< class T>
82  inline constexpr T delta_phi( const T& phi )
83  {
84  if( phi >= M_PI ) return phi - 2*M_PI;
85  else if( phi < -M_PI ) return phi + 2*M_PI;
86  else return phi;
87  }
88 
89  // phi range
90  static constexpr float m_phimin = 0;
91  static constexpr float m_phimax = 2.*M_PI;
92 
93  // TODO: could try to get the r and z range from TPC geometry
94  // r range
95  static constexpr float m_rmin = 20;
96  static constexpr float m_rmax = 78;
97 
98  // z range
99  static constexpr float m_zmin = -105.5;
100  static constexpr float m_zmax = 105.5;
101 
102 }
103 
104 //_____________________________________________________________________
106  SubsysReco( name)
107  , PHParameterInterface(name)
108  , m_matrix_container( new TpcSpaceChargeMatrixContainerv1 )
109 {
111 }
112 
113 //_____________________________________________________________________
115 {
116  m_total_clusters = 0;
118 
121 }
122 
123 //_____________________________________________________________________
125 {
127  m_max_dca = get_double_param( "directlaser_max_dca" );
128  m_max_drphi = get_double_param( "directlaser_max_drphi" );
129  m_max_dz = get_double_param( "directlaser_max_dz" );
130 
131  // print
132  if( Verbosity() )
133  {
134  std::cout
135  << "TpcDirectLaserReconstruction::InitRun\n"
136  << " m_outputfile: " << m_outputfile << "\n"
137  << " m_max_dca: " << m_max_dca << "\n"
138  << " m_max_drphi: " << m_max_drphi << "\n"
139  << " m_max_dz: " << m_max_dz << "\n"
140  << std::endl;
141 
142  // also identify the matrix container
143  m_matrix_container->identify();
144  }
145 
147 }
148 
149 //_____________________________________________________________________
151 {
152  // load nodes
153  const auto res = load_nodes(topNode);
154  if( res != Fun4AllReturnCodes::EVENT_OK ) return res;
155 
156  process_tracks();
158 }
159 
160 //_____________________________________________________________________
162 {
163  // save matrix container in output file
164  if( m_matrix_container )
165  {
166  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
167  outputfile->cd();
168  m_matrix_container->Write( "TpcSpaceChargeMatrixContainer" );
169  }
170 
171  // write evaluation histograms to output
173  {
174  m_histogramfile->cd();
175  for(const auto& o:std::initializer_list<TObject*>({ h_dca_layer, h_deltarphi_layer, h_deltaz_layer, h_entries }))
176  { if( o ) o->Write(); }
177  m_histogramfile->Close();
178  }
179 
180  // print counters
181  std::cout
182  << "TpcDirectLaserReconstruction::End -"
183  << " cluster statistics total: " << m_total_clusters
184  << " accepted: " << m_accepted_clusters << " fraction: "
186  << std::endl;
187 
189 }
190 
191 //___________________________________________________________________________
193 {
194 
195  // DCA cut, to decide whether a cluster should be associated to a given laser track or not
196  set_default_double_param( "directlaser_max_dca", 1.5 );
197 
198 
199 // // residual cuts, used to decide if a given cluster is used to fill SC reconstruction matrices
200 // set_default_double_param( "directlaser_max_drphi", 0.5 );
201 // set_default_double_param( "directlaser_max_dz", 0.5 );
202 
203  set_default_double_param( "directlaser_max_drphi", 2. );
204  set_default_double_param( "directlaser_max_dz", 2. );
205 }
206 
207 //_____________________________________________________________________
209 { m_matrix_container->set_grid_dimensions( phibins, rbins, zbins ); }
210 
211 //_____________________________________________________________________
213 {
214  // acts surface map
215  m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
216  assert( m_surfmaps );
217 
218  // acts geometry
219  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
220  assert( m_tGeometry );
221 
222  // tracks
223  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
224  assert(m_track_map);
225 
226  // hitset container
227  m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
228 
229  // clusters
230  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
231  assert(m_cluster_map);
233 }
234 
235 //_____________________________________________________________________
237 {
238  std::cout << "TpcDirectLaserReconstruction::makeHistograms - writing evaluation histograms to: " << m_histogramfilename << std::endl;
239  m_histogramfile.reset( new TFile(m_histogramfilename.c_str(), "RECREATE") );
240  m_histogramfile->cd();
241 
242  // residuals vs layers
243  h_dca_layer = new TH2F( "dca_layer", ";layer; DCA (cm)", 57, 0, 57, 500, 0, 2 );
244  h_deltarphi_layer = new TH2F( "deltarphi_layer", ";layer; r.#Delta#phi_{track-cluster} (cm)", 57, 0, 57, 2000, -2, 2 );
245  h_deltaz_layer = new TH2F( "deltaz_layer", ";layer; #Deltaz_{track-cluster} (cm)", 57, 0, 57, 400, -2, 2 );
246 
247  // entries vs cell grid
248  /* histogram dimension and axis limits must match that of TpcSpaceChargeMatrixContainer */
249  if( m_matrix_container )
250  {
251  int phibins = 0;
252  int rbins = 0;
253  int zbins = 0;
254  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
255  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 );
256  }
257 }
258 
259 //_____________________________________________________________________
261 {
262  if( !( m_track_map && m_cluster_map ) ) return;
263 
264  // count number of clusters in the TPC
265  for( const auto& [hitsetkey,hitset]:range_adaptor(m_hitsetcontainer->getHitSets()))
266  {
267  if( TrkrDefs::getTrkrId( hitsetkey ) != TrkrDefs::tpcId ) continue;
268 
269  const auto range = m_cluster_map->getClusters(hitsetkey);
270  m_total_clusters += std::distance( range.first, range.second );
271  }
272 
273  // loop over tracks and process
274  for( auto iter = m_track_map->begin(); iter != m_track_map->end(); ++iter )
275  { process_track( iter->second ); }
276 }
277 
278 //_____________________________________________________________________
280 {
281 
282  // get track parameters
283  const TVector3 origin( track->get_x(), track->get_y(), track->get_z() );
284  const TVector3 direction( track->get_px(), track->get_py(), track->get_pz() );
285 
286  if( Verbosity() )
287  { std::cout << "TpcDirectLaserReconstruction::process_track - position: " << origin << " direction: " << direction << std::endl; }
288 
289  // loop over hitsets
290  for( const auto& [hitsetkey,hitset]:range_adaptor(m_hitsetcontainer->getHitSets()))
291  {
292 
293  // only check TPC hitsets
294  if( TrkrDefs::getTrkrId( hitsetkey ) != TrkrDefs::tpcId ) continue;
295 
296  // store layer
297  const auto layer = TrkrDefs::getLayer( hitsetkey );
298 
299  // get corresponding clusters
300  for( const auto& [key,cluster]:range_adaptor(m_cluster_map->getClusters(hitsetkey)))
301  {
302 
303  // get cluster global coordinates
304  const auto global = m_transformer.getGlobalPosition(cluster,m_surfmaps, m_tGeometry);
305 
306  // calculate dca
307  const TVector3 oc( global.x()-origin.x(), global.y()-origin.y(), global.z()-origin.z() );
308  auto t = direction.Dot( oc )/square( direction.Mag() );
309  auto om = direction*t;
310  const auto dca = (oc-om).Mag();
311 
312  // do not associate if dca is too large
313  if( dca > m_max_dca ) continue;
314 
315  // calculate intersection to layer
316  t = line_circle_intersection(origin, direction, get_r( global.x(), global.y() ));
317  if( t < 0 ) continue;
318 
319  // update position on track
320  om = direction*t;
321 
322  // path length
323  const auto pathlength = om.Mag();
324 
325  // get projection to the track
326  const auto projection = origin + om;
327 
328  // create relevant state vector and assign to track
329  SvtxTrackState_v1 state( pathlength );
330  state.set_x( projection.x() );
331  state.set_y( projection.y() );
332  state.set_z( projection.z() );
333 
334  state.set_px( direction.x());
335  state.set_py( direction.y());
336  state.set_pz( direction.z());
337  track->insert_state( &state );
338 
339  // also associate cluster to track
340  track->insert_cluster_key( key );
341 
342  // cluster r, phi and z
343  const auto cluster_r = get_r(global.x(), global.y());
344  const auto cluster_phi = std::atan2(global.y(),global.x());
345  const auto cluster_z = global.z();
346 
347  // cluster errors
348  const auto cluster_rphi_error = cluster->getRPhiError();
349  const auto cluster_z_error = cluster->getZError();
350 
351 // /*
352 // remove clusters with too small errors since they are likely pathological
353 // and have a large contribution to the chisquare
354 // TODO: make these cuts configurable
355 // */
356 // if( cluster_rphi_error < 0.015 ) continue;
357 // if( cluster_z_error < 0.05 ) continue;
358 
359  // track position
360  const auto track_phi = std::atan2( projection.y(), projection.x() );
361  const auto track_z = projection.z();
362 
363  // track angles
364  const auto cosphi( std::cos( track_phi ) );
365  const auto sinphi( std::sin( track_phi ) );
366  const auto track_pphi = -state.get_px()*sinphi + state.get_py()*cosphi;
367  const auto track_pr = state.get_px()*cosphi + state.get_py()*sinphi;
368  const auto track_pz = state.get_pz();
369  const auto talpha = -track_pphi/track_pr;
370  const auto tbeta = -track_pz/track_pr;
371 
372  // sanity check
373  if( std::isnan(talpha) )
374  {
375  std::cout << "TpcDirectLaserReconstruction::process_track - talpha is nan" << std::endl;
376  continue;
377  }
378 
379  if( std::isnan(tbeta) )
380  {
381  std::cout << "TpcDirectLaserReconstruction::process_track - tbeta is nan" << std::endl;
382  continue;
383  }
384 
385  // residuals
386  const auto drp = cluster_r*delta_phi( cluster_phi - track_phi );
387  const auto dz = cluster_z - track_z;
388 
389  // sanity checks
390  if( std::isnan(drp) )
391  {
392  std::cout << "TpcDirectLaserReconstruction::process_track - drp is nan" << std::endl;
393  continue;
394  }
395 
396  if( std::isnan(dz) )
397  {
398  std::cout << "TpcDirectLaserReconstruction::process_track - dz is nan" << std::endl;
399  continue;
400  }
401 
402  if(m_savehistograms)
403  {
404  if(h_dca_layer) h_dca_layer->Fill(layer, dca);
405  if(h_deltarphi_layer) h_deltarphi_layer->Fill(layer, drp);
407  if(h_entries)
408  {
409  auto phi = cluster_phi;
410  while( phi < m_phimin ) phi += 2.*M_PI;
411  while( phi >= m_phimax ) phi -= 2.*M_PI;
412  h_entries->Fill( phi, cluster_r, cluster_z );
413  }
414  }
415 
416 // // check against limits
417 // if( std::abs( drp ) > m_max_drphi ) continue;
418 // if( std::abs( dz ) > m_max_dz ) continue;
419 
420  // residual errors squared
421  const auto erp = square(cluster_rphi_error);
422  const auto ez = square(cluster_z_error);
423 
424  // sanity check
425  if( std::isnan( erp ) )
426  {
427  std::cout << "TpcDirectLaserReconstruction::process_track - erp is nan" << std::endl;
428  continue;
429  }
430 
431  if( std::isnan( ez ) )
432  {
433  std::cout << "TpcDirectLaserReconstruction::process_track - ez is nan" << std::endl;
434  continue;
435  }
436 
437  // get cell
438  const auto i = get_cell_index( global );
439  if( i < 0 )
440  {
441  if( Verbosity() )
442  {
443  std::cout << "TpcDirectLaserReconstruction::process_track - invalid cell index"
444  << " r: " << cluster_r
445  << " phi: " << cluster_phi
446  << " z: " << cluster_z
447  << std::endl;
448  }
449  continue;
450  }
451 
452  // update matrices
453  // see https://indico.bnl.gov/event/7440/contributions/43328/attachments/31334/49446/talk.pdf for details
454  m_matrix_container->add_to_lhs(i, 0, 0, 1./erp );
455  m_matrix_container->add_to_lhs(i, 0, 1, 0 );
456  m_matrix_container->add_to_lhs(i, 0, 2, talpha/erp );
457 
458  m_matrix_container->add_to_lhs(i, 1, 0, 0 );
459  m_matrix_container->add_to_lhs(i, 1, 1, 1./ez );
460  m_matrix_container->add_to_lhs(i, 1, 2, tbeta/ez );
461 
462  m_matrix_container->add_to_lhs(i, 2, 0, talpha/erp );
463  m_matrix_container->add_to_lhs(i, 2, 1, tbeta/ez );
464  m_matrix_container->add_to_lhs(i, 2, 2, square(talpha)/erp + square(tbeta)/ez );
465 
466  m_matrix_container->add_to_rhs(i, 0, drp/erp );
467  m_matrix_container->add_to_rhs(i, 1, dz/ez );
468  m_matrix_container->add_to_rhs(i, 2, talpha*drp/erp + tbeta*dz/ez );
469 
470  // update entries in cell
471  m_matrix_container->add_to_entries(i);
472 
473  // increment number of accepted clusters
475 
476  }
477 
478  }
479 
480 }
481 
482 //_____________________________________________________________________
484 {
485  // get grid dimensions from matrix container
486  int phibins = 0;
487  int rbins = 0;
488  int zbins = 0;
489  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
490 
491  // phi
492  // bound check
493  float phi = std::atan2( global.y(), global.x() );
494  while( phi < m_phimin ) phi += 2.*M_PI;
495  while( phi >= m_phimax ) phi -= 2.*M_PI;
496  int iphi = phibins*(phi-m_phimin)/(m_phimax-m_phimin);
497 
498  // radius
499  const float r = get_r( global.x(), global.y() );
500  if( r < m_rmin || r >= m_rmax ) return -1;
501  int ir = rbins*(r-m_rmin)/(m_rmax-m_rmin);
502 
503  // z
504  const float z = global.z();
505  if( z < m_zmin || z >= m_zmax ) return -1;
506  int iz = zbins*(z-m_zmin)/(m_zmax-m_zmin);
507 
508  return m_matrix_container->get_cell_index( iphi, ir, iz );
509 }