ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcResiduals.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcResiduals.cc
1 #include "PHTpcResiduals.h"
3 
6 #include <phool/getClass.h>
7 #include <phool/phool.h>
8 #include <phool/PHDataNode.h>
9 #include <phool/PHNode.h>
10 #include <phool/PHNodeIterator.h>
11 #include <phool/PHObject.h>
12 #include <phool/PHTimer.h>
13 
14 #include <trackbase/TrkrCluster.h>
20 
22 
24 #include <Acts/Geometry/GeometryIdentifier.hpp>
30 
31 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
32 #include <ActsExamples/EventData/TrkrClusterSourceLink.hpp>
33 
34 
35 #include <cmath>
36 #include <TFile.h>
37 #include <TTree.h>
38 #include <TH1.h>
39 #include <TH2.h>
40 
41 #include <iostream>
42 #include <sstream>
43 
44 namespace
45 {
46 
47  // square
48  template<class T> inline constexpr T square( const T& x ) { return x*x; }
49 
50  // radius
51  template<class T> T get_r( const T& x, const T& y ) { return std::sqrt( square(x) + square(y) ); }
52 
53  template<class T> inline constexpr T deltaPhi(const T& phi)
54  {
55  if (phi > M_PI)
56  return phi - 2. * M_PI;
57  else if (phi <= -M_PI)
58  return phi + 2.* M_PI;
59  else
60  return phi;
61  }
62 }
63 
65  : SubsysReco(name)
66  , m_matrix_container( new TpcSpaceChargeMatrixContainerv1 )
67 {}
68 
70 {
73 }
74 
76 {
79 
82 
84 }
85 
87 {
88  if(Verbosity() > 1)
89  std::cout <<"Starting PHTpcResiduals event "
90  << m_event << std::endl;
91 
92  if(m_event % 1000 == 0)
93  std::cout << "PHTpcResiduals processed " << m_event
94  << " events" << std::endl;
95 
96  int returnVal = processTracks(topNode);
97 
98  if(Verbosity() > 1)
99  std::cout <<"Finished PHTpcResiduals event "
100  << m_event << std::endl;
101 
102  m_event++;
103 
104  return returnVal;
105 }
106 
108 {
109  std::cout << "PHTpcResiduals::End - writing matrices to " << m_outputfile << std::endl;
110  if(Verbosity() > 0)
111  { std::cout << "PHTpcResiduals::End - Number of bad SL propagations " << m_nBadProps << std::endl; }
112 
113  // save matrix container in output file
114  if( m_matrix_container )
115  {
116  std::unique_ptr<TFile> outputfile( TFile::Open( m_outputfile.c_str(), "RECREATE" ) );
117  outputfile->cd();
118  m_matrix_container->Write( "TpcSpaceChargeMatrixContainer" );
119  }
120 
121  // save histograms
123  {
124  m_histogramfile->cd();
125  for( const auto o:std::initializer_list<TObject*>({ h_rphiResid, h_zResid, h_etaResidLayer, h_zResidLayer, h_etaResid, h_index, h_alpha, h_beta, h_deltarphi_layer, h_deltaz_layer, residTup }) )
126  { if( o ) o->Write(); }
127  m_histogramfile->Close();
128  }
130 }
131 
132 
134 {
135 
136  if( Verbosity() )
137  { std::cout << "PHTpcResiduals::processTracks - proto track size " << m_trackMap->size() <<std::endl; }
138 
139  for(auto &[trackKey, track] : *m_trackMap)
140  {
141  if(Verbosity() > 1)
142  std::cout << "Processing track key " << trackKey
143  << std::endl;
144  if(checkTrack(track))
146  }
147 
149 }
150 
152 {
153 
154  if(track->get_pt() < 0.5)
155  return false;
156 
157  if(Verbosity() > 2)
158  std::cout << "Track has pt " << track->get_pt() << std::endl;
159 
160  int nMvtxHits = 0;
161  int nInttHits = 0;
162  int nMMHits = 0;
163 
164  for (SvtxTrack::ConstClusterKeyIter clusIter = track->begin_cluster_keys();
165  clusIter != track->end_cluster_keys();
166  ++clusIter)
167  {
168  auto key = *clusIter;
169  auto trkrId = TrkrDefs::getTrkrId(key);
170  if(trkrId == TrkrDefs::TrkrId::mvtxId)
171  nMvtxHits++;
172  else if(trkrId == TrkrDefs::TrkrId::inttId)
173  nInttHits++;
174  else if(trkrId == TrkrDefs::TrkrId::micromegasId)
175  nMMHits++;
176  }
177 
178  if(Verbosity() > 2)
179  std::cout << "Number of mvtx/intt/MM hits "
180  << nMvtxHits << "/" << nInttHits << "/"
181  << nMMHits << std::endl;
182 
183  // Require at least 2 hits in each detector
184  if( m_useMicromegas )
185  {
186  if(nMvtxHits<2 || nInttHits<2 || nMMHits<2)
187  return false;
188  } else {
189  if(nMvtxHits<2 || nInttHits<2 )
190  return false;
191  }
192 
193  return true;
194 
195 }
196 
198 {
199  auto key = cluster->getClusKey();
200  auto subsurfkey = cluster->getSubSurfKey();
201 
202  // Make a safety check for clusters that couldn't be attached
203  // to a surface
204  auto surf = getSurface(key, subsurfkey);
205  if(!surf)
206  return SourceLink();
207 
208  Acts::BoundVector loc = Acts::BoundVector::Zero();
211 
212  Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
221 
222  SourceLink sl(key, surf, loc, cov);
223 
224  return sl;
225 }
226 
227 //___________________________________________________________________________________
229 {
230  const auto trkrid = TrkrDefs::getTrkrId(cluskey);
231  const auto hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
232 
233  switch( trkrid )
234  {
236  case TrkrDefs::TrkrId::tpcId: return getTpcSurface(hitsetkey, surfkey);
239  {
241  }
242  }
243 
244  // unreachable
245  return nullptr;
246 
247 }
248 
249 //___________________________________________________________________________________
251 {
252  auto surfMap = m_surfMaps->siliconSurfaceMap;
253  auto iter = surfMap.find(hitsetkey);
254  if(iter != surfMap.end())
255  {
256  return iter->second;
257  }
258 
259  // If it can't be found, return nullptr
260  return nullptr;
261 
262 }
263 
264 //___________________________________________________________________________________
266 {
267  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
268  const auto iter = m_surfMaps->tpcSurfaceMap.find(layer);
269  if(iter != m_surfMaps->tpcSurfaceMap.end())
270  {
271  auto surfvec = iter->second;
272  return surfvec.at(surfkey);
273  }
274 
275  // If it can't be found, return nullptr to skip this cluster
276  return nullptr;
277 }
278 
279 //___________________________________________________________________________________
281 {
282  const auto iter = m_surfMaps->mmSurfaceMap.find( hitsetkey );
283  return (iter == m_surfMaps->mmSurfaceMap.end()) ? nullptr:iter->second;
284 }
285 
286 //___________________________________________________________________________________
288 {
289  Acts::Vector3D momentum(track->get_px(),
290  track->get_py(),
291  track->get_pz());
292  double trackQ = track->get_charge() * Acts::UnitConstants::e;
293  double p = track->get_p();
294 
296  for(int i =0; i<6; i++)
297  {
298  for(int j =0; j<6; j++)
299  {
300  cov(i,j) = track->get_acts_covariance(i, j);
301  }
302  }
304  track->get_y() * Acts::UnitConstants::cm,
305  track->get_z() * Acts::UnitConstants::cm);
306 
307  const auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(position);
308  const auto actsFourPos = Acts::Vector4D(position(0), position(1), position(2), 10 * Acts::UnitConstants::ns);
310  actsFourPos, momentum,
311  p, trackQ, cov);
312 
313  return seed;
314 
315 }
316 
318 {
319 
320  if(Verbosity() > 1)
321  std::cout << "Propagating silicon+MM fit params momentum: "
322  << track->get_p() << " and position "
323  << track->get_x() << ", " << track->get_y()
324  << ", " << track->get_z() << " cm "
325  << std::endl;
326 
327  auto trackParams = makeTrackParams(track);
328 
329  int initNBadProps = m_nBadProps;
330  for (SvtxTrack::ConstClusterKeyIter clusIter = track->begin_cluster_keys();
331  clusIter != track->end_cluster_keys();
332  ++clusIter)
333  {
334  auto cluskey = *clusIter;
335 
336  // only propagate to tpc surfaces
338  continue;;
339 
340  auto cluster = m_clusterContainer->findCluster(cluskey);
341 
342  auto sl = makeSourceLink(cluster);
343 
344  auto result = propagateTrackState(trackParams, sl);
345 
346  if(result.ok())
347  {
348  auto pathLength = (*result).first;
349  auto trackStateParams = std::move(*(*result).second);
350  if(Verbosity() > 1)
351  {
352  std::cout << "PHTpcResiduals::processTrack -"
353  << " path length: " << pathLength
354  << " track momentum : "
355  << trackParams.momentum()
356  << " propagator momentum : "
357  << trackStateParams.momentum()
358  << std::endl;
359  }
360 
361  addTrackState( track, pathLength, trackStateParams );
362  calculateTpcResiduals(trackStateParams, cluster);
363 
364  } else {
365 
366  m_nBadProps++;
367  continue;
368 
369  }
370  }
371 
372  if(m_nBadProps > initNBadProps && Verbosity() > 1)
373  std::cout << "Starting track params position/momentum: "
374  << trackParams.position(m_tGeometry->geoContext).transpose()
375  << std::endl << trackParams.momentum().transpose()
376  << std::endl
377  << "Track params phi/eta "
378  << std::atan2(trackParams.momentum().y(),
379  trackParams.momentum().x())
380  << " and "
381  << std::atanh(trackParams.momentum().z() /
382  trackParams.momentum().norm())
383  << std::endl;
384 
385 }
386 
388  const Acts::BoundTrackParameters& params,
389  const SourceLink& sl)
390 {
391  /*
392  std::cout << "Propagating to geo id " << sl.referenceSurface().geometryId() << std::endl;
393  if(sl.referenceSurface().associatedDetectorElement() != nullptr)
394  std::cout << " which has associated detector element " << sl.referenceSurface().associatedDetectorElement()->thickness() << std::endl;
395  */
396  return std::visit([params, sl, this]
397  (auto && inputField) -> ExtrapolationResult {
398  using InputMagneticField =
399  typename std::decay_t<decltype(inputField)>::element_type;
400  using MagneticField = Acts::SharedBField<InputMagneticField>;
403 
404  MagneticField field(inputField);
405  Stepper stepper(field);
406  Propagator propagator(stepper);
407 
409  if(Verbosity() > 10)
410  logLevel = Acts::Logging::VERBOSE;
411 
412  auto logger = Acts::getDefaultLogger("PHTpcResiduals", logLevel);
413 
416  Acts::LoggerWrapper{*logger});
417 
418  auto result = propagator.propagate(params, sl.referenceSurface(),
419  options);
420 
421 
422  if(result.ok())
423  {
424  // return both path length and extrapolated parameters
425  return std::make_pair( (*result).pathLength/Acts::UnitConstants::cm, std::move((*result).endParameters) );
426  } else {
427  return result.error();
428  }
429  },
430  std::move(m_tGeometry->magField));
431 
432 }
433 
435 {
436 
437  /* this is essentially a copy of the code from trackbase_historic/ActsTransformations::fillSvtxTrackStates */
438 
439  // create track state
440  SvtxTrackState_v1 state( pathlength );
441 
442  // save global position
443  const auto global = params.position(m_tGeometry->geoContext);
444  state.set_x(global.x() / Acts::UnitConstants::cm);
445  state.set_y(global.y() / Acts::UnitConstants::cm);
446  state.set_z(global.z() / Acts::UnitConstants::cm);
447 
448  // save momentum
449  const auto momentum = params.momentum();
450  state.set_px(momentum.x());
451  state.set_py(momentum.y());
452  state.set_pz(momentum.z());
453 
454  // covariance
455  ActsTransformations transformer;
456  const auto globalCov = transformer.rotateActsCovToSvtxTrack(params, m_tGeometry->geoContext);
457  for (int i = 0; i < 6; ++i)
458  for (int j = 0; j < 6; ++j)
459  { state.set_error(i, j, globalCov(i,j)); }
460 
461  track->insert_state(&state);
462 }
463 
465  const Acts::BoundTrackParameters &params,
466  TrkrCluster* cluster)
467 {
468 
469  cluskey = cluster->getClusKey();
470  // Get all the relevant information for residual calculation
471  ActsTransformations transformer;
472  const auto globClusPos = transformer.getGlobalPosition(cluster, m_surfMaps, m_tGeometry);
473  clusR = std::sqrt(square(globClusPos(0)) + square(globClusPos(1)));
474  clusPhi = std::atan2(globClusPos(1), globClusPos(0));
475  clusZ = globClusPos(2);
476 
477  clusRPhiErr = cluster->getRPhiError();
478  clusZErr = cluster->getZError();
479 
480  if(Verbosity() > 3)
481  {
482  std::cout << "cluster key is " << cluskey <<std::endl;
483  std::cout << "Cluster r phi and z " << clusR << " "
484  << clusPhi << "+/-" << clusRPhiErr
485  <<" and " << clusZ << "+/-" << clusZErr << std::endl;
486  }
487 
488  if(clusRPhiErr < 0.015)
489  return;
490  if(clusZErr < 0.05)
491  return;
492 
493  const auto globalStatePos = params.position(m_tGeometry->geoContext);
494  const auto globalStateMom = params.momentum();
495  const auto globalStateCov = *params.covariance();
496 
497  stateRPhiErr = sqrt(globalStateCov(Acts::eBoundLoc0,
500  stateZErr = sqrt(globalStateCov(Acts::eBoundLoc1,
503 
504  stateZ = globalStatePos.z() / Acts::UnitConstants::cm;
505 
506  const auto globStateX = globalStatePos.x() / Acts::UnitConstants::cm;
507  const auto globStateY = globalStatePos.y() / Acts::UnitConstants::cm;
508  const auto globStateZ = stateZ;
509 
510  stateR = std::sqrt(square(globStateX) + square(globStateY));
511 
512  const auto dr = clusR - stateR;
513  const auto trackDrDt = (globStateX * globalStateMom(0) +
514  globStateY * globalStateMom(1)) / stateR;
515  const auto trackDxDr = globalStateMom(0) / trackDrDt;
516  const auto trackDyDr = globalStateMom(1) / trackDrDt;
517  const auto trackDzDr = globalStateMom(2) / trackDrDt;
518 
519  const auto trackX = globStateX + dr * trackDxDr;
520  const auto trackY = globStateY + dr * trackDyDr;
521  const auto trackZ = globStateZ + dr * trackDzDr;
522 
523  if(Verbosity() > 2)
524  std::cout << "State Calculations: " << stateR << ", "
525  << dr << ", " << trackDrDt << ", " << trackDxDr
526  << ", " << trackDyDr << ", " << trackDzDr
527  <<" , " << trackX << ", " << trackY << ", "
528  << trackZ << std::endl;
529 
530  statePhi = std::atan2(trackY, trackX);
531  stateZ = trackZ;
532 
533  if(Verbosity() > 3)
534  std::cout << "State r phi and z "
535  << stateR
536  << " " << statePhi << "+/-" << stateRPhiErr
537  << " and " << stateZ << "+/-" << stateZErr << std::endl;
538 
539  const auto erp = square(clusRPhiErr);
540  const auto ez = square(clusZErr);
541 
542  const auto dPhi = clusPhi - statePhi;
543 
544  // Calculate residuals
545  drphi = clusR * deltaPhi(dPhi);
546  dz = clusZ - stateZ;
547 
548  if(Verbosity() > 3)
549  std::cout << "TPC residuals " << drphi << " " << dz << std::endl;
550 
551  const auto trackEta
552  = std::atanh(params.momentum().z() / params.absoluteMomentum());
553  const auto clusEta = std::atanh(clusZ / std::sqrt(
554  square(globClusPos(0)) +
555  square(globClusPos(1)) +
556  square(globClusPos(2))));
557 
558  const auto trackPPhi = -params.momentum()(0) * std::sin(statePhi) +
559  params.momentum()(1) * std::cos(statePhi);
560  const auto trackPR = params.momentum()(0) * std::cos(statePhi) +
561  params.momentum()(1) * std::sin(statePhi);
562  const auto trackPZ = params.momentum()(2);
563  const auto trackAlpha = -trackPPhi / trackPR;
564  const auto trackBeta = -trackPZ / trackPR;
565 
566  if(Verbosity() > 3)
567  std::cout << "Track angles " << trackPPhi << ", " << trackPR
568  << ", " << trackPZ << ", " << trackAlpha << ", " << trackBeta
569  << std::endl;
570 
571  tanBeta = trackBeta;
572  tanAlpha = trackAlpha;
573 
574  // get cell index
575  const auto index = getCell(globClusPos);
576  if(Verbosity() > 3)
577  { std::cout << "Bin index found is " << index << std::endl; }
578 
579  if(index < 0 ) return;
580 
581  if( m_savehistograms )
582  {
583  h_index->Fill(index);
584  h_alpha->Fill(tanAlpha, drphi);
585  h_beta->Fill(tanBeta, dz);
586  h_rphiResid->Fill(clusR , drphi);
587  h_zResid->Fill(stateZ , dz);
588  h_etaResid->Fill(trackEta, clusEta - trackEta);
589  h_zResidLayer->Fill(clusR , dz);
590  h_etaResidLayer->Fill(clusR , clusEta - trackEta);
591 
592  const auto layer = TrkrDefs::getLayer(cluster->getClusKey());
593  h_deltarphi_layer->Fill( layer, drphi );
594  h_deltaz_layer->Fill( layer, dz );
595 
596  residTup->Fill();
597  }
598 
599  // check track angles and residuals agains cuts
600  if(std::abs(trackAlpha) > m_maxTAlpha
602  return;
603 
604  if(std::abs(trackBeta) > m_maxTBeta
606  return;
607 
608  // Fill distortion matrices
609  m_matrix_container->add_to_lhs(index, 0, 0, 1./erp );
610  m_matrix_container->add_to_lhs(index, 0, 1, 0 );
611  m_matrix_container->add_to_lhs(index, 0, 2, trackAlpha/erp );
612 
613  m_matrix_container->add_to_lhs(index, 1, 0, 0 );
614  m_matrix_container->add_to_lhs(index, 1, 1, 1./ez );
615  m_matrix_container->add_to_lhs(index, 1, 2, trackBeta/ez );
616 
617  m_matrix_container->add_to_lhs(index, 2, 0, trackAlpha/erp );
618  m_matrix_container->add_to_lhs(index, 2, 1, trackBeta/ez );
619  m_matrix_container->add_to_lhs(index, 2, 2, square(trackAlpha)/erp + square(trackBeta)/ez );
620 
621  m_matrix_container->add_to_rhs(index, 0, drphi/erp );
622  m_matrix_container->add_to_rhs(index, 1, dz/ez );
623  m_matrix_container->add_to_rhs(index, 2, trackAlpha*drphi/erp + trackBeta*dz/ez );
624 
625  // update entries in cell
626  m_matrix_container->add_to_entries(index);
627 
628  return;
629 }
630 
632 {
633 
634  // get grid dimensions from matrix container
635  int phibins = 0;
636  int rbins = 0;
637  int zbins = 0;
638  m_matrix_container->get_grid_dimensions( phibins, rbins, zbins );
639 
640  // phi
641  float phi = std::atan2(loc(1), loc(0));
642  while( phi < m_phiMin ) phi += 2.*M_PI;
643  while( phi >= m_phiMax ) phi -= 2.*M_PI;
644  const int iphi = phibins * (phi - m_phiMin) / (m_phiMax - m_phiMin);
645 
646  // r
647  const auto r = get_r( loc(0), loc(1) );
648  if( r < m_rMin || r >= m_rMax ) return -1;
649  const int ir = rbins * (r - m_rMin) / (m_rMax - m_rMin);
650 
651  // z
652  const auto z = loc(2);
653  if( z < m_zMin || z >= m_zMax ) return -1;
654  const int iz = zbins * (z - m_zMin) / (m_zMax - m_zMin);
655 
656  // get index from matrix container
657  return m_matrix_container->get_cell_index( iphi, ir, iz );
658 
659 }
660 
662 {
663 
665 }
666 
668 {
669 
670  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
671  "ActsSurfaceMaps");
672  if(!m_surfMaps)
673  {
674  std::cout << PHWHERE << "No Acts surface maps, exiting"
675  << std::endl;
677  }
678 
679  m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,
680  "TRKR_CLUSTER");
681  if(!m_clusterContainer)
682  {
683  std::cout << PHWHERE << "No TRKR_CLUSTER node on node tree. Exiting."
684  << std::endl;
686  }
687 
688  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
689  if(!m_tGeometry)
690  {
691  std::cout << "ActsTrackingGeometry not on node tree. Exiting."
692  << std::endl;
693 
695  }
696 
697  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxSiliconMMTrackMap");
698 
699  if (!m_trackMap)
700  {
701  std::cout << PHWHERE << "SvtxSiliconMMTrackMap not on node tree. Exiting."
702  << std::endl;
704  }
705 
706 
708 }
709 
711 {
712 
713  m_histogramfile.reset( new TFile(m_histogramfilename.c_str(), "RECREATE") );
714  m_histogramfile->cd();
715 
716  const auto total_bins = m_matrix_container->get_grid_size();
717  h_beta = new TH2F("betadz",";tan#beta; #Deltaz [cm]",100,-0.5,0.5,100,-0.5,0.5);
718  h_alpha = new TH2F("alphardphi",";tan#alpha; r#Delta#phi [cm]", 100,-0.5,0.5,100,-0.5,0.5);
719  h_index = new TH1F("index",";index",total_bins, 0, total_bins);
720  h_rphiResid = new TH2F("rphiResid", ";r [cm]; #Deltar#phi [cm]",
721  60, 20, 80, 500, -2, 2);
722  h_zResid = new TH2F("zResid", ";z [cm]; #Deltaz [cm]",
723  200, -100, 100, 1000, -2, 2);
724  h_etaResid = new TH2F("etaResid", ";#eta;#Delta#eta",
725  20, -1, 1, 500, -0.2, 0.2);
726  h_etaResidLayer = new TH2F("etaResidLayer", ";r [cm]; #Delta#eta",
727  60, 20, 80, 500, -0.2, 0.2);
728  h_zResidLayer = new TH2F("zResidLayer", ";r [cm]; #Deltaz [cm]",
729  60, 20, 80, 1000, -2, 2);
730 
731  h_deltarphi_layer = new TH2F( "deltarphi_layer", ";layer; r.#Delta#phi_{track-cluster} (cm)", 57, 0, 57, 500, -2, 2 );
732  h_deltaz_layer = new TH2F( "deltaz_layer", ";layer; #Deltaz_{track-cluster} (cm)", 57, 0, 57, 100, -2, 2 );
733 
734  residTup = new TTree("residTree","tpc residual info");
735  residTup->Branch("tanAlpha",&tanAlpha,"tanAlpha/D");
736  residTup->Branch("tanBeta",&tanBeta,"tanBeta/D");
737  residTup->Branch("drphi",&drphi,"drphi/D");
738  residTup->Branch("dz",&dz,"dz/D");
739 // residTup->Branch("cell",&cell,"cell/I");
740  residTup->Branch("clusR",&clusR,"clusR/D");
741  residTup->Branch("clusPhi",&clusPhi,"clusPhi/D");
742  residTup->Branch("clusZ",&clusZ,"clusZ/D");
743  residTup->Branch("statePhi",&statePhi,"statePhi/D");
744  residTup->Branch("stateZ",&stateZ,"stateZ/D");
745  residTup->Branch("stateR",&stateR,"stateR/D");
746  residTup->Branch("stateRPhiErr",&stateRPhiErr,"stateRPhiErr/D");
747  residTup->Branch("stateZErr",&stateZErr,"stateZErr/D");
748  residTup->Branch("clusRPhiErr",&clusRPhiErr,"clusRPhiErr/D");
749  residTup->Branch("clusZErr",&clusZErr,"clusZErr/D");
750 // residTup->Branch("ir",&ir,"ir/I");
751 // residTup->Branch("iz",&iz,"iz/I");
752 // residTup->Branch("iphi",&iphi,"iphi/I");
753  residTup->Branch("cluskey",&cluskey,"cluskey/l");
754  residTup->Branch("event",&m_event,"event/I");
755 
756 }
757 
758 void PHTpcResiduals::setGridDimensions(const int phiBins, const int rBins, const int zBins)
759 { m_matrix_container->set_grid_dimensions( phiBins, rBins, zBins ); }