Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HelicalFitter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HelicalFitter.cc
1 #include "HelicalFitter.h"
2 
3 #include "Mille.h"
4 
6 #include <trackbase/TrkrDefs.h> // for cluskey, getTrkrId, tpcId
7 #include <trackbase/TpcDefs.h>
8 #include <trackbase/MvtxDefs.h>
9 #include <trackbase/InttDefs.h>
10 #include <trackbase/TrkrClusterv3.h>
11 #include <trackbase/TrkrClusterContainerv4.h>
12 #include <trackbase/TrkrClusterContainer.h>
13 #include <trackbase/TrkrClusterCrossingAssoc.h>
14 #include <trackbase/alignmentTransformationContainer.h>
15 
16 #include <trackbase_historic/TrackSeed_v1.h>
17 #include <trackbase_historic/TrackSeedContainer_v1.h>
18 #include <trackbase_historic/SvtxTrackSeed_v1.h>
19 #include <trackbase_historic/SvtxTrack_v4.h>
20 #include <trackbase_historic/SvtxTrackMap_v2.h>
21 #include <trackbase_historic/SvtxAlignmentState_v1.h>
22 #include <trackbase_historic/SvtxAlignmentStateMap_v1.h>
23 #include <trackbase_historic/SvtxTrackState_v1.h>
24 
25 #include <globalvertex/SvtxVertex.h>
26 #include <globalvertex/SvtxVertexMap.h>
27 
28 #include <Acts/Surfaces/PerigeeSurface.hpp>
29 
30 #include <g4main/PHG4Hit.h> // for PHG4Hit
31 #include <g4main/PHG4Particle.h> // for PHG4Particle
32 #include <g4main/PHG4HitDefs.h> // for keytype
33 
34 #include <fun4all/Fun4AllReturnCodes.h>
35 
36 #include <phool/phool.h>
37 #include <phool/PHCompositeNode.h>
38 #include <phool/getClass.h>
39 
40 #include <TF1.h>
41 #include <TNtuple.h>
42 #include <TFile.h>
43 
44 #include <climits> // for UINT_MAX
45 #include <iostream> // for operator<<, basic_ostream
46 #include <cmath> // for fabs, sqrt
47 #include <set> // for _Rb_tree_const_iterator
48 #include <utility> // for pair
49 #include <memory>
50 
51 using namespace std;
52 
53 //____________________________________________________________________________..
54 HelicalFitter::HelicalFitter(const std::string &name):
55  SubsysReco(name)
56  , PHParameterInterface(name)
57  , _mille(nullptr)
58 {
59  InitializeParameters();
60 
61  vertexPosition(0) = 0;
62  vertexPosition(1) = 0;
63 
64  vtx_sigma(0) = 0.01;
65  vtx_sigma(1) = 0.01;
66 }
67 
68 //____________________________________________________________________________..
69 HelicalFitter::~HelicalFitter()
70 {
71 
72 }
73 
74 //____________________________________________________________________________..
75 int HelicalFitter::InitRun(PHCompositeNode *topNode)
76 {
77  UpdateParametersWithMacro();
78 
79  int ret = GetNodes(topNode);
80  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
81 
82  ret = CreateNodes(topNode);
83  if(ret != Fun4AllReturnCodes::EVENT_OK) return ret;
84 
85  // Instantiate Mille and open output data file
86  if(test_output)
87  {
88  _mille = new Mille(data_outfilename.c_str(), false); // write text in data files, rather than binary, for debugging only
89  }
90  else
91  {
92  _mille = new Mille(data_outfilename.c_str());
93  }
94 
95  // Write the steering file here, and add the data file path to it
96  std::ofstream steering_file(steering_outfilename);
97  steering_file << data_outfilename << std::endl;
98  steering_file.close();
99 
100  if(make_ntuple)
101  {
102  //fout = new TFile("HF_ntuple.root","recreate");
103  fout = new TFile(ntuple_outfilename.c_str(),"recreate");
104  ntp = new TNtuple("ntp","HF ntuple","event:trkid:layer:nsilicon:ntpc:nclus:trkrid:sector:side:subsurf:phi:glbl0:glbl1:glbl2:glbl3:glbl4:glbl5:sensx:sensy:sensz:normx:normy:normz:sensxideal:sensyideal:senszideal:normxideal:normyideal:normzideal:xglobideal:yglobideal:zglobideal:R:X0:Y0:Zs:Z0:xglob:yglob:zglob:xfit:yfit:zfit:pcax:pcay:pcaz:tangx:tangy:tangz:X:Y:fitX:fitY:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdalpha:dXdbeta:dXdgamma:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdalpha:dYdbeta:dYdgamma:dYdx:dYdy:dYdz");
105 
106  track_ntp = new TNtuple("track_ntp","HF track ntuple","track_id:residual_x:residual_y:residualxsigma:residualysigma:dXdR:dXdX0:dXdY0:dXdZs:dXdZ0:dXdx:dXdy:dXdz:dYdR:dYdX0:dYdY0:dYdZs:dYdZ0:dYdx:dYdy:dYdz:xvtx:yvtx:zvtx:event_zvtx:track_phi:perigee_phi");
107 
108  }
109 
110 
111  // print grouping setup to log file:
112  std::cout << "HelicalFitter::InitRun: Surface groupings are mvtx " << mvtx_grp << " intt " << intt_grp << " tpc " << tpc_grp << " mms " << mms_grp << std::endl;
113  std::cout << " possible groupings are:" << std::endl
114  << " mvtx "
115  << AlignmentDefs::mvtxGrp::snsr << " "
116  << AlignmentDefs::mvtxGrp::stv << " "
117  << AlignmentDefs::mvtxGrp::mvtxlyr << " "
118  << AlignmentDefs::mvtxGrp::clamshl << " " << std::endl
119  << " intt "
120  << AlignmentDefs::inttGrp::chp << " "
121  << AlignmentDefs::inttGrp::lad << " "
122  << AlignmentDefs::inttGrp::inttlyr << " "
123  << AlignmentDefs::inttGrp::inttbrl << " " << std::endl
124  << " tpc "
125  << AlignmentDefs::tpcGrp::htst << " "
126  << AlignmentDefs::tpcGrp::sctr << " "
127  << AlignmentDefs::tpcGrp::tp << " " << std::endl
128  << " mms "
129  << AlignmentDefs::mmsGrp::tl << " "
130  << AlignmentDefs::mmsGrp::mm << " " << std::endl;
131 
132  event=-1;
133 
134  return ret;
135 }
136 
137 //_____________________________________________________________________
138 void HelicalFitter::SetDefaultParameters()
139 {
140 
141 
142  return;
143 }
144 
145 //____________________________________________________________________________..
146 int HelicalFitter::process_event(PHCompositeNode*)
147 {
148  // _track_map_tpc contains the TPC seed track stubs
149  // _track_map_silicon contains the silicon seed track stubs
150  // _svtx_seed_map contains the combined silicon and tpc track seeds
151 
152  event++;
153 
154  if(Verbosity() > 0)
155  cout << PHWHERE
156  << " TPC seed map size " << _track_map_tpc->size()
157  << " Silicon seed map size " << _track_map_silicon->size()
158  << endl;
159 
160  if(_track_map_silicon->size() == 0 && _track_map_tpc->size() == 0)
161  return Fun4AllReturnCodes::EVENT_OK;
162 
163  // Decide whether we want to make a helical fit for silicon or TPC
164  unsigned int maxtracks = 0;
165  unsigned int nsilicon = 0;
166  unsigned int ntpc = 0;
167  unsigned int nclus = 0;
168  std::vector<std::vector<Acts::Vector3>> cumulative_global_vec;
169  std::vector<std::vector<TrkrDefs::cluskey>> cumulative_cluskey_vec;
170  std::vector<std::vector<float>> cumulative_fitpars_vec;
171  std::vector<Acts::Vector3> cumulative_vertex;
172  std::vector<TrackSeed> cumulative_someseed;
173  std::vector<SvtxTrack_v4> cumulative_newTrack;
174 
175 
176  if(fittpc) { maxtracks = _track_map_tpc->size(); }
177  if(fitsilicon) { maxtracks = _track_map_silicon->size(); }
178  for(unsigned int trackid = 0; trackid < maxtracks; ++trackid)
179  {
180  TrackSeed* tracklet = nullptr;
181  if(fitsilicon) { tracklet = _track_map_silicon->get(trackid); }
182  else if(fittpc) { tracklet = _track_map_tpc->get(trackid); }
183  if(!tracklet) { continue; }
184 
185  std::vector<Acts::Vector3> global_vec;
186  std::vector<TrkrDefs::cluskey> cluskey_vec;
187 
188  // Get a vector of cluster keys from the tracklet
189  getTrackletClusterList(tracklet, cluskey_vec);
190  // store cluster global positions in a vector global_vec and cluskey_vec
191  TrackFitUtils::getTrackletClusters(_tGeometry, _cluster_map, global_vec, cluskey_vec);
192 
193  correctTpcGlobalPositions( global_vec, cluskey_vec);
194 
195  std::vector<float> fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec); // do helical fit
196 
197  if(fitpars.size() == 0) continue; // discard this track, not enough clusters to fit
198 
199  if(Verbosity() > 1)
200  { std::cout << " Track " << trackid << " radius " << fitpars[0] << " X0 " << fitpars[1]<< " Y0 " << fitpars[2]
201  << " zslope " << fitpars[3] << " Z0 " << fitpars[4] << std::endl; }
202 
204  SvtxTrack_v4 newTrack;
205  newTrack.set_id(trackid);
206  if(fitsilicon) { newTrack.set_silicon_seed(tracklet); }
207  else if(fittpc) { newTrack.set_tpc_seed(tracklet); }
208 
209  // if a full track is requested, get the silicon clusters too and refit
210  if(fittpc && fitfulltrack)
211  {
212  // this associates silicon clusters and adds them to the vectors
213  ntpc = cluskey_vec.size();
214  nsilicon = TrackFitUtils::addClusters(fitpars, dca_cut, _tGeometry, _cluster_map, global_vec, cluskey_vec,0,6);
215  if(nsilicon < 3) continue; // discard this TPC seed, did not get a good match to silicon
216  auto trackseed = std::make_unique<TrackSeed_v1>();
217  for(auto& ckey : cluskey_vec)
218  {
219  if(TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::mvtxId or
220  TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::inttId)
221  {
222  trackseed->insert_cluster_key(ckey);
223  }
224  }
225 
226  newTrack.set_silicon_seed(trackseed.get());
227 
228  // fit the full track now
229  fitpars.clear();
230  fitpars = TrackFitUtils::fitClusters(global_vec, cluskey_vec); // do helical fit
231  if(fitpars.size() == 0) continue; // discard this track, fit failed
232 
233  if(Verbosity() > 1)
234  { std::cout << " Full track " << trackid << " radius " << fitpars[0] << " X0 " << fitpars[1]<< " Y0 " << fitpars[2]
235  << " zslope " << fitpars[3] << " Z0 " << fitpars[4] << std::endl; }
236  }
237  else if(fitsilicon)
238  {
239  nsilicon = cluskey_vec.size();
240  }
241  else if(fittpc && !fitfulltrack)
242  {
243  ntpc = cluskey_vec.size();
244  }
245 
246  Acts::Vector3 beamline(0,0,0);
247  Acts::Vector2 pca2d = TrackFitUtils::get_circle_point_pca(fitpars[0], fitpars[1], fitpars[2], beamline);
248  Acts::Vector3 track_vtx (pca2d(0),pca2d(1),fitpars[4]);
249 
250  newTrack.set_crossing(tracklet->get_crossing());
251  newTrack.set_id(trackid);
252 
255 
256  TrackSeed_v1 someseed;
257  for(auto& ckey : cluskey_vec)
258  { someseed.insert_cluster_key(ckey); }
259  someseed.set_qOverR(tracklet->get_charge() / fitpars[0]);
260 
261  someseed.set_X0(fitpars[1]);
262  someseed.set_Y0(fitpars[2]);
263  someseed.set_Z0(fitpars[4]);
264  someseed.set_slope(fitpars[3]);
265 
266  newTrack.set_x(someseed.get_x());
267  newTrack.set_y(someseed.get_y());
268  newTrack.set_z(someseed.get_z());
269  newTrack.set_px(someseed.get_px(_cluster_map,_tGeometry));
270  newTrack.set_py(someseed.get_py(_cluster_map,_tGeometry));
271  newTrack.set_pz(someseed.get_pz());
272  newTrack.set_charge(tracklet->get_charge());
273 
274  nclus = ntpc+nsilicon;
275 
276  // some basic track quality requirements
277  if(fittpc && ntpc < 35)
278  {
279  if(Verbosity() > 1) { std::cout << " reject this track, ntpc = " << ntpc << std::endl; }
280  continue;
281  }
282  if((fitsilicon || fitfulltrack) && nsilicon < 4)
283  {
284  if(Verbosity() > 1) { std::cout << " reject this track, nsilicon = " << nsilicon << std::endl; }
285  continue;
286  }
287 
288  cumulative_global_vec.push_back(global_vec);
289  cumulative_cluskey_vec.push_back(cluskey_vec);
290  cumulative_vertex.push_back(track_vtx);
291  cumulative_fitpars_vec.push_back(fitpars);
292  cumulative_someseed.push_back(someseed);
293  cumulative_newTrack.push_back(newTrack);
294  }
295 
296  //terminate loop over tracks
297  //Collect fitpars for each track by intializing array of size maxtracks and populaating thorughout the loop
298  //Then start new loop over tracks and for each track go over clsutaer
299  // make vector of global_vecs
300  float xsum = 0;
301  float ysum = 0;
302  float zsum = 0;
303  unsigned int accepted_tracks = cumulative_fitpars_vec.size();
304 
305  //std::cout<<"accepted_tracks: " << accepted_tracks << std::endl;
306  for(unsigned int trackid = 0; trackid < accepted_tracks; ++trackid)
307  {
308  xsum += cumulative_vertex[trackid][0];
309  ysum += cumulative_vertex[trackid][1];
310  zsum += cumulative_vertex[trackid][2];
311  }
312  Acts::Vector3 averageVertex (xsum/accepted_tracks,ysum/accepted_tracks,zsum/accepted_tracks);
313 
314  for(unsigned int trackid = 0; trackid < accepted_tracks; ++trackid)
315  {
316  auto global_vec = cumulative_global_vec[trackid];
317  auto cluskey_vec = cumulative_cluskey_vec[trackid];
318  auto fitpars = cumulative_fitpars_vec[trackid];
319  auto someseed = cumulative_someseed[trackid];
320  auto newTrack = cumulative_newTrack[trackid];
321  //std::cout << "trackid " << trackid << " get id " <<newTrack.get_id()<< std::endl;
322  SvtxAlignmentStateMap::StateVec statevec;
323 
324  // get the residuals and derivatives for all clusters
325  for(unsigned int ivec=0;ivec<global_vec.size(); ++ivec)
326  {
327  auto global = global_vec[ivec];
328  auto cluskey = cluskey_vec[ivec];
329  auto cluster = _cluster_map->findCluster(cluskey);
330  if(!cluster) { continue;}
331 
332  unsigned int trkrid = TrkrDefs::getTrkrId(cluskey);
333 
334  // What we need now is to find the point on the surface at which the helix would intersect
335  // If we have that point, we can transform the fit back to local coords
336  // we have fitpars for the helix, and the cluster key - from which we get the surface
337 
338  Surface surf = _tGeometry->maps().getSurface(cluskey, cluster);
339  Acts::Vector3 helix_pca(0,0,0);
340  Acts::Vector3 helix_tangent(0,0,0);
341  Acts::Vector3 fitpoint = get_helix_surface_intersection(surf, fitpars, global, helix_pca, helix_tangent);
342 
343  // fitpoint is the point where the helical fit intersects the plane of the surface
344  // Now transform the helix fitpoint to local coordinates to compare with cluster local coordinates
345  Acts::Vector3 fitpoint_local = surf->transform(_tGeometry->geometry().getGeoContext()).inverse() * (fitpoint * Acts::UnitConstants::cm);
346 
347  fitpoint_local /= Acts::UnitConstants::cm;
348 
349  auto xloc = cluster->getLocalX(); // in cm
350  auto zloc = cluster->getLocalY();
351 
352  if(trkrid == TrkrDefs::tpcId) { zloc = convertTimeToZ(cluskey, cluster); }
353  //float yloc = 0.0; // Because the fitpoint is on the surface, y will always be zero in local coordinates
354 
355  Acts::Vector2 residual(xloc - fitpoint_local(0), zloc - fitpoint_local(1));
356 
357  unsigned int layer = TrkrDefs::getLayer(cluskey_vec[ivec]);
358  float phi = atan2(global(1), global(0));
359 
360  SvtxTrackState_v1 svtxstate(fitpoint.norm());
361  svtxstate.set_x(fitpoint(0));
362  svtxstate.set_y(fitpoint(1));
363  svtxstate.set_z(fitpoint(2));
364  auto tangent = get_helix_tangent(fitpars, global);
365  svtxstate.set_px(someseed.get_p() * tangent.second.x());
366  svtxstate.set_py(someseed.get_p() * tangent.second.y());
367  svtxstate.set_pz(someseed.get_p() * tangent.second.z());
368  newTrack.insert_state(&svtxstate);
369 
370  if(Verbosity() > 1) {
371  Acts::Vector3 loc_check = surf->transform(_tGeometry->geometry().getGeoContext()).inverse() * (global * Acts::UnitConstants::cm);
372  loc_check /= Acts::UnitConstants::cm;
373  std::cout << " layer " << layer << std::endl
374  << " cluster global " << global(0) << " " << global(1) << " " << global(2) << std::endl
375  << " fitpoint " << fitpoint(0) << " " << fitpoint(1) << " " << fitpoint(2) << std::endl
376  << " fitpoint_local " << fitpoint_local(0) << " " << fitpoint_local(1) << " " << fitpoint_local(2) << std::endl
377  << " cluster local x " << cluster->getLocalX() << " cluster local y " << cluster->getLocalY() << std::endl
378  << " cluster global to local x " << loc_check(0) << " local y " << loc_check(1) << " local z " << loc_check(2) << std::endl
379  << " cluster local residual x " << residual(0) << " cluster local residual y " <<residual(1) << std::endl;
380  }
381 
382  if(Verbosity() > 1)
383  {
384  Acts::Transform3 transform = surf->transform(_tGeometry->geometry().getGeoContext());
385  std::cout << "Transform is:" << std::endl;
386  std::cout << transform.matrix() << std::endl;
387  Acts::Vector3 loc_check = surf->transform(_tGeometry->geometry().getGeoContext()).inverse() * (global * Acts::UnitConstants::cm);
388  loc_check /= Acts::UnitConstants::cm;
389  unsigned int sector = TpcDefs::getSectorId(cluskey_vec[ivec]);
390  unsigned int side = TpcDefs::getSide(cluskey_vec[ivec]);
391  std::cout << " layer " << layer << " sector " << sector << " side " << side << " subsurf " << cluster->getSubSurfKey() << std::endl
392  << " cluster global " << global(0) << " " << global(1) << " " << global(2) << std::endl
393  << " fitpoint " << fitpoint(0) << " " << fitpoint(1) << " " << fitpoint(2) << std::endl
394  << " fitpoint_local " << fitpoint_local(0) << " " << fitpoint_local(1) << " " << fitpoint_local(2) << std::endl
395  << " cluster local x " << cluster->getLocalX() << " cluster local y " << cluster->getLocalY() << std::endl
396  << " cluster global to local x " << loc_check(0) << " local y " << loc_check(1) << " local z " << loc_check(2) << std::endl
397  << " cluster local residual x " << residual(0) << " cluster local residual y " <<residual(1) << std::endl;
398  }
399 
400  // need standard deviation of measurements
401  Acts::Vector2 clus_sigma = getClusterError(cluster, cluskey, global);
402  if(isnan(clus_sigma(0)) || isnan(clus_sigma(1))) { continue; }
403 
404  int glbl_label[AlignmentDefs::NGL];
405  if(layer < 3)
406  {
407  AlignmentDefs::getMvtxGlobalLabels(surf, glbl_label, mvtx_grp);
408  }
409  else if(layer > 2 && layer < 7)
410  {
411  AlignmentDefs::getInttGlobalLabels(surf, glbl_label, intt_grp);
412  }
413  else if (layer < 55)
414  {
415  AlignmentDefs::getTpcGlobalLabels(surf, cluskey, glbl_label, tpc_grp);
416  }
417  else
418  {
419  continue;
420  }
421 
422  // These derivatives are for the local parameters
423  float lcl_derivativeX[AlignmentDefs::NLC];
424  float lcl_derivativeY[AlignmentDefs::NLC];
425 
426  getLocalDerivativesXY(surf, global, fitpars, lcl_derivativeX, lcl_derivativeY, layer);
427 
428  // The global derivs dimensions are [alpha/beta/gamma](x/y/z)
429  float glbl_derivativeX[AlignmentDefs::NGL];
430  float glbl_derivativeY[AlignmentDefs::NGL];
431  getGlobalDerivativesXY(surf, global, fitpoint, fitpars, glbl_derivativeX, glbl_derivativeY, layer);
432 
433  auto alignmentstate = std::make_unique<SvtxAlignmentState_v1>();
434  alignmentstate->set_residual(residual);
435  alignmentstate->set_cluster_key(cluskey);
436  SvtxAlignmentState::GlobalMatrix svtxglob =
437  SvtxAlignmentState::GlobalMatrix::Zero();
438  SvtxAlignmentState::LocalMatrix svtxloc =
439  SvtxAlignmentState::LocalMatrix::Zero();
440  for(int i=0; i<AlignmentDefs::NLC; i++)
441  {
442  svtxloc(0,i) = lcl_derivativeX[i];
443  svtxloc(1,i) = lcl_derivativeY[i];
444  }
445  for(int i=0; i<AlignmentDefs::NGL; i++)
446  {
447  svtxglob(0,i) = glbl_derivativeX[i];
448  svtxglob(1,i) = glbl_derivativeY[i];
449  }
450 
451  alignmentstate->set_local_derivative_matrix(svtxloc);
452  alignmentstate->set_global_derivative_matrix(svtxglob);
453 
454  statevec.push_back(alignmentstate.release());
455 
456  for(unsigned int i = 0; i < AlignmentDefs::NGL; ++i)
457  {
458  if(trkrid == TrkrDefs::mvtxId)
459  {
460  // need stave to get clamshell
461  auto stave = MvtxDefs::getStaveId(cluskey_vec[ivec]);
462  auto clamshell = AlignmentDefs::getMvtxClamshell(layer, stave);
463  if( is_layer_param_fixed(layer, i) || is_mvtx_layer_fixed(layer,clamshell) )
464  {
465  glbl_derivativeX[i] = 0;
466  glbl_derivativeY[i] = 0;
467  }
468  }
469 
470  if(trkrid == TrkrDefs::inttId)
471  {
472  if( is_layer_param_fixed(layer, i) || is_intt_layer_fixed(layer) )
473  {
474  glbl_derivativeX[i] = 0;
475  glbl_derivativeY[i] = 0;
476  }
477  }
478 
479 
480  if(trkrid == TrkrDefs::tpcId)
481  {
482  unsigned int sector = TpcDefs::getSectorId(cluskey_vec[ivec]);
483  unsigned int side = TpcDefs::getSide(cluskey_vec[ivec]);
484  if(is_layer_param_fixed(layer, i) || is_tpc_sector_fixed(layer, sector, side))
485  {
486  glbl_derivativeX[i] = 0;
487  glbl_derivativeY[i] = 0;
488  }
489  }
490  }
491 
492  // Add the measurement separately for each coordinate direction to Mille
493  // set the derivatives non-zero only for parameters we want to be optimized
494  // local parameter numbering is arbitrary:
495  float errinf = 1.0;
496 
497  if(_layerMisalignment.find(layer) != _layerMisalignment.end())
498  {
499  errinf = _layerMisalignment.find(layer)->second;
500  }
501  if(make_ntuple)
502  {
503  // get the local parameters using the ideal transforms
504  alignmentTransformationContainer::use_alignment = false;
505  Acts::Vector3 ideal_center = surf->center(_tGeometry->geometry().getGeoContext()) * 0.1;
506  Acts::Vector3 ideal_norm = -surf->normal(_tGeometry->geometry().getGeoContext());
507  Acts::Vector3 ideal_local(xloc, zloc, 0.0); // cm
508  Acts::Vector3 ideal_glob = surf->transform(_tGeometry->geometry().getGeoContext())*(ideal_local * Acts::UnitConstants::cm);
509  ideal_glob /= Acts::UnitConstants::cm;
510  alignmentTransformationContainer::use_alignment = true;
511 
512  Acts::Vector3 sensorCenter = surf->center(_tGeometry->geometry().getGeoContext()) * 0.1; // cm
513  Acts::Vector3 sensorNormal = -surf->normal(_tGeometry->geometry().getGeoContext());
514  unsigned int sector = TpcDefs::getSectorId(cluskey_vec[ivec]);
515  unsigned int side = TpcDefs::getSide(cluskey_vec[ivec]);
516  unsigned int subsurf = cluster->getSubSurfKey();
517  if(layer < 3)
518  {
519  sector = MvtxDefs::getStaveId(cluskey_vec[ivec]);
520  subsurf = MvtxDefs::getChipId(cluskey_vec[ivec]);
521  }
522  else if(layer >2 && layer < 7)
523  {
524  sector = InttDefs::getLadderPhiId(cluskey_vec[ivec]);
525  subsurf = InttDefs::getLadderZId(cluskey_vec[ivec]);
526  }
527  float ntp_data[75] = {
528  (float) event, (float) trackid,
529  (float) layer, (float) nsilicon, (float) ntpc, (float) nclus, (float) trkrid, (float) sector, (float) side,
530  (float) subsurf, phi,
531  (float) glbl_label[0], (float) glbl_label[1], (float) glbl_label[2], (float) glbl_label[3], (float) glbl_label[4], (float) glbl_label[5],
532  (float) sensorCenter(0), (float) sensorCenter(1), (float) sensorCenter(2),
533  (float) sensorNormal(0), (float) sensorNormal(1), (float) sensorNormal(2),
534  (float) ideal_center(0), (float) ideal_center(1), (float) ideal_center(2),
535  (float) ideal_norm(0), (float) ideal_norm(1), (float) ideal_norm(2),
536  (float) ideal_glob(0), (float) ideal_glob(1), (float) ideal_glob(2),
537  (float) fitpars[0], (float) fitpars[1], (float) fitpars[2], (float) fitpars[3], (float) fitpars[4],
538  (float) global(0), (float) global(1), (float) global(2),
539  (float) fitpoint(0), (float) fitpoint(1), (float) fitpoint(2),
540  (float) helix_pca(0), (float) helix_pca(1), (float) helix_pca(2),
541  (float) helix_tangent(0), (float) helix_tangent(1), (float) helix_tangent(2),
542  xloc,zloc, (float) fitpoint_local(0), (float) fitpoint_local(1),
543  lcl_derivativeX[0],lcl_derivativeX[1],lcl_derivativeX[2],lcl_derivativeX[3],lcl_derivativeX[4],
544  glbl_derivativeX[0],glbl_derivativeX[1],glbl_derivativeX[2],glbl_derivativeX[3],glbl_derivativeX[4],glbl_derivativeX[5],
545  lcl_derivativeY[0],lcl_derivativeY[1],lcl_derivativeY[2],lcl_derivativeY[3],lcl_derivativeY[4],
546  glbl_derivativeY[0],glbl_derivativeY[1],glbl_derivativeY[2],glbl_derivativeY[3],glbl_derivativeY[4],glbl_derivativeY[5] };
547 
548  ntp->Fill(ntp_data);
549 
550  if(Verbosity() > 2)
551  {
552  for(int i=0;i<34;++i)
553  {
554  std::cout << ntp_data[i] << " " ;
555  }
556  std::cout << std::endl;
557  }
558  }
559 
560  // add some cluster cuts
561  if(residual(0) > 0.2) continue; // 2 mm cut
562  if(residual(1) > 0.2) continue; // 2 mm cut
563 
564  if( !isnan(residual(0)) && clus_sigma(0) < 1.0) // discards crazy clusters
565  {
566  _mille->mille(AlignmentDefs::NLC,lcl_derivativeX,AlignmentDefs::NGL,glbl_derivativeX,glbl_label,residual(0), errinf*clus_sigma(0));
567  }
568  if(!isnan(residual(1)) && clus_sigma(1) < 1.0 && trkrid != TrkrDefs::inttId)
569  {
570  _mille->mille(AlignmentDefs::NLC, lcl_derivativeY,AlignmentDefs::NGL,glbl_derivativeY,glbl_label,residual(1), errinf*clus_sigma(1));
571  }
572  }
573 
574  m_alignmentmap->insertWithKey(trackid, statevec);
575  m_trackmap->insertWithKey(&newTrack, trackid);
576 
577  //calculate vertex residual with perigee surface
578  Acts::Vector3 event_vtx(0,0,averageVertex(2));
579 
580  // The residual for the vtx case is (event vtx - track vtx)
581  // that is -dca
582  float dca3dxy;
583  float dca3dz;
584  float dca3dxysigma;
585  float dca3dzsigma;
586  get_dca(newTrack,dca3dxy,dca3dz,dca3dxysigma,dca3dzsigma,event_vtx);
587  Acts::Vector2 vtx_residual(-dca3dxy, -dca3dz);
588 
589  float lclvtx_derivativeX[AlignmentDefs::NLC];
590  float lclvtx_derivativeY[AlignmentDefs::NLC];
591  getLocalVtxDerivativesXY(newTrack, event_vtx, fitpars, lclvtx_derivativeX, lclvtx_derivativeY);
592 
593  // The global derivs dimensions are [alpha/beta/gamma](x/y/z)
594  float glblvtx_derivativeX[3];
595  float glblvtx_derivativeY[3];
596  getGlobalVtxDerivativesXY(newTrack, event_vtx, glblvtx_derivativeX, glblvtx_derivativeY);
597 
598  if(use_event_vertex)
599  {
600  if(Verbosity() > 3)
601  {
602  std::cout << "vertex info for track " << trackid << " with charge " << newTrack.get_charge() << std::endl;
603 
604  std::cout << "vertex is " << event_vtx.transpose() << std::endl;
605  std::cout << "vertex residuals " << vtx_residual.transpose()
606  << std::endl;
607  std::cout << "local derivatives " << std::endl;
608  for(int i=0; i<AlignmentDefs::NLC; i++)
609  std::cout << lclvtx_derivativeX[i] << ", ";
610  std::cout << std::endl;
611  for(int i=0; i<AlignmentDefs::NLC; i++)
612  std::cout << lclvtx_derivativeY[i] << ", ";
613  std::cout << "global vtx derivaties " << std::endl;
614  for(int i=0; i<3; i++) std::cout << glblvtx_derivativeX[i] << ", ";
615  std::cout << std::endl;
616  for(int i=0; i<3; i++) std::cout << glblvtx_derivativeY[i] << ", ";
617  }
618 
619  // add some track cuts
620  if(fabs(newTrack.get_z() - event_vtx(2)) > 0.2) continue; // 2 mm cut
621  if(fabs(newTrack.get_x()) > 0.2) continue; // 2 mm cut
622  if(fabs(newTrack.get_y()) > 0.2) continue; // 2 mm cut
623 
624  if(!isnan(vtx_residual(0)))
625  {
626  _mille->mille(AlignmentDefs::NLC,lclvtx_derivativeX,AlignmentDefs::NGLVTX,glblvtx_derivativeX,AlignmentDefs::glbl_vtx_label,vtx_residual(0), vtx_sigma(0));
627  }
628  if(!isnan(vtx_residual(1)))
629  {
630  _mille->mille(AlignmentDefs::NLC,lclvtx_derivativeY,AlignmentDefs::NGLVTX,glblvtx_derivativeY,AlignmentDefs::glbl_vtx_label,vtx_residual(1), vtx_sigma(1));
631  }
632  }
633 
634  if(make_ntuple)
635  {
636  Acts::Vector3 mom(newTrack.get_px(),newTrack.get_py(),newTrack.get_pz());
637  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
638  float perigee_phi = atan2(r(1), r(0));
639  float track_phi = atan2(newTrack.get_py(), newTrack.get_px());
640  // float ntp_data[30] = {(float) trackid,dca3dxy,dca3dz,(float) vtx_sigma(0),(float) vtx_sigma(1),
641  float ntp_data[30] = {(float) trackid,(float) vtx_residual(0),(float) vtx_residual(1),(float) vtx_sigma(0),(float) vtx_sigma(1),
642  lclvtx_derivativeX[0],lclvtx_derivativeX[1],lclvtx_derivativeX[2],lclvtx_derivativeX[3],lclvtx_derivativeX[4],
643  glblvtx_derivativeX[0],glblvtx_derivativeX[1],glblvtx_derivativeX[2],
644  lclvtx_derivativeY[0],lclvtx_derivativeY[1],lclvtx_derivativeY[2],lclvtx_derivativeY[3],lclvtx_derivativeY[4],
645  glblvtx_derivativeY[0],glblvtx_derivativeY[1],glblvtx_derivativeY[2],
646  newTrack.get_x(), newTrack.get_y(), newTrack.get_z(),(float) event_vtx(2),track_phi, perigee_phi};
647 
648  track_ntp->Fill(ntp_data);
649  }
650 
651  if(Verbosity()>1)
652  {
653  std::cout << "vtx_residual xy: " << vtx_residual(0)<< " vtx_residual z: " << vtx_residual(1) << " vtx_sigma xy: " << vtx_sigma(0) << " vtx_sigma z: " << vtx_sigma(1) << std::endl;
654  std::cout << "track_x" << newTrack.get_x()<<"track_y" << newTrack.get_y()<<"track_z" << newTrack.get_z()<<std::endl;
655  }
656 
657 
658  // close out this track
659  _mille->end();
660 
661  } // end loop over tracks
662 
663 
664  return Fun4AllReturnCodes::EVENT_OK;
665 }
666 
667 Acts::Vector3 HelicalFitter::get_helix_surface_intersection(Surface surf, std::vector<float>& fitpars, Acts::Vector3 global)
668 {
669  // we want the point where the helix intersects the plane of the surface
670  // get the plane of the surface
671  Acts::Vector3 sensorCenter = surf->center(_tGeometry->geometry().getGeoContext()) * 0.1; // convert to cm
672  Acts::Vector3 sensorNormal = -surf->normal(_tGeometry->geometry().getGeoContext());
673  sensorNormal /= sensorNormal.norm();
674 
675  // there are analytic solutions for a line-plane intersection.
676  // to use this, need to get the vector tangent to the helix near the measurement and a point on it.
677  std::pair<Acts::Vector3, Acts::Vector3> line = get_helix_tangent(fitpars, global);
678  Acts::Vector3 pca = line.first;
679  Acts::Vector3 tangent = line.second;
680 
681  Acts::Vector3 intersection = get_line_plane_intersection(pca, tangent, sensorCenter, sensorNormal);
682 
683  return intersection;
684 }
685 
686 Acts::Vector3 HelicalFitter::get_helix_surface_intersection(Surface surf, std::vector<float>& fitpars, Acts::Vector3 global, Acts::Vector3& pca, Acts::Vector3& tangent)
687 {
688  // we want the point where the helix intersects the plane of the surface
689 
690  // get the plane of the surface
691  Acts::Vector3 sensorCenter = surf->center(_tGeometry->geometry().getGeoContext()) * 0.1; // convert to cm
692  Acts::Vector3 sensorNormal = -surf->normal(_tGeometry->geometry().getGeoContext());
693  sensorNormal /= sensorNormal.norm();
694 
695  // there are analytic solutions for a line-plane intersection.
696  // to use this, need to get the vector tangent to the helix near the measurement and a point on it.
697  std::pair<Acts::Vector3, Acts::Vector3> line = get_helix_tangent(fitpars, global);
698  pca = line.first;
699  tangent = line.second;
700 
701  Acts::Vector3 intersection = get_line_plane_intersection(pca, tangent, sensorCenter, sensorNormal);
702 
703  return intersection;
704 }
705 
706 
707 Acts::Vector3 HelicalFitter::get_helix_vtx(Acts::Vector3 event_vtx, const std::vector<float>& fitpars)
708 {
709  Acts::Vector2 pca2d = TrackFitUtils::get_circle_point_pca(fitpars[0], fitpars[1], fitpars[2], event_vtx);
710  Acts::Vector3 helix_vtx (pca2d(0),pca2d(1),fitpars[4]);
711 
712  return helix_vtx;
713 }
714 
715 
716 Acts::Vector3 HelicalFitter::get_line_plane_intersection(Acts::Vector3 PCA, Acts::Vector3 tangent, Acts::Vector3 sensor_center, Acts::Vector3 sensor_normal)
717 {
718  // get the intersection of the line made by PCA and tangent with the plane of the sensor
719 
720  // For a point on the line
721  // p = PCA + d * tangent;
722  // for a point on the plane
723  // (p - sensor_center).sensor_normal = 0
724 
725  // The solution is:
726  float d = (sensor_center - PCA).dot(sensor_normal) / tangent.dot(sensor_normal);
727  Acts::Vector3 intersection = PCA + d * tangent;
728 
729  return intersection;
730 }
731 
732 
733 std::pair<Acts::Vector3, Acts::Vector3> HelicalFitter::get_helix_tangent(const std::vector<float>& fitpars, Acts::Vector3 global)
734 {
735  auto pair = TrackFitUtils::get_helix_tangent(fitpars, global);
736  /*
737  save for posterity purposes
738  if(Verbosity() > 2)
739  {
740  // different method for checking:
741  // project the circle PCA vector an additional small amount and find the helix PCA to that point
742 
743  float projection = 0.25; // cm
744  Acts::Vector3 second_point = pca + projection * pca/pca.norm();
745  Acts::Vector2 second_point_pca_circle = TrackFitUtils::get_circle_point_pca(radius, x0, y0, second_point);
746  float second_point_pca_z = second_point_pca_circle.norm() * zslope + z0;
747  Acts::Vector3 second_point_pca2(second_point_pca_circle(0), second_point_pca_circle(1), second_point_pca_z);
748  Acts::Vector3 tangent2 = (second_point_pca2 - pca) / (second_point_pca2 - pca).norm();
749  Acts::Vector3 final_pca2 = getPCALinePoint(global, tangent2, pca);
750 
751  std::cout << " get_helix_tangent: getting tangent at angle_pca: " << angle_pca * 180.0 / M_PI << std::endl
752  << " original first point pca " << pca(0) << " " << pca(1) << " " << pca(2) << std::endl
753  << " original second point pca " << second_point_pca(0) << " " << second_point_pca(1) << " " << second_point_pca(2) << std::endl
754  << " original tangent " << tangent(0) << " " << tangent(1) << " " << tangent(2) << std::endl
755  << " original final pca from line " << final_pca(0) << " " << final_pca(1) << " " << final_pca(2) << std::endl;
756 
757  if(Verbosity() > 3)
758  {
759  std::cout << " Check: 2nd point pca meth 2 "<< second_point_pca2(0)<< " "<< second_point_pca2(1) << " "<< second_point_pca2(2) << std::endl
760  << " check tangent " << tangent2(0) << " " << tangent2(1) << " " << tangent2(2) << std::endl
761  << " check final pca from line " << final_pca2(0) << " " << final_pca2(1) << " " << final_pca2(2)
762  << std::endl;
763  }
764 
765  }
766  */
767 
768 
769  return pair;
770 }
771 
772 int HelicalFitter::End(PHCompositeNode* )
773 {
774  // closes output file in destructor
775  delete _mille;
776 
777  if(make_ntuple)
778  {
779  fout->Write();
780  fout->Close();
781  }
782 
783  return Fun4AllReturnCodes::EVENT_OK;
784 }
785 int HelicalFitter::CreateNodes(PHCompositeNode* topNode)
786 {
787  PHNodeIterator iter(topNode);
788 
789  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
790 
791  if (!dstNode)
792  {
793  std::cerr << "DST node is missing, quitting" << std::endl;
794  throw std::runtime_error("Failed to find DST node in PHActsTrkFitter::createNodes");
795  }
796 
797  PHNodeIterator dstIter(topNode);
798  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
799  if (!svtxNode)
800  {
801  svtxNode = new PHCompositeNode("SVTX");
802  dstNode->addNode(svtxNode);
803  }
804 
805  m_trackmap = findNode::getClass<SvtxTrackMap>(topNode, "HelicalFitterTrackMap");
806  if(!m_trackmap)
807  {
808  m_trackmap = new SvtxTrackMap_v2;
809  PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(m_trackmap,"HelicalFitterTrackMap","PHObject");
810  svtxNode->addNode(node);
811  }
812 
813  m_alignmentmap = findNode::getClass<SvtxAlignmentStateMap>(topNode,"HelicalFitterAlignmentStateMap");
814  if(!m_alignmentmap)
815  {
816  m_alignmentmap = new SvtxAlignmentStateMap_v1;
817  PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(m_alignmentmap,"HelicalFitterAlignmentStateMap","PHObject");
818  svtxNode->addNode(node);
819  }
820 
821  return Fun4AllReturnCodes::EVENT_OK;
822 }
823 int HelicalFitter::GetNodes(PHCompositeNode* topNode)
824 {
825  //---------------------------------
826  // Get additional objects off the Node Tree
827  //---------------------------------
828 
829  _track_map_silicon = findNode::getClass<TrackSeedContainer>(topNode, _silicon_track_map_name);
830  if (!_track_map_silicon)
831  {
832  cerr << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer " << endl;
833  return Fun4AllReturnCodes::ABORTEVENT;
834  }
835 
836  _track_map_tpc = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
837  if (!_track_map_tpc)
838  {
839  cerr << PHWHERE << " ERROR: Can't find " << _track_map_name.c_str() << endl;
840  return Fun4AllReturnCodes::ABORTEVENT;
841  }
842 
843  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
844  if (!_cluster_map)
845  {
846  std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
847  return Fun4AllReturnCodes::ABORTEVENT;
848  }
849 
850  _tGeometry = findNode::getClass<ActsGeometry>(topNode,"ActsGeometry");
851  if(!_tGeometry)
852  {
853  std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
854  return Fun4AllReturnCodes::ABORTEVENT;
855  }
856 
857  return Fun4AllReturnCodes::EVENT_OK;
858 }
859 
860 Acts::Vector3 HelicalFitter::get_helix_pca(std::vector<float>& fitpars, Acts::Vector3 global)
861 {
862  return TrackFitUtils::get_helix_pca(fitpars, global);
863 }
864 
865 
866 Acts::Vector3 HelicalFitter::getPCALinePoint(Acts::Vector3 global, Acts::Vector3 tangent, Acts::Vector3 posref)
867 {
868  // Approximate track with a straight line consisting of the state position posref and the vector (px,py,pz)
869 
870  // The position of the closest point on the line to global is:
871  // posref + projection of difference between the point and posref on the tangent vector
872  Acts::Vector3 pca = posref + ( (global - posref).dot(tangent) ) * tangent;
873 
874  return pca;
875 }
876 
877 
878 float HelicalFitter::convertTimeToZ(TrkrDefs::cluskey cluster_key, TrkrCluster *cluster)
879 {
880  // must convert local Y from cluster average time of arival to local cluster z position
881  double drift_velocity = _tGeometry->get_drift_velocity();
882  double zdriftlength = cluster->getLocalY() * drift_velocity;
883  double surfCenterZ = 52.89; // 52.89 is where G4 thinks the surface center is
884  double zloc = surfCenterZ - zdriftlength; // converts z drift length to local z position in the TPC in north
885  unsigned int side = TpcDefs::getSide(cluster_key);
886  if(side == 0) zloc = -zloc;
887  float z = zloc; // in cm
888 
889  return z;
890 }
891 
892 void HelicalFitter::makeTpcGlobalCorrections(TrkrDefs::cluskey cluster_key, short int crossing, Acts::Vector3& global)
893 {
894  // make all corrections to global position of TPC cluster
895  unsigned int side = TpcDefs::getSide(cluster_key);
896  float z = m_clusterCrossingCorrection.correctZ(global[2], side, crossing);
897  global[2] = z;
898 
899  // apply distortion corrections
900  if(_dcc_static) { global = _distortionCorrection.get_corrected_position( global, _dcc_static ); }
901  if(_dcc_average) { global = _distortionCorrection.get_corrected_position( global, _dcc_average ); }
902  if(_dcc_fluctuation) { global = _distortionCorrection.get_corrected_position( global, _dcc_fluctuation ); }
903 }
904 
905 void HelicalFitter::getTrackletClusters(TrackSeed *tracklet, std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey>& cluskey_vec)
906 {
907  getTrackletClusterList(tracklet, cluskey_vec);
908  // store cluster global positions in a vector
909  TrackFitUtils::getTrackletClusters(_tGeometry, _cluster_map, global_vec, cluskey_vec);
910 }
911 
912 void HelicalFitter::getTrackletClusterList(TrackSeed *tracklet, std::vector<TrkrDefs::cluskey>& cluskey_vec)
913 {
914  for (auto clusIter = tracklet->begin_cluster_keys();
915  clusIter != tracklet->end_cluster_keys();
916  ++clusIter)
917  {
918  auto key = *clusIter;
919  auto cluster = _cluster_map->findCluster(key);
920  if(!cluster)
921  {
922  std::cout << "Failed to get cluster with key " << key << std::endl;
923  continue;
924  }
925 
927  auto surf = _tGeometry->maps().getSurface(key, cluster);
928  if(!surf) { continue; }
929 
930  // drop some bad layers in the TPC completely
931  unsigned int layer = TrkrDefs::getLayer(key);
932  if(layer == 7 || layer == 22 || layer == 23 || layer == 38 || layer == 39) {continue;}
933 
934  cluskey_vec.push_back(key);
935 
936  } // end loop over clusters for this track
937 }
938 
939 std::vector<float> HelicalFitter::fitClusters(std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey> cluskey_vec)
940 {
941  return TrackFitUtils::fitClusters(global_vec, cluskey_vec); // do helical fit
942 }
943 
944 Acts::Vector2 HelicalFitter::getClusterError(TrkrCluster *cluster, TrkrDefs::cluskey cluskey, Acts::Vector3& global)
945 {
946  Acts::Vector2 clus_sigma(0,0);
947 
948  double clusRadius = sqrt(global[0]*global[0] + global[1]*global[1]);
949  auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster,clusRadius,cluskey);
950  double phierror = sqrt(para_errors.first);
951  double zerror = sqrt(para_errors.second);
952  clus_sigma(1) = zerror;
953  clus_sigma(0) = phierror;
954 
955  return clus_sigma;
956 }
957 
958 // new one
959 void HelicalFitter::getLocalDerivativesXY(Surface surf, Acts::Vector3 global, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5], unsigned int layer)
960 {
961  // Calculate the derivatives of the residual wrt the track parameters numerically
962  std::vector<float> temp_fitpars;
963 
964  std::vector<float> fitpars_delta;
965  fitpars_delta.push_back(0.1); // radius, cm
966  fitpars_delta.push_back(0.1); // X0, cm
967  fitpars_delta.push_back(0.1); // Y0, cm
968  fitpars_delta.push_back(0.1); // zslope, cm
969  fitpars_delta.push_back(0.1); // Z0, cm
970 
971  for(unsigned int ip = 0; ip < fitpars.size(); ++ip)
972  {
973  temp_fitpars.push_back(fitpars[ip]);
974  }
975 
976  // calculate projX and projY vectors once for the optimum fit parameters
977  if(Verbosity() > 1) std::cout << "Call get_helix_tangent for best fit fitpars" << std::endl;
978  std::pair<Acts::Vector3, Acts::Vector3> tangent = get_helix_tangent(fitpars, global);
979 
980  Acts::Vector3 projX(0,0,0), projY(0,0,0);
981  get_projectionXY(surf, tangent, projX, projY);
982 
983  Acts::Vector3 intersection = get_helix_surface_intersection(surf, temp_fitpars, global);
984 
985  // loop over the track fit parameters
986  for(unsigned int ip = 0; ip < fitpars.size(); ++ip)
987  {
988  Acts::Vector3 intersection_delta[2];
989  for(int ipm = 0; ipm < 2; ++ipm)
990  {
991  temp_fitpars[ip] = fitpars[ip]; // reset to best fit value
992  float deltapm = pow(-1.0, ipm);
993  temp_fitpars[ip] += deltapm * fitpars_delta[ip];
994 
995  Acts::Vector3 temp_intersection = get_helix_surface_intersection(surf, temp_fitpars, global);
996  intersection_delta[ipm] = temp_intersection - intersection;
997  }
998  Acts::Vector3 average_intersection_delta = (intersection_delta[0] - intersection_delta[1]) / (2 * fitpars_delta[ip]);
999 
1000  if(Verbosity() > 1)
1001  {std::cout << " average_intersection_delta / delta " << average_intersection_delta(0) << " " << average_intersection_delta(1) << " " << average_intersection_delta(2) << std::endl;}
1002 
1003  // calculate the change in fit for X and Y
1004  // - note negative sign from ATLAS paper is dropped here because mille wants the derivative of the fit, not the derivative of the residual
1005  lcl_derivativeX[ip] = average_intersection_delta.dot(projX);
1006  lcl_derivativeY[ip] = average_intersection_delta.dot(projY);
1007  if(Verbosity() > 1)
1008  {std::cout << " layer " << layer << " ip " << ip << " derivativeX " << lcl_derivativeX[ip] << " "
1009  << " derivativeY " << lcl_derivativeY[ip] << std::endl;}
1010 
1011  temp_fitpars[ip] = fitpars[ip];
1012  }
1013 }
1014 
1015 
1016 void HelicalFitter::getLocalVtxDerivativesXY(SvtxTrack& track, Acts::Vector3 event_vtx, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5])
1017 {
1018  // Calculate the derivatives of the residual wrt the track parameters numerically
1019  std::vector<float> temp_fitpars;
1020  Acts::Vector3 track_vtx (track.get_x(),track.get_y(),track.get_z());
1021 
1022  std::vector<float> fitpars_delta;
1023  fitpars_delta.push_back(0.1); // radius, cm
1024  fitpars_delta.push_back(0.1); // X0, cm
1025  fitpars_delta.push_back(0.1); // Y0, cm
1026  fitpars_delta.push_back(0.1); // zslope, cm
1027  fitpars_delta.push_back(0.1); // Z0, cm
1028 
1029  for(unsigned int ip = 0; ip < fitpars.size(); ++ip)
1030  {
1031  temp_fitpars.push_back(fitpars[ip]);
1032  }
1033 
1034  // calculate projX and projY vectors once for the optimum fit parameters
1035  if(Verbosity() > 1) {std::cout << "Call get_helix_tangent for best fit fitpars" << std::endl;}
1036 
1037  Acts::Vector3 perigeeNormal (track.get_px(),track.get_py(),track.get_pz()) ;
1038 
1039  // loop over the track fit parameters
1040  for(unsigned int ip = 0; ip < fitpars.size(); ++ip)
1041  {
1042  Acts::Vector3 localPerturb[2];
1043  Acts::Vector3 paperPerturb[2]; // for local derivative calculation like from the paper
1044 
1045  for(int ipm = 0; ipm < 2; ++ipm)
1046  {
1047  temp_fitpars[ip] = fitpars[ip]; // reset to best fit value
1048  float deltapm = pow(-1.0, ipm);
1049  temp_fitpars[ip] += deltapm * fitpars_delta[ip];
1050 
1051  Acts::Vector3 temp_track_vtx = get_helix_vtx(event_vtx, temp_fitpars); // temporary pca
1052  paperPerturb[ipm] = temp_track_vtx; // for og local derivative calculation
1053 
1054  // old method is next two lines
1055  Acts::Vector3 localtemp_track_vtx = globalvtxToLocalvtx(track, event_vtx, temp_track_vtx);
1056  localPerturb[ipm] = localtemp_track_vtx ;
1057 
1058  if(Verbosity() > 1)
1059  {
1060  std::cout << "vtx local parameter " << ip << " with ipm " << ipm << " deltapm " << deltapm << " :" << std::endl;
1061  std::cout<<" fitpars "<< fitpars[ip]<<" temp_fitpars "<<temp_fitpars[ip]<<std::endl;
1062  std::cout << " localtmp_track_vtx: "<< localtemp_track_vtx<<std::endl;
1063  }
1064  }
1065 
1066  Acts::Vector3 projX(0,0,0), projY(0,0,0);
1067  get_projectionVtxXY(track, event_vtx, projX, projY);
1068 
1069  Acts::Vector3 average_vtxX = (paperPerturb[0] - paperPerturb[1]) / (2 * fitpars_delta[ip]);
1070  Acts::Vector3 average_vtxY = (paperPerturb[0] - paperPerturb[1]) / (2 * fitpars_delta[ip]);
1071 
1072  // average_vtxX and average_vtxY are the numerical results in global coords for d(fit)/d(par)
1073  // The ATLAS paper formula is for the derivative of the residual, which is (measurement - fit = event vertex - track vertex)
1074  // Millepede wants the derivative of the fit, so we drop the minus sign from the paper
1075  lcl_derivativeX[ip] = average_vtxX.dot(projX); //
1076  lcl_derivativeY[ip] = average_vtxY.dot(projY);
1077 
1078  temp_fitpars[ip] = fitpars[ip];
1079  }
1080 }
1081 
1082 void HelicalFitter::getGlobalDerivativesXY(Surface surf, Acts::Vector3 global, Acts::Vector3 fitpoint, const std::vector<float>& fitpars, float glbl_derivativeX[6], float glbl_derivativeY[6], unsigned int layer)
1083 {
1084  Acts::Vector3 unitx(1, 0, 0);
1085  Acts::Vector3 unity(0, 1, 0);
1086  Acts::Vector3 unitz(0, 0, 1);
1087 
1088  // calculate projX and projY vectors once for the optimum fit parameters
1089  std::pair<Acts::Vector3, Acts::Vector3> tangent = get_helix_tangent(fitpars, global); // should this be global, not fitpoint?
1090 
1091  Acts::Vector3 projX(0,0,0), projY(0,0,0);
1092  get_projectionXY(surf, tangent, projX, projY);
1093 
1094  // translations
1095  glbl_derivativeX[3] = unitx.dot(projX);
1096  glbl_derivativeX[4] = unity.dot(projX);
1097  glbl_derivativeX[5] = unitz.dot(projX);
1098 
1099  glbl_derivativeY[3] = unitx.dot(projY);
1100  glbl_derivativeY[4] = unity.dot(projY);
1101  glbl_derivativeY[5] = unitz.dot(projY);
1102 
1103  /*
1104  // note: the global derivative sign should be reversed from the ATLAS paper
1105  // because mille wants the derivative of the fit, while the ATLAS paper gives the derivative of the residual.
1106  // But this sign reversal does NOT work.
1107  // Verified that not reversing the sign here produces the correct sign of the prediction of the residual..
1108  for(unsigned int i = 3; i < 6; ++i)
1109  {
1110  glbl_derivativeX[i] *= -1.0;
1111  glbl_derivativeY[i] *= -1.0;
1112  }
1113  */
1114  // rotations
1115  // need center of sensor to intersection point
1116  Acts::Vector3 sensorCenter = surf->center(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm; // convert to cm
1117  Acts::Vector3 OM = fitpoint - sensorCenter; // this effectively reverses the sign from the ATLAS paper
1118 
1119  glbl_derivativeX[0] = (unitx.cross(OM)).dot(projX);
1120  glbl_derivativeX[1] = (unity.cross(OM)).dot(projX);
1121  glbl_derivativeX[2] = (unitz.cross(OM)).dot(projX);
1122 
1123  glbl_derivativeY[0] = (unitx.cross(OM)).dot(projY);
1124  glbl_derivativeY[1] = (unity.cross(OM)).dot(projY);
1125  glbl_derivativeY[2] = (unitz.cross(OM)).dot(projY);
1126 
1127 
1128  if(Verbosity() > 1)
1129  {
1130  for(int ip=0;ip<6;++ip)
1131  {
1132  std::cout << " layer " << layer << " ip " << ip << " glbl_derivativeX " << glbl_derivativeX[ip] << " "
1133  << " glbl_derivativeY " << glbl_derivativeY[ip] << std::endl;
1134  }
1135  }
1136 
1137 }
1138 
1139 
1140 void HelicalFitter::getGlobalVtxDerivativesXY(SvtxTrack& track, Acts::Vector3 event_vtx, float glbl_derivativeX[3], float glbl_derivativeY[3])
1141 {
1142  Acts::Vector3 unitx(1, 0, 0);
1143  Acts::Vector3 unity(0, 1, 0);
1144  Acts::Vector3 unitz(0, 0, 1);
1145 
1146  Acts::Vector3 track_vtx (track.get_x(),track.get_y(),track.get_z());
1147  Acts::Vector3 mom(track.get_px(),track.get_py(),track.get_pz());
1148 
1149  // calculate projX and projY vectors once for the optimum fit parameters
1150  Acts::Vector3 projX(0,0,0), projY(0,0,0);
1151  get_projectionVtxXY(track, event_vtx, projX, projY);
1152 
1153  // translations
1154  glbl_derivativeX[0] = unitx.dot(projX);
1155  glbl_derivativeX[1] = unity.dot(projX);
1156  glbl_derivativeX[2] = unitz.dot(projX);
1157  glbl_derivativeY[0] = unitx.dot(projY);
1158  glbl_derivativeY[1] = unity.dot(projY);
1159  glbl_derivativeY[2] = unitz.dot(projY);
1160 
1161  // The derivation in the ATLAS paper used above gives the derivative of the residual (= measurement - fit)
1162  // pede wants the derivative of the fit, so we reverse that - valid if our residual is (event vertex - track vertex)
1163 
1164  // Verified that reversing these signs produces the correct sign and magnitude for the prediction of the residual.
1165  // tested this by offsetting the simulated event vertex with zero misalignments. Pede fit reproduced simulated (xvtx, yvtx) within 7%.
1166  // -- test gave zero for zvtx param, since this is determined relative to the measured event z vertex.
1167  for(int i = 0; i<3;++i)
1168  {
1169  glbl_derivativeX[i] *= -1.0;
1170  glbl_derivativeY[i] *= -1.0;
1171  }
1172 
1173 }
1174 
1175 void HelicalFitter::get_projectionXY(Surface surf, std::pair<Acts::Vector3, Acts::Vector3> tangent, Acts::Vector3& projX, Acts::Vector3& projY)
1176 {
1177  // we only need the direction part of the tangent
1178  Acts::Vector3 tanvec = tangent.second;
1179  // get the plane of the surface
1180  Acts::Vector3 sensorCenter = surf->center(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm; // convert to cm
1181  // sensorNormal is the Z vector
1182  Acts::Vector3 Z = -surf->normal(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
1183  // get surface X and Y unit vectors in global frame
1184  // transform Xlocal = 1.0 to global, subtract the surface center, normalize to 1
1185  Acts::Vector3 xloc(1.0,0.0,0.0); //local coord unit vector in x
1186  Acts::Vector3 xglob = surf->transform(_tGeometry->geometry().getGeoContext()) * (xloc * Acts::UnitConstants::cm);
1187  xglob /= Acts::UnitConstants::cm;
1188  Acts::Vector3 yloc(0.0,1.0,0.0);
1189  Acts::Vector3 yglob = surf->transform(_tGeometry->geometry().getGeoContext()) * (yloc * Acts::UnitConstants::cm);
1190  yglob /= Acts::UnitConstants::cm;
1191  Acts::Vector3 X = (xglob-sensorCenter) / (xglob-sensorCenter).norm();
1192  Acts::Vector3 Y = (yglob-sensorCenter) / (yglob-sensorCenter).norm();
1193  // see equation 31 of the ATLAS paper (and discussion) for this
1194  projX = X - (tanvec.dot(X) / tanvec.dot(Z)) * Z;
1195  projY = Y - (tanvec.dot(Y) / tanvec.dot(Z)) * Z;
1196  return;
1197 }
1198 
1199 void HelicalFitter::get_projectionVtxXY(SvtxTrack& track, Acts::Vector3 event_vtx, Acts::Vector3& projX, Acts::Vector3& projY)
1200 {
1201  Acts::Vector3 tanvec(track.get_px(),track.get_py(),track.get_pz());
1202  Acts::Vector3 normal(track.get_px(),track.get_py(),0);
1203 
1204  tanvec /= tanvec.norm();
1205  normal /= normal.norm();
1206 
1207  // get surface X and Y unit vectors in global frame
1208  Acts::Vector3 xloc(1.0,0.0,0.0);
1209  Acts::Vector3 yloc(0.0,0.0,1.0); // local y
1210  Acts::Vector3 xglob = localvtxToGlobalvtx(track, event_vtx, xloc);
1211  Acts::Vector3 yglob = yloc + event_vtx;
1212  Acts::Vector3 X = (xglob-event_vtx) / (xglob-event_vtx).norm(); // local unit vector transformed to global coordinates
1213  Acts::Vector3 Y = (yglob-event_vtx) / (yglob-event_vtx).norm();
1214 
1215  // see equation 31 of the ATLAS paper (and discussion) for this
1216  projX = X - (tanvec.dot(X) / tanvec.dot(normal)) * normal;
1217  projY = Y - (tanvec.dot(Y) / tanvec.dot(normal)) * normal;
1218 
1219  return;
1220 }
1221 
1222 unsigned int HelicalFitter::addSiliconClusters(std::vector<float>& fitpars, std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey>& cluskey_vec)
1223 {
1224 
1225  return TrackFitUtils::addClusters(fitpars, dca_cut, _tGeometry, _cluster_map, global_vec, cluskey_vec,0,6);
1226 }
1227 
1228 bool HelicalFitter::is_intt_layer_fixed(unsigned int layer)
1229 {
1230  bool ret = false;
1231  auto it = fixed_intt_layers.find(layer);
1232  if(it != fixed_intt_layers.end())
1233  ret = true;
1234 
1235  return ret;
1236 }
1237 
1238 bool HelicalFitter::is_mvtx_layer_fixed(unsigned int layer, unsigned int clamshell)
1239 {
1240  bool ret = false;
1241 
1242  std::pair pair = std::make_pair(layer,clamshell);
1243  auto it = fixed_mvtx_layers.find(pair);
1244  if(it != fixed_mvtx_layers.end())
1245  ret = true;
1246 
1247  return ret;
1248 }
1249 
1250 void HelicalFitter::set_intt_layer_fixed(unsigned int layer)
1251 {
1252  fixed_intt_layers.insert(layer);
1253 }
1254 
1255 void HelicalFitter::set_mvtx_layer_fixed(unsigned int layer, unsigned int clamshell)
1256 {
1257  fixed_mvtx_layers.insert(std::make_pair(layer,clamshell));
1258 }
1259 
1260 
1261 bool HelicalFitter::is_layer_param_fixed(unsigned int layer, unsigned int param)
1262 {
1263  bool ret = false;
1264  std::pair<unsigned int, unsigned int> pair = std::make_pair(layer, param);
1265  auto it = fixed_layer_params.find(pair);
1266  if(it != fixed_layer_params.end())
1267  ret = true;
1268 
1269  return ret;
1270 }
1271 
1272 void HelicalFitter::set_layer_param_fixed(unsigned int layer, unsigned int param)
1273 {
1274  std::pair<unsigned int, unsigned int> pair = std::make_pair(layer, param);
1275  fixed_layer_params.insert(pair);
1276 }
1277 
1278 void HelicalFitter::set_tpc_sector_fixed(unsigned int region, unsigned int sector, unsigned int side)
1279 {
1280  // make a combined subsector index
1281  unsigned int subsector = region * 24 + side * 12 + sector;
1282  fixed_sectors.insert(subsector);
1283 }
1284 
1285 bool HelicalFitter::is_tpc_sector_fixed(unsigned int layer, unsigned int sector, unsigned int side)
1286 {
1287  bool ret = false;
1288  unsigned int region = AlignmentDefs::getTpcRegion(layer);
1289  unsigned int subsector = region * 24 + side * 12 + sector;
1290  auto it = fixed_sectors.find(subsector);
1291  if(it != fixed_sectors.end())
1292  ret = true;
1293 
1294  return ret;
1295 }
1296 
1297 void HelicalFitter::correctTpcGlobalPositions(std::vector<Acts::Vector3> global_vec, std::vector<TrkrDefs::cluskey> cluskey_vec)
1298 {
1299  for(unsigned int iclus=0;iclus<cluskey_vec.size();++iclus)
1300  {
1301  auto cluskey = cluskey_vec[iclus];
1302  auto global = global_vec[iclus];
1303  const unsigned int trkrId = TrkrDefs::getTrkrId(cluskey);
1304  if(trkrId == TrkrDefs::tpcId)
1305  {
1306  // have to add corrections for TPC clusters after transformation to global
1307  int crossing = 0; // for now
1308  makeTpcGlobalCorrections(cluskey, crossing, global);
1309  global_vec[iclus] = global;
1310  }
1311  }
1312 }
1313 
1314 float HelicalFitter::getVertexResidual(Acts::Vector3 vtx)
1315 {
1316  float phi = atan2(vtx(1),vtx(0));
1317  float r = vtx(0)/cos(phi);
1318  float test_r = sqrt(vtx(0)*vtx(0)+vtx(1)*vtx(1));
1319 
1320  if(Verbosity() > 1)
1321  {
1322  std::cout << "my method position: " << vtx << std::endl;
1323  std::cout << "r " << r <<" phi: " << phi*180/M_PI<<" test_r"<< test_r << std::endl;
1324  }
1325  return r;
1326 }
1327 
1328 void HelicalFitter::get_dca(SvtxTrack& track,float& dca3dxy, float& dca3dz, float& dca3dxysigma, float& dca3dzsigma, Acts::Vector3 event_vertex)
1329 {
1330  //give trackseed
1331  dca3dxy = NAN;
1332  Acts::Vector3 track_vtx(track.get_x(),track.get_y(),track.get_z());
1333  Acts::Vector3 mom(track.get_px(),track.get_py(),track.get_pz());
1334 
1335  track_vtx -= event_vertex; // difference between track_vertex and event_vtx
1336 
1337  Acts::ActsSquareMatrix<3> posCov;
1338  for(int i = 0; i < 3; ++i)
1339  {
1340  for(int j = 0; j < 3; ++j)
1341  {
1342  posCov(i, j) = track.get_error(i, j);
1343  }
1344  }
1345 
1346  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
1347 
1348  float phi = atan2(r(1), r(0));
1349  Acts::RotationMatrix3 rot;
1350  Acts::RotationMatrix3 rot_T;
1351  phi *= -1;
1352  rot(0,0) = cos(phi);
1353  rot(0,1) = -sin(phi);
1354  rot(0,2) = 0;
1355  rot(1,0) = sin(phi);
1356  rot(1,1) = cos(phi);
1357  rot(1,2) = 0;
1358  rot(2,0) = 0;
1359  rot(2,1) = 0;
1360  rot(2,2) = 1;
1361  rot_T = rot.transpose();
1362 
1363  Acts::Vector3 pos_R = rot * track_vtx;
1364  Acts::ActsSquareMatrix<3> rotCov = rot * posCov * rot_T;
1365  dca3dxy = pos_R(0);
1366  dca3dz = pos_R(2);
1367  dca3dxysigma = sqrt(rotCov(0,0));
1368  dca3dzsigma = sqrt(rotCov(2,2));
1369 
1370  if(Verbosity()>1)
1371  {
1372  std::cout << " momentum X z: "<<r<< " phi: " << phi*180/M_PI << std::endl;
1373  std::cout << "dca3dxy " << dca3dxy << " dca3dz: " << dca3dz << " pos_R(1): "<<pos_R(1)<<" dca3dxysigma " << dca3dxysigma << " dca3dzsigma " << dca3dzsigma << std::endl;
1374  }
1375 
1376 }
1377 
1378 Acts::Vector3 HelicalFitter::globalvtxToLocalvtx(SvtxTrack& track, Acts::Vector3 event_vertex)
1379 {
1380 
1381  Acts::Vector3 track_vtx(track.get_x(),track.get_y(),track.get_z());
1382  Acts::Vector3 mom(track.get_px(),track.get_py(),track.get_pz());
1383  track_vtx -= event_vertex; // difference between track_vertex and event_vtx
1384 
1385  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
1386  float phi = atan2(r(1), r(0));
1387  Acts::RotationMatrix3 rot;
1388  Acts::RotationMatrix3 rot_T;
1389  phi *= -1;
1390  rot(0,0) = cos(phi);
1391  rot(0,1) = -sin(phi);
1392  rot(0,2) = 0;
1393  rot(1,0) = sin(phi);
1394  rot(1,1) = cos(phi);
1395  rot(1,2) = 0;
1396  rot(2,0) = 0;
1397  rot(2,1) = 0;
1398  rot(2,2) = 1;
1399  rot_T = rot.transpose();
1400 
1401  Acts::Vector3 pos_R = rot * track_vtx;
1402 
1403  if(Verbosity()>1)
1404  {
1405  std::cout << " momentum X z: "<<r<< " phi: " << phi*180/M_PI << std::endl;
1406  std::cout << " pos_R(0): "<<pos_R(0)<<" pos_R(1): "<<pos_R(1) << std::endl;
1407  }
1408  return pos_R;
1409 }
1410 
1411 Acts::Vector3 HelicalFitter::globalvtxToLocalvtx(SvtxTrack& track, Acts::Vector3 event_vertex, Acts::Vector3 PCA)
1412 {
1413  Acts::Vector3 mom(track.get_px(),track.get_py(),track.get_pz());
1414  PCA -= event_vertex; // difference between track_vertex and event_vtx
1415 
1416  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
1417  float phi = atan2(r(1), r(0));
1418  Acts::RotationMatrix3 rot;
1419  Acts::RotationMatrix3 rot_T;
1420  phi *= -1;
1421  rot(0,0) = cos(phi);
1422  rot(0,1) = -sin(phi);
1423  rot(0,2) = 0;
1424  rot(1,0) = sin(phi);
1425  rot(1,1) = cos(phi);
1426  rot(1,2) = 0;
1427  rot(2,0) = 0;
1428  rot(2,1) = 0;
1429  rot(2,2) = 1;
1430  rot_T = rot.transpose();
1431 
1432  Acts::Vector3 pos_R = rot * PCA;
1433 
1434  if(Verbosity()>1)
1435  {
1436  std::cout << " momentum X z: "<<r<< " phi: " << phi*180/M_PI << std::endl;
1437  std::cout << " pos_R(0): "<<pos_R(0)<<" pos_R(1): "<<pos_R(1) << std::endl;
1438  }
1439  return pos_R;
1440 }
1441 
1442 
1443 Acts::Vector3 HelicalFitter::localvtxToGlobalvtx(SvtxTrack& track, Acts::Vector3 event_vtx, Acts::Vector3 local)
1444 {
1445 
1446  //Acts::Vector3 track_vtx = local;
1447  Acts::Vector3 mom(track.get_px(),track.get_py(),track.get_pz());
1448  //std::cout << "first pos: " << pos << " mom: " << mom << std::endl;
1449  //local -= event_vertex; // difference between track_vertex and event_vtx
1450 
1451  Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
1452  float phi = atan2(r(1), r(0));
1453  Acts::RotationMatrix3 rot;
1454  Acts::RotationMatrix3 rot_T;
1455  // phi *= -1;
1456  rot(0,0) = cos(phi);
1457  rot(0,1) = -sin(phi);
1458  rot(0,2) = 0;
1459  rot(1,0) = sin(phi);
1460  rot(1,1) = cos(phi);
1461  rot(1,2) = 0;
1462  rot(2,0) = 0;
1463  rot(2,1) = 0;
1464  rot(2,2) = 1;
1465 
1466  rot_T = rot.transpose();
1467 
1468  Acts::Vector3 pos_R = rot * local;
1469  pos_R += event_vtx;
1470  if(Verbosity()>1)
1471  {
1472  std::cout << " momentum X z: "<<r<< " phi: " << phi*180/M_PI << std::endl;
1473  std::cout << " pos_R(0): "<<pos_R(0)<<" pos_R(1): "<<pos_R(1)<<"pos_R(2): "<<pos_R(2) << std::endl;
1474  }
1475  return pos_R;
1476 }
1477