Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SVReco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SVReco.C
1 #include "SVReco.h"
2 /*#include <trackbase_historic/SvtxCluster.h>
3 #include <trackbase_historic/SvtxClusterMap.h>
4 #include <trackbase_historic/SvtxHitMap.h>*/
5 #include <trackbase_historic/SvtxTrack.h>
6 #include <trackbase_historic/SvtxTrack_v1.h>
7 #include <trackbase_historic/SvtxTrackState_v1.h>
8 #include <trackbase_historic/SvtxVertex.h>
9 #include <trackbase_historic/SvtxVertex_v1.h>
10 #include <trackbase_historic/SvtxTrackMap.h>
11 #include <trackbase_historic/SvtxVertexMap.h>
12 
13 #include <trackbase/TrkrClusterContainer.h>
14 #include <trackbase/TrkrClusterv1.h>
15 #include <trackbase/TrkrHitSet.h>
16 #include <trackbase/TrkrHitSetContainer.h>
17 #include <trackbase/TrkrClusterHitAssoc.h>
18 #include <mvtx/MvtxDefs.h>
19 #include <intt/InttDefs.h>
20 
21 #include <g4jets/JetMap.h>
22 
23 #include <fun4all/Fun4AllReturnCodes.h>
24 #include <fun4all/PHTFileServer.h>
25 
26 //#include <g4detectors/PHG4CellContainer.h>
27 #include <g4detectors/PHG4Cell.h>
28 #include <g4detectors/PHG4CylinderGeomContainer.h>
29 #include <mvtx/CylinderGeom_Mvtx.h>
30 #include <intt/CylinderGeomIntt.h>
31 
32 #include <g4main/PHG4Hit.h>
33 #include <g4main/PHG4HitContainer.h>
34 #include <g4main/PHG4TruthInfoContainer.h>
35 #include <g4main/PHG4Particle.h>
36 
37 #include <phgenfit/Fitter.h>
38 #include <phgenfit/PlanarMeasurement.h>
39 #include <phgenfit/Track.h>
40 #include <phgenfit/SpacepointMeasurement.h>
41 
42 #include <phool/getClass.h>
43 #include <phool/phool.h>
44 #include <phool/PHCompositeNode.h>
45 #include <phool/PHIODataNode.h>
46 #include <phool/PHNodeIterator.h>
47 
48 #include <phgeom/PHGeomUtility.h>
49 #include <phfield/PHFieldUtility.h>
50 
51 #include <GenFit/FieldManager.h>
52 #include <GenFit/GFRaveVertex.h>
53 #include <GenFit/GFRaveVertexFactory.h>
54 #include <GenFit/MeasuredStateOnPlane.h>
55 #include <GenFit/RKTrackRep.h>
56 #include <GenFit/StateOnPlane.h>
57 #include <GenFit/Track.h>
58 #include <GenFit/KalmanFitterInfo.h>
59 
60 //#include <HFJetTruthGeneration/HFJetDefs.h>
61 
62 #include <TClonesArray.h>
63 #include <TMatrixDSym.h>
64 #include <TTree.h>
65 #include <TVector3.h>
66 #include <TRandom3.h>
67 #include <TRotation.h>
68 
69 #include <iostream>
70 #include <utility>
71 #include <memory>
72 
73 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
74 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
75 #define LogWarning(exp) std::cout<<"WARNING: "<<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
76 
77 using namespace std;
78 
79 /*Rave
80 #include <rave/Version.h>
81 //#include <rave/Track.h>
82 #include <rave/VertexFactory.h>
83 #include <rave/ConstantMagneticField.h>
84 */
85 //GenFit
86 //#include <GenFit/GFRaveConverters.h>
87 
88 
89 /*
90  * Constructor
91  */
92 SVReco::SVReco(const string &name) :
93  _mag_field_file_name("/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"),
94  _mag_field_re_scaling_factor(1.4 / 1.5), //what is this and why?
95  _reverse_mag_field(false),
96  _fitter(NULL),
97  _track_fitting_alg_name("DafRef"),
98  _n_maps_layer(3),
99  _n_intt_layer(4),
100  _primary_pid_guess(11), //for the tracks
101  _cut_Ncluster(false),
102  _cut_min_pT(0.1),
103  _cut_dca(5.0), //probably in mm
104  _cut_chi2_ndf(5),
105  _use_ladder_geom(false),
106  _vertex_finder(NULL),
107  _vertexing_method("avf-smoothing:1"), /*need a list of these and their proper domains*/
108  _clustermap(NULL),
109  _trackmap(NULL),
110  _vertexmap(NULL),
111  _do_eval(false),
112  _verbosity(10),
113  _do_evt_display(false)
114 {
115 
116 }
117 
118 
119 int SVReco::InitEvent(PHCompositeNode *topNode) {
120  cout<<"getting SVR nodes"<<endl;
121  GetNodes(topNode);
122  //cout<<"got vertexing nodes"<<endl;
124  for(auto p : _main_rf_phgf_tracks) delete p;
125  _main_rf_phgf_tracks.clear();
126 
127  _svtxtrk_gftrk_map.clear();
128  _svtxtrk_wt_map.clear();
129  _svtxtrk_id.clear();
130 
131  //iterate over all tracks to find priary vertex and make rave/genfit objects
132  cout<<"start SVR track loop"<<endl;
133  for (SvtxTrackMap::Iter iter = _trackmap->begin(); iter != _trackmap->end();
134  ++iter) {
135  SvtxTrack* svtx_track = iter->second;
136  // do track cuts
137  if (!svtx_track || svtx_track->get_ndf()<40 || svtx_track->get_pt()<_cut_min_pT ||
138  (svtx_track->get_chisq()/svtx_track->get_ndf())>_cut_chi2_ndf ||
139  fabs(svtx_track->get_dca3d_xy())>_cut_dca || fabs(svtx_track->get_dca3d_z())>_cut_dca)
140  continue;
141 
142  int n_MVTX = 0, n_INTT = 0, n_TPC = 0;
143  //cout<<"Keys:";
144  for (SvtxTrack::ConstClusterKeyIter iter2 = svtx_track->begin_cluster_keys(); iter2!=svtx_track->end_cluster_keys(); iter2++) {
145  TrkrDefs::cluskey cluster_key = *iter2;
146  //cout<<cluster_key<<',';
147  //count where the hits are
148  float layer = (float) TrkrDefs::getLayer(cluster_key);
149  if (layer<_n_maps_layer) n_MVTX++;
150  else if (layer<_n_maps_layer+_n_intt_layer) n_INTT++;
151  else n_TPC++;
152  }//cluster loop
153  //cluster cuts
154  //cout<<"\n cluster loop with n_MVTX="<<n_MVTX<<" n_INTT="<<n_INTT<<" and nTPC="<<n_TPC<<endl;
155  if ( _cut_Ncluster && (n_MVTX<2 || n_INTT<2 || n_TPC<25) ){
156  continue;
157  }
158  //cout << (svtx_track->get_chisq()/svtx_track->get_ndf()) << ", " << n_TPC << ", " << svtx_track->get_pt() << endl;
159  //cout << svtx_track->get_ndf() << ", " << svtx_track->size_clusters() << endl;
160  //cout<<"making genfit"<<endl;
161  PHGenFit::Track* rf_phgf_track = MakeGenFitTrack(svtx_track); //convert SvtxTrack to GenFit Track
162  //cout<<"made genfit"<<endl;
163 
164  if (rf_phgf_track) {
165  //svtx_track->identify();
166  //make a map to connect SvtxTracks to their respective GenFit Tracks
167  _svtxtrk_id.push_back(svtx_track->get_id());
168  _svtxtrk_gftrk_map[svtx_track->get_id()] = _main_rf_phgf_tracks.size();
169  _main_rf_phgf_tracks.push_back(rf_phgf_track); //to be used by findSecondaryVerticies
170  }
171  }
172  cout<<"exit track loop ntracks="<<_main_rf_phgf_tracks.size()<<'\n';
173 
174  return Fun4AllReturnCodes::EVENT_OK;
175 }
176 
177 PHGenFit::Track* SVReco::getPHGFTrack(SvtxTrack* svtxtrk){
178  if(svtxtrk)return _main_rf_phgf_tracks[_svtxtrk_gftrk_map[svtxtrk->get_id()]];
179  else return NULL;
180 }
181 
182 
183 int SVReco::InitRun(PHCompositeNode *topNode) {
184 
185  CreateNodes(topNode);
186 
187  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
188  PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
189 
190  _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, field,
191  _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
192 
193  if (!_fitter) {
194  cerr << PHWHERE << endl;
195  return Fun4AllReturnCodes::ABORTRUN;
196  }
197 
198  _vertex_finder = new genfit::GFRaveVertexFactory(_verbosity);
199  _vertex_finder->setMethod(_vertexing_method.data());
200 
201  if (!_vertex_finder) {
202  std::cerr<< PHWHERE<<" bad run init no SVR"<<endl;
203  return Fun4AllReturnCodes::ABORTRUN;
204  }
205 
206  return Fun4AllReturnCodes::EVENT_OK;
207 }
208 
209 //note that this might only work for conversion like events
210 genfit::GFRaveVertex* SVReco::findSecondaryVertex(SvtxTrack* track1, SvtxTrack* track2) {
211  //_vertex_finder->setMethod("avr-smoothing:1");
212  _vertex_finder->setMethod("avr");
213  //_vertex_finder->setMethod("avf-smoothing:1");
214  //_vertex_finder->setMethod("kalman");
215  vector<genfit::GFRaveVertex*> rave_vertices_conversion;
216  vector<genfit::Track*> rf_gf_tracks_conversion;
217 
218  if (_svtxtrk_gftrk_map.find(track1->get_id())!=_svtxtrk_gftrk_map.end()&&
219  _svtxtrk_gftrk_map.find(track2->get_id())!=_svtxtrk_gftrk_map.end())
220  {
221  unsigned int trk_index = _svtxtrk_gftrk_map[track1->get_id()];
222  PHGenFit::Track* rf_phgf_track = _main_rf_phgf_tracks[trk_index];
223  rf_gf_tracks_conversion.push_back(rf_phgf_track->getGenFitTrack());
224  //check the genfit track is working
225  /*cout<<"Track Comparison Original:\n";
226  track1->identify();*/
227  //printGenFitTrackKinematics(rf_phgf_track);
228 
229  trk_index = _svtxtrk_gftrk_map[track2->get_id()];
230  rf_phgf_track = _main_rf_phgf_tracks[trk_index];
231  // track2->identify();
232  // printGenFitTrackKinematics(rf_phgf_track);
233  rf_gf_tracks_conversion.push_back(rf_phgf_track->getGenFitTrack());
234  }
235  if (rf_gf_tracks_conversion.size()>1){
236  try{
237  _vertex_finder->findVertices(&rave_vertices_conversion, rf_gf_tracks_conversion);
238  }catch (...){
239  std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
240  }
241  /* if(TMath::Abs(rave_vertices_conversion[0]->getChi2()-1.)>.3){
242  cout<<"PRINTING:\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\n";
243  rave_vertices_conversion[0]->Print();
244  track1->identify();
245  track2->identify();
246  cout<<"\n\n\n\n\n";
247  }*/
248  //if a vertex is found return ownership
249  //TODO check rave_verticies_conversion for mem leak
250  if(rave_vertices_conversion.size()>0) return rave_vertices_conversion[0];
251  else return NULL;
252  }//more than 1 track
253  else return NULL;
254 }
255 
256 SVReco::~SVReco(){
257  cout<<PHWHERE<<"delete"<<endl;
258  if(_fitter){
259  delete _fitter;
260  _fitter=NULL;
261  }
262  cout<<PHWHERE<<"delete"<<endl;
263  if(_vertex_finder){
264  delete _vertex_finder;
265  _vertex_finder=NULL;
266  }
267  cout<<PHWHERE<<"delete"<<endl;
268  for (std::vector<PHGenFit::Track*>::iterator i = _main_rf_phgf_tracks.begin(); i != _main_rf_phgf_tracks.end(); ++i)
269  {
270  if(*i)
271  delete *i;
272  *i=NULL;
273  }
274 }
275 
276 //prints the track using the trival verex for point extrapolation
277 void SVReco::printGenFitTrackKinematics(PHGenFit::Track* track){
278 
279  TVector3 vertex_position(0, 0, 0);
280 
281  std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = NULL;
282  try {
283  gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane> (track->extrapolateToPoint(vertex_position));
284  } catch (...) {
285  if (_verbosity >= 2)
286  LogWarning("extrapolateToPoint failed!");
287  }
288  TVector3 mom = gf_state_vertex_ca->getMom();
289  TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
290 
291  cout << "OUT Ex: " << sqrt(cov[0][0]) << ", Ey: " << sqrt(cov[1][1]) << ", Ez: " << sqrt(cov[2][2]) << endl;
292  cout << "OUT Px: " << mom.X() << ", Py: " << mom.Y() << ", Pz: " << mom.Z() << endl;
293 }
294 
295 
296 
297 //should be deprecated
298 void SVReco::reset_eval_variables(){
299  gf_prim_vtx_ntrk = rv_prim_vtx_ntrk = 0;
300  for (int i=0; i<3; i++){
301  gf_prim_vtx[i] = gf_prim_vtx_err[i] = -999;
302  rv_prim_vtx[i] = rv_prim_vtx_err[i] = -999;
303  }//i
304 
305  rv_sv_njets = 0;
306 
307  for (int ijet=0; ijet<10; ijet++){
308  rv_sv_jet_id[ijet] = -999;
309  rv_sv_jet_prop[ijet][0] = rv_sv_jet_prop[ijet][1] = -999;
310  rv_sv_jet_pT[ijet] = -999;
311  rv_sv_jet_px[ijet] = rv_sv_jet_py[ijet] = rv_sv_jet_pz[ijet] = -999;
312 
313  rv_sv_pT00_nvtx[ijet] = rv_sv_pT05_nvtx[ijet] = rv_sv_pT10_nvtx[ijet] = rv_sv_pT15_nvtx[ijet] = 0;
314 
315  for (int ivtx=0; ivtx<30; ivtx++){
316  rv_sv_pT00_vtx_ntrk[ijet][ivtx] = rv_sv_pT05_vtx_ntrk[ijet][ivtx] = rv_sv_pT10_vtx_ntrk[ijet][ivtx] = rv_sv_pT15_vtx_ntrk[ijet][ivtx] = 0;
317  rv_sv_pT00_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT05_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT10_vtx_ntrk_good[ijet][ivtx] = rv_sv_pT15_vtx_ntrk_good[ijet][ivtx] = 0;
318  rv_sv_pT00_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT05_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT10_vtx_ntrk_good_pv[ijet][ivtx] = rv_sv_pT15_vtx_ntrk_good_pv[ijet][ivtx] = 0;
319 
320  rv_sv_pT00_vtx_mass[ijet][ivtx] = rv_sv_pT00_vtx_mass_corr[ijet][ivtx] = rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
321  rv_sv_pT05_vtx_mass[ijet][ivtx] = rv_sv_pT05_vtx_mass_corr[ijet][ivtx] = rv_sv_pT05_vtx_pT[ijet][ivtx] = -999;
322  rv_sv_pT10_vtx_mass[ijet][ivtx] = rv_sv_pT10_vtx_mass_corr[ijet][ivtx] = rv_sv_pT10_vtx_pT[ijet][ivtx] = -999;
323  rv_sv_pT15_vtx_mass[ijet][ivtx] = rv_sv_pT15_vtx_mass_corr[ijet][ivtx] = rv_sv_pT15_vtx_pT[ijet][ivtx] = -999;
324 
325  rv_sv_pT00_vtx_x[ijet][ivtx] = rv_sv_pT00_vtx_y[ijet][ivtx] = rv_sv_pT00_vtx_z[ijet][ivtx] = -999;
326  rv_sv_pT00_vtx_ex[ijet][ivtx] = rv_sv_pT00_vtx_ey[ijet][ivtx] = rv_sv_pT00_vtx_ez[ijet][ivtx] = -999;
327  rv_sv_pT05_vtx_x[ijet][ivtx] = rv_sv_pT05_vtx_y[ijet][ivtx] = rv_sv_pT05_vtx_z[ijet][ivtx] = -999;
328  rv_sv_pT05_vtx_ex[ijet][ivtx] = rv_sv_pT05_vtx_ey[ijet][ivtx] = rv_sv_pT05_vtx_ez[ijet][ivtx] = -999;
329  rv_sv_pT10_vtx_x[ijet][ivtx] = rv_sv_pT10_vtx_y[ijet][ivtx] = rv_sv_pT10_vtx_z[ijet][ivtx] = -999;
330  rv_sv_pT10_vtx_ex[ijet][ivtx] = rv_sv_pT10_vtx_ey[ijet][ivtx] = rv_sv_pT10_vtx_ez[ijet][ivtx] = -999;
331  rv_sv_pT15_vtx_x[ijet][ivtx] = rv_sv_pT15_vtx_y[ijet][ivtx] = rv_sv_pT15_vtx_z[ijet][ivtx] = -999;
332  rv_sv_pT15_vtx_ex[ijet][ivtx] = rv_sv_pT15_vtx_ey[ijet][ivtx] = rv_sv_pT15_vtx_ez[ijet][ivtx] = -999;
333 
334  rv_sv_pT00_vtx_jet_theta[ijet][ivtx] = rv_sv_pT00_vtx_pT[ijet][ivtx] = -999;
335  rv_sv_pT05_vtx_jet_theta[ijet][ivtx] = rv_sv_pT05_vtx_pT[ijet][ivtx] = -999;
336  rv_sv_pT10_vtx_jet_theta[ijet][ivtx] = rv_sv_pT10_vtx_pT[ijet][ivtx] = -999;
337  rv_sv_pT15_vtx_jet_theta[ijet][ivtx] = rv_sv_pT15_vtx_pT[ijet][ivtx] = -999;
338 
339  rv_sv_pT00_vtx_chi2[ijet][ivtx] = rv_sv_pT00_vtx_ndf[ijet][ivtx] = -999;
340  rv_sv_pT05_vtx_chi2[ijet][ivtx] = rv_sv_pT05_vtx_ndf[ijet][ivtx] = -999;
341  rv_sv_pT10_vtx_chi2[ijet][ivtx] = rv_sv_pT10_vtx_ndf[ijet][ivtx] = -999;
342  rv_sv_pT15_vtx_chi2[ijet][ivtx] = rv_sv_pT15_vtx_ndf[ijet][ivtx] = -999;
343 
344  }//ivtx
345  }//ijet
346 
347  return;
348 }
349 
350 int SVReco::CreateNodes(PHCompositeNode *topNode){
351 
352  return Fun4AllReturnCodes::EVENT_OK;
353 }
354 
355 /*
356  * GetNodes():
357  * Get all the all the required nodes off the node tree
358  */
359 int SVReco::GetNodes(PHCompositeNode * topNode){
360  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
361  if (!_clustermap){
362  cout << PHWHERE << " TRKR_CLUSTERS node not found on node tree"
363  << endl;
364  return Fun4AllReturnCodes::ABORTEVENT;
365  }
366  // Input Svtx Tracks
367  _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
368  if (!_trackmap){
369  cout << PHWHERE << " SvtxTrackMap node not found on node tree"
370  << endl;
371  return Fun4AllReturnCodes::ABORTEVENT;
372  }
373  // Input Svtx Vertices
374  _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
375  if (!_vertexmap){
376  cout << PHWHERE << " SvtxVertexrMap node not found on node tree"
377  << endl;
378  return Fun4AllReturnCodes::ABORTEVENT;
379  }
380 
381  _geom_container_intt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
382  if (!_geom_container_intt)
383  {
384  cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
385  << endl;
386  return Fun4AllReturnCodes::ABORTEVENT;
387  }
388  _geom_container_maps = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
389  if (!_geom_container_maps)
390  {
391  cout << PHWHERE << "CYLINDERGEOM_MVTX node not found on node tree"
392  << endl;
393  return Fun4AllReturnCodes::ABORTEVENT;
394  }
395  if(!_vertex_finder){
396  std::cerr<< PHWHERE<<" bad run init no SVR"<<endl;
397  return Fun4AllReturnCodes::ABORTRUN;
398  }
399 
400 
401  return Fun4AllReturnCodes::EVENT_OK;
402 }
403 
404 //Edited version original from @sh-lim
405 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack){
406  if (!intrack){
407  cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
408  return NULL;
409  }
410 
411  if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
412  cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
413  return NULL;
414  }
415 
416  TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
417  TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
418  TMatrixDSym seed_cov(6);
419  for (int i=0; i<6; i++){
420  for (int j=0; j<6; j++){
421  seed_cov[i][j] = intrack->get_error(i,j);
422  }
423  }
424 
425  // Create measurements
426  std::vector<PHGenFit::Measurement*> measurements;
427 
428  for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
429  // unsigned int cluster_id = *iter;
430  TrkrCluster* cluster = _clustermap->findCluster(*iter);
431  if (!cluster) {
432  LogError("No cluster Found!");
433  continue;
434  }
435  float x = cluster->getPosition(0);
436  float y = cluster->getPosition(1);
437  float radius = sqrt(x*x+y*y);
438  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
439  seed_mom.SetPhi(pos.Phi());
440  seed_mom.SetTheta(pos.Theta());
441 
442  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
443  //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
444 
445  if (_use_ladder_geom){ //I don't understand this bool
446  unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
447  unsigned int layer = TrkrDefs::getLayer(*iter);
448  if (trkrid == TrkrDefs::mvtxId) {
449  int stave_index = MvtxDefs::getStaveId(*iter);
450  int chip_index = MvtxDefs::getChipId(*iter);
451 
452  double ladder_location[3] = { 0.0, 0.0, 0.0 };
453  //not exactly sure where the cylinder geoms are currently objectified. check this
454  CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
455  // returns the center of the sensor in world coordinates - used to get the ladder phi location
456  geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
457  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
458  n.RotateZ(geom->get_stave_phi_tilt());
459  }
460  else if (trkrid == TrkrDefs::inttId) {
461  //this may bug but it looks ok for now
462  CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
463  double hit_location[3] = { 0.0, 0.0, 0.0 };
464  geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
465 
466  n.SetXYZ(hit_location[0], hit_location[1], 0);
467  n.RotateZ(geom->get_strip_phi_tilt());
468  }
469  }//if use_ladder_geom
470  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
471  measurements.push_back(meas);
472  }//cluster loop
473  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
474  PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
475  track->addMeasurements(measurements);
476 
477  if (_fitter->processTrack(track, false) != 0) {
478  if (_verbosity >= 1)
479  LogWarning("Track fitting failed");
480  return NULL;
481  }
482 
483  track->getGenFitTrack()->setMcTrackId(intrack->get_id());
484 
485  return track;
486 }
487 
488 //inspired by PHG4TrackKalmanFitter
489 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack, const SvtxVertex* invertex){
490  if (!intrack){
491  cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
492  return NULL;
493  }
494 
495  if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
496  cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
497  return NULL;
498  }
499 
500  // Create measurements
501  std::vector<PHGenFit::Measurement*> measurements;
502 
503  //create space point measurement from vtx
504  if (invertex and invertex->size_tracks() > 1) {
505  TVector3 pos(invertex->get_x(), invertex->get_y(), invertex->get_z());
506  TMatrixDSym cov(3);
507  cov.Zero();
508  bool is_vertex_cov_sane = true;
509  cout<<"Making covarience for vtx measurement"<<endl;
510  for (unsigned int i = 0; i < 3; i++)
511  for (unsigned int j = 0; j < 3; j++) {
512  cov(i, j) = invertex->get_error(i, j);
513  }
514 
515  cout<<"Made covarience"<<endl;
516  if (is_vertex_cov_sane) {
517  PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
518  pos, cov);
519  measurements.push_back(meas);
520  }
521 
522  //convert SvtxTrack to matricies
523  TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
524  TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
525  TMatrixDSym seed_cov(6);
526  for (int i=0; i<6; i++){
527  for (int j=0; j<6; j++){
528  seed_cov[i][j] = intrack->get_error(i,j);
529  }
530  }
531  cout<<"Making track cluster measurments"<<endl;
532  //make measurements from the track clusters
533  for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
534  // unsigned int cluster_id = *iter;
535  TrkrCluster* cluster = _clustermap->findCluster(*iter);
536  if (!cluster) {
537  LogError("No cluster Found!");
538  continue;
539  }
540  float x = cluster->getPosition(0);
541  float y = cluster->getPosition(1);
542  float radius = sqrt(x*x+y*y);
543  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
544  seed_mom.SetPhi(pos.Phi());
545  seed_mom.SetTheta(pos.Theta());
546 
547  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
548  //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
549 
550  if (_use_ladder_geom){ //I don't understand this bool
551  unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
552  unsigned int layer = TrkrDefs::getLayer(*iter);
553  if (trkrid == TrkrDefs::mvtxId) {
554  int stave_index = MvtxDefs::getStaveId(*iter);
555  int chip_index = MvtxDefs::getChipId(*iter);
556 
557  double ladder_location[3] = { 0.0, 0.0, 0.0 };
558  //not exactly sure where the cylinder geoms are currently objectified. check this
559  CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
560  // returns the center of the sensor in world coordinates - used to get the ladder phi location
561  geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
562  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
563  n.RotateZ(geom->get_stave_phi_tilt());
564  }
565  else if (trkrid == TrkrDefs::inttId) {
566  //this may bug but it looks ok for now
567  CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
568  double hit_location[3] = { 0.0, 0.0, 0.0 };
569  geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
570 
571  n.SetXYZ(hit_location[0], hit_location[1], 0);
572  n.RotateZ(geom->get_strip_phi_tilt());
573  }
574  }//if use_ladder_geom
575  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
576  measurements.push_back(meas);
577  }//cluster loop
578  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
579  PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
580  track->addMeasurements(measurements);
581 
582  if (_fitter->processTrack(track, false) != 0) {
583  if (_verbosity >= 1)
584  LogWarning("Track fitting failed");
585  return NULL;
586  }
587  track->getGenFitTrack()->setMcTrackId(intrack->get_id());
588  return track;
589  }//valid vtx
590  else{
591  cerr<<PHWHERE<<": invalid vertex"<<endl;
592  return NULL;
593  }
594 
595 }
596 
597 //inspired by PHG4TrackKalmanFitter
598 PHGenFit::Track* SVReco::MakeGenFitTrack(const SvtxTrack* intrack, const genfit::GFRaveVertex* invertex){
599  if (!intrack){
600  cerr << PHWHERE << " Input SvtxTrack is NULL!" << endl;
601  return NULL;
602  }
603 
604  if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
605  cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
606  return NULL;
607  }
608 
609  // Create measurements
610  std::vector<PHGenFit::Measurement*> measurements;
611 
612  //create space point measurement from vtx
613  //TODO check that getNTracks is properly initialized
614  if (invertex and invertex->getNTracks() > 1) {
615  TVector3 pos=invertex->getPos();
616  TMatrixDSym cov = invertex->getCov();
617  PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
618  pos, cov);
619  measurements.push_back(meas);
620 
621  //convert SvtxTrack to matricies
622  TVector3 seed_pos(intrack->get_x(), intrack->get_y(), intrack->get_z());
623  TVector3 seed_mom(intrack->get_px(), intrack->get_py(), intrack->get_pz()); //mom stands for momentum
624  TMatrixDSym seed_cov(6);
625  for (int i=0; i<6; i++){
626  for (int j=0; j<6; j++){
627  seed_cov[i][j] = intrack->get_error(i,j);
628  }
629  }
630  cout<<"Making track cluster measurments"<<endl;
631  //make measurements from the track clusters
632  for (auto iter = intrack->begin_cluster_keys(); iter != intrack->end_cluster_keys(); ++iter){
633  // unsigned int cluster_id = *iter;
634  TrkrCluster* cluster = _clustermap->findCluster(*iter);
635  if (!cluster) {
636  LogError("No cluster Found!");
637  continue;
638  }
639  float x = cluster->getPosition(0);
640  float y = cluster->getPosition(1);
641  float radius = sqrt(x*x+y*y);
642  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
643  seed_mom.SetPhi(pos.Phi());
644  seed_mom.SetTheta(pos.Theta());
645 
646  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
647  //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
648 
649  if (_use_ladder_geom){ //I don't understand this bool
650  unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
651  unsigned int layer = TrkrDefs::getLayer(*iter);
652  if (trkrid == TrkrDefs::mvtxId) {
653  int stave_index = MvtxDefs::getStaveId(*iter);
654  int chip_index = MvtxDefs::getChipId(*iter);
655 
656  double ladder_location[3] = { 0.0, 0.0, 0.0 };
657  //not exactly sure where the cylinder geoms are currently objectified. check this
658  CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
659  // returns the center of the sensor in world coordinates - used to get the ladder phi location
660  geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
661  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
662  n.RotateZ(geom->get_stave_phi_tilt());
663  }
664  else if (trkrid == TrkrDefs::inttId) {
665  //this may bug but it looks ok for now
666  CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
667  double hit_location[3] = { 0.0, 0.0, 0.0 };
668  geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
669 
670  n.SetXYZ(hit_location[0], hit_location[1], 0);
671  n.RotateZ(geom->get_strip_phi_tilt());
672  }
673  }//if use_ladder_geom
674  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
675  measurements.push_back(meas);
676  }//cluster loop
677  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
678  PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
679  track->addMeasurements(measurements);
680 
681  if (_fitter->processTrack(track, false) != 0) {
682  if (_verbosity >= 1)
683  LogWarning("Track fitting failed");
684  return NULL;
685  }
686  track->getGenFitTrack()->setMcTrackId(intrack->get_id());
687  return track;
688  }//valid vtx
689  else{
690  cerr<<PHWHERE<<": invalid vertex"<<endl;
691  return NULL;
692  }
693 
694 }
695 
696 /*inspired by PHG4TrackKalmanFitter
697 PHGenFit::Track* SVReco::MakeGenFitTrack(const PHGenFit::Track* intrack, const genfit::GFRaveVertex* invertex){
698  if (!intrack){
699  cerr << PHWHERE << " Input PHGenFit::Track is NULL!" << endl;
700  return NULL;
701  }
702 
703  if (_use_ladder_geom and !_geom_container_intt and !_geom_container_maps) {
704  cout << PHWHERE << "No PHG4CylinderGeomContainer found!" << endl;
705  return NULL;
706  }
707 
708  // Create measurements
709  std::vector<PHGenFit::Measurement*> measurements;
710 
711  //create space point measurement from vtx
712  //TODO check that getNTracks is properly initialized
713  if (invertex and invertex->getNTracks() > 1) {
714  TVector3 pos=invertex->getPos();
715  TMatrixDSym cov = invertex->getCov();
716  PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(
717  pos, cov);
718  measurements.push_back(meas);
719 
720  //convert SvtxTrack to matricies
721  TVector3 seed_pos = intrack->getGenFitTrack()->getPos();
722  TVector3 seed_mom = intrack->get_mom(); //mom stands for momentum
723  TMatrixDSym seed_cov = intrack->getGenFitTrack()->getCov();
724 
725  cout<<"Making track cluster measurments"<<endl;
726  //make measurements from the track clusters
727  for (auto iter = intrack->get_cluster_keys().begin(); iter != intrack->get_cluster_keys().end(); ++iter){
728  // unsigned int cluster_id = *iter;
729  TrkrCluster* cluster = _clustermap->findCluster(*iter);
730  if (!cluster) {
731  LogError("No cluster Found!");
732  continue;
733  }
734  float x = cluster->getPosition(0);
735  float y = cluster->getPosition(1);
736  float radius = sqrt(x*x+y*y);
737  TVector3 pos(cluster->getPosition(0), cluster->getPosition(1), cluster->getPosition(2));
738  seed_mom.SetPhi(pos.Phi());
739  seed_mom.SetTheta(pos.Theta());
740 
741  TVector3 n(cluster->getPosition(0), cluster->getPosition(1), 0);
742  //cout<<"Cluster with {"<<cluster->getPosition(0)<<','<<cluster->getPosition(0)<<"}\n";
743 
744  if (_use_ladder_geom){ //I don't understand this bool
745  unsigned int trkrid = TrkrDefs::getTrkrId(*iter);
746  unsigned int layer = TrkrDefs::getLayer(*iter);
747  if (trkrid == TrkrDefs::mvtxId) {
748  int stave_index = MvtxDefs::getStaveId(*iter);
749  int chip_index = MvtxDefs::getChipId(*iter);
750 
751  double ladder_location[3] = { 0.0, 0.0, 0.0 };
752  //not exactly sure where the cylinder geoms are currently objectified. check this
753  CylinderGeom_Mvtx *geom = (CylinderGeom_Mvtx*) _geom_container_maps->GetLayerGeom(layer);
754  // returns the center of the sensor in world coordinates - used to get the ladder phi location
755  geom->find_sensor_center(stave_index, 0, 0, chip_index, ladder_location);//the mvtx module and half stave are 0
756  n.SetXYZ(ladder_location[0], ladder_location[1], 0);
757  n.RotateZ(geom->get_stave_phi_tilt());
758  }
759  else if (trkrid == TrkrDefs::inttId) {
760  //this may bug but it looks ok for now
761  CylinderGeomIntt* geom = (CylinderGeomIntt*) _geom_container_intt->GetLayerGeom(layer);
762  double hit_location[3] = { 0.0, 0.0, 0.0 };
763  geom->find_segment_center(InttDefs::getLadderZId(*iter),InttDefs::getLadderPhiId(*iter), hit_location);
764 
765  n.SetXYZ(hit_location[0], hit_location[1], 0);
766  n.RotateZ(geom->get_strip_phi_tilt());
767  }
768  }//if use_ladder_geom
769  PHGenFit::Measurement* meas = new PHGenFit::PlanarMeasurement(pos, n,radius*cluster->getPhiError(), cluster->getZError());
770  measurements.push_back(meas);
771  }//cluster loop
772  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(_primary_pid_guess);
773  PHGenFit::Track* track(new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov));
774  track->addMeasurements(measurements);
775 
776  if (_fitter->processTrack(track, false) != 0) {
777  if (_verbosity >= 1)
778  LogWarning("Track fitting failed");
779  return NULL;
780  }
781  track->getGenFitTrack()->setMcTrackId(intrack->get_id());
782  return track;
783  }//valid vtx
784  else{
785  cerr<<PHWHERE<<": invalid vertex"<<endl;
786  return NULL;
787  }
788 
789 }*/
790 
791 //From PHG4TrackKalmanFitter
792 SvtxVertex* SVReco::GFRVVtxToSvtxVertex(genfit::GFRaveVertex* rave_vtx)const{
793  SvtxVertex* svtx_vtx= new SvtxVertex_v1();
794 
795  svtx_vtx->set_chisq(rave_vtx->getChi2());
796  svtx_vtx->set_ndof(rave_vtx->getNdf());
797  svtx_vtx->set_position(0, rave_vtx->getPos().X());
798  svtx_vtx->set_position(1, rave_vtx->getPos().Y());
799  svtx_vtx->set_position(2, rave_vtx->getPos().Z());
800 
801  for (int i = 0; i < 3; i++)
802  for (int j = 0; j < 3; j++)
803  svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
804 
805  for (unsigned int i = 0; i < rave_vtx->getNTracks(); i++) {
806  //TODO Assume id's are sync'ed between _trackmap_refit and gf_tracks, need to change?
807  const genfit::Track* rave_track =
808  rave_vtx->getParameters(i)->getTrack();
809  for (unsigned int j = 0; j < _main_rf_phgf_tracks.size(); j++) {
810  if (rave_track == _main_rf_phgf_tracks[j]->getGenFitTrack())
811  svtx_vtx->insert_track(j);
812  }
813  }
814  return svtx_vtx;
815 }
816 
817 /*
818  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
819  */
820 void SVReco::FillVertexMap(
821  const std::vector<genfit::GFRaveVertex*>& rave_vertices,
822  const std::vector<genfit::Track*>& gf_tracks){
823 
824  for (unsigned int ivtx=0; ivtx<rave_vertices.size(); ++ivtx){
825  genfit::GFRaveVertex* rave_vtx = rave_vertices[ivtx];
826  rv_prim_vtx_ntrk = rave_vtx->getNTracks();
827 
828  //TVector3 vertex_position(rv_prim_vtx[0], rv_prim_vtx[1], rv_prim_vtx[2]);
829 
830  //cout << "N TRK gf: " << gf_tracks.size() << ", rv: " << rv_prim_vtx_ntrk << endl;
831 
832  for (int itrk=0; itrk< rv_prim_vtx_ntrk; itrk++){
833 
834  TVector3 rvtrk_mom = rave_vtx->getParameters(itrk)->getMom();
835  float rvtrk_w = rave_vtx->getParameters(itrk)->getWeight();
836 
837  unsigned int rvtrk_mc_id = rave_vtx->getParameters(itrk)->getTrack()->getMcTrackId();
838  _svtxtrk_wt_map[rvtrk_mc_id] = rvtrk_w;
839 
840  //cout << "w: " << rvtrk_w << ", mc id: " << rvtrk_mc_id << endl;
841  /*
842  SvtxTrack* svtx_track = _trackmap->get(rvtrk_mc_id);
843 
844  cout << "rave trk, px: " << rvtrk_mom.Px() << ", py: " << rvtrk_mom.Py() << ", pz: " << rvtrk_mom.Pz() << endl;
845  cout << "svtx trk, px: " << svtx_track->get_px() << ", py: " << svtx_track->get_py() << ", pz: " << svtx_track->get_pz() << endl;
846  */
847 
848  /*
849  for (SvtxTrackMap::ConstIter iter3=_trackmap->begin(); iter3!=_trackmap->end(); iter3++)
850  {
851  SvtxTrack* svtx_track = iter3->second;
852 
853  if ( fabs((svtx_track->get_px()-rvtrk_mom.Px())/svtx_track->get_px())<0.10
854  && fabs((svtx_track->get_py()-rvtrk_mom.Py())/svtx_track->get_py())<0.10
855  && fabs((svtx_track->get_pz()-rvtrk_mom.Pz())/svtx_track->get_pz())<0.10 ){
856  cout << "rave trk, px: " << rvtrk_mom.Px() << ", py: " << rvtrk_mom.Py() << ", pz: " << rvtrk_mom.Pz() << endl;
857  //cout << "ggff trk, px: " << gftrk_mom.Px() << ", py: " << gftrk_mom.Py() << ", pz: " << gftrk_mom.Pz() << endl;
858  cout << "svtx trk, px: " << svtx_track->get_px() << ", py: " << svtx_track->get_py() << ", pz: " << svtx_track->get_pz() << endl;
859  }
860  }//iter3
861  */
862 
863  //unsigned int trk_id = svtxtrk_id[itrk];
864 
865  /*
866  cout << "W: " << w_trk
867  << ", id: " << rave_vtx->getParameters(itrk)->GetUniqueID()
868  << ", px: " << rave_vtx->getParameters(itrk)->getMom().Px()
869  << ", py: " << rave_vtx->getParameters(itrk)->getMom().Py()
870  << ", pz: " << rave_vtx->getParameters(itrk)->getMom().Pz()
871  << endl;
872  */
873 
874  //if (svtxtrk_gftrk_map.find(svtx_track->get_id())!=svtxtrk_gftrk_map.end()){
875  //}
876 
877  //TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom();
878  }//primvtx_tracks loop
879  }//rave verticies loop
880  return;
881 }
882 
883 int SVReco::GetSVMass_mom(
884  const genfit::GFRaveVertex* rave_vtx,
885  float & vtx_mass,
886  float & vtx_px,
887  float & vtx_py,
888  float & vtx_pz,
889  int & ntrk_good_pv
890  ){
891 
892  float sum_E = 0, sum_px = 0, sum_py = 0, sum_pz = 0;
893 
894  int N_good = 0, N_good_pv = 0;
895 
896  for (unsigned int itrk=0; itrk<rave_vtx->getNTracks(); itrk++){
897  TVector3 mom_trk = rave_vtx->getParameters(itrk)->getMom();
898 
899  double w_trk = rave_vtx->getParameters(itrk)->getWeight();
900 
901  sum_px += mom_trk.X();
902  sum_py += mom_trk.Y();
903  sum_pz += mom_trk.Z();
904  sum_E += sqrt(mom_trk.Mag2() + 0.140*0.140);
905 
906  //cout << "W: " << w_trk << endl;
907  if ( w_trk>0.7 ){
908  N_good++;
909 
910  unsigned int rvtrk_mc_id = rave_vtx->getParameters(itrk)->getTrack()->getMcTrackId();
911  //cout << "mc_id: " << rvtrk_mc_id << ", wt: " << svtxtrk_wt_map[rvtrk_mc_id] << endl;
912  if ( _svtxtrk_wt_map[rvtrk_mc_id]>0.7 ){
913  N_good_pv++;
914  }
915  }//
916 
917  }//itrk
918 
919  vtx_mass = sqrt(sum_E*sum_E - sum_px*sum_px - sum_py*sum_py - sum_pz*sum_pz);
920  vtx_px = sum_px;
921  vtx_py = sum_py;
922  vtx_pz = sum_pz;
923 
924  ntrk_good_pv = N_good_pv;
925 
926  //cout << "Mass: " << vtx_mass << ", pT: " << vtx_pT << endl;
927  return N_good;
928 }
929 
930 //Seems deprecated
931 void SVReco::FillSVMap(
932  const std::vector<genfit::GFRaveVertex*>& rave_vertices,
933  const std::vector<genfit::Track*>& gf_tracks){
934 
935  return;
936 }
937 
938 PHGenFit::Track* SVReco::refitTrack(SvtxVertex* vtx, SvtxTrack* svtxtrk){
939  PHGenFit::Track* gftrk = MakeGenFitTrack(svtxtrk,vtx);
940  if(gftrk) {
941  cout<<"good genfit refit"<<endl;
942  //MakeSvtxTrack(svtxtrk,gftrk,vtx);
943  return gftrk;
944  }
945  else {
946  cout<<"No refit possible"<<endl;
947  return NULL;
948  }
949 }
950 
951 PHGenFit::Track* SVReco::refitTrack(genfit::GFRaveVertex* vtx, SvtxTrack* svtxtrk){
952  PHGenFit::Track* gftrk = MakeGenFitTrack(svtxtrk,vtx);
953  if(gftrk) {
954  cout<<"good genfit refit"<<endl;
955  //MakeSvtxTrack(svtxtrk,gftrk,vtx);
956  return gftrk;
957  }
958  else {
959  cout<<"No refit possible"<<endl;
960  return NULL;
961  }
962 }
963 
964 /*PHGenFit::Track* SVReco::refitTrack(genfit::GFRaveVertex* vtx, PHGenFit::Track* phgf_track){
965  PHGenFit::Track* gftrk = MakeGenFitTrack(phgf_track,vtx);
966  if(gftrk) {
967  cout<<"good genfit refit"<<endl;
968  //MakeSvtxTrack(svtxtrk,gftrk,vtx);
969  return gftrk;
970  }
971  else {
972  cout<<"No refit possible"<<endl;
973  return NULL;
974  }
975 }*/
976 
977 
978 /*FIXME this code is broken I have made zero attempt to find out why
979 * I decided not to use the SvtxTrack after they are refit
980 * may need to make the phgf_track pointer shared
981 */
982 std::shared_ptr<SvtxTrack> SVReco::MakeSvtxTrack(const SvtxTrack* svtx_track,
983  const PHGenFit::Track* phgf_track, const SvtxVertex* vertex) {
984 
985  double chi2 = phgf_track->get_chi2();
986  double ndf = phgf_track->get_ndf();
987 
988  TVector3 vertex_position(0, 0, 0);
989  TMatrixF vertex_cov(3,3);
990  double dvr2 = 0;
991  double dvz2 = 0;
992 
993  if (vertex) {
994  vertex_position.SetXYZ(vertex->get_x(), vertex->get_y(),
995  vertex->get_z());
996  dvr2 = vertex->get_error(0, 0) + vertex->get_error(1, 1);
997  dvz2 = vertex->get_error(2, 2);
998 
999  for(int i=0;i<3;i++)
1000  for(int j=0;j<3;j++)
1001  vertex_cov[i][j] = vertex->get_error(i,j);
1002  }
1003  else{
1004  cerr<<PHWHERE<<": No vertex to make SvtxTrack"<<endl;
1005  }
1006 
1007  //genfit::MeasuredStateOnPlane* gf_state_beam_line_ca = nullptr;
1008  std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_beam_line_ca = nullptr;
1009  try {
1010  gf_state_beam_line_ca = std::shared_ptr<genfit::MeasuredStateOnPlane>(phgf_track->extrapolateToLine(vertex_position,
1011  TVector3(0., 0., 1.)));
1012  } catch (...) {
1013  if (_verbosity >= 2)
1014  LogWarning("extrapolateToLine failed!");
1015  }
1016  if(!gf_state_beam_line_ca) return nullptr;
1017 
1025  double u = gf_state_beam_line_ca->getState()[3];
1026  double v = gf_state_beam_line_ca->getState()[4];
1027 
1028  double du2 = gf_state_beam_line_ca->getCov()[3][3];
1029  double dv2 = gf_state_beam_line_ca->getCov()[4][4];
1030 
1031  //delete gf_state_beam_line_ca;
1032 
1033  //const SvtxTrack_v1* temp_track = static_cast<const SvtxTrack_v1*> (svtx_track);
1034  // SvtxTrack_v1* out_track = new SvtxTrack_v1(
1035  // *static_cast<const SvtxTrack_v1*>(svtx_track));
1036  std::shared_ptr < SvtxTrack_v1 > out_track = std::shared_ptr < SvtxTrack_v1
1037  > (new SvtxTrack_v1(*static_cast<const SvtxTrack_v1*>(svtx_track)));
1038 
1039  out_track->set_dca2d(u);
1040  out_track->set_dca2d_error(sqrt(du2 + dvr2));
1041 
1042  std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state_vertex_ca = nullptr;
1043  try {
1044  gf_state_vertex_ca = std::shared_ptr < genfit::MeasuredStateOnPlane
1045  > (phgf_track->extrapolateToPoint(vertex_position));
1046  } catch (...) {
1047  if (_verbosity >= 2)
1048  LogWarning("extrapolateToPoint failed!");
1049  }
1050  if (!gf_state_vertex_ca) {
1051  //delete out_track;
1052  return nullptr;
1053  }
1054 
1055  TVector3 mom = gf_state_vertex_ca->getMom();
1056  TVector3 pos = gf_state_vertex_ca->getPos();
1057  TMatrixDSym cov = gf_state_vertex_ca->get6DCov();
1058 
1059  // genfit::MeasuredStateOnPlane* gf_state_vertex_ca =
1060  // phgf_track->extrapolateToLine(vertex_position,
1061  // TVector3(0., 0., 1.));
1062 
1063  u = gf_state_vertex_ca->getState()[3];
1064  v = gf_state_vertex_ca->getState()[4];
1065 
1066  du2 = gf_state_vertex_ca->getCov()[3][3];
1067  dv2 = gf_state_vertex_ca->getCov()[4][4];
1068 
1069 
1070  double dca3d = sqrt(u * u + v * v);
1071  double dca3d_error = sqrt(du2 + dv2 + dvr2 + dvz2);
1072 
1073  out_track->set_dca(dca3d);
1074  out_track->set_dca_error(dca3d_error);
1075 
1081  /*
1082  // Rotate from u,v,n to r: n X Z, Z': n X r, n using 5D state/cov
1083  // commented on 2017-10-09
1084  TMatrixF pos_in(3,1);
1085  TMatrixF cov_in(3,3);
1086  pos_in[0][0] = gf_state_vertex_ca->getState()[3];
1087  pos_in[1][0] = gf_state_vertex_ca->getState()[4];
1088  pos_in[2][0] = 0.;
1089  cov_in[0][0] = gf_state_vertex_ca->getCov()[3][3];
1090  cov_in[0][1] = gf_state_vertex_ca->getCov()[3][4];
1091  cov_in[0][2] = 0.;
1092  cov_in[1][0] = gf_state_vertex_ca->getCov()[4][3];
1093  cov_in[1][1] = gf_state_vertex_ca->getCov()[4][4];
1094  cov_in[1][2] = 0.;
1095  cov_in[2][0] = 0.;
1096  cov_in[2][1] = 0.;
1097  cov_in[2][2] = 0.;
1098  TMatrixF pos_out(3,1);
1099  TMatrixF cov_out(3,3);
1100  TVector3 vu = gf_state_vertex_ca->getPlane().get()->getU();
1101  TVector3 vv = gf_state_vertex_ca->getPlane().get()->getV();
1102  TVector3 vn = vu.Cross(vv);
1103  pos_cov_uvn_to_rz(vu, vv, vn, pos_in, cov_in, pos_out, cov_out);
1105  TMatrixF vertex_cov_out(3,3);
1106  get_vertex_error_uvn(vu,vv,vn, vertex_cov, vertex_cov_out);
1107  float dca3d_xy = pos_out[0][0];
1108  float dca3d_z = pos_out[1][0];
1109  float dca3d_xy_error = sqrt(cov_out[0][0] + vertex_cov_out[0][0]);
1110  float dca3d_z_error = sqrt(cov_out[1][1] + vertex_cov_out[1][1]);
1111  //Begin DEBUG
1112  // LogDebug("rotation debug---------- ");
1113  // gf_state_vertex_ca->Print();
1114  // LogDebug("dca rotation---------- ");
1115  // pos_out = pos_in;
1116  // cov_out = cov_in;
1117  // pos_in.Print();
1118  // cov_in.Print();
1119  // pos_out.Print();
1120  // cov_out.Print();
1121  // cout
1122  // <<"dca3d_xy: "<<dca3d_xy <<" +- "<<dca3d_xy_error*dca3d_xy_error
1123  // <<"; dca3d_z: "<<dca3d_z<<" +- "<< dca3d_z_error*dca3d_z_error
1124  // <<"\n";
1125  // gf_state_vertex_ca->get6DCov().Print();
1126  // LogDebug("vertex rotation---------- ");
1127  // vertex_position.Print();
1128  // vertex_cov.Print();
1129  // vertex_cov_out.Print();
1130  //End DEBUG
1131  */
1132 
1133  //
1134  // in: X, Y, Z; out; r: n X Z, Z X r, Z
1135 
1136  float dca3d_xy = NAN;
1137  float dca3d_z = NAN;
1138  float dca3d_xy_error = NAN;
1139  float dca3d_z_error = NAN;
1140 
1141  try{
1142  TMatrixF pos_in(3,1);
1143  TMatrixF cov_in(3,3);
1144  TMatrixF pos_out(3,1);
1145  TMatrixF cov_out(3,3);
1146 
1147  TVectorD state6(6); // pos(3), mom(3)
1148  TMatrixDSym cov6(6,6); //
1149 
1150  gf_state_vertex_ca->get6DStateCov(state6, cov6);
1151 
1152  TVector3 vn(state6[3], state6[4], state6[5]);
1153 
1154  // mean of two multivariate gaussians Pos - Vertex
1155  pos_in[0][0] = state6[0] - vertex_position.X();
1156  pos_in[1][0] = state6[1] - vertex_position.Y();
1157  pos_in[2][0] = state6[2] - vertex_position.Z();
1158 
1159 
1160  for(int i=0;i<3;++i){
1161  for(int j=0;j<3;++j){
1162  cov_in[i][j] = cov6[i][j] + vertex_cov[i][j];
1163  }
1164  }
1165 
1166  dca3d_xy = pos_out[0][0];
1167  dca3d_z = pos_out[2][0];
1168  dca3d_xy_error = sqrt(cov_out[0][0]);
1169  dca3d_z_error = sqrt(cov_out[2][2]);
1170 
1171 #ifdef _DEBUG_
1172  cout<<__LINE__<<": Vertex: ----------------"<<endl;
1173  vertex_position.Print();
1174  vertex_cov.Print();
1175 
1176  cout<<__LINE__<<": State: ----------------"<<endl;
1177  state6.Print();
1178  cov6.Print();
1179 
1180  cout<<__LINE__<<": Mean: ----------------"<<endl;
1181  pos_in.Print();
1182  cout<<"===>"<<endl;
1183  pos_out.Print();
1184 
1185  cout<<__LINE__<<": Cov: ----------------"<<endl;
1186  cov_in.Print();
1187  cout<<"===>"<<endl;
1188  cov_out.Print();
1189 
1190  cout<<endl;
1191 #endif
1192 
1193  } catch (...) {
1194  if (_verbosity > 0)
1195  LogWarning("DCA calculationfailed!");
1196  }
1197 
1198  out_track->set_dca3d_xy(dca3d_xy);
1199  out_track->set_dca3d_z(dca3d_z);
1200  out_track->set_dca3d_xy_error(dca3d_xy_error);
1201  out_track->set_dca3d_z_error(dca3d_z_error);
1202 
1203  //if(gf_state_vertex_ca) delete gf_state_vertex_ca;
1204 
1205  out_track->set_chisq(chi2);
1206  out_track->set_ndf(ndf);
1207  out_track->set_charge(phgf_track->get_charge());
1208 
1209  out_track->set_px(mom.Px());
1210  out_track->set_py(mom.Py());
1211  out_track->set_pz(mom.Pz());
1212 
1213  out_track->set_x(pos.X());
1214  out_track->set_y(pos.Y());
1215  out_track->set_z(pos.Z());
1216 
1217  for (int i = 0; i < 6; i++) {
1218  for (int j = i; j < 6; j++) {
1219  out_track->set_error(i, j, cov[i][j]);
1220  }
1221  }
1222 
1223  // for (SvtxTrack::ConstClusterIter iter = svtx_track->begin_clusters();
1224  // iter != svtx_track->end_clusters(); ++iter) {
1225  // unsigned int cluster_id = *iter;
1226  // SvtxCluster* cluster = _clustermap->get(cluster_id);
1227  // if (!cluster) {
1228  // LogError("No cluster Found!");
1229  // continue;
1230  // }
1231  // //cluster->identify(); //DEBUG
1232  //
1233  // //unsigned int l = cluster->get_layer();
1234  //
1235  // TVector3 pos(cluster->get_x(), cluster->get_y(), cluster->get_z());
1236  //
1237  // double radius = pos.Pt();
1238  //
1239  // std::shared_ptr<genfit::MeasuredStateOnPlane> gf_state = nullptr;
1240  // try {
1241  // gf_state = std::shared_ptr < genfit::MeasuredStateOnPlane
1242  // > (phgf_track->extrapolateToCylinder(radius,
1243  // TVector3(0, 0, 0), TVector3(0, 0, 1), 0));
1244  // } catch (...) {
1245  // if (Verbosity() >= 2)
1246  // LogWarning("Exrapolation failed!");
1247  // }
1248  // if (!gf_state) {
1249  // if (Verbosity() > 1)
1250  // LogWarning("Exrapolation failed!");
1251  // continue;
1252  // }
1253  //
1254  // //SvtxTrackState* state = new SvtxTrackState_v1(radius);
1255  // std::shared_ptr<SvtxTrackState> state = std::shared_ptr<SvtxTrackState> (new SvtxTrackState_v1(radius));
1256  // state->set_x(gf_state->getPos().x());
1257  // state->set_y(gf_state->getPos().y());
1258  // state->set_z(gf_state->getPos().z());
1259  //
1260  // state->set_px(gf_state->getMom().x());
1261  // state->set_py(gf_state->getMom().y());
1262  // state->set_pz(gf_state->getMom().z());
1263  //
1264  // //gf_state->getCov().Print();
1265  //
1266  // for (int i = 0; i < 6; i++) {
1267  // for (int j = i; j < 6; j++) {
1268  // state->set_error(i, j, gf_state->get6DCov()[i][j]);
1269  // }
1270  // }
1271  //
1272  // out_track->insert_state(state.get());
1273  //
1274  //#ifdef _DEBUG_
1275  // cout
1276  // <<__LINE__
1277  // <<": " << radius <<" => "
1278  // <<sqrt(state->get_x()*state->get_x() + state->get_y()*state->get_y())
1279  // <<endl;
1280  //#endif
1281  // }
1282 
1283 #ifdef _DEBUG_
1284  cout << __LINE__ << endl;
1285 #endif
1286 
1287  const genfit::Track *gftrack = phgf_track->getGenFitTrack();
1288  const genfit::AbsTrackRep *rep = gftrack->getCardinalRep();
1289  for(unsigned int id = 0; id< gftrack->getNumPointsWithMeasurement();++id) {
1290  genfit::TrackPoint *trpoint = gftrack->getPointWithMeasurementAndFitterInfo(id, gftrack->getCardinalRep());
1291 
1292  if(!trpoint) {
1293  if (_verbosity > 1)
1294  LogWarning("!trpoint");
1295  continue;
1296  }
1297 
1298  genfit::KalmanFitterInfo* kfi = static_cast<genfit::KalmanFitterInfo*>( trpoint->getFitterInfo(rep) );
1299  if(!kfi) {
1300  if (_verbosity > 1)
1301  LogWarning("!kfi");
1302  continue;
1303  }
1304 
1305  std::shared_ptr<const genfit::MeasuredStateOnPlane> gf_state = nullptr;
1306  try {
1307  //gf_state = std::shared_ptr <genfit::MeasuredStateOnPlane> (const_cast<genfit::MeasuredStateOnPlane*> (&(kfi->getFittedState(true))));
1308  const genfit::MeasuredStateOnPlane* temp_state = &(kfi->getFittedState(true));
1309  gf_state = std::shared_ptr <genfit::MeasuredStateOnPlane> (new genfit::MeasuredStateOnPlane(*temp_state));
1310  } catch (...) {
1311  if (_verbosity > 1)
1312  LogWarning("Exrapolation failed!");
1313  }
1314  if (!gf_state) {
1315  if (_verbosity > 1)
1316  LogWarning("Exrapolation failed!");
1317  continue;
1318  }
1319  genfit::MeasuredStateOnPlane temp;
1320  float pathlength = -phgf_track->extrapolateToPoint(temp,vertex_position,id);
1321 
1322  std::shared_ptr<SvtxTrackState> state = std::shared_ptr<SvtxTrackState> (new SvtxTrackState_v1(pathlength));
1323  state->set_x(gf_state->getPos().x());
1324  state->set_y(gf_state->getPos().y());
1325  state->set_z(gf_state->getPos().z());
1326 
1327  state->set_px(gf_state->getMom().x());
1328  state->set_py(gf_state->getMom().y());
1329  state->set_pz(gf_state->getMom().z());
1330 
1331  //gf_state->getCov().Print();
1332 
1333  for (int i = 0; i < 6; i++) {
1334  for (int j = i; j < 6; j++) {
1335  state->set_error(i, j, gf_state->get6DCov()[i][j]);
1336  }
1337  }
1338 
1339  out_track->insert_state(state.get());
1340 
1341 #ifdef _DEBUG_
1342  cout
1343  <<__LINE__
1344  <<": " << id
1345  <<": " << pathlength <<" => "
1346  <<sqrt(state->get_x()*state->get_x() + state->get_y()*state->get_y())
1347  <<endl;
1348 #endif
1349  }
1350 
1351  return out_track;
1352 }