Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackResiduals.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackResiduals.cc
1 
2 #include "TrackResiduals.h"
3 
4 #include <trackbase/ClusterErrorPara.h>
5 #include <trackbase/InttDefs.h>
6 #include <trackbase/MvtxDefs.h>
7 #include <trackbase/TpcDefs.h>
8 #include <trackbase/TrackFitUtils.h>
9 #include <trackbase/TrkrCluster.h>
10 #include <trackbase/TrkrClusterContainer.h>
11 #include <trackbase/TrkrHit.h>
12 #include <trackbase/TrkrHitSet.h>
13 #include <trackbase/TrkrHitSetContainer.h>
14 
15 #include <g4detectors/PHG4CylinderGeomContainer.h>
16 #include <g4detectors/PHG4TpcCylinderGeom.h>
17 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
18 
19 #include <intt/CylinderGeomIntt.h>
20 
21 #include <micromegas/CylinderGeomMicromegas.h>
22 #include <micromegas/MicromegasDefs.h>
23 
24 #include <mvtx/CylinderGeom_Mvtx.h>
25 
26 #include <trackbase_historic/ActsTransformations.h>
27 #include <trackbase_historic/SvtxAlignmentState.h>
28 #include <trackbase_historic/SvtxAlignmentStateMap.h>
29 #include <trackbase_historic/SvtxTrack.h>
30 #include <trackbase_historic/SvtxTrackMap.h>
31 #include <trackbase_historic/TrackSeed.h>
32 
33 #include <ffarawobjects/Gl1RawHit.h>
34 
35 #include <globalvertex/GlobalVertex.h>
36 #include <globalvertex/GlobalVertexMap.h>
37 #include <globalvertex/SvtxVertex.h>
38 #include <globalvertex/SvtxVertexMap.h>
39 
40 #include <fun4all/Fun4AllReturnCodes.h>
41 
42 #include <phool/PHCompositeNode.h>
43 #include <phool/PHNodeIterator.h>
44 #include <phool/getClass.h>
45 
46 namespace
47 {
48  template <class T>
49  inline T square(const T& t)
50  {
51  return t * t;
52  }
53  template <class T>
54  inline T r(const T& x, const T& y)
55  {
56  return std::sqrt(square(x) + square(y));
57  }
58 
59  std::vector<TrkrDefs::cluskey> get_cluster_keys(SvtxTrack* track)
60  {
61  std::vector<TrkrDefs::cluskey> out;
62  for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
63  {
64  if (seed)
65  {
66  std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
67  }
68  }
69 
70  return out;
71  }
72 } // namespace
73 
74 //____________________________________________________________________________..
75 TrackResiduals::TrackResiduals(const std::string& name)
76  : SubsysReco(name)
77 {
78 }
79 
80 //____________________________________________________________________________..
81 TrackResiduals::~TrackResiduals()
82 {
83 }
84 
85 //____________________________________________________________________________..
86 int TrackResiduals::Init(PHCompositeNode*)
87 {
88  return Fun4AllReturnCodes::EVENT_OK;
89 }
90 
91 //____________________________________________________________________________..
92 int TrackResiduals::InitRun(PHCompositeNode*)
93 {
94  m_outfile = new TFile(m_outfileName.c_str(), "RECREATE");
95  createBranches();
96 
97  return Fun4AllReturnCodes::EVENT_OK;
98 }
99 void TrackResiduals::clearClusterStateVectors()
100 {
101  m_cluskeys.clear();
102  m_idealsurfcenterx.clear();
103  m_idealsurfcentery.clear();
104  m_idealsurfcenterz.clear();
105  m_idealsurfnormx.clear();
106  m_idealsurfnormy.clear();
107  m_idealsurfnormz.clear();
108  m_missurfcenterx.clear();
109  m_missurfcentery.clear();
110  m_missurfcenterz.clear();
111  m_missurfnormx.clear();
112  m_missurfnormy.clear();
113  m_missurfnormz.clear();
114  m_clusgxideal.clear();
115  m_clusgyideal.clear();
116  m_clusgzideal.clear();
117  m_missurfalpha.clear();
118  m_missurfbeta.clear();
119  m_missurfgamma.clear();
120  m_idealsurfalpha.clear();
121  m_idealsurfbeta.clear();
122  m_idealsurfgamma.clear();
123 
124  m_statelxglobderivdx.clear();
125  m_statelxglobderivdy.clear();
126  m_statelxglobderivdz.clear();
127  m_statelxglobderivdalpha.clear();
128  m_statelxglobderivdbeta.clear();
129  m_statelxglobderivdgamma.clear();
130 
131  m_statelxlocderivd0.clear();
132  m_statelxlocderivz0.clear();
133  m_statelxlocderivphi.clear();
134  m_statelxlocderivtheta.clear();
135  m_statelxlocderivqop.clear();
136 
137  m_statelzglobderivdx.clear();
138  m_statelzglobderivdy.clear();
139  m_statelzglobderivdz.clear();
140  m_statelzglobderivdalpha.clear();
141  m_statelzglobderivdbeta.clear();
142  m_statelzglobderivdgamma.clear();
143 
144  m_statelzlocderivd0.clear();
145  m_statelzlocderivz0.clear();
146  m_statelzlocderivphi.clear();
147  m_statelzlocderivtheta.clear();
148  m_statelzlocderivqop.clear();
149 
150  m_clusedge.clear();
151  m_clusoverlap.clear();
152  m_cluslx.clear();
153  m_cluslz.clear();
154  m_cluselx.clear();
155  m_cluselz.clear();
156  m_clusgx.clear();
157  m_clusgy.clear();
158  m_clusgz.clear();
159  m_cluslayer.clear();
160  m_clussize.clear();
161  m_clushitsetkey.clear();
162 
163  m_statelx.clear();
164  m_statelz.clear();
165  m_stateelx.clear();
166  m_stateelz.clear();
167  m_stategx.clear();
168  m_stategy.clear();
169  m_stategz.clear();
170  m_statepx.clear();
171  m_statepy.clear();
172  m_statepz.clear();
173 }
174 //____________________________________________________________________________..
175 int TrackResiduals::process_event(PHCompositeNode* topNode)
176 {
177  auto trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
178  auto clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
179  auto geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
180  auto vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
181  auto alignmentmap = findNode::getClass<SvtxAlignmentStateMap>(topNode, m_alignmentMapName);
182  auto hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
183  auto tpcGeom =
184  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
185  auto mvtxGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
186  auto inttGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
187  auto mmGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
188  if (!mmGeom)
189  {
190  mmGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS");
191  }
192  if (!trackmap or !clustermap or !geometry or !hitmap)
193  {
194  std::cout << "Missing node, can't continue" << std::endl;
195  return Fun4AllReturnCodes::ABORTEVENT;
196  }
197  auto gl1 = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");
198  if (gl1)
199  {
200  m_bco = gl1->get_bco();
201  auto lbshift = m_bco << 24;
202  m_bcotr = lbshift >> 24;
203  }
204  else
205  {
206  m_bco = std::numeric_limits<uint64_t>::quiet_NaN();
207  m_bcotr = std::numeric_limits<uint64_t>::quiet_NaN();
208  }
209  if (Verbosity() > 1)
210  {
211  std::cout << "Track map size is " << trackmap->size() << std::endl;
212  }
213 
214  if (m_doHits)
215  {
216  fillHitTree(hitmap, geometry, tpcGeom, mvtxGeom, inttGeom, mmGeom);
217  }
218 
219  if (m_doClusters)
220  {
221  fillClusterTree(clustermap, geometry);
222  }
223 
224  for (const auto& [key, track] : *trackmap)
225  {
226  if (!track)
227  {
228  continue;
229  }
230 
231  m_trackid = key;
232  m_crossing = track->get_crossing();
233  m_px = track->get_px();
234  m_py = track->get_py();
235  m_pz = track->get_pz();
236  m_pt = std::sqrt(square(m_px) + square(m_py));
237  m_eta = atanh(m_pz / std::sqrt(square(m_pt) + square(m_pz)));
238  m_phi = atan2(m_py, m_px);
239  float CVxx = track->get_error(3, 3);
240  float CVxy = track->get_error(3, 4);
241  float CVyy = track->get_error(4, 4);
242  m_deltapt = std::sqrt((CVxx * square(m_px) + 2 * CVxy * m_px * m_py + CVyy * square(m_py)) / (square(m_px) + square(m_py)));
243 
244  m_charge = track->get_charge();
245  m_quality = track->get_quality();
246  m_chisq = track->get_chisq();
247  m_ndf = track->get_ndf();
248  m_nmaps = 0;
249  m_nintt = 0;
250  m_ntpc = 0;
251  m_nmms = 0;
252  m_vertexid = track->get_vertex_id();
253  if (vertexmap)
254  {
255  auto vertexit = vertexmap->find(m_vertexid);
256  if (vertexit != vertexmap->end())
257  {
258  auto vertex = vertexit->second;
259  m_vx = vertex->get_x();
260  m_vy = vertex->get_y();
261  m_vz = vertex->get_z();
262  }
263  }
264  m_pcax = track->get_x();
265  m_pcay = track->get_y();
266  m_pcaz = track->get_z();
267 
268  clearClusterStateVectors();
269  if (Verbosity() > 1)
270  {
271  std::cout << "Track " << key << " has cluster/states"
272  << std::endl;
273  }
274 
275  if (!m_doAlignment)
276  {
277  std::vector<TrkrDefs::cluskey> keys;
278  for (const auto& ckey : get_cluster_keys(track))
279  {
280  if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::tpcId)
281  {
282  keys.push_back(ckey);
283  }
284  }
285  if (m_zeroField)
286  {
287  lineFitClusters(keys, geometry, clustermap);
288  }
289  for (const auto& ckey : get_cluster_keys(track))
290  {
291  fillClusterBranches(ckey, track, topNode);
292  }
293  }
294  m_nhits = m_nmaps + m_nintt + m_ntpc + m_nmms;
295 
296  if (m_doAlignment)
297  {
299  clearClusterStateVectors();
300 
301  if (alignmentmap and alignmentmap->find(key) != alignmentmap->end())
302  {
303  auto& statevec = alignmentmap->find(key)->second;
304 
305  for (auto& state : statevec)
306  {
307  auto ckey = state->get_cluster_key();
308 
309  fillClusterBranches(ckey, track, topNode);
310 
311  auto& globderivs = state->get_global_derivative_matrix();
312  auto& locderivs = state->get_local_derivative_matrix();
313 
314  m_statelxglobderivdalpha.push_back(globderivs(0, 0));
315  m_statelxglobderivdbeta.push_back(globderivs(0, 1));
316  m_statelxglobderivdgamma.push_back(globderivs(0, 2));
317  m_statelxglobderivdx.push_back(globderivs(0, 3));
318  m_statelxglobderivdy.push_back(globderivs(0, 4));
319  m_statelxglobderivdz.push_back(globderivs(0, 5));
320 
321  m_statelzglobderivdalpha.push_back(globderivs(1, 0));
322  m_statelzglobderivdbeta.push_back(globderivs(1, 1));
323  m_statelzglobderivdgamma.push_back(globderivs(1, 2));
324  m_statelzglobderivdx.push_back(globderivs(1, 3));
325  m_statelzglobderivdy.push_back(globderivs(1, 4));
326  m_statelzglobderivdz.push_back(globderivs(1, 5));
327 
328  m_statelxlocderivd0.push_back(locderivs(0, 0));
329  m_statelxlocderivz0.push_back(locderivs(0, 1));
330  m_statelxlocderivphi.push_back(locderivs(0, 2));
331  m_statelxlocderivtheta.push_back(locderivs(0, 3));
332  m_statelxlocderivqop.push_back(locderivs(0, 4));
333 
334  m_statelzlocderivd0.push_back(locderivs(1, 0));
335  m_statelzlocderivz0.push_back(locderivs(1, 1));
336  m_statelzlocderivphi.push_back(locderivs(1, 2));
337  m_statelzlocderivtheta.push_back(locderivs(1, 3));
338  m_statelzlocderivqop.push_back(locderivs(1, 4));
339  }
340  }
341  }
342  m_tree->Fill();
343  }
344 
345  m_event++;
346  return Fun4AllReturnCodes::EVENT_OK;
347 }
348 
349 float TrackResiduals::convertTimeToZ(ActsGeometry* geometry, TrkrDefs::cluskey cluster_key, TrkrCluster* cluster)
350 {
351  // must convert local Y from cluster average time of arival to local cluster z position
352  double drift_velocity = geometry->get_drift_velocity();
353  double zdriftlength = cluster->getLocalY() * drift_velocity;
354  double surfCenterZ = 52.89; // 52.89 is where G4 thinks the surface center is
355  double zloc = surfCenterZ - zdriftlength; // converts z drift length to local z position in the TPC in north
356  unsigned int side = TpcDefs::getSide(cluster_key);
357  if (side == 0) zloc = -zloc;
358  float z = zloc; // in cm
359 
360  return z;
361 }
362 void TrackResiduals::lineFitClusters(std::vector<TrkrDefs::cluskey>& keys,
363  ActsGeometry* geometry,
364  TrkrClusterContainer* clusters)
365 {
366  std::vector<Acts::Vector3> clusPos;
367  TrackFitUtils::getTrackletClusters(geometry, clusters,
368  clusPos, keys);
369  TrackFitUtils::position_vector_t xypoints, rzpoints;
370  for (auto& pos : clusPos)
371  {
372  xypoints.push_back(std::make_pair(pos.x(), pos.y()));
373  float clusr = r(pos.x(), pos.y());
374  if (pos.y() < 0) clusr *= -1;
375  rzpoints.push_back(std::make_pair(pos.z(), clusr));
376  }
377 
378  auto xyparams = TrackFitUtils::line_fit(xypoints);
379  auto rzparams = TrackFitUtils::line_fit(rzpoints);
380  m_xyint = std::get<1>(xyparams);
381  m_xyslope = std::get<0>(xyparams);
382  m_rzint = std::get<1>(rzparams);
383  m_rzslope = std::get<0>(rzparams);
384 }
385 
386 void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters,
387  ActsGeometry* geometry)
388 {
389  for (auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId,
390  TrkrDefs::TrkrId::tpcId, TrkrDefs::TrkrId::micromegasId})
391  {
392  for (const auto& hitsetkey : clusters->getHitSetKeys(det))
393  {
394  m_scluslayer = TrkrDefs::getLayer(hitsetkey);
395  auto range = clusters->getClusters(hitsetkey);
396  for (auto iter = range.first; iter != range.second; ++iter)
397  {
398  auto key = iter->first;
399  auto cluster = clusters->findCluster(key);
400  auto glob = geometry->getGlobalPosition(key, cluster);
401  m_sclusgx = glob.x();
402  m_sclusgy = glob.y();
403  m_sclusgz = glob.z();
404  m_sclusgr = r(m_sclusgx, m_sclusgy);
405  m_sclusphi = atan2(glob.y(), glob.x());
406  m_scluseta = acos(glob.z() / std::sqrt(square(glob.x()) + square(glob.y()) + square(glob.z())));
407  m_adc = cluster->getAdc();
408  m_clusmaxadc = cluster->getMaxAdc();
409  m_scluslx = cluster->getLocalX();
410  m_scluslz = cluster->getLocalY();
411  auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster, m_sclusgr, key);
412  m_phisize = cluster->getPhiSize();
413  m_zsize = cluster->getZSize();
414  m_scluselx = std::sqrt(para_errors.first);
415  m_scluselz = std::sqrt(para_errors.second);
416 
418  switch (det)
419  {
420  case TrkrDefs::TrkrId::mvtxId:
421  m_staveid = MvtxDefs::getStaveId(key);
422  m_chipid = MvtxDefs::getChipId(key);
423  m_strobeid = MvtxDefs::getStrobeId(key);
424 
425  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
426  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
427  m_timebucket = std::numeric_limits<int>::quiet_NaN();
428  m_clussector = std::numeric_limits<int>::quiet_NaN();
429  m_side = std::numeric_limits<int>::quiet_NaN();
430  m_segtype = std::numeric_limits<int>::quiet_NaN();
431  m_tileid = std::numeric_limits<int>::quiet_NaN();
432  break;
433  case TrkrDefs::TrkrId::inttId:
434  m_ladderzid = InttDefs::getLadderZId(key);
435  m_ladderphiid = InttDefs::getLadderPhiId(key);
436  m_timebucket = InttDefs::getTimeBucketId(key);
437 
438  m_staveid = std::numeric_limits<int>::quiet_NaN();
439  m_chipid = std::numeric_limits<int>::quiet_NaN();
440  m_strobeid = std::numeric_limits<int>::quiet_NaN();
441  m_clussector = std::numeric_limits<int>::quiet_NaN();
442  m_side = std::numeric_limits<int>::quiet_NaN();
443  m_segtype = std::numeric_limits<int>::quiet_NaN();
444  m_tileid = std::numeric_limits<int>::quiet_NaN();
445  break;
446  case TrkrDefs::TrkrId::tpcId:
447  m_clussector = TpcDefs::getSectorId(key);
448  m_side = TpcDefs::getSide(key);
449 
450  m_staveid = std::numeric_limits<int>::quiet_NaN();
451  m_chipid = std::numeric_limits<int>::quiet_NaN();
452  m_strobeid = std::numeric_limits<int>::quiet_NaN();
453  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
454  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
455  m_timebucket = std::numeric_limits<int>::quiet_NaN();
456  m_segtype = std::numeric_limits<int>::quiet_NaN();
457  m_tileid = std::numeric_limits<int>::quiet_NaN();
458  break;
459  case TrkrDefs::TrkrId::micromegasId:
460  m_segtype = (int) MicromegasDefs::getSegmentationType(key);
461  m_tileid = MicromegasDefs::getTileId(key);
462 
463  m_staveid = std::numeric_limits<int>::quiet_NaN();
464  m_chipid = std::numeric_limits<int>::quiet_NaN();
465  m_strobeid = std::numeric_limits<int>::quiet_NaN();
466  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
467  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
468  m_timebucket = std::numeric_limits<int>::quiet_NaN();
469  m_clussector = std::numeric_limits<int>::quiet_NaN();
470  m_side = std::numeric_limits<int>::quiet_NaN();
471  break;
472  default:
473  break;
474  }
475 
476  m_clustree->Fill();
477  }
478  }
479  }
480 }
481 
482 //____________________________________________________________________________..
483 int TrackResiduals::End(PHCompositeNode*)
484 {
485  m_outfile->cd();
486  m_tree->Write();
487  if (m_doClusters)
488  {
489  m_clustree->Write();
490  }
491  if (m_doHits)
492  {
493  m_hittree->Write();
494  }
495  m_outfile->Close();
496 
497  return Fun4AllReturnCodes::EVENT_OK;
498 }
499 void TrackResiduals::fillHitTree(TrkrHitSetContainer* hitmap,
500  ActsGeometry* geometry,
501  PHG4TpcCylinderGeomContainer* tpcGeom,
502  PHG4CylinderGeomContainer* mvtxGeom,
503  PHG4CylinderGeomContainer* inttGeom,
504  PHG4CylinderGeomContainer* mmGeom)
505 {
506  if (!tpcGeom or !mvtxGeom or !inttGeom or !mmGeom)
507  {
508  std::cout << PHWHERE << "missing hit map, can't continue with hit tree"
509  << std::endl;
510  return;
511  }
512  TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
513  for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
514  hitsetiter != all_hitsets.second;
515  ++hitsetiter)
516  {
517  m_hitsetkey = hitsetiter->first;
518  TrkrHitSet* hitset = hitsetiter->second;
519 
520  m_hitlayer = TrkrDefs::getLayer(m_hitsetkey);
521  auto det = TrkrDefs::getTrkrId(m_hitsetkey);
523  switch (det)
524  {
525  case TrkrDefs::TrkrId::mvtxId:
526  {
527  m_staveid = MvtxDefs::getStaveId(m_hitsetkey);
528  m_chipid = MvtxDefs::getChipId(m_hitsetkey);
529  m_strobeid = MvtxDefs::getStrobeId(m_hitsetkey);
530 
531  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
532  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
533  m_timebucket = std::numeric_limits<int>::quiet_NaN();
534  m_sector = std::numeric_limits<int>::quiet_NaN();
535  m_side = std::numeric_limits<int>::quiet_NaN();
536  m_segtype = std::numeric_limits<int>::quiet_NaN();
537  m_tileid = std::numeric_limits<int>::quiet_NaN();
538  break;
539  }
540  case TrkrDefs::TrkrId::inttId:
541  {
542  m_ladderzid = InttDefs::getLadderZId(m_hitsetkey);
543  m_ladderphiid = InttDefs::getLadderPhiId(m_hitsetkey);
544  m_timebucket = InttDefs::getTimeBucketId(m_hitsetkey);
545 
546  m_staveid = std::numeric_limits<int>::quiet_NaN();
547  m_chipid = std::numeric_limits<int>::quiet_NaN();
548  m_strobeid = std::numeric_limits<int>::quiet_NaN();
549  m_sector = std::numeric_limits<int>::quiet_NaN();
550  m_side = std::numeric_limits<int>::quiet_NaN();
551  m_segtype = std::numeric_limits<int>::quiet_NaN();
552  m_tileid = std::numeric_limits<int>::quiet_NaN();
553  break;
554  }
555  case TrkrDefs::TrkrId::tpcId:
556  {
557  m_sector = TpcDefs::getSectorId(m_hitsetkey);
558  m_side = TpcDefs::getSide(m_hitsetkey);
559 
560  m_staveid = std::numeric_limits<int>::quiet_NaN();
561  m_chipid = std::numeric_limits<int>::quiet_NaN();
562  m_strobeid = std::numeric_limits<int>::quiet_NaN();
563  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
564  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
565  m_timebucket = std::numeric_limits<int>::quiet_NaN();
566  m_segtype = std::numeric_limits<int>::quiet_NaN();
567  m_tileid = std::numeric_limits<int>::quiet_NaN();
568 
569  break;
570  }
571  case TrkrDefs::TrkrId::micromegasId:
572  {
573  m_segtype = (int) MicromegasDefs::getSegmentationType(m_hitsetkey);
574  m_tileid = MicromegasDefs::getTileId(m_hitsetkey);
575 
576  m_staveid = std::numeric_limits<int>::quiet_NaN();
577  m_chipid = std::numeric_limits<int>::quiet_NaN();
578  m_strobeid = std::numeric_limits<int>::quiet_NaN();
579  m_ladderzid = std::numeric_limits<int>::quiet_NaN();
580  m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
581  m_timebucket = std::numeric_limits<int>::quiet_NaN();
582  m_sector = std::numeric_limits<int>::quiet_NaN();
583  m_side = std::numeric_limits<int>::quiet_NaN();
584  break;
585  }
586  default:
587  break;
588  }
589 
590  // Got all stave/ladder/sector/tile info, now get the actual hit info
591  auto hitrangei = hitset->getHits();
592  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
593  hitr != hitrangei.second;
594  ++hitr)
595  {
596  auto hitkey = hitr->first;
597  auto hit = hitr->second;
598  m_adc = hit->getAdc();
599 
600  switch (det)
601  {
602  case TrkrDefs::TrkrId::mvtxId:
603  {
604  m_row = MvtxDefs::getRow(hitkey);
605  m_col = MvtxDefs::getCol(hitkey);
606  auto layergeom = dynamic_cast<CylinderGeom_Mvtx*>(mvtxGeom->GetLayerGeom(m_hitlayer));
607  auto local_coords = layergeom->get_local_coords_from_pixel(m_row, m_col);
608  TVector2 local;
609  local.SetX(local_coords.X());
610  local.SetY(local_coords.Z());
611  auto surf = geometry->maps().getSiliconSurface(m_hitsetkey);
612  auto glob = layergeom->get_world_from_local_coords(surf, geometry, local);
613  m_hitgx = glob.X();
614  m_hitgy = glob.Y();
615  m_hitgz = glob.Z();
616 
617  m_segtype = std::numeric_limits<int>::quiet_NaN();
618  m_tileid = std::numeric_limits<int>::quiet_NaN();
619  m_strip = std::numeric_limits<int>::quiet_NaN();
620  m_hitpad = std::numeric_limits<int>::quiet_NaN();
621  m_hittbin = std::numeric_limits<int>::quiet_NaN();
622  break;
623  }
624  case TrkrDefs::TrkrId::inttId:
625  {
626  m_row = InttDefs::getRow(hitkey);
627  m_col = InttDefs::getCol(hitkey);
628  auto geom = dynamic_cast<CylinderGeomIntt*>(inttGeom->GetLayerGeom(m_hitlayer));
629  double local_hit_loc[3] = {0, 0, 0};
630  geom->find_strip_center_localcoords(m_ladderzid, m_row, m_col, local_hit_loc);
631  auto surf = geometry->maps().getSiliconSurface(m_hitsetkey);
632  TVector2 local;
633  local.SetX(local_hit_loc[1]);
634  local.SetY(local_hit_loc[2]);
635  auto glob = geom->get_world_from_local_coords(surf, geometry, local);
636  m_hitgx = glob.X();
637  m_hitgy = glob.Y();
638  m_hitgz = glob.Z();
639  m_segtype = std::numeric_limits<int>::quiet_NaN();
640  m_tileid = std::numeric_limits<int>::quiet_NaN();
641  m_strip = std::numeric_limits<int>::quiet_NaN();
642  m_hitpad = std::numeric_limits<int>::quiet_NaN();
643  m_hittbin = std::numeric_limits<int>::quiet_NaN();
644  break;
645  }
646  case TrkrDefs::TrkrId::tpcId:
647  {
648  m_row = std::numeric_limits<int>::quiet_NaN();
649  m_col = std::numeric_limits<int>::quiet_NaN();
650  m_segtype = std::numeric_limits<int>::quiet_NaN();
651  m_tileid = std::numeric_limits<int>::quiet_NaN();
652  m_strip = std::numeric_limits<int>::quiet_NaN();
653 
654  m_hitpad = TpcDefs::getPad(hitkey);
655  m_hittbin = TpcDefs::getTBin(hitkey);
656 
657  auto geoLayer = tpcGeom->GetLayerCellGeom(m_hitlayer);
658  auto phi = geoLayer->get_phicenter(m_hitpad);
659  auto radius = geoLayer->get_radius();
660  float AdcClockPeriod = 53.0; // ns (?)
661  double zdriftlength = m_hittbin * geometry->get_drift_velocity() * AdcClockPeriod;
662  unsigned short NTBins = (unsigned short) geoLayer->get_zbins();
663  double tdriftmax = AdcClockPeriod * NTBins / 2.0;
664  m_hitgz = (tdriftmax * geometry->get_drift_velocity()) - zdriftlength;
665  if (m_side == 0)
666  {
667  m_hitgz *= -1;
668  }
669  m_hitgx = radius * std::cos(phi);
670  m_hitgy = radius * std::sin(phi);
671  break;
672  }
673  case TrkrDefs::TrkrId::micromegasId:
674  {
675  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(mmGeom->GetLayerGeom(m_hitlayer));
676  m_strip = MicromegasDefs::getStrip(hitkey);
677  const auto global_coord = layergeom->get_world_coordinates(m_tileid, geometry, m_strip);
678  m_hitgx = global_coord.X();
679  m_hitgy = global_coord.Y();
680  m_hitgz = global_coord.Z();
681  m_row = std::numeric_limits<int>::quiet_NaN();
682  m_col = std::numeric_limits<int>::quiet_NaN();
683  m_segtype = std::numeric_limits<int>::quiet_NaN();
684  m_tileid = std::numeric_limits<int>::quiet_NaN();
685  m_hitpad = std::numeric_limits<int>::quiet_NaN();
686  m_hittbin = std::numeric_limits<int>::quiet_NaN();
687  }
688  default:
689  break;
690  }
691 
692  m_hittree->Fill();
693  }
694  }
695 }
696 
697 void TrackResiduals::fillClusterBranches(TrkrDefs::cluskey ckey, SvtxTrack* track,
698  PHCompositeNode* topNode)
699 {
700  auto clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
701  auto geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
702 
703  ActsTransformations transformer;
704  TrkrCluster* cluster = clustermap->findCluster(ckey);
705  switch (TrkrDefs::getTrkrId(ckey))
706  {
707  case TrkrDefs::mvtxId:
708  m_nmaps++;
709  break;
710  case TrkrDefs::inttId:
711  m_nintt++;
712  break;
713  case TrkrDefs::tpcId:
714  m_ntpc++;
715  break;
716  case TrkrDefs::micromegasId:
717  m_nmms++;
718  break;
719  }
720 
721  Acts::Vector3 clusglob = geometry->getGlobalPosition(ckey, cluster);
722 
723  SvtxTrackState* state = nullptr;
724 
725  for (auto state_iter = track->begin_states();
726  state_iter != track->end_states();
727  ++state_iter)
728  {
729  SvtxTrackState* tstate = state_iter->second;
730  auto stateckey = tstate->get_cluskey();
731  if (stateckey == ckey)
732  {
733  state = tstate;
734  break;
735  }
736  }
737 
738  m_cluskeys.push_back(ckey);
739 
741  m_cluslx.push_back(cluster->getLocalX());
742  m_clusedge.push_back(cluster->getEdge());
743  m_clusoverlap.push_back(cluster->getOverlap());
744  float clusz = cluster->getLocalY();
745 
746  if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::tpcId)
747  {
748  clusz = convertTimeToZ(geometry, ckey, cluster);
749  }
750 
751  m_cluslz.push_back(clusz);
752  float clusr = r(clusglob.x(), clusglob.y());
753  auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster,
754  clusr, ckey);
755  m_cluselx.push_back(sqrt(para_errors.first));
756  m_cluselz.push_back(sqrt(para_errors.second));
757  m_clusgx.push_back(clusglob.x());
758  m_clusgy.push_back(clusglob.y());
759  m_clusgz.push_back(clusglob.z());
760  m_cluslayer.push_back(TrkrDefs::getLayer(ckey));
761  m_clussize.push_back(cluster->getPhiSize() + cluster->getZSize());
762  m_clushitsetkey.push_back(TrkrDefs::getHitSetKeyFromClusKey(ckey));
763 
764  if (Verbosity() > 1)
765  {
766  std::cout << "Track state/clus in layer "
767  << (unsigned int) TrkrDefs::getLayer(ckey) << " with pos "
768  << clusglob.transpose() << std::endl;
769  }
770  if (!state)
771  {
772  if (m_zeroField)
773  {
774  fillStatesWithLineFit(ckey, cluster, geometry);
775  }
779  m_idealsurfalpha.push_back(NAN);
780  m_idealsurfbeta.push_back(NAN);
781  m_idealsurfgamma.push_back(NAN);
782  m_missurfalpha.push_back(NAN);
783  m_missurfbeta.push_back(NAN);
784  m_missurfgamma.push_back(NAN);
785  m_idealsurfcenterx.push_back(NAN);
786  m_idealsurfcentery.push_back(NAN);
787  m_idealsurfcenterz.push_back(NAN);
788  m_idealsurfnormx.push_back(NAN);
789  m_idealsurfnormy.push_back(NAN);
790  m_idealsurfnormz.push_back(NAN);
791  m_missurfcenterx.push_back(NAN);
792  m_missurfcentery.push_back(NAN);
793  m_missurfcenterz.push_back(NAN);
794  m_missurfnormx.push_back(NAN);
795  m_missurfnormy.push_back(NAN);
796  m_missurfnormz.push_back(NAN);
797  m_clusgxideal.push_back(NAN);
798  m_clusgyideal.push_back(NAN);
799  m_clusgzideal.push_back(NAN);
800  m_statepx.push_back(NAN);
801  m_statepy.push_back(NAN);
802  m_statepz.push_back(NAN);
803  m_statepl.push_back(NAN);
804  return;
805  }
806 
807  auto surf = geometry->maps().getSurface(ckey, cluster);
808  Acts::Vector3 stateglob(state->get_x(), state->get_y(), state->get_z());
809  Acts::Vector2 stateloc;
810  auto misaligncenter = surf->center(geometry->geometry().getGeoContext());
811  auto misalignnorm = -1 * surf->normal(geometry->geometry().getGeoContext());
812  auto misrot = surf->transform(geometry->geometry().getGeoContext()).rotation();
813  auto result = surf->globalToLocal(geometry->geometry().getGeoContext(),
814  stateglob * Acts::UnitConstants::cm,
815  misalignnorm);
816 
817  float mgamma = atan2(-misrot(1, 0), misrot(0, 0));
818  float mbeta = -asin(misrot(0, 1));
819  float malpha = atan2(misrot(1, 1), misrot(2, 1));
820 
822  alignmentTransformationContainer::use_alignment = false;
823  auto idealcenter = surf->center(geometry->geometry().getGeoContext());
824  auto idealnorm = -1 * surf->normal(geometry->geometry().getGeoContext());
825  Acts::Vector3 ideal_local(cluster->getLocalX(), clusz, 0.0);
826  Acts::Vector3 ideal_glob = surf->transform(geometry->geometry().getGeoContext()) * (ideal_local * Acts::UnitConstants::cm);
827  auto idealrot = surf->transform(geometry->geometry().getGeoContext()).rotation();
828 
837  float igamma = atan2(-idealrot(1, 0), idealrot(0, 0));
838  float ibeta = -asin(idealrot(0, 1));
839  float ialpha = atan2(idealrot(1, 1), idealrot(2, 1));
840 
841  alignmentTransformationContainer::use_alignment = true;
842 
843  idealcenter /= Acts::UnitConstants::cm;
844  misaligncenter /= Acts::UnitConstants::cm;
845  ideal_glob /= Acts::UnitConstants::cm;
846 
847  m_idealsurfalpha.push_back(ialpha);
848  m_idealsurfbeta.push_back(ibeta);
849  m_idealsurfgamma.push_back(igamma);
850  m_missurfalpha.push_back(malpha);
851  m_missurfbeta.push_back(mbeta);
852  m_missurfgamma.push_back(mgamma);
853 
854  m_idealsurfcenterx.push_back(idealcenter.x());
855  m_idealsurfcentery.push_back(idealcenter.y());
856  m_idealsurfcenterz.push_back(idealcenter.z());
857  m_idealsurfnormx.push_back(idealnorm.x());
858  m_idealsurfnormy.push_back(idealnorm.y());
859  m_idealsurfnormz.push_back(idealnorm.z());
860  m_missurfcenterx.push_back(misaligncenter.x());
861  m_missurfcentery.push_back(misaligncenter.y());
862  m_missurfcenterz.push_back(misaligncenter.z());
863  m_missurfnormx.push_back(misalignnorm.x());
864  m_missurfnormy.push_back(misalignnorm.y());
865  m_missurfnormz.push_back(misalignnorm.z());
866  m_clusgxideal.push_back(ideal_glob.x());
867  m_clusgyideal.push_back(ideal_glob.y());
868  m_clusgzideal.push_back(ideal_glob.z());
869 
870  if (result.ok())
871  {
872  stateloc = result.value() / Acts::UnitConstants::cm;
873  }
874  else
875  {
877  Acts::Vector3 loct = surf->transform(geometry->geometry().getGeoContext()).inverse() * (stateglob * Acts::UnitConstants::cm);
878  loct /= Acts::UnitConstants::cm;
879  stateloc(0) = loct(0);
880  stateloc(1) = loct(1);
881  }
882 
883  const auto actscov =
884  transformer.rotateSvtxTrackCovToActs(state);
885 
886  m_statelx.push_back(stateloc(0));
887  m_statelz.push_back(stateloc(1));
888  m_stateelx.push_back(std::sqrt(actscov(Acts::eBoundLoc0, Acts::eBoundLoc0)) / Acts::UnitConstants::cm);
889  m_stateelz.push_back(std::sqrt(actscov(Acts::eBoundLoc1, Acts::eBoundLoc1)) / Acts::UnitConstants::cm);
890  m_stategx.push_back(state->get_x());
891  m_stategy.push_back(state->get_y());
892  m_stategz.push_back(state->get_z());
893  m_statepx.push_back(state->get_px());
894  m_statepy.push_back(state->get_py());
895  m_statepz.push_back(state->get_pz());
896  m_statepl.push_back(state->get_pathlength());
897 }
898 void TrackResiduals::fillStatesWithLineFit(const TrkrDefs::cluskey& key,
899  TrkrCluster* cluster, ActsGeometry* geometry)
900 {
901  auto surf = geometry->maps().getSurface(key, cluster);
902 
906  float x1 = -1;
907  float x2 = 5;
908  float y1 = m_xyslope * x1 + m_xyint;
909  float y2 = m_xyslope * x2 + m_xyint;
910 
913  float r1 = r(x1, y1);
914  float r2 = r(x2, y2);
915  if (y1 < 0) r1 *= -1;
916  if (y2 < 0) r2 *= -1;
917  float z1 = (r1 - m_rzint) / m_rzslope;
918  float z2 = (r2 - m_rzint) / m_rzslope;
919  Acts::Vector3 v1(x1, y1, z1), v2(x2, y2, z2);
920 
921  Acts::Vector3 surfcenter = surf->center(geometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
922  Acts::Vector3 surfnorm = surf->normal(geometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
923 
924  Acts::Vector3 u = v2 - v1;
925  float dot = surfnorm.dot(u);
926  if (abs(dot) > 1e-6)
927  {
928  Acts::Vector3 w = v1 - surfcenter;
929  float fac = -surfnorm.dot(w) / dot;
930  u *= fac;
931  Acts::Vector3 intersection = v1 + u;
932 
933  auto locstateres = surf->globalToLocal(geometry->geometry().getGeoContext(),
934  intersection * Acts::UnitConstants::cm,
935  surfnorm);
936  if (locstateres.ok())
937  {
938  Acts::Vector2 loc = locstateres.value() / Acts::UnitConstants::cm;
939  m_statelx.push_back(loc(0));
940  m_statelz.push_back(loc(1));
941  }
942  else
943  {
944  Acts::Vector3 loct = surf->transform(geometry->geometry().getGeoContext()).inverse() * (intersection * Acts::UnitConstants::cm);
945  loct /= Acts::UnitConstants::cm;
946  m_statelx.push_back(loct(0));
947  m_statelz.push_back(loct(1));
948  }
949  m_stategx.push_back(intersection.x());
950  m_stategy.push_back(intersection.y());
951  m_stategz.push_back(intersection.z());
952  }
953  else
954  {
957  m_statelx.push_back(NAN);
958  m_statelz.push_back(NAN);
959  m_stategx.push_back(NAN);
960  m_stategy.push_back(NAN);
961  m_stategz.push_back(NAN);
962  }
963 }
964 void TrackResiduals::createBranches()
965 {
966  m_hittree = new TTree("hittree", "A tree with all hits");
967  m_hittree->Branch("event", &m_event, "m_event/I");
968  m_hittree->Branch("gl1bco", &m_bco, "m_bco/l");
969  m_hittree->Branch("trbco", &m_bcotr, "m_bcotr/l");
970  m_hittree->Branch("hitsetkey", &m_hitsetkey, "m_hitsetkey/i");
971  m_hittree->Branch("gx", &m_hitgx, "m_hitgx/F");
972  m_hittree->Branch("gy", &m_hitgy, "m_hitgy/F");
973  m_hittree->Branch("gz", &m_hitgz, "m_hitgz/F");
974  m_hittree->Branch("layer", &m_hitlayer, "m_hitlayer/I");
975  m_hittree->Branch("sector", &m_sector, "m_sector/I");
976  m_hittree->Branch("side", &m_side, "m_side/I");
977  m_hittree->Branch("stave", &m_staveid, "m_staveid/I");
978  m_hittree->Branch("chip", &m_chipid, "m_chipid/I");
979  m_hittree->Branch("strobe", &m_strobeid, "m_strobeid/I");
980  m_hittree->Branch("ladderz", &m_ladderzid, "m_ladderzid/I");
981  m_hittree->Branch("ladderphi", m_ladderphiid, "m_ladderphiid/I");
982  m_hittree->Branch("timebucket", &m_timebucket, "m_timebucket/I");
983  m_hittree->Branch("pad", &m_hitpad, "m_hitpad/I");
984  m_hittree->Branch("tbin", &m_hittbin, "m_hittbin/I");
985  m_hittree->Branch("col", &m_col, "m_col/I");
986  m_hittree->Branch("row", &m_row, "m_row/I");
987  m_hittree->Branch("segtype", &m_segtype, "m_segtype/I");
988  m_hittree->Branch("tile", &m_tileid, "m_tileid/I");
989  m_hittree->Branch("strip", &m_strip, "m_strip/I");
990  m_hittree->Branch("adc", &m_adc, "m_adc/F");
991 
992  m_clustree = new TTree("clustertree", "A tree with all clusters");
993  m_clustree->Branch("event", &m_event, "m_event/I");
994  m_clustree->Branch("gl1bco", &m_bco, "m_bco/l");
995  m_clustree->Branch("trbco", &m_bcotr, "m_bcotr/l");
996  m_clustree->Branch("lx", &m_scluslx, "m_scluslx/F");
997  m_clustree->Branch("lz", &m_scluslz, "m_scluslz/F");
998  m_clustree->Branch("gx", &m_sclusgx, "m_sclusgx/F");
999  m_clustree->Branch("gy", &m_sclusgy, "m_sclusgy/F");
1000  m_clustree->Branch("gz", &m_sclusgz, "m_sclusgz/F");
1001  m_clustree->Branch("r", &m_sclusgr, "m_sclusgr/F");
1002  m_clustree->Branch("phi", &m_sclusphi, "m_sclusphi/F");
1003  m_clustree->Branch("eta", &m_scluseta, "m_scluseta/F");
1004  m_clustree->Branch("adc", &m_adc, "m_adc/F");
1005  m_clustree->Branch("phisize", &m_phisize, "m_phisize/I");
1006  m_clustree->Branch("zsize", &m_zsize, "m_zsize/I");
1007  m_clustree->Branch("layer", &m_scluslayer, "m_scluslayer/I");
1008  m_clustree->Branch("erphi", &m_scluselx, "m_scluselx/F");
1009  m_clustree->Branch("ez", &m_scluselz, "m_scluselz/F");
1010  m_clustree->Branch("maxadc", &m_clusmaxadc, "m_clusmaxadc/F");
1011  m_clustree->Branch("sector", &m_clussector, "m_clussector/I");
1012  m_clustree->Branch("side", &m_side, "m_side/I");
1013  m_clustree->Branch("stave", &m_staveid, "m_staveid/I");
1014  m_clustree->Branch("chip", &m_chipid, "m_chipid/I");
1015  m_clustree->Branch("strobe", &m_strobeid, "m_strobeid/I");
1016  m_clustree->Branch("ladderz", &m_ladderzid, "m_ladderzid/I");
1017  m_clustree->Branch("ladderphi", m_ladderphiid, "m_ladderphiid/I");
1018  m_clustree->Branch("timebucket", &m_timebucket, "m_timebucket/I");
1019  m_clustree->Branch("segtype", &m_segtype, "m_segtype/I");
1020  m_clustree->Branch("tile", &m_tileid, "m_tileid/I");
1021 
1022  m_tree = new TTree("residualtree", "A tree with track, cluster, and state info");
1023  m_tree->Branch("event", &m_event, "m_event/I");
1024  m_tree->Branch("trackid", &m_trackid, "m_trackid/I");
1025  m_tree->Branch("gl1bco", &m_bco, "m_bco/l");
1026  m_tree->Branch("trbco", &m_bcotr, "m_bcotr/l");
1027  m_tree->Branch("crossing", &m_crossing, "m_crossing/I");
1028  m_tree->Branch("px", &m_px, "m_px/F");
1029  m_tree->Branch("py", &m_py, "m_py/F");
1030  m_tree->Branch("pz", &m_pz, "m_pz/F");
1031  m_tree->Branch("pt", &m_pt, "m_pt/F");
1032  m_tree->Branch("eta", &m_eta, "m_eta/F");
1033  m_tree->Branch("phi", &m_phi, "m_phi/F");
1034  m_tree->Branch("deltapt", &m_deltapt, "m_deltapt/F");
1035  m_tree->Branch("charge", &m_charge, "m_charge/I");
1036  m_tree->Branch("quality", &m_quality, "m_quality/F");
1037  m_tree->Branch("ndf", &m_ndf, "m_ndf/F");
1038  m_tree->Branch("nhits", &m_nhits, "m_nhits/I");
1039  m_tree->Branch("nmaps", &m_nmaps, "m_nmaps/I");
1040  m_tree->Branch("nintt", &m_nintt, "m_nintt/I");
1041  m_tree->Branch("ntpc", &m_ntpc, "m_ntpc/I");
1042  m_tree->Branch("nmms", &m_nmms, "m_nmms/I");
1043  m_tree->Branch("vertexid", &m_vertexid, "m_vertexid/I");
1044  m_tree->Branch("vx", &m_vx, "m_vx/F");
1045  m_tree->Branch("vy", &m_vy, "m_vy/F");
1046  m_tree->Branch("vz", &m_vz, "m_vz/F");
1047  m_tree->Branch("pcax", &m_pcax, "m_pcax/F");
1048  m_tree->Branch("pcay", &m_pcay, "m_pcay/F");
1049  m_tree->Branch("pcaz", &m_pcaz, "m_pcaz/F");
1050  m_tree->Branch("rzslope", &m_rzslope, "m_rzslope/F");
1051  m_tree->Branch("xyslope", &m_xyslope, "m_xyslope/F");
1052  m_tree->Branch("rzint", &m_rzint, "m_rzint/F");
1053  m_tree->Branch("xyint", &m_xyint, "m_xyint/F");
1054 
1055  m_tree->Branch("cluskeys", &m_cluskeys);
1056  m_tree->Branch("clusedge", &m_clusedge);
1057  m_tree->Branch("clusoverlap", &m_clusoverlap);
1058  m_tree->Branch("cluslx", &m_cluslx);
1059  m_tree->Branch("cluslz", &m_cluslz);
1060  m_tree->Branch("cluselx", &m_cluselx);
1061  m_tree->Branch("cluselz", &m_cluselz);
1062  m_tree->Branch("clusgx", &m_clusgx);
1063  m_tree->Branch("clusgy", &m_clusgy);
1064  m_tree->Branch("clusgz", &m_clusgz);
1065  m_tree->Branch("cluslayer", &m_cluslayer);
1066  m_tree->Branch("clussize", &m_clussize);
1067  m_tree->Branch("clushitsetkey", &m_clushitsetkey);
1068  m_tree->Branch("idealsurfcenterx", &m_idealsurfcenterx);
1069  m_tree->Branch("idealsurfcentery", &m_idealsurfcentery);
1070  m_tree->Branch("idealsurfcenterz", &m_idealsurfcenterz);
1071  m_tree->Branch("idealsurfnormx", &m_idealsurfnormx);
1072  m_tree->Branch("idealsurfnormy", &m_idealsurfnormy);
1073  m_tree->Branch("idealsurfnormz", &m_idealsurfnormz);
1074  m_tree->Branch("missurfcenterx", &m_missurfcenterx);
1075  m_tree->Branch("missurfcentery", &m_missurfcentery);
1076  m_tree->Branch("missurfcenterz", &m_missurfcenterz);
1077  m_tree->Branch("missurfnormx", &m_missurfnormx);
1078  m_tree->Branch("missurfnormy", &m_missurfnormy);
1079  m_tree->Branch("missurfnormz", &m_missurfnormz);
1080  m_tree->Branch("clusgxideal", &m_clusgxideal);
1081  m_tree->Branch("clusgyideal", &m_clusgyideal);
1082  m_tree->Branch("clusgzideal", &m_clusgzideal);
1083  m_tree->Branch("missurfalpha", &m_missurfalpha);
1084  m_tree->Branch("missurfbeta", &m_missurfbeta);
1085  m_tree->Branch("missurfgamma", &m_missurfgamma);
1086  m_tree->Branch("idealsurfalpha", &m_idealsurfalpha);
1087  m_tree->Branch("idealsurfbeta", &m_idealsurfbeta);
1088  m_tree->Branch("idealsurfgamma", &m_idealsurfgamma);
1089 
1090  m_tree->Branch("statelx", &m_statelx);
1091  m_tree->Branch("statelz", &m_statelz);
1092  m_tree->Branch("stateelx", &m_stateelx);
1093  m_tree->Branch("stateelz", &m_stateelz);
1094  m_tree->Branch("stategx", &m_stategx);
1095  m_tree->Branch("stategy", &m_stategy);
1096  m_tree->Branch("stategz", &m_stategz);
1097  m_tree->Branch("statepx", &m_statepx);
1098  m_tree->Branch("statepy", &m_statepy);
1099  m_tree->Branch("statepz", &m_statepz);
1100  m_tree->Branch("statepl", &m_statepl);
1101 
1102  m_tree->Branch("statelxglobderivdx", &m_statelxglobderivdx);
1103  m_tree->Branch("statelxglobderivdy", &m_statelxglobderivdy);
1104  m_tree->Branch("statelxglobderivdz", &m_statelxglobderivdz);
1105  m_tree->Branch("statelxglobderivdalpha", &m_statelxglobderivdalpha);
1106  m_tree->Branch("statelxglobderivdbeta", &m_statelxglobderivdbeta);
1107  m_tree->Branch("statelxglobderivdgamma", &m_statelxglobderivdgamma);
1108 
1109  m_tree->Branch("statelxlocderivd0", &m_statelxlocderivd0);
1110  m_tree->Branch("statelxlocderivz0", &m_statelxlocderivz0);
1111  m_tree->Branch("statelxlocderivphi", &m_statelxlocderivphi);
1112  m_tree->Branch("statelxlocderivtheta", &m_statelxlocderivtheta);
1113  m_tree->Branch("statelxlocderivqop", &m_statelxlocderivqop);
1114 
1115  m_tree->Branch("statelzglobderivdx", &m_statelzglobderivdx);
1116  m_tree->Branch("statelzglobderivdy", &m_statelzglobderivdy);
1117  m_tree->Branch("statelzglobderivdz", &m_statelzglobderivdz);
1118  m_tree->Branch("statelzglobderivdalpha", &m_statelzglobderivdalpha);
1119  m_tree->Branch("statelzglobderivdbeta", &m_statelzglobderivdbeta);
1120  m_tree->Branch("statelzglobderivdgamma", &m_statelzglobderivdgamma);
1121 
1122  m_tree->Branch("statelzlocderivd0", &m_statelzlocderivd0);
1123  m_tree->Branch("statelzlocderivz0", &m_statelzlocderivz0);
1124  m_tree->Branch("statelzlocderivphi", &m_statelzlocderivphi);
1125  m_tree->Branch("statelzlocderivtheta", &m_statelzlocderivtheta);
1126  m_tree->Branch("statelzlocderivqop", &m_statelzlocderivqop);
1127 }