EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test_with_vertex.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file test_with_vertex.cpp
1 #include <iostream>
2 #include "TFile.h"
3 #include "TTree.h"
4 #include "HelixHough.h"
5 #include "HelixRange.h"
6 #include "HelixResolution.h"
7 #include "VertexFitFunc.h"
8 #include "NewtonMinimizerGradHessian.h"
9 #include <sys/time.h>
10 #include <Eigen/LU>
11 #include <Eigen/Core>
12 #include <math.h>
13 #include "FourHitSeedFinder/FourHitSeedFinder.h"
14 
15 using namespace std;
16 using namespace FitNewton;
17 using namespace Eigen;
18 
19 int main(int argc, char** argv)
20 {
21  bool print_truth = false;
22  bool print_reco = false;
23 
24  TFile infile(argv[1]);
25  TTree* etree=0;
26  TTree* ttree=0;
27  infile.GetObject("events", etree);
28  etree->SetBranchAddress("tracklist", &ttree);
29  unsigned int track=0;
30  unsigned int nhits=0;
31  int layer[12];
32  double x_hits[12];
33  double y_hits[12];
34  double z_hits[12];
35  double trk_kappa=0.;
36  double trk_d=0.;
37  double trk_phi=0.;
38  double trk_dzdl=0.;
39  double trk_z0=0.;
40  unsigned char charge=0;
41 
42 
43  int nlayers = 4;
44  vector<float> radii;
45  radii.assign(nlayers,0.);
46  radii[0]=2.5;
47  radii[1]=5.0;
48  radii[2]=10.0;
49  radii[3]=14.0;
50  vector<float> smear_xy_layer;smear_xy_layer.assign(nlayers,0);
51  vector<float> smear_z_layer;smear_z_layer.assign(nlayers,0);
52  float sqrt_12 = sqrt(12.);
53  smear_xy_layer[0] = (50.0e-4/sqrt_12);
54  smear_z_layer[0] = (425.0e-4/sqrt_12);
55  smear_xy_layer[1] = (50.0e-4/sqrt_12);
56  smear_z_layer[1] = (425.0e-4/sqrt_12);
57  smear_xy_layer[2] = (80.0e-4/sqrt_12);
58  smear_z_layer[2] = (1000.0e-4/sqrt_12);
59  smear_xy_layer[3] = (80.0e-4/sqrt_12);
60  smear_z_layer[3] = (1000.0e-4/sqrt_12);
61 
62  //phi,d,kappa,dzdl,z0
63  HelixResolution min_res(0.01, 0.5, 0.002, 0.01, 2.);
64  HelixResolution max_res(0.001, 0.5, 0.001, 0.01, 2.);
65  HelixRange top_range(0., 1.*M_PI, -0.0025, 0.0025, 0., 0.03, -0.9, 0.9, -0.005, 0.005);
66  FourHitSeedFinder tracker(radii, 4, 1, 2, 4, 1, min_res, max_res, top_range);
67  tracker.setLayerResolution(smear_xy_layer, smear_z_layer);
68  tracker.setVertexResolution(0.005, 0.01);
69  tracker.setUsingVertex(true);
70  tracker.setChi2Cut(3.0);
71  tracker.setPrintTimings(true);
72  unsigned int max_hits = 5;
73 
74 
75 
76 
77 
78  for(unsigned int ev=0;ev<etree->GetEntries();ev++)
79  {
80  etree->GetEntry(ev);
81  ttree->SetBranchAddress("track", &track);
82  ttree->SetBranchAddress("nhits", &nhits);
83  ttree->SetBranchAddress("x_hits", &x_hits);
84  ttree->SetBranchAddress("y_hits", &y_hits);
85  ttree->SetBranchAddress("z_hits", &z_hits);
86  ttree->SetBranchAddress("layer", &layer);
87  ttree->SetBranchAddress("trk_kappa", &trk_kappa);
88  ttree->SetBranchAddress("trk_d", &trk_d);
89  ttree->SetBranchAddress("trk_phi", &trk_phi);
90  ttree->SetBranchAddress("trk_dzdl", &trk_dzdl);
91  ttree->SetBranchAddress("trk_z0", &trk_z0);
92  ttree->SetBranchAddress("charge", &charge);
93 
94  cout<<"event "<<ev<<":"<<endl<<endl;
95 
96  vector<SimpleHit3D> hits; // 3d space point + hit index
97  vector<SimpleTrack3D> tracks; // list of space points _and_ helix parameters
98  unsigned int index=0;
99  cout<<"total MC tracks = "<<ttree->GetEntries()<<endl;
100  cout<<endl;
101 
102  for(unsigned int trk=0;trk<ttree->GetEntries();trk++)
103  {
104  ttree->GetEntry(trk);
105  if(print_truth==true)
106  {
107  cout<<"truth track : "<<endl;
108  cout<<"phi = "<<trk_phi<<endl;
109  cout<<"d = "<<trk_d<<endl;
110  cout<<"kappa = "<<trk_kappa<<endl;
111  cout<<"dzdl = "<<trk_dzdl<<endl;
112  cout<<"z0 = "<<trk_z0<<endl;
113  cout<<endl;
114  }
115  for(unsigned int hit=0;hit<nhits;hit++)
116  {
117  float phi = atan2(y_hits[hit], x_hits[hit]);
118  float xy_error = smear_xy_layer[layer[hit]]*sqrt_12*0.5;
119  float x_error = fabs(xy_error*sin(phi));
120  float y_error = fabs(xy_error*cos(phi));
121  float z_error = smear_z_layer[layer[hit]]*sqrt_12*0.5;
122 
123  hits.push_back(SimpleHit3D(x_hits[hit],x_error, y_hits[hit],y_error, z_hits[hit],z_error, index, layer[hit]));
124  index++;
125  }
126  }
127 
128  double z_avg = 0.;
129  double z_weight = 0.;
130  for(unsigned int ht=0;ht<hits.size();++ht)
131  {
132  double temp = 1./(hits[ht].dz);
133  z_avg += (hits[ht].z)*temp;
134  z_weight += temp;
135  }
136  z_avg/=z_weight;
137  cout<<"z average = "<<z_avg<<endl;
138 
139  HelixResolution min_res_init(0.01, 10.5, 0.01, 0.02, 0.01);
140  HelixResolution max_res_init(0.001, 0.5, 0.01, 0.02, 0.01);
141  HelixRange top_range_init(0., 1.*M_PI, -0.1, 0.1, 0., 0.005, -0.9, 0.9, -0.1 + z_avg, 0.1 + z_avg);
142  FourHitSeedFinder tracker_init(radii, 2, 1, 2, 4, 2, min_res_init, max_res_init, top_range_init);
143  tracker_init.setLayerResolution(smear_xy_layer, smear_z_layer);
144  tracker_init.setVertexResolution(0.001, 0.005);
145  tracker_init.setUsingVertex(false);
146  tracker_init.setChi2Cut(3.0);
147  tracker_init.setPrintTimings(true);
148  unsigned int max_hits_init = 8;
149  unsigned int maxtracks = 160;
150 
151 
152  maxtracks = 15*(hits.size())/1000;
153  if(maxtracks < 40){maxtracks=40;}
154 
155  timeval t1,t2;
156  double time1=0.;
157  double time2=0.;
158  gettimeofday(&t1, NULL);
159  tracker_init.findHelices(hits, 4, max_hits_init, tracks, maxtracks);
160  gettimeofday(&t2, NULL);
161  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec)/1000000.);
162  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec)/1000000.);
163  cout<<"initial tracking time = "<<(time2 - time1)<<endl;
164 
165  unsigned int ngood = 0;
166  unsigned int n_trk = tracks.size();
167  for(unsigned int trk=0;trk<n_trk;++trk)
168  {
169  unsigned int high = 0;
170  unsigned int low = (1<<31);
171  unsigned int n_hit = tracks[trk].hits.size();
172  for(unsigned int ht=0;ht<n_hit;++ht)
173  {
174  tracks[trk].hits[ht].index > high ? high = tracks[trk].hits[ht].index : high = high;
175  tracks[trk].hits[ht].index < low ? low = tracks[trk].hits[ht].index : low = low;
176  }
177  if( (high - low) <= 3 ){ngood++;}
178  }
179 
180  cout<<"tracks found before vertex finding = "<<tracks.size()<<endl;
181  cout<<"good tracks found before vertex finding = "<<ngood<<endl;
182 
183  VertexFitFunc vfit;
184  vfit.setTracks(&tracks);
185 
186  // very general minimization routine (good for other things too)
187  NewtonMinimizerGradHessian minimizer;
188  minimizer.setFunction(&vfit);
189 
190  VectorXd start_point = VectorXd::Zero(3); // input initial guess
191  VectorXd min_point = VectorXd::Zero(3); // output vector for minimize method below
192 
193  //first use a large sigma for expo-dca^2 calculation (same in x,y,z)
194  vfit.setFixedPar(0, 3.0);
195  // -> guess
196  // <- best fit
197  // -> relative tolerance
198  // -> max number of iterations
199  // -> absolute tolerance
200  // return true if successful
201  minimizer.minimize(start_point, min_point, 1.0e-9, 16, 1.0e-15);
202  //then continue with a smaller sigma
203  vfit.setFixedPar(0, 0.3);
204  start_point = min_point;
205  minimizer.minimize(start_point, min_point, 1.0e-9, 16, 1.0e-15);
206  vfit.setFixedPar(0, 0.03);
207  start_point = min_point;
208  minimizer.minimize(start_point, min_point, 1.0e-9, 16, 1.0e-15);
209 
210 
211  // store output vertex spatial point
212  vector<double> vertex;
213  vertex.assign(3,0.);
214  vertex[0] = min_point(0);
215  vertex[1] = min_point(1);
216  vertex[2] = min_point(2);
217 
218 // vertex[0] = 0.;
219 // vertex[1] = 0.;
220 // vertex[2] = 0.;
221 
222  cout << "Found Event vertex = " << vertex[0] << " " << vertex[1] << " " << vertex[2] << endl << endl;
223 
224  tracks.clear();
225 
226 
227  for(unsigned int ht=0;ht<hits.size();++ht)
228  {
229  hits[ht].x -= vertex[0];
230  hits[ht].y -= vertex[1];
231  hits[ht].z -= vertex[2];
232  }
233 
234  gettimeofday(&t1, NULL);
235  // tracker.findHelices(hits, 4, max_hits, tracks);
236  tracker.findHelices(hits, 4, max_hits, tracks);
237  gettimeofday(&t2, NULL);
238  time1 = ((double)(t1.tv_sec) + (double)(t1.tv_usec)/1000000.);
239  time2 = ((double)(t2.tv_sec) + (double)(t2.tv_usec)/1000000.);
240  cout<<"primary tracking time = "<<(time2 - time1)<<endl;
241 
242  for(unsigned int trk=0;trk<tracks.size();++trk)
243  {
244  for(unsigned int ht=0;ht<tracks[trk].hits.size();++ht)
245  {
246  tracks[trk].hits[ht].x += vertex[0];
247  tracks[trk].hits[ht].y += vertex[1];
248  tracks[trk].hits[ht].z += vertex[2];
249  }
250  }
251 
252 
253 
254  ngood = 0;
255  n_trk = tracks.size();
256  for(unsigned int trk=0;trk<n_trk;++trk)
257  {
258  unsigned int high = 0;
259  unsigned int low = (1<<31);
260  unsigned int n_hit = tracks[trk].hits.size();
261  for(unsigned int ht=0;ht<n_hit;++ht)
262  {
263  tracks[trk].hits[ht].index > high ? high = tracks[trk].hits[ht].index : high = high;
264  tracks[trk].hits[ht].index < low ? low = tracks[trk].hits[ht].index : low = low;
265  }
266  if( (high - low) <= 3 ){ngood++;}
267  }
268 
269  if(print_reco==true)
270  {
271  for(unsigned int trk=0;trk<n_trk;++trk)
272  {
273  cout<<"track "<<trk<<" : "<<endl;
274  unsigned int n_hit = tracks[trk].hits.size();
275  for(unsigned int ht=0;ht<n_hit;++ht)
276  {
277  cout<<"hit "<<ht<<" : ";
278  cout<<tracks[trk].hits[ht].x<<" ";
279  cout<<tracks[trk].hits[ht].y<<" ";
280  cout<<tracks[trk].hits[ht].z<<" ";
281  cout<<tracks[trk].hits[ht].index<<endl;
282  }
283  cout<<"phi = "<<tracks[trk].phi<<endl;
284  cout<<"d = "<<tracks[trk].d<<endl;
285  cout<<"kappa = "<<tracks[trk].kappa<<endl;
286  cout<<"dzdl = "<<tracks[trk].dzdl<<endl;
287  cout<<"z0 = "<<tracks[trk].z0<<endl;
288  cout<<endl;
289  }
290  }
291 
292 
293  cout<<n_trk<<" tracks found"<<endl;
294  cout<<ngood<<" good tracks found"<<endl;
295  cout<<endl<<endl;
296  }
297 
298 
299  return 0;
300 
301 }