ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcSpaceChargeReconstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcSpaceChargeReconstruction.cc
1 
10 
14 #include <phool/getClass.h>
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHNodeIterator.h>
17 #include <trackbase/TrkrCluster.h>
21 
22 #include <TFile.h>
23 
24 #include <cassert>
25 #include <memory>
26 
27 namespace
28 {
29 
31  template<class T> inline constexpr T square( const T& x ) { return x*x; }
32 
34  template< class T>
35  inline constexpr T delta_phi( const T& phi )
36  {
37  if( phi >= M_PI ) return phi - 2*M_PI;
38  else if( phi < -M_PI ) return phi + 2*M_PI;
39  else return phi;
40  }
41 
43  template<class T> T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
44 
46  template<int type>
47  int get_clusters( SvtxTrack* track )
48  {
49  return std::count_if( track->begin_cluster_keys(), track->end_cluster_keys(),
50  []( const TrkrDefs::cluskey& key ) { return TrkrDefs::getTrkrId(key) == type; } );
51  }
52 
53  // phi range
54  static constexpr float m_phimin = 0;
55  static constexpr float m_phimax = 2.*M_PI;
56 
57  // TODO: could try to get the r and z range from TPC geometry
58  // r range
59  static constexpr float m_rmin = 20;
60  static constexpr float m_rmax = 78;
61 
62  // z range
63  static constexpr float m_zmin = -105.5;
64  static constexpr float m_zmax = 105.5;
65 
66 }
67 
68 //_____________________________________________________________________
70  SubsysReco( name)
71  , PHParameterInterface(name)
72  , m_matrix_container( new TpcSpaceChargeMatrixContainerv1 )
73 {
75 }
76 
77 //_____________________________________________________________________
79 { m_matrix_container->set_grid_dimensions( phibins, rbins, zbins ); }
80 
81 //_____________________________________________________________________
83 {
84  // reset counters
85  m_total_tracks = 0;
87 
88  m_total_clusters = 0;
90 
92 }
93 
94 //_____________________________________________________________________
96 {
97 
98  // load parameters
100  m_max_talpha = get_double_param( "spacecharge_max_talpha" );
101  m_max_drphi = get_double_param( "spacecharge_max_drphi" );
102  m_max_tbeta = get_double_param( "spacecharge_max_tbeta" );
103  m_max_dz = get_double_param( "spacecharge_max_dz" );
104 
105  // print
106  if( Verbosity() )
107  {
108  std::cout
109  << "TpcSpaceChargeReconstruction::InitRun\n"
110  << " m_outputfile: " << m_outputfile << "\n"
111  << " m_use_micromegas: " << std::boolalpha << m_use_micromegas << "\n"
112  << " m_max_talpha: " << m_max_talpha << "\n"
113  << " m_max_drphi: " << m_max_drphi << "\n"
114  << " m_max_tbeta: " << m_max_tbeta << "\n"
115  << " m_max_dz: " << m_max_dz << "\n"
116  << std::endl;
117 
118  // also identify the matrix container
119  m_matrix_container->identify();
120  }
121 
123 }
124 
125 //_____________________________________________________________________
127 {
128  // load nodes
129  const auto res = load_nodes(topNode);
130  if( res != Fun4AllReturnCodes::EVENT_OK ) return res;
131 
132  process_tracks();
134 }
135 
136 //_____________________________________________________________________
138 {
139 
140  // save matrix container in output file
141  if( m_matrix_container )
142  {
143  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
144  outputfile->cd();
145  m_matrix_container->Write( "TpcSpaceChargeMatrixContainer" );
146  }
147 
148  // print counters
149  std::cout
150  << "TpcSpaceChargeReconstruction::End -"
151  << " track statistics total: " << m_total_tracks
152  << " accepted: " << m_accepted_tracks
153  << " fraction: " << 100.*m_accepted_tracks/m_total_tracks << "%"
154  << std::endl;
155 
156  std::cout
157  << "TpcSpaceChargeReconstruction::End -"
158  << " cluster statistics total: " << m_total_clusters
159  << " accepted: " << m_accepted_clusters << " fraction: "
161  << std::endl;
162 
164 }
165 
166 //___________________________________________________________________________
168 {
169  // residual cuts
170  set_default_double_param( "spacecharge_max_talpha", 0.6 );
171  set_default_double_param( "spacecharge_max_drphi", 0.5 );
172  set_default_double_param( "spacecharge_max_tbeta", 1.5 );
173  set_default_double_param( "spacecharge_max_dz", 0.5 );
174 }
175 
176 //_____________________________________________________________________
178 {
179  // get necessary nodes
180  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
181  assert(m_track_map);
182 
183  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
184  assert(m_cluster_map);
186 }
187 
188 //_____________________________________________________________________
190 {
191  if( !( m_track_map && m_cluster_map ) ) return;
192  for( auto iter = m_track_map->begin(); iter != m_track_map->end(); ++iter )
193  {
194  ++m_total_tracks;
195  if( accept_track( iter->second ) )
196  {
198  process_track( iter->second );
199  }
200  }
201 }
202 
203 //_____________________________________________________________________
205 {
206  // ignore tracks whose transverse momentum is too small
207  const auto pt = std::sqrt( square( track->get_px() ) + square( track->get_py() ) );
208  if( pt < 0.5 ) return false;
209 
210  // ignore tracks with too few mvtx, intt and micromegas hits
211  if( get_clusters<TrkrDefs::mvtxId>(track) < 2 ) return false;
212  if( get_clusters<TrkrDefs::inttId>(track) < 2 ) return false;
213  if( m_use_micromegas && get_clusters<TrkrDefs::micromegasId>(track) < 2 ) return false;
214 
215  // all tests passed
216  return true;
217 }
218 
219 //_____________________________________________________________________
221 {
222 
223  // running track state
224  auto state_iter = track->begin_states();
225 
226  // loop over clusters
227  for( auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
228  {
229 
230  const auto& cluster_key = *key_iter;
231  auto cluster = m_cluster_map->findCluster( cluster_key );
232  if( !cluster )
233  {
234  std::cout << PHWHERE << " unable to find cluster for key " << cluster_key << std::endl;
235  continue;
236  }
237 
239 
240  // make sure
241  const auto detId = TrkrDefs::getTrkrId(cluster_key);
242  if(detId != TrkrDefs::tpcId) continue;
243 
244  // cluster r, phi and z
245  const auto cluster_r = get_r( cluster->getX(), cluster->getY() );
246  const auto cluster_phi = std::atan2( cluster->getY(), cluster->getX() );
247  const auto cluster_z = cluster->getZ();
248 
249  // cluster errors
250  const auto cluster_rphi_error = cluster->getRPhiError();
251  const auto cluster_z_error = cluster->getZError();
252 
253  /*
254  remove clusters with too small errors since they are likely pathological
255  and have a large contribution to the chisquare
256  TODO: make these cuts configurable
257  */
258  if( cluster_rphi_error < 0.015 ) continue;
259  if( cluster_z_error < 0.05 ) continue;
260 
261  // find track state that is the closest to cluster
262  /* this assumes that both clusters and states are sorted along r */
263  float dr_min = -1;
264  for( auto iter = state_iter; iter != track->end_states(); ++iter )
265  {
266  const auto dr = std::abs( cluster_r - get_r( iter->second->get_x(), iter->second->get_y() ) );
267  if( dr_min < 0 || dr < dr_min )
268  {
269  state_iter = iter;
270  dr_min = dr;
271  } else break;
272  }
273 
274  // get relevant track state
275  const auto state = state_iter->second;
276 
277  // extrapolate track parameters to the cluster r
278  const auto track_r = get_r( state->get_x(), state->get_y() );
279  const auto dr = cluster_r - track_r;
280  const auto track_drdt = (state->get_x()*state->get_px() + state->get_y()*state->get_py())/track_r;
281  const auto track_dxdr = state->get_px()/track_drdt;
282  const auto track_dydr = state->get_py()/track_drdt;
283  const auto track_dzdr = state->get_pz()/track_drdt;
284 
285  // store state position
286  const auto track_x = state->get_x() + dr*track_dxdr;
287  const auto track_y = state->get_y() + dr*track_dydr;
288  const auto track_z = state->get_z() + dr*track_dzdr;
289  const auto track_phi = std::atan2( track_y, track_x );
290 
291  // get track angles
292  const auto cosphi( std::cos( track_phi ) );
293  const auto sinphi( std::sin( track_phi ) );
294  const auto track_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
295  const auto track_pr = state->get_px()*cosphi + state->get_py()*sinphi;
296  const auto track_pz = state->get_pz();
297  const auto talpha = -track_pphi/track_pr;
298  const auto tbeta = -track_pz/track_pr;
299 
300  // sanity check
301  if( std::isnan(talpha) )
302  {
303  std::cout << "TpcSpaceChargeReconstruction::process_track - talpha is nan" << std::endl;
304  continue;
305  }
306 
307  if( std::isnan(tbeta) )
308  {
309  std::cout << "TpcSpaceChargeReconstruction::process_track - tbeta is nan" << std::endl;
310  continue;
311  }
312 
313  // check against limits
314  if( std::abs( talpha ) > m_max_talpha ) continue;
315  if( std::abs( tbeta ) > m_max_tbeta ) continue;
316 
317  // track errors
318  const auto track_rphi_error = state->get_rphi_error();
319  const auto track_z_error = state->get_z_error();
320 
321  // residuals
322  const auto drp = cluster_r*delta_phi( cluster_phi - track_phi );
323  const auto dz = cluster_z - track_z;
324 
325  // sanity checks
326  if( std::isnan(drp) )
327  {
328  std::cout << "TpcSpaceChargeReconstruction::process_track - drp is nan" << std::endl;
329  continue;
330  }
331 
332  if( std::isnan(dz) )
333  {
334  std::cout << "TpcSpaceChargeReconstruction::process_track - dz is nan" << std::endl;
335  continue;
336  }
337 
338  // check against limits
339  if( std::abs( drp ) > m_max_drphi ) continue;
340  if( std::abs( dz ) > m_max_dz ) continue;
341 
342  // residual errors squared
343  const auto erp = square(track_rphi_error) + square(cluster_rphi_error);
344  const auto ez = square(track_z_error) + square(cluster_z_error);
345 
346  // sanity check
347  // TODO: check whether this happens and fix upstream
348  if( std::isnan( erp ) )
349  {
350  std::cout << "TpcSpaceChargeReconstruction::process_track - erp is nan" << std::endl;
351  continue;
352  }
353 
354  if( std::isnan( ez ) )
355  {
356  std::cout << "TpcSpaceChargeReconstruction::process_track - ez is nan" << std::endl;
357  continue;
358  }
359 
360  // get cell
361  const auto i = get_cell_index( cluster );
362  if( i < 0 )
363  {
364  std::cout << "TpcSpaceChargeReconstruction::process_track - invalid cell index" << std::endl;
365  continue;
366  }
367 
368  // update matrices
369  // see https://indico.bnl.gov/event/7440/contributions/43328/attachments/31334/49446/talk.pdf for details
370  m_matrix_container->add_to_lhs(i, 0, 0, 1./erp );
371  m_matrix_container->add_to_lhs(i, 0, 1, 0 );
372  m_matrix_container->add_to_lhs(i, 0, 2, talpha/erp );
373 
374  m_matrix_container->add_to_lhs(i, 1, 0, 0 );
375  m_matrix_container->add_to_lhs(i, 1, 1, 1./ez );
376  m_matrix_container->add_to_lhs(i, 1, 2, tbeta/ez );
377 
378  m_matrix_container->add_to_lhs(i, 2, 0, talpha/erp );
379  m_matrix_container->add_to_lhs(i, 2, 1, tbeta/ez );
380  m_matrix_container->add_to_lhs(i, 2, 2, square(talpha)/erp + square(tbeta)/ez );
381 
382  m_matrix_container->add_to_rhs(i, 0, drp/erp );
383  m_matrix_container->add_to_rhs(i, 1, dz/ez );
384  m_matrix_container->add_to_rhs(i, 2, talpha*drp/erp + tbeta*dz/ez );
385 
386  // update entries in cell
387  m_matrix_container->add_to_entries(i);
388 
389  // increment number of accepted clusters
391 
392  }
393 
394 }
395 
396 //_____________________________________________________________________
398 {
399  // get grid dimensions from matrix container
400  int phibins = 0;
401  int rbins = 0;
402  int zbins = 0;
403  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
404 
405  // phi
406  // bound check
407  float phi = std::atan2( cluster->getY(), cluster->getX() );
408  while( phi < m_phimin ) phi += 2.*M_PI;
409  while( phi >= m_phimax ) phi -= 2.*M_PI;
410  int iphi = phibins*(phi-m_phimin)/(m_phimax-m_phimin);
411 
412  // radius
413  const float r = get_r( cluster->getX(), cluster->getY() );
414  if( r < m_rmin || r >= m_rmax ) return -1;
415  int ir = rbins*(r-m_rmin)/(m_rmax-m_rmin);
416 
417  // z
418  const float z = cluster->getZ();
419  if( z < m_zmin || z >= m_zmax ) return -1;
420  int iz = zbins*(z-m_zmin)/(m_zmax-m_zmin);
421 
422  return m_matrix_container->get_cell_index( iphi, ir, iz );
423 }