Class Reference for E1039 Core & Analysis Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KalmanDSTrk.cxx
Go to the documentation of this file.
1 
10 //#include "KalmanFitter.h"
11 #include "KalmanDSTrk.h"
12 #include "TriggerRoad.h"
13 #include "PatternDBUtil.h"
14 
15 
16 #include <phfield/PHField.h>
17 #include <phool/PHTimer.h>
18 #include <fun4all/Fun4AllBase.h>
19 
20 #include <iostream>
21 #include <fstream>
22 #include <algorithm>
23 #include <cmath>
24 #include <limits>
25 
26 #include <TCanvas.h>
27 #include <TGraphErrors.h>
28 #include <TBox.h>
29 #include <TMatrixD.h>
30 #include <TFile.h>
31 #include <TNtuple.h>
32 
33 //#define _DEBUG_ON
34 #define _DEBUG_YUHW_
35 
36 KalmanDSTrk::KalmanDSTrk(
37  const PHField* field,
38  const TGeoManager *geom,
39  bool enable_KF,
40  int DS_level,
41  const std::string sim_db_name,
42  const std::string pattern_db_name
43  )
44 : verbosity(10)
45 ,_enable_KF(enable_KF)
46 ,_DS_level(DS_level)
47 ,_sim_db_name(sim_db_name)
48 ,_pattern_db_name(pattern_db_name)
49 , _pattern_db(nullptr)
50 {
51  using namespace std;
52 
53 #ifdef _DEBUG_ON
54  cout << "Initialization of KalmanDSTrk ..." << endl;
55  cout << "========================================" << endl;
56 #endif
57 
58  _timers.insert(std::make_pair<std::string, PHTimer*>("build_db", new PHTimer("build_db")));
59  _timers.insert(std::make_pair<std::string, PHTimer*>("load_db", new PHTimer("load_db")));
60 
61  _timers.insert(std::make_pair<std::string, PHTimer*>("event", new PHTimer("event")));
62  _timers.insert(std::make_pair<std::string, PHTimer*>("st2", new PHTimer("st2")));
63  _timers.insert(std::make_pair<std::string, PHTimer*>("search_db_2", new PHTimer("search_db_2")));
64  _timers.insert(std::make_pair<std::string, PHTimer*>("st3", new PHTimer("st3")));
65 
66  _timers.insert(std::make_pair<std::string, PHTimer*>("st23", new PHTimer("st23")));
67  _timers.insert(std::make_pair<std::string, PHTimer*>("search_db_23", new PHTimer("search_db_23")));
68  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_fit1", new PHTimer("st23_fit1")));
69  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_prop", new PHTimer("st23_prop")));
70  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_fit2", new PHTimer("st23_fit2")));
71  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_hodo", new PHTimer("st23_hodo")));
72  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_lr40", new PHTimer("st23_lr40")));
73  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_lr150", new PHTimer("st23_lr150")));
74  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_rm_hits", new PHTimer("st23_rm_hits")));
75  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_acc", new PHTimer("st23_acc")));
76  _timers.insert(std::make_pair<std::string, PHTimer*>("st23_reduce", new PHTimer("st23_reduce")));
77 
78  _timers.insert(std::make_pair<std::string, PHTimer*>("global", new PHTimer("global")));
79  _timers.insert(std::make_pair<std::string, PHTimer*>("search_db_glb", new PHTimer("search_db_glb")));
80  _timers.insert(std::make_pair<std::string, PHTimer*>("global_st1", new PHTimer("global_st1")));
81  _timers.insert(std::make_pair<std::string, PHTimer*>("global_link", new PHTimer("global_link")));
82  _timers.insert(std::make_pair<std::string, PHTimer*>("global_kalman", new PHTimer("global_kalman")));
83 
84  _timers.insert(std::make_pair<std::string, PHTimer*>("kalman", new PHTimer("kalman")));
85 
86  //Initialize jobOpts service
87  p_jobOptsSvc = JobOptsSvc::instance();
88 
89  //Initialize Kalman fitter
90  if(_enable_KF)
91  {
92  kmfitter = new KalmanFitter(field, geom);
93  kmfitter->setControlParameter(50, 0.001);
94  }
95 
96  //Load Patter DB from file
97 
98 
99  if(_DS_level > KalmanDSTrk::NO_DS) {
100 
101  if(_pattern_db_name != "") {
102  std::cout <<"KalmanDSTrk::KalmanDSTrk: load DB from pattern db: "<< _pattern_db_name << std::endl;
103  _timers["load_db"]->restart();
104  _pattern_db = PatternDBUtil::LoadPatternDB(_pattern_db_name);
105  _timers["load_db"]->stop();
106  } else if(_sim_db_name != "") {
107  std::cout <<"KalmanDSTrk::KalmanDSTrk: load DB from sim db: "<< _sim_db_name << std::endl;
108  _timers["build_db"]->restart();
109  _pattern_db = new PatternDB();
110  PatternDBUtil::BuildPatternDB(_sim_db_name, "PatternDB_tmp.root", *_pattern_db);
111  _timers["build_db"]->stop();
112  } else {
113  std::cout <<"KalmanDSTrk::KalmanDSTrk: no sim or pattern DB" << std::endl;
114  }
115 
116  if(_pattern_db) {
117  std::cout <<"KalmanDSTrk::KalmanDSTrk: DB loaded. St23 size: "<< _pattern_db->St23.size() << std::endl;
118  std::cout << "================================================================" << std::endl;
119  std::cout << "Build DB "<<_timers["build_db"]->get_accumulated_time()/1000. << " sec" <<std::endl;
120  std::cout << "Load DB "<<_timers["load_db"]->get_accumulated_time()/1000. << " sec" <<std::endl;
121  std::cout << "================================================================" << std::endl;
122  } else {
123  std::cout <<"KalmanDSTrk::KalmanDSTrk: DB NOT loaded - _DS_level set to 0 " << std::endl;
124  _DS_level = 0;
125  }
126  }
127 
128  //Initialize minuit minimizer
129  minimizer[0] = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Simplex");
130  minimizer[1] = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
131  fcn = ROOT::Math::Functor(&tracklet_curr, &Tracklet::Eval, p_jobOptsSvc->m_enableKMag ? 5 : 4);
132 
133  for(int i = 0; i < 2; ++i)
134  {
135  minimizer[i]->SetMaxFunctionCalls(1000000);
136  minimizer[i]->SetMaxIterations(100);
137  minimizer[i]->SetTolerance(1E-2);
138  minimizer[i]->SetFunction(fcn);
139  minimizer[i]->SetPrintLevel(0);
140  }
141 
142  //Minimize ROOT output
143  extern Int_t gErrorIgnoreLevel;
144  gErrorIgnoreLevel = 9999;
145 
146  //Initialize geometry service
147  p_geomSvc = GeomSvc::instance();
148 #ifdef _DEBUG_ON
149  p_geomSvc->printTable();
150  p_geomSvc->printWirePosition();
151  p_geomSvc->printAlignPar();
152 #endif
153 
154  //Initialize plane angles for all planes
155  for(int i = 1; i <= nChamberPlanes; ++i)
156  {
157  costheta_plane[i] = p_geomSvc->getCostheta(i);
158  sintheta_plane[i] = p_geomSvc->getSintheta(i);
159  }
160 
161  //Initialize hodoscope IDs
162  detectorIDs_mask[0] = p_geomSvc->getDetectorIDs("H1");
163  detectorIDs_mask[1] = p_geomSvc->getDetectorIDs("H2");
164  detectorIDs_mask[2] = p_geomSvc->getDetectorIDs("H3");
165  detectorIDs_mask[3] = p_geomSvc->getDetectorIDs("H4");
166  detectorIDs_maskX[0] = p_geomSvc->getDetectorIDs("H1[TB]");
167  detectorIDs_maskX[1] = p_geomSvc->getDetectorIDs("H2[TB]");
168  detectorIDs_maskX[2] = p_geomSvc->getDetectorIDs("H3[TB]");
169  detectorIDs_maskX[3] = p_geomSvc->getDetectorIDs("H4[TB]");
170  detectorIDs_maskY[0] = p_geomSvc->getDetectorIDs("H1[LR]");
171  detectorIDs_maskY[1] = p_geomSvc->getDetectorIDs("H2[LR]");
172  detectorIDs_maskY[2] = p_geomSvc->getDetectorIDs("H4Y1[LR]");
173  detectorIDs_maskY[3] = p_geomSvc->getDetectorIDs("H4Y2[LR]");
174  detectorIDs_muidHodoAid[0] = p_geomSvc->getDetectorIDs("H4[TB]");
175  detectorIDs_muidHodoAid[1] = p_geomSvc->getDetectorIDs("H4Y");
176 
177  //Register masking stations for tracklets in station-0/1, 2, 3+/-
178  stationIDs_mask[0].push_back(1);
179  stationIDs_mask[1].push_back(1);
180  stationIDs_mask[2].push_back(2);
181  stationIDs_mask[3].push_back(3);
182  stationIDs_mask[4].push_back(3);
183 
184  //Masking stations for back partial
185  stationIDs_mask[5].push_back(2);
186  stationIDs_mask[5].push_back(3);
187  stationIDs_mask[5].push_back(4);
188 
189  //Masking stations for global track
190  stationIDs_mask[6].push_back(1);
191  stationIDs_mask[6].push_back(2);
192  stationIDs_mask[6].push_back(3);
193  stationIDs_mask[6].push_back(4);
194 
195  //prop. tube IDs for mu id
196  detectorIDs_muid[0][0] = p_geomSvc->getDetectorID("P1X1");
197  detectorIDs_muid[0][1] = p_geomSvc->getDetectorID("P1X2");
198  detectorIDs_muid[0][2] = p_geomSvc->getDetectorID("P2X1");
199  detectorIDs_muid[0][3] = p_geomSvc->getDetectorID("P2X2");
200  detectorIDs_muid[1][0] = p_geomSvc->getDetectorID("P1Y1");
201  detectorIDs_muid[1][1] = p_geomSvc->getDetectorID("P1Y2");
202  detectorIDs_muid[1][2] = p_geomSvc->getDetectorID("P2Y1");
203  detectorIDs_muid[1][3] = p_geomSvc->getDetectorID("P2Y2");
204 
205  //Reference z_ref for mu id
206  z_ref_muid[0][0] = MUID_Z_REF;
207  z_ref_muid[0][1] = MUID_Z_REF;
208  z_ref_muid[0][2] = 0.5*(p_geomSvc->getPlanePosition(detectorIDs_muid[0][0]) + p_geomSvc->getPlanePosition(detectorIDs_muid[0][1]));
209  z_ref_muid[0][3] = z_ref_muid[0][2];
210 
211  z_ref_muid[1][0] = MUID_Z_REF;
212  z_ref_muid[1][1] = MUID_Z_REF;
213  z_ref_muid[1][2] = 0.5*(p_geomSvc->getPlanePosition(detectorIDs_muid[1][0]) + p_geomSvc->getPlanePosition(detectorIDs_muid[1][1]));
214  z_ref_muid[1][3] = z_ref_muid[1][2];
215 
216  //Initialize masking window sizes, with optimized contingency
217  for(int i = nChamberPlanes+1; i <= nChamberPlanes+nHodoPlanes+nPropPlanes; i++)
218  {
219  double factor = 0.;
220  if(i > nChamberPlanes && i <= nChamberPlanes+4) factor = 0.25; //for station-1 hodo
221  if(i > nChamberPlanes+4 && i <= nChamberPlanes+8) factor = 0.2; //for station-2 hodo
222  if(i > nChamberPlanes+8 && i <= nChamberPlanes+10) factor = 0.15; //for station-3 hodo
223  if(i > nChamberPlanes+10 && i <= nChamberPlanes+nHodoPlanes) factor = 0.; //for station-4 hodo
224  if(i > nChamberPlanes+nHodoPlanes) factor = 0.15; //for station-4 proptube
225 
226  z_mask[i-nChamberPlanes-1] = p_geomSvc->getPlanePosition(i);
227  for(int j = 1; j <= p_geomSvc->getPlaneNElements(i); j++)
228  {
229  double x_min, x_max, y_min, y_max;
230  p_geomSvc->get2DBoxSize(i, j, x_min, x_max, y_min, y_max);
231 
232  x_min -= (factor*(x_max - x_min));
233  x_max += (factor*(x_max - x_min));
234  y_min -= (factor*(y_max - y_min));
235  y_max += (factor*(y_max - y_min));
236 
237  x_mask_min[i-nChamberPlanes-1][j-1] = x_min;
238  x_mask_max[i-nChamberPlanes-1][j-1] = x_max;
239  y_mask_min[i-nChamberPlanes-1][j-1] = y_min;
240  y_mask_max[i-nChamberPlanes-1][j-1] = y_max;
241  }
242  }
243 
244 #ifdef _DEBUG_ON
245  cout << "========================" << endl;
246  cout << "Hodo. masking settings: " << endl;
247  for(int i = 0; i < 4; i++)
248  {
249  cout << "For station " << i+1 << endl;
250  for(std::vector<int>::iterator iter = detectorIDs_mask[i].begin(); iter != detectorIDs_mask[i].end(); ++iter) cout << "All: " << *iter << endl;
251  for(std::vector<int>::iterator iter = detectorIDs_maskX[i].begin(); iter != detectorIDs_maskX[i].end(); ++iter) cout << "X: " << *iter << endl;
252  for(std::vector<int>::iterator iter = detectorIDs_maskY[i].begin(); iter != detectorIDs_maskY[i].end(); ++iter) cout << "Y: " << *iter << endl;
253  }
254 
255  for(int i = 0; i < nStations; ++i)
256  {
257  std::cout << "Masking stations for tracklets with stationID = " << i + 1 << ": " << std::endl;
258  for(std::vector<int>::iterator iter = stationIDs_mask[i].begin(); iter != stationIDs_mask[i].end(); ++iter)
259  {
260  std::cout << *iter << " ";
261  }
262  std::cout << std::endl;
263  }
264 #endif
265 
266  //Initialize super stationIDs
267  for(int i = 0; i < nChamberPlanes/6+2; i++) superIDs[i].clear();
268  superIDs[0].push_back((p_geomSvc->getDetectorIDs("D0X")[0] + 1)/2);
269  superIDs[0].push_back((p_geomSvc->getDetectorIDs("D0U")[0] + 1)/2);
270  superIDs[0].push_back((p_geomSvc->getDetectorIDs("D0V")[0] + 1)/2);
271  superIDs[1].push_back((p_geomSvc->getDetectorIDs("D1X")[0] + 1)/2);
272  superIDs[1].push_back((p_geomSvc->getDetectorIDs("D1U")[0] + 1)/2);
273  superIDs[1].push_back((p_geomSvc->getDetectorIDs("D1V")[0] + 1)/2);
274  superIDs[2].push_back((p_geomSvc->getDetectorIDs("D2X")[0] + 1)/2);
275  superIDs[2].push_back((p_geomSvc->getDetectorIDs("D2U")[0] + 1)/2);
276  superIDs[2].push_back((p_geomSvc->getDetectorIDs("D2V")[0] + 1)/2);
277  superIDs[3].push_back((p_geomSvc->getDetectorIDs("D3pX")[0] + 1)/2);
278  superIDs[3].push_back((p_geomSvc->getDetectorIDs("D3pU")[0] + 1)/2);
279  superIDs[3].push_back((p_geomSvc->getDetectorIDs("D3pV")[0] + 1)/2);
280  superIDs[4].push_back((p_geomSvc->getDetectorIDs("D3mX")[0] + 1)/2);
281  superIDs[4].push_back((p_geomSvc->getDetectorIDs("D3mU")[0] + 1)/2);
282  superIDs[4].push_back((p_geomSvc->getDetectorIDs("D3mV")[0] + 1)/2);
283 
284  superIDs[5].push_back((p_geomSvc->getDetectorIDs("P1X")[0] + 1)/2);
285  superIDs[5].push_back((p_geomSvc->getDetectorIDs("P2X")[0] + 1)/2);
286  superIDs[6].push_back((p_geomSvc->getDetectorIDs("P1Y")[0] + 1)/2);
287  superIDs[6].push_back((p_geomSvc->getDetectorIDs("P2Y")[0] + 1)/2);
288 
289 #ifdef _DEBUG_ON
290  cout << "=============" << endl;
291  cout << "Chamber IDs: " << endl;
292  TString stereoNames[3] = {"X", "U", "V"};
293  for(int i = 0; i < nChamberPlanes/6; i++)
294  {
295  for(int j = 0; j < 3; j++) cout << i << " " << stereoNames[j].Data() << ": " << superIDs[i][j] << endl;
296  }
297 
298  cout << "Proptube IDs: " << endl;
299  for(int i = nChamberPlanes/6; i < nChamberPlanes/6+2; i++)
300  {
301  for(int j = 0; j < 2; j++) cout << i << " " << j << ": " << superIDs[i][j] << endl;
302  }
303 
304  //Initialize widow sizes for X-U matching and z positions of all chambers
305  cout << "======================" << endl;
306  cout << "U plane window sizes: " << endl;
307 #endif
308 
309  double u_factor[] = {5., 5., 5., 15., 15.};
310  for(int i = 0; i < nChamberPlanes/6; i++)
311  {
312  int xID = 2*superIDs[i][0] - 1;
313  int uID = 2*superIDs[i][1] - 1;
314  int vID = 2*superIDs[i][2] - 1;
315  double spacing = p_geomSvc->getPlaneSpacing(uID);
316  double x_span = p_geomSvc->getPlaneScaleY(uID);
317 
318  z_plane_x[i] = 0.5*(p_geomSvc->getPlanePosition(xID) + p_geomSvc->getPlanePosition(xID+1));
319  z_plane_u[i] = 0.5*(p_geomSvc->getPlanePosition(uID) + p_geomSvc->getPlanePosition(uID+1));
320  z_plane_v[i] = 0.5*(p_geomSvc->getPlanePosition(vID) + p_geomSvc->getPlanePosition(vID+1));
321 
322  u_costheta[i] = costheta_plane[uID];
323  u_sintheta[i] = sintheta_plane[uID];
324 
325  //u_win[i] = fabs(0.5*x_span/(spacing/sintheta_plane[uID])) + 2.*spacing + u_factor[i];
326  u_win[i] = fabs(0.5*x_span*sintheta_plane[uID]) + TX_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_costheta[i]) + TY_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_sintheta[i]) + 2.*spacing + u_factor[i];
327 
328 #ifdef _DEBUG_ON
329  cout << "Station " << i << ": " << xID << " " << uID << " " << vID << " " << u_win[i] << endl;
330 
331  cout
332  << "x_span = " << x_span << endl
333  << "sintheta_plane = " << sintheta_plane[uID] << endl
334  << "u_sintheta = " << u_sintheta[i] << endl
335  << "z_plane_u[i] - z_plane_x[i] = " << z_plane_u[i] - z_plane_x[i] << endl
336  << "fabs(0.5*x_span*sintheta_plane[uID]) = " << fabs(0.5*x_span*sintheta_plane[uID]) << endl
337  << "TX_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_costheta[i]) = " << TX_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_costheta[i]) << endl
338  << "TY_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_sintheta[i]) = " << TY_MAX*fabs((z_plane_u[i] - z_plane_x[i])*u_sintheta[i]) << endl
339  << "2.*spacing = " << 2.*spacing << endl
340  << "u_factor[i] = " << u_factor[i] << endl;
341 #endif
342  }
343 
344  //Initialize Z positions and maximum parameters of all planes
345  for(int i = 1; i <= nChamberPlanes; i++)
346  {
347  z_plane[i] = p_geomSvc->getPlanePosition(i);
348  slope_max[i] = costheta_plane[i]*TX_MAX + sintheta_plane[i]*TY_MAX;
349  intersection_max[i] = costheta_plane[i]*X0_MAX + sintheta_plane[i]*Y0_MAX;
350 
351 #ifdef COARSE_MODE
352  resol_plane[i] = 3.*p_geomSvc->getPlaneSpacing(i)/sqrt(12.);
353 #else
354  if(i <= 6)
355  {
356  resol_plane[i] = p_jobOptsSvc->m_st0_reject;
357  }
358  else if(i <= 12)
359  {
360  resol_plane[i] = p_jobOptsSvc->m_st1_reject;
361  }
362  else if(i <= 18)
363  {
364  resol_plane[i] = p_jobOptsSvc->m_st2_reject;
365  }
366  else if(i <= 24)
367  {
368  resol_plane[i] = p_jobOptsSvc->m_st3p_reject;
369  }
370  else
371  {
372  resol_plane[i] = p_jobOptsSvc->m_st3m_reject;
373  }
374 #endif
375  spacing_plane[i] = p_geomSvc->getPlaneSpacing(i);
376  }
377 
378 #ifdef _DEBUG_ON
379  cout << "======================================" << endl;
380  cout << "Maximum local slope and intersection: " << endl;
381 #endif
382  for(int i = 1; i <= nChamberPlanes/2; i++)
383  {
384  double d_slope = (p_geomSvc->getPlaneResolution(2*i - 1) + p_geomSvc->getPlaneResolution(2*i))/(z_plane[2*i] - z_plane[2*i-1]);
385  double d_intersection = d_slope*z_plane[2*i];
386 
387  slope_max[2*i-1] += d_slope;
388  intersection_max[2*i-1] += d_intersection;
389  slope_max[2*i] += d_slope;
390  intersection_max[2*i] += d_intersection;
391 
392 #ifdef _DEBUG_ON
393  cout << "Super plane " << i << ": " << slope_max[2*i-1] << " " << intersection_max[2*i-1] << endl;
394 #endif
395  }
396 
397  //Initialize sagitta ratios, index 0, 1, 2 are for X, U, V, this is the incrementing order of plane type
398  s_detectorID[0] = p_geomSvc->getDetectorID("D2X");
399  s_detectorID[1] = p_geomSvc->getDetectorID("D2Up");
400  s_detectorID[2] = p_geomSvc->getDetectorID("D2Vp");
401 
402  _ana_mode = true;
403  if(_ana_mode) {
404  _fana = TFile::Open("ktracker_ana.root","recreate");
405  _tana_Event = new TNtuple("Event", "Event",
406  "NSt2:NSt3:NSt23:NSt1:NGlobal:NKalman:"
407  "TSt2:TSt3:TSt23:TGKalman:TGlobal:TKalman:"
408  "TEvent"
409  );
410  _tana_St1 = new TNtuple("St1", "St1", "in:DS:out");
411  _tana_St2 = new TNtuple("St2", "St2", "in:DS:out");
412  _tana_St3 = new TNtuple("St3", "St3", "in:DS:out");
413  _tana_St23 = new TNtuple("St23", "St23", "in:DS:out:time");
414  _tana_St123 = new TNtuple("St123","St123", "in:DS:out:time");
415  }
416 
417  _event = 0;
418 }
419 
420 KalmanDSTrk::~KalmanDSTrk()
421 {
422  if(_enable_KF) delete kmfitter;
423  delete minimizer[0];
424  delete minimizer[1];
425 
426  if(_ana_mode) {
427  _fana->cd();
428  _tana_Event->Write();
429  _tana_St1->Write();
430  _tana_St2->Write();
431  _tana_St3->Write();
432  _tana_St23->Write();
433  _tana_St123->Write();
434  _fana->Close();
435  }
436 }
437 
438 void KalmanDSTrk::setRawEventDebug(SRawEvent* event_input)
439 {
440  rawEvent = event_input;
441  hitAll = event_input->getAllHits();
442 }
443 
444 
445 int KalmanDSTrk::setRawEvent(SRawEvent* event_input)
446 {
447  // reset all timers
448  for(auto iter=_timers.begin(); iter!=_timers.end();++iter) {
449  iter->second->reset();
450  }
451 
452  _timers["event"]->restart();
453  int ret = this->setRawEventWorker(event_input);
454  _timers["event"]->stop();
455 
456 
457  if(_ana_mode) {
458  float data[] = {
459  static_cast<float>(trackletsInSt[1].size()), // St2
460  static_cast<float>(trackletsInSt[2].size()), // St3
461  static_cast<float>(trackletsInSt[3].size()), // St23
462  static_cast<float>(trackletsInSt[0].size()), // St1
463  static_cast<float>(trackletsInSt[4].size()), // Global
464  static_cast<float>(stracks.size()),
465  static_cast<float>(_timers["st2"]->get_accumulated_time()/1000.),
466  static_cast<float>(_timers["st3"]->get_accumulated_time()/1000.),
467  static_cast<float>(_timers["st23"]->get_accumulated_time()/1000.),
468  static_cast<float>(_timers["global_kalman"]->get_accumulated_time()/1000.),
469  static_cast<float>(_timers["global"]->get_accumulated_time()/1000.),
470  static_cast<float>(_timers["kalman"]->get_accumulated_time()/1000.),
471  static_cast<float>(_timers["event"]->get_accumulated_time()/1000.),
472  };
473 
474  _tana_Event->Fill(data);
475  }
476 
477  ++_event;
478 
479  return ret;
480 }
481 
482 int KalmanDSTrk::setRawEventWorker(SRawEvent* event_input)
483 {
484  //Initialize tracklet lists
485  for(int i = 0; i < 5; i++) trackletsInSt[i].clear();
486  stracks.clear();
487 
488  //pre-tracking cuts
489  rawEvent = event_input;
490  if(!acceptEvent(rawEvent)) return TFEXIT_FAIL_MULTIPLICITY;
491  hitAll = event_input->getAllHits();
492 #ifdef _DEBUG_ON
493  for(std::vector<Hit>::iterator iter = hitAll.begin(); iter != hitAll.end(); ++iter) iter->print();
494 #endif
495 
496  //Initialize hodo and masking IDs
497  for(int i = 0; i < 4; i++)
498  {
499  //std::cout << "For station " << i << std::endl;
500  hitIDs_mask[i].clear();
501  if(p_jobOptsSvc->m_mcMode || rawEvent->isFPGATriggered())
502  {
503  hitIDs_mask[i] = rawEvent->getHitsIndexInDetectors(detectorIDs_maskX[i]);
504  }
505  else
506  {
507  hitIDs_mask[i] = rawEvent->getHitsIndexInDetectors(detectorIDs_maskY[i]);
508  }
509 
510  //for(std::list<int>::iterator iter = hitIDs_mask[i].begin(); iter != hitIDs_mask[i].end(); ++iter) std::cout << *iter << " " << hitAll[*iter].detectorID << " === ";
511  //std::cout << std::endl;
512  }
513 
514  //Initialize prop. tube IDs
515  for(int i = 0; i < 2; ++i)
516  {
517  for(int j = 0; j < 4; ++j)
518  {
519  hitIDs_muid[i][j].clear();
520  hitIDs_muid[i][j] = rawEvent->getHitsIndexInDetector(detectorIDs_muid[i][j]);
521  }
522  hitIDs_muidHodoAid[i].clear();
523  hitIDs_muidHodoAid[i] = rawEvent->getHitsIndexInDetectors(detectorIDs_muidHodoAid[i]);
524  }
525 
526 #ifndef ALIGNMENT_MODE
527 #ifdef _DEBUG_ON
528  LogInfo("ALIGNMENT_MODE OFF");
529 #endif
530  buildPropSegments();
531  if(propSegs[0].empty() || propSegs[1].empty())
532  {
533 #ifdef _DEBUG_ON
534  LogInfo("Failed in prop tube segment building: " << propSegs[0].size() << ", " << propSegs[1].size());
535 #endif
536  //return TFEXIT_FAIL_ROUGH_MUONID;
537  }
538 #endif
539 
540  //Build tracklets in station 2, 3+, 3-
541  _timers["st2"]->restart();
542  buildTrackletsInStation(3, 1); //3 for station-2, 1 for list position 1
543  _timers["st2"]->stop();
544 
545  if(verbosity >= 2) {
546  LogInfo("NTracklets in St2: " << trackletsInSt[1].size());
547  }
548 
549  if(trackletsInSt[1].empty())
550  {
551 #ifdef _DEBUG_ON
552  LogInfo("Failed in tracklet build at station 2");
553 #endif
554  return TFEXIT_FAIL_ST2_TRACKLET;
555  }
556 
557  _timers["st3"]->restart();
558  buildTrackletsInStation(4, 2); //4 for station-3+
559  buildTrackletsInStation(5, 2); //5 for station-3-
560  _timers["st3"]->stop();
561 
562  if(verbosity >= 2) {
563  LogInfo("NTracklets in St3: " << trackletsInSt[2].size());
564  }
565 
566  if(trackletsInSt[2].empty())
567  {
568 #ifdef _DEBUG_ON
569  LogInfo("Failed in tracklet build at station 3");
570 #endif
571  return TFEXIT_FAIL_ST3_TRACKLET;
572  }
573 
574  //Build back partial tracks in station 2, 3+ and 3-
575  buildBackPartialTracks();
576 
577  if(verbosity >= 2) {
578  LogInfo("NTracklets St2+St3: " << trackletsInSt[3].size());
579  }
580 
581  if(trackletsInSt[3].empty())
582  {
583 #ifdef _DEBUG_ON
584  LogInfo("Failed in connecting station-2 and 3 tracks");
585 #endif
586  return TFEXIT_FAIL_BACKPARTIAL;
587  }
588 
589  //Connect tracklets in station 2/3 and station 1 to form global tracks
590  buildGlobalTracks();
591 
592  if(verbosity >= 2) {
593  LogInfo("NTracklets Global: " << trackletsInSt[4].size());
594  }
595 
596 #ifdef _DEBUG_ON
597  for(int i = 0; i < 2; ++i)
598  {
599  std::cout << "=======================================================================================" << std::endl;
600  LogInfo("Prop tube segments in " << (i == 0 ? "X-Z" : "Y-Z"));
601  for(std::list<PropSegment>::iterator seg = propSegs[i].begin(); seg != propSegs[i].end(); ++seg)
602  {
603  seg->print();
604  }
605  std::cout << "=======================================================================================" << std::endl;
606  }
607 
608  for(int i = 0; i <= 4; i++)
609  {
610  std::cout << "=======================================================================================" << std::endl;
611  LogInfo("Final tracklets in station: " << i+1 << " is " << trackletsInSt[i].size());
612  for(std::list<Tracklet>::iterator tracklet = trackletsInSt[i].begin(); tracklet != trackletsInSt[i].end(); ++tracklet)
613  {
614  tracklet->print();
615  }
616  std::cout << "=======================================================================================" << std::endl;
617  }
618 #endif
619 
620  if(trackletsInSt[4].empty()) return TFEXIT_FAIL_GLOABL;
621  if(!_enable_KF) return TFEXIT_SUCCESS;
622 
623  //Build kalman tracks
624  _timers["kalman"]->restart();
625  for(std::list<Tracklet>::iterator tracklet = trackletsInSt[4].begin(); tracklet != trackletsInSt[4].end(); ++tracklet)
626  {
627  SRecTrack strack = processOneTracklet(*tracklet);
628  stracks.push_back(strack);
629  }
630  _timers["kalman"]->stop();
631 
632 #ifdef _DEBUG_ON
633  LogInfo(stracks.size() << " final tracks:");
634  for(std::list<SRecTrack>::iterator strack = stracks.begin(); strack != stracks.end(); ++strack)
635  {
636  strack->print();
637  }
638 #endif
639 
640  return TFEXIT_SUCCESS;
641 }
642 
643 bool KalmanDSTrk::acceptEvent(SRawEvent* rawEvent)
644 {
645 
646  if(Verbosity() >= Fun4AllBase::VERBOSITY_A_LOT) {
647  LogInfo("D0: " << rawEvent->getNHitsInD0());
648  LogInfo("D1: " << rawEvent->getNHitsInD1());
649  LogInfo("D2: " << rawEvent->getNHitsInD2());
650  LogInfo("D3p: " << rawEvent->getNHitsInD3p());
651  LogInfo("D3m: " << rawEvent->getNHitsInD3m());
652  LogInfo("H1: " << rawEvent->getNHitsInDetectors(detectorIDs_maskX[0]));
653  LogInfo("H2: " << rawEvent->getNHitsInDetectors(detectorIDs_maskX[1]));
654  LogInfo("H3: " << rawEvent->getNHitsInDetectors(detectorIDs_maskX[2]));
655  LogInfo("H4: " << rawEvent->getNHitsInDetectors(detectorIDs_maskX[3]));
656  LogInfo("Prop:" << rawEvent->getNPropHitsAll());
657  LogInfo("NRoadsPos: " << rawEvent->getNRoadsPos());
658  LogInfo("NRoadsNeg: " << rawEvent->getNRoadsNeg());
659  }
660 
661  if(rawEvent->getNHitsInD0() > 350) return false;
662  if(rawEvent->getNHitsInD1() > 350) return false; // 31% - 58 hit per plane
663  if(rawEvent->getNHitsInD2() > 170) return false; // 23% - 28 hit per plane
664  if(rawEvent->getNHitsInD3p() > 140) return false;// 18% - 23 hit per plane
665  if(rawEvent->getNHitsInD3m() > 140) return false;// 18% - 23 hit per plane
666 
667  /*
668  if(rawEvent->getNHitsInDetectors(detectorIDs_maskX[0]) > 15) return false;
669  if(rawEvent->getNHitsInDetectors(detectorIDs_maskX[1]) > 10) return false;
670  if(rawEvent->getNHitsInDetectors(detectorIDs_maskX[2]) > 10) return false;
671  if(rawEvent->getNHitsInDetectors(detectorIDs_maskX[3]) > 10) return false;
672  if(rawEvent->getNPropHitsAll() > 300) return false;
673 
674  if(rawEvent->getNRoadsPos() > 5) return false;
675  if(rawEvent->getNRoadsNeg() > 5) return false;
676  */
677 
678  return true;
679 }
680 
681 void KalmanDSTrk::buildBackPartialTracks()
682 {
683  _timers["st23"]->restart();
684 #ifndef ALIGNMENT_MODE
685  //Temporary container for a simple chisq fit
686  int nHitsX2, nHitsX3;
687  double z_fit[4], x_fit[4];
688  double a, b;
689 #endif
690 
691 #ifdef _DEBUG_YUHW_
692  std::map<std::string, int> counter = {
693  {"in", 0},
694  {"DS", 0},
695  {"out", 0}
696  };
697 #endif
698 
699  for(std::list<Tracklet>::iterator tracklet3 = trackletsInSt[2].begin(); tracklet3 != trackletsInSt[2].end(); ++tracklet3)
700  {
701 #ifndef ALIGNMENT_MODE
702  //Extract the X hits only from station-3 tracks
703  nHitsX3 = 0;
704  for(std::list<SignedHit>::iterator ptr_hit = tracklet3->hits.begin(); ptr_hit != tracklet3->hits.end(); ++ptr_hit)
705  {
706  if(ptr_hit->hit.index < 0) continue;
707  if(p_geomSvc->getPlaneType(ptr_hit->hit.detectorID) == 1)
708  {
709  z_fit[nHitsX3] = z_plane[ptr_hit->hit.detectorID];
710  x_fit[nHitsX3] = ptr_hit->hit.pos;
711  ++nHitsX3;
712  }
713  }
714 #endif
715  Tracklet tracklet_best;
716  for(std::list<Tracklet>::iterator tracklet2 = trackletsInSt[1].begin(); tracklet2 != trackletsInSt[1].end(); ++tracklet2)
717  {
718  if(fabs(tracklet2->tx - tracklet3->tx) > 0.15 || fabs(tracklet2->ty - tracklet3->ty) > 0.1) continue;
719 #ifdef _DEBUG_YUHW_
720  counter["in"]++;
721 #endif
722  // Pattern dictionary search
723  if(_DS_level >= KalmanDSTrk::ST23_DS) {
724  _timers["search_db_23"]->restart();
725  bool matched = false;
726 
727  auto key2 = PatternDBUtil::GetTrackletKey(*tracklet2, PatternDB::DC2);
728  auto key3p = PatternDBUtil::GetTrackletKey(*tracklet3, PatternDB::DC3p);
729  auto key3m = PatternDBUtil::GetTrackletKey(*tracklet3, PatternDB::DC3m);
730 
731  PartTrackKey key23;
732  if(key2!=PatternDB::ERR_KEY and key3p!=PatternDB::ERR_KEY) {key23 = PartTrackKey(key2, key3p);}
733  else if(key2!=PatternDB::ERR_KEY and key3m!=PatternDB::ERR_KEY) {key23 = PartTrackKey(key2, key3m);}
734  else continue;
735 
736  //LogInfo(key23);
737  if(_pattern_db->St23.find(key23)!=_pattern_db->St23.end()) matched = true;
738 
739  _timers["search_db_23"]->stop();
740  if(!matched) {
741  if(Verbosity() > 20) {
742  LogInfo("St23 Pattern NOT Found!");
743  }
744  continue;
745  }
746 
747 #ifdef _DEBUG_YUHW_
748  counter["DS"]++;
749 #endif
750  }
751 
752 #ifndef ALIGNMENT_MODE
753  //Extract the X hits from station-2 tracke
754  nHitsX2 = nHitsX3;
755  for(std::list<SignedHit>::iterator ptr_hit = tracklet2->hits.begin(); ptr_hit != tracklet2->hits.end(); ++ptr_hit)
756  {
757  if(ptr_hit->hit.index < 0) continue;
758  if(p_geomSvc->getPlaneType(ptr_hit->hit.detectorID) == 1)
759  {
760  z_fit[nHitsX2] = z_plane[ptr_hit->hit.detectorID];
761  x_fit[nHitsX2] = ptr_hit->hit.pos;
762  ++nHitsX2;
763  }
764  }
765 
766  if(verbosity>2) _timers["st23_fit1"]->restart();
767  //Apply a simple linear fit to get rough estimation of X-Z slope and intersection
768  chi2fit(nHitsX2, z_fit, x_fit, a, b);
769  if(verbosity>2) _timers["st23_fit1"]->stop();
770  if(fabs(a) > 2.*TX_MAX || fabs(b) > 2.*X0_MAX) continue;
771 
772  if(verbosity>2) _timers["st23_prop"]->restart();
773  //Project to proportional tubes to see if there is enough
774  int nPropHits = 0;
775  for(int i = 0; i < 4; ++i)
776  {
777  double x_exp = a*z_mask[detectorIDs_muid[0][i] - nChamberPlanes - 1] + b;
778  for(std::list<int>::iterator iter = hitIDs_muid[0][i].begin(); iter != hitIDs_muid[0][i].end(); ++iter)
779  {
780  if(fabs(hitAll[*iter].pos - x_exp) < 5.08)
781  {
782  ++nPropHits;
783  break;
784  }
785  }
786  if(nPropHits > 0) break;
787  }
788  if(verbosity>2) _timers["st23_prop"]->stop();
789  if(nPropHits == 0) continue;
790 #endif
791 
792  Tracklet tracklet_23 = (*tracklet2) + (*tracklet3);
793 #ifdef _DEBUG_ON
794  LogInfo("Using following two tracklets:");
795  tracklet2->print();
796  tracklet3->print();
797  LogInfo("Yield this combination:");
798  tracklet_23.print();
799 #endif
800  if(verbosity>2) _timers["st23_fit2"]->restart();
801  fitTracklet(tracklet_23);
802  if(verbosity>2) _timers["st23_fit2"]->stop();
803 
804  if(tracklet_23.chisq > 9000.)
805  {
806 #ifdef _DEBUG_ON
807  tracklet_23.print();
808  LogInfo("Impossible combination!");
809 #endif
810  continue;
811  }
812 
813  if(verbosity>2) _timers["st23_hodo"]->restart();
814  if(!hodoMask(tracklet_23))
815  {
816  if(verbosity>2) _timers["st23_hodo"]->stop();
817 #ifdef _DEBUG_ON
818  LogInfo("Hodomasking failed!");
819 #endif
820  continue;
821  }
822 #ifdef _DEBUG_ON
823  LogInfo("Hodomasking Scucess!");
824 #endif
825 
826 #ifndef COARSE_MODE
827  if(verbosity>2) _timers["st23_lr40"]->restart();
828  resolveLeftRight(tracklet_23, 40.);
829  if(verbosity>2) _timers["st23_lr40"]->stop();
830  if(verbosity>2) _timers["st23_lr150"]->restart();
831  resolveLeftRight(tracklet_23, 150.);
832  if(verbosity>2) _timers["st23_lr150"]->stop();
833 #endif
834  if(verbosity>2) _timers["st23_rm_hits"]->restart();
836  removeBadHits(tracklet_23);
837  if(verbosity>2) _timers["st23_rm_hits"]->stop();
838 
839 #ifdef _DEBUG_ON
840  LogInfo("New tracklet: ");
841  tracklet_23.print();
842 
843  LogInfo("Current best:");
844  tracklet_best.print();
845 
846  LogInfo("Comparison: " << (tracklet_23 < tracklet_best));
847  LogInfo("Quality: " << acceptTracklet(tracklet_23));
848 #endif
849 
850  if(verbosity>2) _timers["st23_acc"]->restart();
851  //If current tracklet is better than the best tracklet up-to-now
852  if(acceptTracklet(tracklet_23) && tracklet_23 < tracklet_best)
853  {
854  tracklet_best = tracklet_23;
855  }
856 #ifdef _DEBUG_YUHW_
857  counter["out"]++;
858 #endif
859  if(verbosity>2) _timers["st23_acc"]->stop();
860 #ifdef _DEBUG_ON
861  else
862  {
863  LogInfo("Rejected!!");
864  }
865 #endif
866  }
867 
868  if(verbosity>2) _timers["st23_reduce"]->restart();
869  if(tracklet_best.isValid()) trackletsInSt[3].push_back(tracklet_best);
870  if(verbosity>2) _timers["st23_reduce"]->stop();
871  }
872 
873  reduceTrackletList(trackletsInSt[3]);
874  trackletsInSt[3].sort();
875 
876  _timers["st23"]->stop();
877 
878 #ifdef _DEBUG_YUHW_
879  if(Verbosity() >= Fun4AllBase::VERBOSITY_A_LOT) {
880  LogInfo("");
881  std::cout
882  << counter["in"] << " => "
883  << counter["DS"] << " => "
884  << counter["out"] << std::endl;
885  }
886 
887  if(_ana_mode) {
888  float tana_data[] = {
889  static_cast<float>(counter["in"]),
890  static_cast<float>(counter["DS"]),
891  static_cast<float>(counter["out"]),
892  static_cast<float>(_timers["st23"]->get_accumulated_time()/1000.)
893  };
894  _tana_St23->Fill(tana_data);
895  }
896 #endif
897 }
898 
899 void KalmanDSTrk::buildGlobalTracks()
900 {
901  _timers["global"]->restart();
902 
903 #ifdef _DEBUG_YUHW_
904  std::map<std::string, int> counter = {
905  {"in", 0},
906  {"DS", 0},
907  {"out", 0}
908  };
909 #endif
910 
911  double pos_exp[3], window[3];
912  for(std::list<Tracklet>::iterator tracklet23 = trackletsInSt[3].begin(); tracklet23 != trackletsInSt[3].end(); ++tracklet23)
913  {
914  Tracklet tracklet_best[2];
915  for(int i = 0; i < 2; ++i) //for two station-1 chambers
916  {
917  trackletsInSt[0].clear();
918 
919  //Calculate the window in station 1
920  if(p_jobOptsSvc->m_enableKMag)
921  {
922  getSagittaWindowsInSt1(*tracklet23, pos_exp, window, i+1);
923  }
924  else
925  {
926  getExtrapoWindowsInSt1(*tracklet23, pos_exp, window, i+1);
927  }
928 
929 #ifdef _DEBUG_ON
930  LogInfo("Using this back partial: ");
931  tracklet23->print();
932  for(int j = 0; j < 3; j++) LogInfo("Extrapo: " << pos_exp[j] << " " << window[j]);
933 #endif
934 
935  _timers["global_st1"]->restart();
936  buildTrackletsInStation(i+1, 0, pos_exp, window);
937  _timers["global_st1"]->stop();
938 
939  _timers["global_link"]->restart();
940  Tracklet tracklet_best_prob, tracklet_best_vtx;
941  for(std::list<Tracklet>::iterator tracklet1 = trackletsInSt[0].begin(); tracklet1 != trackletsInSt[0].end(); ++tracklet1)
942  {
943 #ifdef _DEBUG_ON
944  LogInfo("With this station 1 track:");
945  tracklet1->print();
946 #endif
947 
948 #ifdef _DEBUG_YUHW_
949  counter["in"]++;
950 #endif
951  // Pattern dictionary search
952  if(_DS_level >= KalmanDSTrk::ST123_DS) {
953  _timers["search_db_glb"]->restart();
954  bool matched = false;
955 
956  auto key1 = PatternDBUtil::GetTrackletKey(*tracklet1, PatternDB::DC1);
957  auto key2 = PatternDBUtil::GetTrackletKey(*tracklet23, PatternDB::DC2);
958  auto key3p = PatternDBUtil::GetTrackletKey(*tracklet23, PatternDB::DC3p);
959  auto key3m = PatternDBUtil::GetTrackletKey(*tracklet23, PatternDB::DC3m);
960 
961  GlobTrackKey key123;
962  if(
963  key1!=PatternDB::ERR_KEY
964  and key2!=PatternDB::ERR_KEY
965  and key3p!=PatternDB::ERR_KEY) {
966  key123 = GlobTrackKey(key1, key2, key3p);
967  } else if (
968  key1!=PatternDB::ERR_KEY
969  and key2!=PatternDB::ERR_KEY
970  and key3m!=PatternDB::ERR_KEY) {
971  key123 = GlobTrackKey(key1, key2, key3m);
972  } else continue;
973 
974  //std::cout << key123;
975  if(_pattern_db->St123.find(key123)!=_pattern_db->St123.end()) matched = true;
976 
977  _timers["search_db_glb"]->stop();
978 
979 
980  if(Verbosity() > 20) {
981  std::cout << key123;
982  if(matched) {
983  std::cout<< "St123 Pattern Found!";
984  } else {
985  std::cout<< "St123 Pattern NOT Found!";
986  }
987  std::cout << std::endl;
988  }
989 
990  if(!matched) {
991  continue;
992  }
993 
994 #ifdef _DEBUG_YUHW_
995  counter["DS"]++;
996 #endif
997  }
998 
999  Tracklet tracklet_global = (*tracklet23) * (*tracklet1);
1000  fitTracklet(tracklet_global);
1001  if(!hodoMask(tracklet_global)) continue;
1002 
1003 #ifndef COARSE_MODE
1004  resolveLeftRight(tracklet_global, 75.);
1006  resolveLeftRight(tracklet_global, 150.);
1007  resolveSingleLeftRight(tracklet_global);
1008 #endif
1009  removeBadHits(tracklet_global);
1011 
1012  //Most basic cuts
1013  if(!acceptTracklet(tracklet_global)) continue;
1014 
1015  //Get the tracklets that has the best prob
1016  if(tracklet_global < tracklet_best_prob) tracklet_best_prob = tracklet_global;
1017 
1018 #if !defined(ALIGNMENT_MODE) && defined(_ENABLE_KF)
1019  _timers["global_kalman"]->restart();
1021  SRecTrack recTrack = processOneTracklet(tracklet_global);
1022  _timers["global_kalman"]->stop();
1023  tracklet_global.chisq_vtx = recTrack.getChisqVertex();
1024 
1025  if(recTrack.isValid() && tracklet_global.chisq_vtx < tracklet_best_vtx.chisq_vtx) tracklet_best_vtx = tracklet_global;
1026 #endif
1027 
1028 #ifdef _DEBUG_ON
1029  LogInfo("New tracklet: ");
1030  tracklet_global.print();
1031 
1032  LogInfo("Current best by prob:");
1033  tracklet_best_prob.print();
1034 
1035  LogInfo("Comparison I: " << (tracklet_global < tracklet_best_prob));
1036  LogInfo("Quality I : " << acceptTracklet(tracklet_global));
1037 
1038 #if !defined(ALIGNMENT_MODE) && defined(_ENABLE_KF)
1039  LogInfo("Current best by vtx:");
1040  tracklet_best_vtx.print();
1041 
1042  LogInfo("Comparison II: " << (tracklet_global.chisq_vtx < tracklet_best_vtx.chisq_vtx));
1043  LogInfo("Quality II : " << recTrack.isValid());
1044 #endif
1045 #endif
1046  }
1047  _timers["global_link"]->stop();
1048 
1049 
1050 #if !defined(ALIGNMENT_MODE) && defined(_ENABLE_KF)
1051  //The selection logic is, prefer the tracks with best p-value, as long as it's not low-pz
1052  if(tracklet_best_prob.isValid() && 1./tracklet_best_prob.invP > 18.)
1053  {
1054  tracklet_best[i] = tracklet_best_prob;
1055  }
1056  else if(tracklet_best_vtx.isValid()) //otherwise select the one with best vertex chisq, TODO: maybe add a z-vtx constraint
1057  {
1058  tracklet_best[i] = tracklet_best_vtx;
1059  }
1060  else if(tracklet_best_prob.isValid()) //then fall back to the default only choice
1061  {
1062  tracklet_best[i] = tracklet_best_prob;
1063  }
1064 #else
1065  if(tracklet_best_prob.isValid()) //then fall back to the default only choice
1066  {
1067  tracklet_best[i] = tracklet_best_prob;
1068  }
1069 #endif
1070  }
1071 
1072  //Merge the tracklets from two stations if necessary
1073  Tracklet tracklet_merge;
1074  if(fabs(tracklet_best[0].getMomentum() - tracklet_best[1].getMomentum())/tracklet_best[0].getMomentum() < MERGE_THRES)
1075  {
1076  //Merge the track and re-fit
1077  tracklet_merge = tracklet_best[0].merge(tracklet_best[1]);
1078  fitTracklet(tracklet_merge);
1079 
1080 #ifdef _DEBUG_ON
1081  LogInfo("Merging two track candidates with momentum: " << tracklet_best[0].getMomentum() << " " << tracklet_best[1].getMomentum());
1082  LogInfo("tracklet_best_1:"); tracklet_best[0].print();
1083  LogInfo("tracklet_best_2:"); tracklet_best[1].print();
1084  LogInfo("tracklet_merge:"); tracklet_merge.print();
1085 #endif
1086  }
1087 
1088  if(tracklet_merge.isValid() && tracklet_merge < tracklet_best[0] && tracklet_merge < tracklet_best[1])
1089  {
1090 #ifdef _DEBUG_ON
1091  LogInfo("Choose merged tracklet");
1092 #endif
1093 #ifdef _DEBUG_YUHW_
1094  counter["out"]++;
1095 #endif
1096  trackletsInSt[4].push_back(tracklet_merge);
1097  }
1098  else if(tracklet_best[0].isValid() && tracklet_best[0] < tracklet_best[1])
1099  {
1100 #ifdef _DEBUG_ON
1101  LogInfo("Choose tracklet with station-0");
1102 #endif
1103 #ifdef _DEBUG_YUHW_
1104  counter["out"]++;
1105 #endif
1106  trackletsInSt[4].push_back(tracklet_best[0]);
1107  }
1108  else if(tracklet_best[1].isValid())
1109  {
1110 #ifdef _DEBUG_ON
1111  LogInfo("Choose tracklet with station-1");
1112 #endif
1113 #ifdef _DEBUG_YUHW_
1114  counter["out"]++;
1115 #endif
1116  trackletsInSt[4].push_back(tracklet_best[1]);
1117  }
1118  }
1119 
1120  trackletsInSt[4].sort();
1121 
1122  _timers["global"]->stop();
1123 
1124 #ifdef _DEBUG_YUHW_
1125  if(Verbosity() >= Fun4AllBase::VERBOSITY_A_LOT) {
1126  LogInfo("");
1127  std::cout
1128  << counter["in"] << " => "
1129  << counter["DS"] << " => "
1130  << counter["out"] << std::endl;
1131  }
1132 
1133  if(_ana_mode) {
1134  float tana_data[] = {
1135  static_cast<float>(counter["in"]),
1136  static_cast<float>(counter["DS"]),
1137  static_cast<float>(counter["out"]),
1138  static_cast<float>(_timers["global"]->get_accumulated_time()/1000.)
1139  };
1140  _tana_St123->Fill(tana_data);
1141  }
1142 #endif
1143 }
1144 
1145 void KalmanDSTrk::resolveLeftRight(Tracklet& tracklet, double threshold)
1146 {
1147 #ifdef _DEBUG_ON
1148  LogInfo("Left right for this track..");
1149  tracklet.print();
1150 #endif
1151 
1152  //Check if the track has been updated
1153  bool isUpdated = false;
1154 
1155  //Four possibilities
1156  int possibility[4][2] = {{1, 1}, {1, -1}, {-1, 1}, {-1, -1}};
1157 
1158  //Total number of hit pairs in this tracklet
1159  int nPairs = tracklet.hits.size()/2;
1160 
1161  int nResolved = 0;
1162  std::list<SignedHit>::iterator hit1 = tracklet.hits.begin();
1163  std::list<SignedHit>::iterator hit2 = tracklet.hits.begin();
1164  ++hit2;
1165  while(true)
1166  {
1167 #ifdef _DEBUG_ON
1168  LogInfo(hit1->hit.index << " " << hit2->sign << " === " << hit2->hit.index << " " << hit2->sign);
1169  int detectorID1 = hit1->hit.detectorID;
1170  int detectorID2 = hit2->hit.detectorID;
1171  LogInfo("Hit1: " << tracklet.getExpPositionX(z_plane[detectorID1])*costheta_plane[detectorID1] + tracklet.getExpPositionY(z_plane[detectorID1])*sintheta_plane[detectorID1] << " " << hit1->hit.pos + hit1->hit.driftDistance << " " << hit1->hit.pos - hit1->hit.driftDistance);
1172  LogInfo("Hit2: " << tracklet.getExpPositionX(z_plane[detectorID2])*costheta_plane[detectorID2] + tracklet.getExpPositionY(z_plane[detectorID2])*sintheta_plane[detectorID2] << " " << hit2->hit.pos + hit2->hit.driftDistance << " " << hit2->hit.pos - hit2->hit.driftDistance);
1173 #endif
1174 
1175  if(hit1->hit.index > 0 && hit2->hit.index > 0 && hit1->sign*hit2->sign == 0)
1176  {
1177  int index_min = -1;
1178  double pull_min = 1E6;
1179  for(int i = 0; i < 4; i++)
1180  {
1181  double slope_local = (hit1->pos(possibility[i][0]) - hit2->pos(possibility[i][1]))/(z_plane[hit1->hit.detectorID] - z_plane[hit2->hit.detectorID]);
1182  double inter_local = hit1->pos(possibility[i][0]) - slope_local*z_plane[hit1->hit.detectorID];
1183 
1184  if(fabs(slope_local) > slope_max[hit1->hit.detectorID] || fabs(inter_local) > intersection_max[hit1->hit.detectorID]) continue;
1185 
1186  double tx, ty, x0, y0;
1187  double err_tx, err_ty, err_x0, err_y0;
1188  if(tracklet.stationID == 6 && hit1->hit.detectorID <= 6)
1189  {
1190  tracklet.getXZInfoInSt1(tx, x0);
1191  tracklet.getXZErrorInSt1(err_tx, err_x0);
1192  }
1193  else
1194  {
1195  tx = tracklet.tx;
1196  x0 = tracklet.x0;
1197  err_tx = tracklet.err_tx;
1198  err_x0 = tracklet.err_x0;
1199  }
1200  ty = tracklet.ty;
1201  y0 = tracklet.y0;
1202  err_ty = tracklet.err_ty;
1203  err_y0 = tracklet.err_y0;
1204 
1205  double slope_exp = costheta_plane[hit1->hit.detectorID]*tx + sintheta_plane[hit1->hit.detectorID]*ty;
1206  double err_slope = fabs(costheta_plane[hit1->hit.detectorID]*err_tx) + fabs(sintheta_plane[hit2->hit.detectorID]*err_ty);
1207  double inter_exp = costheta_plane[hit1->hit.detectorID]*x0 + sintheta_plane[hit1->hit.detectorID]*y0;
1208  double err_inter = fabs(costheta_plane[hit1->hit.detectorID]*err_x0) + fabs(sintheta_plane[hit2->hit.detectorID]*err_y0);
1209 
1210  double pull = sqrt((slope_exp - slope_local)*(slope_exp - slope_local)/err_slope/err_slope + (inter_exp - inter_local)*(inter_exp - inter_local)/err_inter/err_inter);
1211  if(pull < pull_min)
1212  {
1213  index_min = i;
1214  pull_min = pull;
1215  }
1216 
1217 #ifdef _DEBUG_ON
1218  LogInfo(hit1->hit.detectorID << ": " << i << " " << possibility[i][0] << " " << possibility[i][1]);
1219  LogInfo(tx << " " << x0 << " " << ty << " " << y0);
1220  LogInfo("Slope: " << slope_local << " " << slope_exp << " " << err_slope);
1221  LogInfo("Intersection: " << inter_local << " " << inter_exp << " " << err_inter);
1222  LogInfo("Current: " << pull << " " << index_min << " " << pull_min);
1223 #endif
1224  }
1225 
1226  //LogInfo("Final: " << index_min << " " << pull_min);
1227  if(index_min >= 0 && pull_min < threshold)//((tracklet.stationID == 5 && pull_min < 25.) || (tracklet.stationID == 6 && pull_min < 100.)))
1228  {
1229  //Origin imp.
1230  hit1->sign = possibility[index_min][0];
1231  hit2->sign = possibility[index_min][1];
1232  isUpdated = true;
1233  }
1234  }
1235 
1236  ++nResolved;
1237  if(nResolved >= nPairs) break;
1238 
1239  ++hit1;
1240  ++hit1;
1241  ++hit2;
1242  ++hit2;
1243  }
1244 #ifdef _DEBUG_ON
1245  tracklet.print();
1246 #endif
1247  if(isUpdated) fitTracklet(tracklet);
1248 #ifdef _DEBUG_ON
1249  tracklet.print();
1250 #endif
1251 }
1252 
1253 void KalmanDSTrk::resolveSingleLeftRight(Tracklet& tracklet)
1254 {
1255 #ifdef _DEBUG_ON
1256  LogInfo("Single left right for this track..");
1257  tracklet.print();
1258 #endif
1259 
1260  //Check if the track has been updated
1261  bool isUpdated = false;
1262  for(std::list<SignedHit>::iterator hit_sign = tracklet.hits.begin(); hit_sign != tracklet.hits.end(); ++hit_sign)
1263  {
1264  if(hit_sign->hit.index < 0 || hit_sign->sign != 0) continue;
1265 
1266  int detectorID = hit_sign->hit.detectorID;
1267  double pos_exp = tracklet.getExpPositionX(z_plane[detectorID])*costheta_plane[detectorID] + tracklet.getExpPositionY(z_plane[detectorID])*sintheta_plane[detectorID];
1268  hit_sign->sign = pos_exp > hit_sign->hit.pos ? 1 : -1;
1269 
1270  isUpdated = true;
1271  }
1272 
1273  if(isUpdated) fitTracklet(tracklet);
1274 }
1275 
1276 void KalmanDSTrk::removeBadHits(Tracklet& tracklet)
1277 {
1278 #ifdef _DEBUG_ON
1279  LogInfo("Removing hits for this track..");
1280  tracklet.calcChisq();
1281  tracklet.print();
1282 #endif
1283 
1284  //Check if the track has beed updated
1285  int signflipflag[nChamberPlanes];
1286  for(int i = 0; i < nChamberPlanes; ++i) signflipflag[i] = 0;
1287 
1288  bool isUpdated = true;
1289  while(isUpdated)
1290  {
1291  isUpdated = false;
1292  tracklet.calcChisq();
1293 
1294  SignedHit* hit_remove = NULL;
1295  SignedHit* hit_neighbour = NULL;
1296  double res_remove1 = -1.;
1297  double res_remove2 = -1.;
1298  for(std::list<SignedHit>::iterator hit_sign = tracklet.hits.begin(); hit_sign != tracklet.hits.end(); ++hit_sign)
1299  {
1300  if(hit_sign->hit.index < 0) continue;
1301 
1302  int detectorID = hit_sign->hit.detectorID;
1303  double res_curr = fabs(tracklet.residual[detectorID-1]);
1304  if(res_remove1 < res_curr)
1305  {
1306  res_remove1 = res_curr;
1307  res_remove2 = fabs(tracklet.residual[detectorID-1] - 2.*hit_sign->sign*hit_sign->hit.driftDistance);
1308  hit_remove = &(*hit_sign);
1309 
1310  std::list<SignedHit>::iterator iter = hit_sign;
1311  hit_neighbour = detectorID % 2 == 0 ? &(*(--iter)) : &(*(++iter));
1312  }
1313  }
1314  if(hit_remove == NULL) continue;
1315  if(hit_remove->sign == 0 && tracklet.isValid()) continue; //if sign is undecided, and chisq is OKay, then pass
1316 
1317  double cut = hit_remove->sign == 0 ? hit_remove->hit.driftDistance + resol_plane[hit_remove->hit.detectorID] : resol_plane[hit_remove->hit.detectorID];
1318  if(res_remove1 > cut)
1319  {
1320 #ifdef _DEBUG_ON
1321  LogInfo("Dropping this hit: " << res_remove1 << " " << res_remove2 << " " << signflipflag[hit_remove->hit.detectorID-1] << " " << cut);
1322  hit_remove->hit.print();
1323  hit_neighbour->hit.print();
1324 #endif
1325 
1326  //can only be changed less than twice
1327  if(res_remove2 < cut && signflipflag[hit_remove->hit.detectorID-1] < 2)
1328  {
1329  hit_remove->sign = -hit_remove->sign;
1330  hit_neighbour->sign = 0;
1331  ++signflipflag[hit_remove->hit.detectorID-1];
1332 #ifdef _DEBUG_ON
1333  LogInfo("Only changing the sign.");
1334 #endif
1335  }
1336  else
1337  {
1338  //Set the index of the hit to be removed to -1 so it's not used anymore
1339  //also set the sign assignment of the neighbour hit to 0 (i.e. undecided)
1340  hit_remove->hit.index = -1;
1341  hit_neighbour->sign = 0;
1342  int planeType = p_geomSvc->getPlaneType(hit_remove->hit.detectorID);
1343  if(planeType == 1)
1344  {
1345  --tracklet.nXHits;
1346  }
1347  else if(planeType == 2)
1348  {
1349  --tracklet.nUHits;
1350  }
1351  else
1352  {
1353  --tracklet.nVHits;
1354  }
1355 
1356  //If both hit pairs are not included, the track can be rejected
1357  if(hit_neighbour->hit.index < 0)
1358  {
1359 #ifdef _DEBUG_ON
1360  LogInfo("Both hits in a view are missing! Will exit the bad hit removal...");
1361 #endif
1362  return;
1363  }
1364  }
1365  isUpdated = true;
1366  }
1367 
1368  if(isUpdated)
1369  {
1370  fitTracklet(tracklet);
1371  resolveSingleLeftRight(tracklet);
1372  }
1373  }
1374 }
1375 
1376 void KalmanDSTrk::resolveLeftRight(SRawEvent::hit_pair hpair, int& LR1, int& LR2)
1377 {
1378  LR1 = 0;
1379  LR2 = 0;
1380 
1381  //If either hit is missing, no left-right can be assigned
1382  if(hpair.first < 0 || hpair.second < 0)
1383  {
1384  return;
1385  }
1386 
1387  int possibility[4][2] = {{1, 1}, {1, -1}, {-1, 1}, {-1, -1}};
1388  int nResolved = 0;
1389  for(int i = 0; i < 4; i++)
1390  {
1391  if(nResolved > 1) break;
1392 
1393  int hitID1 = hpair.first;
1394  int hitID2 = hpair.second;
1395  double slope_local = (hitAll[hitID1].pos + possibility[i][0]*hitAll[hitID1].driftDistance - hitAll[hitID2].pos - possibility[i][1]*hitAll[hitID2].driftDistance)/(z_plane[hitAll[hitID1].detectorID] - z_plane[hitAll[hitID2].detectorID]);
1396  double intersection_local = hitAll[hitID1].pos + possibility[i][0]*hitAll[hitID1].driftDistance - slope_local*z_plane[hitAll[hitID1].detectorID];
1397 
1398  //LogInfo(i << " " << nResolved << " " << slope_local << " " << intersection_local);
1399  if(fabs(slope_local) < slope_max[hitAll[hitID1].detectorID] && fabs(intersection_local) < intersection_max[hitAll[hitID1].detectorID])
1400  {
1401  nResolved++;
1402  LR1 = possibility[i][0];
1403  LR2 = possibility[i][1];
1404  }
1405  }
1406 
1407  if(nResolved > 1)
1408  {
1409  LR1 = 0;
1410  LR2 = 0;
1411  }
1412 
1413  //LogInfo("Final: " << LR1 << " " << LR2);
1414 }
1415 
1416 void KalmanDSTrk::buildTrackletsInStation(int stationID, int listID, double* pos_exp, double* window)
1417 {
1418 #ifdef _DEBUG_ON
1419  LogInfo("Building tracklets in station " << stationID);
1420 #endif
1421 
1422  //actuall ID of the tracklet lists
1423  int sID = stationID - 1;
1424 
1425  //Extract the X, U, V hit pairs
1426  std::list<SRawEvent::hit_pair> pairs_X, pairs_U, pairs_V;
1427  if(pos_exp == NULL)
1428  {
1429  pairs_X = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][0]);
1430  pairs_U = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][1]);
1431  pairs_V = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][2]);
1432  }
1433  else
1434  {
1435  //Note that in pos_exp[], index 0 stands for X, index 1 stands for U, index 2 stands for V
1436  pairs_X = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][0], pos_exp[0], window[0]);
1437  pairs_U = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][1], pos_exp[1], window[1]);
1438  pairs_V = rawEvent->getPartialHitPairsInSuperDetector(superIDs[sID][2], pos_exp[2], window[2]);
1439  }
1440 
1441 #ifdef _DEBUG_ON_
1442  LogInfo("Hit pairs in this event: ");
1443  for(std::list<SRawEvent::hit_pair>::iterator iter = pairs_X.begin(); iter != pairs_X.end(); ++iter) LogInfo("X :" << iter->first << " " << iter->second << " " << hitAll[iter->first].index << " " << (iter->second < 0 ? -1 : hitAll[iter->second].index));
1444  for(std::list<SRawEvent::hit_pair>::iterator iter = pairs_U.begin(); iter != pairs_U.end(); ++iter) LogInfo("U :" << iter->first << " " << iter->second << " " << hitAll[iter->first].index << " " << (iter->second < 0 ? -1 : hitAll[iter->second].index));
1445  for(std::list<SRawEvent::hit_pair>::iterator iter = pairs_V.begin(); iter != pairs_V.end(); ++iter) LogInfo("V :" << iter->first << " " << iter->second << " " << hitAll[iter->first].index << " " << (iter->second < 0 ? -1 : hitAll[iter->second].index));
1446 #endif
1447 
1448  if(pairs_X.empty() || pairs_U.empty() || pairs_V.empty())
1449  {
1450 #ifdef _DEBUG_ON
1451  LogInfo("Not all view has hits in station " << stationID);
1452 #endif
1453  return;
1454  }
1455 
1456 #ifdef _DEBUG_YUHW_
1457  std::map<std::string, int> counter = {
1458  {"in", 0},
1459  {"DS", 0},
1460  {"out", 0}
1461  };
1462 #endif
1463 
1464  //X-U combination first, then add V pairs
1465  for(std::list<SRawEvent::hit_pair>::iterator xiter = pairs_X.begin(); xiter != pairs_X.end(); ++xiter)
1466  {
1467  //U projections from X plane
1468  double x_pos = xiter->second >= 0 ? 0.5*(hitAll[xiter->first].pos + hitAll[xiter->second].pos) : hitAll[xiter->first].pos;
1469  double u_min = x_pos*u_costheta[sID] - u_win[sID];
1470  double u_max = u_min + 2.*u_win[sID];
1471 
1472 #ifdef _DEBUG_ON
1473  LogInfo("Trying X hits " << xiter->first << " " << xiter->second << " " << hitAll[xiter->first].elementID << " at " << x_pos);
1474  LogInfo("U plane window:" << u_min << " " << u_max);
1475 #endif
1476  for(std::list<SRawEvent::hit_pair>::iterator uiter = pairs_U.begin(); uiter != pairs_U.end(); ++uiter)
1477  {
1478  double u_pos = uiter->second >= 0 ? 0.5*(hitAll[uiter->first].pos + hitAll[uiter->second].pos) : hitAll[uiter->first].pos;
1479 #ifdef _DEBUG_ON
1480  LogInfo("Trying U hits " << uiter->first << " " << uiter->second << " " << hitAll[uiter->first].elementID << " at " << u_pos);
1481 #endif
1482  if(u_pos < u_min || u_pos > u_max) continue;
1483 
1484  //V projections from X and U plane
1485  double z_x = xiter->second >= 0 ? z_plane_x[sID] : z_plane[hitAll[xiter->first].detectorID];
1486  double z_u = uiter->second >= 0 ? z_plane_u[sID] : z_plane[hitAll[uiter->first].detectorID];
1487  double z_v = z_plane_v[sID];
1488  double v_win1 = spacing_plane[hitAll[uiter->first].detectorID]*2.*u_costheta[sID];
1489  double v_win2 = fabs((z_u + z_v - 2.*z_x)*u_costheta[sID]*TX_MAX);
1490  double v_win3 = fabs((z_v - z_u)*u_sintheta[sID]*TY_MAX);
1491  double v_win = v_win1 + v_win2 + v_win3 + 2.*spacing_plane[hitAll[uiter->first].detectorID];
1492  double v_min = 2*x_pos*u_costheta[sID] - u_pos - v_win;
1493  double v_max = v_min + 2.*v_win;
1494 
1495 #ifdef _DEBUG_ON
1496  LogInfo("V plane window:" << v_min << " " << v_max);
1497 
1498  std::cout
1499  << "v_win1 = " << v_win1 << std::endl
1500  << "v_win2 = " << v_win2 << std::endl
1501  << "v_win3 = " << v_win3 << std::endl
1502  << "2.*spacing_plane[hitAll[uiter->first].detectorID] = " << 2.*spacing_plane[hitAll[uiter->first].detectorID] << std::endl;
1503 #endif
1504  for(std::list<SRawEvent::hit_pair>::iterator viter = pairs_V.begin(); viter != pairs_V.end(); ++viter)
1505  {
1506 #ifdef _DEBUG_YUHW_
1507  counter["in"]++;
1508 #endif
1509  if(_DS_level >= KalmanDSTrk::IN_ST_DS) {
1510  if(stationID==3)
1511  _timers["search_db_2"]->restart();
1512 
1513  PatternDB::STATION station = PatternDB::ERROR_STATION;
1514  if(stationID==1 or stationID==2) station = PatternDB::DC1;
1515  if(stationID==3) station = PatternDB::DC2;
1516  if(stationID==4) station = PatternDB::DC3p;
1517  if(stationID==5) station = PatternDB::DC3m;
1518  bool matched = false;
1519 
1520  std::vector< std::pair<unsigned int, unsigned int> > det_elem_pairs;
1521 
1522  det_elem_pairs.push_back({hitAll[xiter->first].detectorID, hitAll[xiter->first].elementID});
1523  det_elem_pairs.push_back({hitAll[xiter->second].detectorID, hitAll[xiter->second].elementID});
1524  det_elem_pairs.push_back({hitAll[uiter->first].detectorID, hitAll[uiter->first].elementID});
1525  det_elem_pairs.push_back({hitAll[uiter->second].detectorID, hitAll[uiter->second].elementID});
1526  det_elem_pairs.push_back({hitAll[viter->first].detectorID, hitAll[viter->first].elementID});
1527  det_elem_pairs.push_back({hitAll[viter->second].detectorID, hitAll[viter->second].elementID});
1528 
1529  TrackletKey key = PatternDBUtil::GetTrackletKey(det_elem_pairs, station);
1530 
1531  if(station == PatternDB::DC1
1532  and _pattern_db->St1.find(key)!=_pattern_db->St1.end()) matched = true;
1533  else if(station == PatternDB::DC2
1534  and _pattern_db->St2.find(key)!=_pattern_db->St2.end()) matched = true;
1535  else if( ( station == PatternDB::DC3p or station == PatternDB::DC3m )
1536  and _pattern_db->St3.find(key)!=_pattern_db->St3.end()) matched = true;
1537 
1538  if(stationID==3)
1539  _timers["search_db_2"]->stop();
1540 
1541  if(Verbosity() > 20) {
1542  std::cout
1543  << hitAll[xiter->first].index << ", "
1544  << hitAll[xiter->second].index << ", "
1545  << hitAll[uiter->first].index << ", "
1546  << hitAll[uiter->second].index << ", "
1547  << hitAll[viter->first].index << ", "
1548  << hitAll[viter->second].index
1549  << std::endl;
1550 
1551  std::cout << key;
1552  if(matched) {
1553  std::cout<< "St" << station << " Pattern Found!";
1554  } else {
1555  std::cout<< "St" << station << " Pattern NOT Found!";
1556  }
1557  std::cout << std::endl;
1558  }
1559 
1560  if(!matched) {
1561  continue;
1562  }
1563 
1564 #ifdef _DEBUG_YUHW_
1565  counter["DS"]++;
1566 #endif
1567  }
1568 
1569  double v_pos = viter->second >= 0 ? 0.5*(hitAll[viter->first].pos + hitAll[viter->second].pos) : hitAll[viter->first].pos;
1570 #ifdef _DEBUG_ON
1571  LogInfo("Trying V hits " << viter->first << " " << viter->second << " " << hitAll[viter->first].elementID << " at " << v_pos);
1572 #endif
1573  if(v_pos < v_min || v_pos > v_max) continue;
1574 
1575  //Now add the tracklet
1576  int LR1 = 0;
1577  int LR2 = 0;
1578  Tracklet tracklet_new;
1579  tracklet_new.stationID = stationID;
1580 
1581  //resolveLeftRight(*xiter, LR1, LR2);
1582  if(xiter->first >= 0)
1583  {
1584  tracklet_new.hits.push_back(SignedHit(hitAll[xiter->first], LR1));
1585  tracklet_new.nXHits++;
1586  }
1587  if(xiter->second >= 0)
1588  {
1589  tracklet_new.hits.push_back(SignedHit(hitAll[xiter->second], LR2));
1590  tracklet_new.nXHits++;
1591  }
1592 
1593  //resolveLeftRight(*uiter, LR1, LR2);
1594  if(uiter->first >= 0)
1595  {
1596  tracklet_new.hits.push_back(SignedHit(hitAll[uiter->first], LR1));
1597  tracklet_new.nUHits++;
1598  }
1599  if(uiter->second >= 0)
1600  {
1601  tracklet_new.hits.push_back(SignedHit(hitAll[uiter->second], LR2));
1602  tracklet_new.nUHits++;
1603  }
1604 
1605  //resolveLeftRight(*viter, LR1, LR2);
1606  if(viter->first >= 0)
1607  {
1608  tracklet_new.hits.push_back(SignedHit(hitAll[viter->first], LR1));
1609  tracklet_new.nVHits++;
1610  }
1611  if(viter->second >= 0)
1612  {
1613  tracklet_new.hits.push_back(SignedHit(hitAll[viter->second], LR2));
1614  tracklet_new.nVHits++;
1615  }
1616 
1617  tracklet_new.sortHits();
1618  if(!tracklet_new.isValid())
1619  {
1620  fitTracklet(tracklet_new);
1621  }
1622  else
1623  {
1624  continue;
1625  }
1626 #ifdef _DEBUG_ON
1627  if(Verbosity()>=Fun4AllBase::VERBOSITY_A_LOT) {
1628  tracklet_new.print();
1629  }
1630 #endif
1631  if(acceptTracklet(tracklet_new))
1632  {
1633  trackletsInSt[listID].push_back(tracklet_new);
1634 #ifdef _DEBUG_YUHW_
1635  counter["out"]++;
1636 #endif
1637  }
1638 #ifdef _DEBUG_ON
1639  else
1640  {
1641  LogInfo("Rejected!!!");
1642  }
1643 #endif
1644  }
1645  }
1646  }
1647 
1648 #ifdef _DEBUG_YUHW_
1649  if(Verbosity() >= Fun4AllBase::VERBOSITY_A_LOT) {
1650  LogInfo("");
1651  std::cout
1652  << counter["in"] << " => "
1653  << counter["DS"] << " => "
1654  << counter["out"] << std::endl;
1655  }
1656 
1657  if(_ana_mode) {
1658  float tana_data[] = {
1659  static_cast<float>(counter["in"]),
1660  static_cast<float>(counter["DS"]),
1661  static_cast<float>(counter["out"])
1662  };
1663  if(stationID==2) _tana_St1->Fill(tana_data);
1664  if(stationID==3) _tana_St2->Fill(tana_data);
1665  if(stationID==4||stationID==5) _tana_St3->Fill(tana_data);
1666  }
1667 #endif
1668 
1669  //Reduce the tracklet list and add dummy hits
1670  //reduceTrackletList(trackletsInSt[listID]);
1671  for(std::list<Tracklet>::iterator iter = trackletsInSt[listID].begin(); iter != trackletsInSt[listID].end(); ++iter)
1672  {
1673  iter->addDummyHits();
1674  }
1675 
1676  //Only retain the best 200 tracklets if exceeded
1677  if(trackletsInSt[listID].size() > 200)
1678  {
1679  trackletsInSt[listID].sort();
1680  trackletsInSt[listID].resize(200);
1681  }
1682 }
1683 
1684 bool KalmanDSTrk::acceptTracklet(Tracklet& tracklet)
1685 {
1686  //Tracklet itself is okay with enough hits (4-out-of-6) and small chi square
1687  if(!tracklet.isValid())
1688  {
1689 #ifdef _DEBUG_ON
1690  LogInfo("Failed in quality check!");
1691 #endif
1692  return false;
1693  }
1694 
1695  //Hodoscope masking requirement
1696  if(!hodoMask(tracklet)) return false;
1697 
1698  //For back partials, require projection inside KMAG, and muon id in prop. tubes
1699  if(tracklet.stationID > nStations-2)
1700  {
1701  if(!p_geomSvc->isInKMAG(tracklet.getExpPositionX(Z_KMAG_BEND), tracklet.getExpPositionY(Z_KMAG_BEND))) return false;
1702 #ifndef ALIGNMENT_MODE
1703  if(!(muonID_comp(tracklet) || muonID_search(tracklet))) return false;
1704 #endif
1705  }
1706 
1707  //If everything is fine ...
1708 #ifdef _DEBUG_ON
1709  LogInfo("AcceptTracklet!!!");
1710 #endif
1711  return true;
1712 }
1713 
1714 bool KalmanDSTrk::hodoMask(Tracklet& tracklet)
1715 {
1716  //LogInfo(tracklet.stationID);
1717  int nHodoHits = 0;
1718  for(std::vector<int>::iterator stationID = stationIDs_mask[tracklet.stationID-1].begin(); stationID != stationIDs_mask[tracklet.stationID-1].end(); ++stationID)
1719  {
1720  bool masked = false;
1721  for(std::list<int>::iterator iter = hitIDs_mask[*stationID-1].begin(); iter != hitIDs_mask[*stationID-1].end(); ++iter)
1722  {
1723  int detectorID = hitAll[*iter].detectorID;
1724  int elementID = hitAll[*iter].elementID;
1725 
1726  int idx1 = detectorID - nChamberPlanes - 1;
1727  int idx2 = elementID - 1;
1728 
1729  double factor = tracklet.stationID == nChamberPlanes/6-2 ? 5. : 3.; //special for station-2, based on real data tuning
1730  double xfudge = tracklet.stationID < nStations-1 ? 0.5*(x_mask_max[idx1][idx2] - x_mask_min[idx1][idx2]) : 0.15*(x_mask_max[idx1][idx2] - x_mask_min[idx1][idx2]);
1731  double z_hodo = z_mask[idx1];
1732  double x_hodo = tracklet.getExpPositionX(z_hodo);
1733  double y_hodo = tracklet.getExpPositionY(z_hodo);
1734  double err_x = factor*tracklet.getExpPosErrorX(z_hodo) + xfudge;
1735  double err_y = factor*tracklet.getExpPosErrorY(z_hodo);
1736 
1737  double x_min = x_mask_min[idx1][idx2] - err_x;
1738  double x_max = x_mask_max[idx1][idx2] + err_x;
1739  double y_min = y_mask_min[idx1][idx2] - err_y;
1740  double y_max = y_mask_max[idx1][idx2] + err_y;
1741 
1742 #ifdef _DEBUG_ON
1743  LogInfo(*iter);
1744  hitAll[*iter].print();
1745  LogInfo(nHodoHits << "/" << stationIDs_mask[tracklet.stationID-1].size() << ": " << z_hodo << " " << x_hodo << " +/- " << err_x << " " << y_hodo << " +/-" << err_y << " : " << x_min << " " << x_max << " " << y_min << " " << y_max);
1746 #endif
1747  if(x_hodo > x_min && x_hodo < x_max && y_hodo > y_min && y_hodo < y_max)
1748  {
1749  nHodoHits++;
1750  masked = true;
1751 
1752  break;
1753  }
1754  }
1755 
1756  if(!masked) return false;
1757  }
1758 
1759 #ifdef _DEBUG_ON
1760  LogInfo(tracklet.stationID << " " << nHodoHits << " " << stationIDs_mask[tracklet.stationID-1].size());
1761 #endif
1762  return true;
1763 }
1764 
1765 bool KalmanDSTrk::muonID_search(Tracklet& tracklet)
1766 {
1767  //Set the cut value on multiple scattering
1768  //multiple scattering: sigma = 0.0136*sqrt(L/L0)*(1. + 0.038*ln(L/L0))/P, L = 1m, L0 = 1.76cm
1769  double cut = 0.03;
1770  if(tracklet.stationID == nStations)
1771  {
1772  double cut_the = MUID_THE_P0*tracklet.invP;
1773  double cut_emp = MUID_EMP_P0 + MUID_EMP_P1/tracklet.invP + MUID_EMP_P2/tracklet.invP/tracklet.invP;
1774  cut = MUID_REJECT*(cut_the > cut_emp ? cut_the : cut_emp);
1775  }
1776 
1777  double slope[2] = {tracklet.tx, tracklet.ty};
1778  double pos_absorb[2] = {tracklet.getExpPositionX(MUID_Z_REF), tracklet.getExpPositionY(MUID_Z_REF)};
1779  PropSegment* segs[2] = {&(tracklet.seg_x), &(tracklet.seg_y)};
1780  for(int i = 0; i < 2; ++i)
1781  {
1782  //this shorting circuting can only be done to X-Z, Y-Z needs more complicated thing
1783  //if(i == 0 && segs[i]->getNHits() > 2 && segs[i]->isValid() && fabs(slope[i] - segs[i]->a) < cut) continue;
1784 
1785  segs[i]->init();
1786  for(int j = 0; j < 4; ++j)
1787  {
1788  int index = detectorIDs_muid[i][j] - nChamberPlanes - 1;
1789  double pos_ref = j < 2 ? pos_absorb[i] : segs[i]->getPosRef(pos_absorb[i] + slope[i]*(z_ref_muid[i][j] - MUID_Z_REF));
1790  double pos_exp = slope[i]*(z_mask[index] - z_ref_muid[i][j]) + pos_ref;
1791 
1792  if(!p_geomSvc->isInPlane(detectorIDs_muid[i][j], tracklet.getExpPositionX(z_mask[index]), tracklet.getExpPositionY(z_mask[index]))) continue;
1793 
1794  double win_tight = cut*(z_mask[index] - z_ref_muid[i][j]);
1795  win_tight = win_tight > 2.54 ? win_tight : 2.54;
1796  double win_loose = win_tight*2;
1797  double dist_min = 1E6;
1798  for(std::list<int>::iterator iter = hitIDs_muid[i][j].begin(); iter != hitIDs_muid[i][j].end(); ++iter)
1799  {
1800  double pos = hitAll[*iter].pos;
1801  double dist = pos - pos_exp;
1802  if(dist < -win_loose) continue;
1803  if(dist > win_loose) break;
1804 
1805  double dist_l = fabs(pos - hitAll[*iter].driftDistance - pos_exp);
1806  double dist_r = fabs(pos + hitAll[*iter].driftDistance - pos_exp);
1807  dist = dist_l < dist_r ? dist_l : dist_r;
1808  if(dist < dist_min)
1809  {
1810  dist_min = dist;
1811  if(dist < win_tight)
1812  {
1813  segs[i]->hits[j].hit = hitAll[*iter];
1814  segs[i]->hits[j].sign = fabs(pos - hitAll[*iter].driftDistance - pos_exp) < fabs(pos + hitAll[*iter].driftDistance - pos_exp) ? -1 : 1;
1815  }
1816  }
1817  }
1818  }
1819  segs[i]->fit();
1820 
1821  //this shorting circuting can only be done to X-Z, Y-Z needs more complicated thing
1822  //if(i == 0 && !(segs[i]->isValid() && fabs(slope[i] - segs[i]->a) < cut)) return false;
1823  }
1824 
1825  muonID_hodoAid(tracklet);
1826  if(segs[0]->getNHits() + segs[1]->getNHits() >= 5)
1827  {
1828  return true;
1829  }
1830  else if(segs[1]->getNHits() == 1 || segs[1]->getNPlanes() == 1)
1831  {
1832  return segs[1]->nHodoHits >= 2;
1833  }
1834  return false;
1835 }
1836 
1837 bool KalmanDSTrk::muonID_comp(Tracklet& tracklet)
1838 {
1839  //Set the cut value on multiple scattering
1840  //multiple scattering: sigma = 0.0136*sqrt(L/L0)*(1. + 0.038*ln(L/L0))/P, L = 1m, L0 = 1.76cm
1841  double cut = 0.03;
1842  if(tracklet.stationID == nStations)
1843  {
1844  double cut_the = MUID_THE_P0*tracklet.invP;
1845  double cut_emp = MUID_EMP_P0 + MUID_EMP_P1/tracklet.invP + MUID_EMP_P2/tracklet.invP/tracklet.invP;
1846  cut = MUID_REJECT*(cut_the > cut_emp ? cut_the : cut_emp);
1847  }
1848 #ifdef _DEBUG_ON
1849  LogInfo("Muon ID cut is: " << cut << " rad.");
1850 #endif
1851 
1852  double slope[2] = {tracklet.tx, tracklet.ty};
1853  PropSegment* segs[2] = {&(tracklet.seg_x), &(tracklet.seg_y)};
1854 
1855  for(int i = 0; i < 2; ++i)
1856  {
1857 #ifdef _DEBUG_ON
1858  if(i == 0) LogInfo("Working in X-Z:");
1859  if(i == 1) LogInfo("Working in Y-Z:");
1860 #endif
1861 
1862  double pos_ref = i == 0 ? tracklet.getExpPositionX(MUID_Z_REF) : tracklet.getExpPositionY(MUID_Z_REF);
1863  if(segs[i]->getNHits() > 2 && segs[i]->isValid() && fabs(slope[i] - segs[i]->a) < cut && fabs(segs[i]->getExpPosition(MUID_Z_REF) - pos_ref) < MUID_R_CUT)
1864  {
1865 #ifdef _DEBUG_ON
1866  LogInfo("Muon ID are already avaiable!");
1867 #endif
1868  continue;
1869  }
1870 
1871  for(std::list<PropSegment>::iterator iter = propSegs[i].begin(); iter != propSegs[i].end(); ++iter)
1872  {
1873 #ifdef _DEBUG_ON
1874  LogInfo("Testing this prop segment, with ref pos = " << pos_ref << ", slope_ref = " << slope[i]);
1875  iter->print();
1876 #endif
1877  if(fabs(iter->a - slope[i]) < cut && fabs(iter->getExpPosition(MUID_Z_REF) - pos_ref) < MUID_R_CUT)
1878  {
1879  *(segs[i]) = *iter;
1880 #ifdef _DEBUG_ON
1881  LogInfo("Accepted!");
1882 #endif
1883  break;
1884  }
1885  }
1886 
1887  if(!segs[i]->isValid()) return false;
1888  }
1889 
1890  if(segs[0]->getNHits() + segs[1]->getNHits() < 5) return false;
1891  return true;
1892 }
1893 
1894 bool KalmanDSTrk::muonID_hodoAid(Tracklet& tracklet)
1895 {
1896  double win = 0.03;
1897  double factor = 5.;
1898  if(tracklet.stationID == nStations)
1899  {
1900  double win_the = MUID_THE_P0*tracklet.invP;
1901  double win_emp = MUID_EMP_P0 + MUID_EMP_P1/tracklet.invP + MUID_EMP_P2/tracklet.invP/tracklet.invP;
1902  win = MUID_REJECT*(win_the > win_emp ? win_the : win_emp);
1903  factor = 3.;
1904  }
1905 
1906  PropSegment* segs[2] = {&(tracklet.seg_x), &(tracklet.seg_y)};
1907  for(int i = 0; i < 2; ++i)
1908  {
1909  segs[i]->nHodoHits = 0;
1910  for(std::list<int>::iterator iter = hitIDs_muidHodoAid[i].begin(); iter != hitIDs_muidHodoAid[i].end(); ++iter)
1911  {
1912  int detectorID = hitAll[*iter].detectorID;
1913  int elementID = hitAll[*iter].elementID;
1914 
1915  int idx1 = detectorID - nChamberPlanes - 1;
1916  int idx2 = elementID - 1;
1917 
1918  double z_hodo = z_mask[idx1];
1919  double x_hodo = tracklet.getExpPositionX(z_hodo);
1920  double y_hodo = tracklet.getExpPositionY(z_hodo);
1921  double err_x = factor*tracklet.getExpPosErrorX(z_hodo) + win*(z_hodo - MUID_Z_REF);
1922  double err_y = factor*tracklet.getExpPosErrorY(z_hodo) + win*(z_hodo - MUID_Z_REF);
1923 
1924  err_x = err_x/(x_mask_max[idx1][idx2] - x_mask_min[idx1][idx2]) > 0.25 ? 0.25*err_x/(x_mask_max[idx1][idx2] - x_mask_min[idx1][idx2]) : err_x;
1925  err_y = err_y/(y_mask_max[idx1][idx2] - y_mask_min[idx1][idx2]) > 0.25 ? 0.25*err_y/(y_mask_max[idx1][idx2] - y_mask_min[idx1][idx2]) : err_y;
1926 
1927  double x_min = x_mask_min[idx1][idx2] - err_x;
1928  double x_max = x_mask_max[idx1][idx2] + err_x;
1929  double y_min = y_mask_min[idx1][idx2] - err_y;
1930  double y_max = y_mask_max[idx1][idx2] + err_y;
1931 
1932  if(x_hodo > x_min && x_hodo < x_max && y_hodo > y_min && y_hodo < y_max)
1933  {
1934  segs[i]->hodoHits[segs[i]->nHodoHits++] = hitAll[*iter];
1935  if(segs[i]->nHodoHits > 4) break;
1936  }
1937  }
1938  }
1939 
1940  return true;
1941 }
1942 
1943 void KalmanDSTrk::buildPropSegments()
1944 {
1945 #ifdef _DEBUG_ON
1946  LogInfo("Building prop. tube segments");
1947 #endif
1948 
1949  for(int i = 0; i < 2; ++i)
1950  {
1951  propSegs[i].clear();
1952 
1953  //note for prop tubes superID index starts from 4
1954  std::list<SRawEvent::hit_pair> pairs_forward = rawEvent->getPartialHitPairsInSuperDetector(superIDs[i+5][0]);
1955  std::list<SRawEvent::hit_pair> pairs_backward = rawEvent->getPartialHitPairsInSuperDetector(superIDs[i+5][1]);
1956 
1957 #ifdef _DEBUG_ON
1958  std::cout << "superID: " << superIDs[i+5][0] << ", " << superIDs[i+5][1] << std::endl;
1959  for(std::list<SRawEvent::hit_pair>::iterator iter = pairs_forward.begin(); iter != pairs_forward.end(); ++iter)
1960  LogInfo("Forward: " << iter->first << " " << iter->second << " " << hitAll[iter->first].index << " " << (iter->second < 0 ? -1 : hitAll[iter->second].index));
1961  for(std::list<SRawEvent::hit_pair>::iterator iter = pairs_backward.begin(); iter != pairs_backward.end(); ++iter)
1962  LogInfo("Backward: " << iter->first << " " << iter->second << " " << hitAll[iter->first].index << " " << (iter->second < 0 ? -1 : hitAll[iter->second].index));
1963 #endif
1964 
1965  for(std::list<SRawEvent::hit_pair>::iterator fiter = pairs_forward.begin(); fiter != pairs_forward.end(); ++fiter)
1966  {
1967 #ifdef _DEBUG_ON
1968  LogInfo("Trying forward pair " << fiter->first << " " << fiter->second);
1969 #endif
1970  for(std::list<SRawEvent::hit_pair>::iterator biter = pairs_backward.begin(); biter != pairs_backward.end(); ++biter)
1971  {
1972 #ifdef _DEBUG_ON
1973  LogInfo("Trying backward pair " << biter->first << " " << biter->second);
1974 #endif
1975 
1976  PropSegment seg;
1977 
1978  //Note that the backward plane comes as the first in pair
1979  if(fiter->first >= 0) seg.hits[1] = SignedHit(hitAll[fiter->first], 0);
1980  if(fiter->second >= 0) seg.hits[0] = SignedHit(hitAll[fiter->second], 0);
1981  if(biter->first >= 0) seg.hits[3] = SignedHit(hitAll[biter->first], 0);
1982  if(biter->second >= 0) seg.hits[2] = SignedHit(hitAll[biter->second], 0);
1983 
1984 #ifdef _DEBUG_ON
1985  seg.print();
1986 #endif
1987  seg.fit();
1988 #ifdef _DEBUG_ON
1989  seg.print();
1990 #endif
1991 
1992  if(seg.isValid())
1993  {
1994  propSegs[i].push_back(seg);
1995  }
1996 #ifdef _DEBUG_ON
1997  else
1998  {
1999  LogInfo("Rejected!");
2000  }
2001 #endif
2002  }
2003  }
2004  }
2005 }
2006 
2007 
2008 int KalmanDSTrk::fitTracklet(Tracklet& tracklet)
2009 {
2010  tracklet_curr = tracklet;
2011 
2012  //idx = 0, using simplex; idx = 1 using migrad
2013  int idx = 1;
2014 #ifdef _ENABLE_MULTI_MINI
2015  if(tracklet.stationID < nStations-1) idx = 0;
2016 #endif
2017 
2018  minimizer[idx]->SetLimitedVariable(0, "tx", tracklet.tx, 0.001, -TX_MAX, TX_MAX);
2019  minimizer[idx]->SetLimitedVariable(1, "ty", tracklet.ty, 0.001, -TY_MAX, TY_MAX);
2020  minimizer[idx]->SetLimitedVariable(2, "x0", tracklet.x0, 0.1, -X0_MAX, X0_MAX);
2021  minimizer[idx]->SetLimitedVariable(3, "y0", tracklet.y0, 0.1, -Y0_MAX, Y0_MAX);
2022  if(p_jobOptsSvc->m_enableKMag)
2023  {
2024  minimizer[idx]->SetLimitedVariable(4, "invP", tracklet.invP, 0.001*tracklet.invP, INVP_MIN, INVP_MAX);
2025  }
2026  minimizer[idx]->Minimize();
2027 
2028  tracklet.tx = minimizer[idx]->X()[0];
2029  tracklet.ty = minimizer[idx]->X()[1];
2030  tracklet.x0 = minimizer[idx]->X()[2];
2031  tracklet.y0 = minimizer[idx]->X()[3];
2032 
2033  tracklet.err_tx = minimizer[idx]->Errors()[0];
2034  tracklet.err_ty = minimizer[idx]->Errors()[1];
2035  tracklet.err_x0 = minimizer[idx]->Errors()[2];
2036  tracklet.err_y0 = minimizer[idx]->Errors()[3];
2037 
2038  if(p_jobOptsSvc->m_enableKMag && tracklet.stationID == nStations)
2039  {
2040  tracklet.invP = minimizer[idx]->X()[4];
2041  tracklet.err_invP = minimizer[idx]->Errors()[4];
2042  }
2043 
2044  tracklet.chisq = minimizer[idx]->MinValue();
2045 
2046  int status = minimizer[idx]->Status();
2047  return status;
2048 }
2049 
2050 int KalmanDSTrk::reduceTrackletList(std::list<Tracklet>& tracklets)
2051 {
2052  std::list<Tracklet> targetList;
2053 
2054  tracklets.sort();
2055  while(!tracklets.empty())
2056  {
2057  targetList.push_back(tracklets.front());
2058  tracklets.pop_front();
2059 
2060 #ifdef _DEBUG_ON_LEVEL_2
2061  LogInfo("Current best tracklet in reduce");
2062  targetList.back().print();
2063 #endif
2064 
2065  for(std::list<Tracklet>::iterator iter = tracklets.begin(); iter != tracklets.end(); )
2066  {
2067  if(iter->similarity(targetList.back()))
2068  {
2069 #ifdef _DEBUG_ON_LEVEL_2
2070  LogInfo("Removing this tracklet: ");
2071  iter->print();
2072 #endif
2073  iter = tracklets.erase(iter);
2074  continue;
2075  }
2076  else
2077  {
2078  ++iter;
2079  }
2080  }
2081  }
2082 
2083  tracklets.assign(targetList.begin(), targetList.end());
2084  return 0;
2085 }
2086 
2087 void KalmanDSTrk::getExtrapoWindowsInSt1(Tracklet& tracklet, double* pos_exp, double* window, int st1ID)
2088 {
2089  if(tracklet.stationID != nStations-1)
2090  {
2091  for(int i = 0; i < 3; i++)
2092  {
2093  pos_exp[i] = 9999.;
2094  window[i] = 0.;
2095  }
2096  return;
2097  }
2098 
2099  for(int i = 0; i < 3; i++)
2100  {
2101  int detectorID = (st1ID-1)*6 + 2*i + 2;
2102  int idx = p_geomSvc->getPlaneType(detectorID) - 1;
2103 
2104  double z_st1 = z_plane[detectorID];
2105  double x_st1 = tracklet.getExpPositionX(z_st1);
2106  double y_st1 = tracklet.getExpPositionY(z_st1);
2107  double err_x = tracklet.getExpPosErrorX(z_st1);
2108  double err_y = tracklet.getExpPosErrorY(z_st1);
2109 
2110  pos_exp[idx] = p_geomSvc->getUinStereoPlane(detectorID, x_st1, y_st1);
2111  window[idx] = 5.*(fabs(costheta_plane[detectorID]*err_x) + fabs(sintheta_plane[detectorID]*err_y));
2112  }
2113 }
2114 
2115 void KalmanDSTrk::getSagittaWindowsInSt1(Tracklet& tracklet, double* pos_exp, double* window, int st1ID)
2116 {
2117  if(tracklet.stationID != nStations-1)
2118  {
2119  for(int i = 0; i < 3; i++)
2120  {
2121  pos_exp[i] = 9999.;
2122  window[i] = 0.;
2123  }
2124  return;
2125  }
2126 
2127  double z_st3 = z_plane[tracklet.hits.back().hit.detectorID];
2128  double x_st3 = tracklet.getExpPositionX(z_st3);
2129  double y_st3 = tracklet.getExpPositionY(z_st3);
2130 
2131  //For U, X, and V planes
2132  for(int i = 0; i < 3; i++)
2133  {
2134  int detectorID = (st1ID-1)*6 + 2*i + 2;
2135  int idx = p_geomSvc->getPlaneType(detectorID) - 1;
2136 
2137  if(!(idx >= 0 && idx <3)) continue;
2138 
2139  double pos_st3 = p_geomSvc->getUinStereoPlane(s_detectorID[idx], x_st3, y_st3);
2140 
2141  double z_st1 = z_plane[detectorID];
2142  double z_st2 = z_plane[s_detectorID[idx]];
2143  double x_st2 = tracklet.getExpPositionX(z_st2);
2144  double y_st2 = tracklet.getExpPositionY(z_st2);
2145  double pos_st2 = p_geomSvc->getUinStereoPlane(s_detectorID[idx], x_st2, y_st2);
2146 
2147  double s2_target = pos_st2 - pos_st3*(z_st2 - Z_TARGET)/(z_st3 - Z_TARGET);
2148  double s2_dump = pos_st2 - pos_st3*(z_st2 - Z_DUMP)/(z_st3 - Z_DUMP);
2149 
2150  double pos_exp_target = SAGITTA_TARGET_CENTER*s2_target + pos_st3*(z_st1 - Z_TARGET)/(z_st3 - Z_TARGET);
2151  double pos_exp_dump = SAGITTA_DUMP_CENTER*s2_dump + pos_st3*(z_st1 - Z_DUMP)/(z_st3 - Z_DUMP);
2152  double win_target = fabs(s2_target*SAGITTA_TARGET_WIN);
2153  double win_dump = fabs(s2_dump*SAGITTA_DUMP_WIN);
2154 
2155  double p_min = std::min(pos_exp_target - win_target, pos_exp_dump - win_dump);
2156  double p_max = std::max(pos_exp_target + win_target, pos_exp_dump + win_dump);
2157 
2158  pos_exp[idx] = 0.5*(p_max + p_min);
2159  window[idx] = 0.5*(p_max - p_min);
2160  }
2161 }
2162 
2163 void KalmanDSTrk::printAtDetectorBack(int stationID, std::string outputFileName)
2164 {
2165  TCanvas c1;
2166 
2167  std::vector<double> x, y, dx, dy;
2168  for(std::list<Tracklet>::iterator iter = trackletsInSt[stationID].begin(); iter != trackletsInSt[stationID].end(); ++iter)
2169  {
2170  double z = p_geomSvc->getPlanePosition(iter->stationID*6);
2171  x.push_back(iter->getExpPositionX(z));
2172  y.push_back(iter->getExpPositionY(z));
2173  dx.push_back(iter->getExpPosErrorX(z));
2174  dy.push_back(iter->getExpPosErrorY(z));
2175  }
2176 
2177  TGraphErrors gr(x.size(), &x[0], &y[0], &dx[0], &dy[0]);
2178  gr.SetMarkerStyle(8);
2179 
2180  //Add detector frames
2181  std::vector<double> x_f, y_f, dx_f, dy_f;
2182  x_f.push_back(p_geomSvc->getPlaneCenterX(stationID*6 + 6));
2183  y_f.push_back(p_geomSvc->getPlaneCenterY(stationID*6 + 6));
2184  dx_f.push_back(p_geomSvc->getPlaneScaleX(stationID*6 + 6)*0.5);
2185  dy_f.push_back(p_geomSvc->getPlaneScaleY(stationID*6 + 6)*0.5);
2186 
2187  if(stationID == 2)
2188  {
2189  x_f.push_back(p_geomSvc->getPlaneCenterX(stationID*6 + 12));
2190  y_f.push_back(p_geomSvc->getPlaneCenterY(stationID*6 + 12));
2191  dx_f.push_back(p_geomSvc->getPlaneScaleX(stationID*6 + 12)*0.5);
2192  dy_f.push_back(p_geomSvc->getPlaneScaleY(stationID*6 + 12)*0.5);
2193  }
2194 
2195  TGraphErrors gr_frame(x_f.size(), &x_f[0], &y_f[0], &dx_f[0], &dy_f[0]);
2196  gr_frame.SetLineColor(kRed);
2197  gr_frame.SetLineWidth(2);
2198  gr_frame.SetFillColor(15);
2199 
2200  c1.cd();
2201  gr_frame.Draw("A2[]");
2202  gr.Draw("Psame");
2203 
2204  c1.SaveAs(outputFileName.c_str());
2205 }
2206 
2207 SRecTrack KalmanDSTrk::processOneTracklet(Tracklet& tracklet)
2208 {
2209  //tracklet.print();
2210  KalmanTrack kmtrk;
2211 
2212  //Set the whole hit and node list
2213  for(std::list<SignedHit>::iterator iter = tracklet.hits.begin(); iter != tracklet.hits.end(); ++iter)
2214  {
2215  if(iter->hit.index < 0) continue;
2216 
2217  Node node_add(*iter);
2218  kmtrk.getNodeList().push_back(node_add);
2219  kmtrk.getHitIndexList().push_back(iter->sign*iter->hit.index);
2220  }
2221 
2222  //Set initial state
2223  TrkPar trkpar_curr;
2224  trkpar_curr._z = p_geomSvc->getPlanePosition(kmtrk.getNodeList().back().getHit().detectorID);
2225  //FIXME Debug Testing: sign reverse
2226  trkpar_curr._state_kf[0][0] = -tracklet.getCharge()*tracklet.invP/sqrt(1. + tracklet.tx*tracklet.tx + tracklet.ty*tracklet.ty);
2227  trkpar_curr._state_kf[1][0] = tracklet.tx;
2228  trkpar_curr._state_kf[2][0] = tracklet.ty;
2229  trkpar_curr._state_kf[3][0] = tracklet.getExpPositionX(trkpar_curr._z);
2230  trkpar_curr._state_kf[4][0] = tracklet.getExpPositionY(trkpar_curr._z);
2231 
2232  trkpar_curr._covar_kf.Zero();
2233  trkpar_curr._covar_kf[0][0] = 0.001;//1E6*tracklet.err_invP*tracklet.err_invP;
2234  trkpar_curr._covar_kf[1][1] = 0.01;//1E6*tracklet.err_tx*tracklet.err_tx;
2235  trkpar_curr._covar_kf[2][2] = 0.01;//1E6*tracklet.err_ty*tracklet.err_ty;
2236  trkpar_curr._covar_kf[3][3] = 100;//1E6*tracklet.getExpPosErrorX(trkpar_curr._z)*tracklet.getExpPosErrorX(trkpar_curr._z);
2237  trkpar_curr._covar_kf[4][4] = 100;//1E6*tracklet.getExpPosErrorY(trkpar_curr._z)*tracklet.getExpPosErrorY(trkpar_curr._z);
2238 
2239  kmtrk.setCurrTrkpar(trkpar_curr);
2240  kmtrk.getNodeList().back().getPredicted() = trkpar_curr;
2241 
2242  //Fit the track first with possibily a few nodes unresolved
2243  if(!fitTrack(kmtrk))
2244  {
2245 #ifdef _DEBUG_ON
2246  LogInfo("!fitTrack(kmtrk) - try flip charge");
2247 #endif
2248  trkpar_curr._state_kf[0][0] *= -1.;
2249  kmtrk.setCurrTrkpar(trkpar_curr);
2250  kmtrk.getNodeList().back().getPredicted() = trkpar_curr;
2251  if(!fitTrack(kmtrk))
2252  {
2253 #ifdef _DEBUG_ON
2254  LogInfo("!fitTrack(kmtrk) - failed flip charge also");
2255 #endif
2256  SRecTrack strack = tracklet.getSRecTrack();
2257  strack.setKalmanStatus(-1);
2258  return strack;
2259  }
2260  }
2261 
2262  if(!kmtrk.isValid()) {
2263 #ifdef _DEBUG_ON
2264  LogInfo("!kmtrk.isValid() Chi2 = " << kmtrk.getChisq() << " - try flip charge");
2265 #endif
2266  trkpar_curr._state_kf[0][0] *= -1.;
2267  kmtrk.setCurrTrkpar(trkpar_curr);
2268  kmtrk.getNodeList().back().getPredicted() = trkpar_curr;
2269  if(!fitTrack(kmtrk))
2270  {
2271 #ifdef _DEBUG_ON
2272  LogInfo("!fitTrack(kmtrk) - failed flip charge also");
2273 #endif
2274  SRecTrack strack = tracklet.getSRecTrack();
2275  strack.setKalmanStatus(-1);
2276  return strack;
2277  }
2278 
2279 #ifdef _DEBUG_ON
2280  LogInfo("Chi2 after flip charge: " << kmtrk.getChisq());
2281  if(kmtrk.isValid()) {
2282  LogInfo("flip charge worked!");
2283  }
2284 #endif
2285  }
2286 
2287 #ifdef _DEBUG_ON
2288  LogInfo("kmtrk.print()");
2289  kmtrk.print();
2290  LogInfo("kmtrk.printNodes()");
2291  kmtrk.printNodes();
2292 #endif
2293 
2294  //Resolve left-right based on the current solution, re-fit if anything changed
2295  //resolveLeftRight(kmtrk);
2296  if(kmtrk.isValid())
2297  {
2298  SRecTrack strack = kmtrk.getSRecTrack();
2299 
2300  //Set trigger road ID
2301  TriggerRoad road(tracklet);
2302  strack.setTriggerRoad(road.getRoadID());
2303 
2304  //Set prop tube slopes
2305  strack.setNHitsInPT(tracklet.seg_x.getNHits(), tracklet.seg_y.getNHits());
2306  strack.setPTSlope(tracklet.seg_x.a, tracklet.seg_y.a);
2307 
2308  strack.setKalmanStatus(1);
2309 
2310  return strack;
2311  }
2312  else
2313  {
2314 #ifdef _DEBUG_ON
2315  LogInfo("!kmtrk.isValid()");
2316 #endif
2317  SRecTrack strack = tracklet.getSRecTrack();
2318  strack.setKalmanStatus(-1);
2319 
2320  return strack;
2321  }
2322 }
2323 
2324 bool KalmanDSTrk::fitTrack(KalmanTrack& kmtrk)
2325 {
2326  if(kmtrk.getNodeList().empty()) return false;
2327 
2328  if(kmfitter->processOneTrack(kmtrk) == 0)
2329  {
2330  return false;
2331  }
2332  kmfitter->updateTrack(kmtrk);
2333 
2334  return true;
2335 }
2336 
2337 void KalmanDSTrk::resolveLeftRight(KalmanTrack& kmtrk)
2338 {
2339  bool isUpdated = false;
2340 
2341  std::list<int>::iterator hitID = kmtrk.getHitIndexList().begin();
2342  for(std::list<Node>::iterator node = kmtrk.getNodeList().begin(); node != kmtrk.getNodeList().end(); )
2343  {
2344  if(*hitID == 0)
2345  {
2346  double x_hit = node->getSmoothed().get_x();
2347  double y_hit = node->getSmoothed().get_y();
2348  double pos_hit = p_geomSvc->getUinStereoPlane(node->getHit().detectorID, x_hit, y_hit);
2349 
2350  int sign = 0;
2351  if(pos_hit > node->getHit().pos)
2352  {
2353  sign = 1;
2354  }
2355  else
2356  {
2357  sign = -1;
2358  }
2359 
2360  //update the node list
2361  TMatrixD m(1, 1), dm(1, 1);
2362  m[0][0] = node->getHit().pos + sign*node->getHit().driftDistance;
2363  dm[0][0] = p_geomSvc->getPlaneResolution(node->getHit().detectorID)*p_geomSvc->getPlaneResolution(node->getHit().detectorID);
2364  node->setMeasurement(m, dm);
2365  *hitID = sign*node->getHit().index;
2366 
2367  isUpdated = true;
2368  }
2369 
2370  ++node;
2371  ++hitID;
2372  }
2373 
2374  if(isUpdated) fitTrack(kmtrk);
2375 }
2376 
2377 void KalmanDSTrk::printTimers() {
2378  std::cout
2379  <<"KalmanDSTrk::printTimers: event: " << _event
2380  << ": " << _timers["event"]->get_accumulated_time()/1000. << " sec" << std::endl;
2381  std::cout << "================================================================" << std::endl;
2382  std::cout << "Tracklet St2 "<<_timers["st2"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2383  std::cout << " Search DB "<<_timers["search_db_2"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2384 
2385  std::cout << "Tracklet St3 "<<_timers["st3"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2386 
2387  std::cout << "Tracklet St23 "<<_timers["st23"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2388  std::cout << " Search DB "<<_timers["search_db_23"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2389  std::cout << " st23_fit1 "<<_timers["st23_fit1"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2390  std::cout << " st23_prop "<<_timers["st23_prop"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2391  std::cout << " st23_fit2 "<<_timers["st23_fit2"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2392  std::cout << " st23_hodo "<<_timers["st23_hodo"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2393  std::cout << " st23_lr40 "<<_timers["st23_lr40"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2394  std::cout << " st23_lr150 "<<_timers["st23_lr150"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2395  std::cout << " st23_rm_hits "<<_timers["st23_rm_hits"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2396  std::cout << " st23_acc "<<_timers["st23_acc"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2397  std::cout << " st23_reduce "<<_timers["st23_reduce"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2398  std::cout << " st23_prop "<<_timers["st23_prop"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2399 
2400  std::cout << "Tracklet Global "<<_timers["global"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2401  std::cout << " Global St1 "<<_timers["global_st1"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2402  std::cout << " Search DB "<<_timers["search_db_glb"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2403  std::cout << " Global Link "<<_timers["global_link"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2404  std::cout << " Link Kalman "<<_timers["global_kalman"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2405  std::cout << "Tracklet Kalman "<<_timers["kalman"]->get_accumulated_time()/1000. << " sec" <<std::endl;
2406  std::cout << "================================================================" << std::endl;
2407 }
2408 
2409 void KalmanDSTrk::chi2fit(int n, double x[], double y[], double& a, double& b)
2410 {
2411  double sum = 0.;
2412  double sx = 0.;
2413  double sy = 0.;
2414  double sxx = 0.;
2415  double syy = 0.;
2416  double sxy = 0.;
2417 
2418  for(int i = 0; i < n; ++i)
2419  {
2420  ++sum;
2421  sx += x[i];
2422  sy += y[i];
2423  sxx += (x[i]*x[i]);
2424  syy += (y[i]*y[i]);
2425  sxy += (x[i]*y[i]);
2426  }
2427 
2428  double det = sum*sxx - sx*sx;
2429  if(fabs(det) < 1E-20)
2430  {
2431  a = 0.;
2432  b = 0.;
2433 
2434  return;
2435  }
2436 
2437  a = (sum*sxy - sx*sy)/det;
2438  b = (sy*sxx - sxy*sx)/det;
2439 }
2440 
2441 
2442 
2443 
2444 
2445 
2446 
2447 
2448 
2449 
2450 
KalmanDSTrk::removeBadHits
void removeBadHits(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1276
KalmanTrack::setCurrTrkpar
void setCurrTrkpar(TrkPar &_trkpar)
set the current track parameter
Definition: KalmanTrack.h:77
Hit::driftDistance
Float_t driftDistance
Definition: SRawEvent.h:81
PatternDB
PatternDB interface objects.
Definition: PatternDB.h:117
SRecTrack::getChisqVertex
Double_t getChisqVertex()
Definition: SRecEvent.h:182
KalmanDSTrk::getExtrapoWindowsInSt1
void getExtrapoWindowsInSt1(Tracklet &tracklet, double *pos_exp, double *window, int st1ID)
Definition: KalmanDSTrk.cxx:2087
SRawEvent::getAllHits
std::vector< Hit > & getAllHits()
Definition: SRawEvent.h:141
PropSegment::hits
SignedHit hits[4]
Definition: FastTracklet.h:118
SRawEvent::getNHitsInD3p
Int_t getNHitsInD3p()
Definition: SRawEvent.cxx:471
KalmanTrack::getSRecTrack
SRecTrack getSRecTrack()
Output to SRecTrack.
Definition: KalmanTrack.cxx:801
KalmanDSTrk::setRawEvent
int setRawEvent(SRawEvent *event_input)
Definition: KalmanDSTrk.cxx:445
Tracklet::calcChisq
double calcChisq()
Definition: FastTracklet.cxx:934
Tracklet::getXZErrorInSt1
void getXZErrorInSt1(double &err_tx_st1, double &err_x0_st1)
Definition: FastTracklet.cxx:772
Hit::detectorID
Short_t detectorID
Definition: SRawEvent.h:78
KalmanDSTrk::muonID_search
bool muonID_search(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1765
nPropPlanes
#define nPropPlanes
Definition: GlobalConsts.h:8
Tracklet::getExpPosErrorX
double getExpPosErrorX(double z) const
Definition: FastTracklet.cxx:637
KalmanTrack::print
void print()
Debugging print.
Definition: KalmanTrack.cxx:820
TFEXIT_FAIL_MULTIPLICITY
#define TFEXIT_FAIL_MULTIPLICITY
Definition: GlobalConsts.h:19
Tracklet::Eval
double Eval(const double *par)
Definition: FastTracklet.cxx:990
Tracklet::residual
double residual[nChamberPlanes]
Definition: FastTracklet.h:247
GeomSvc::getPlaneCenterX
double getPlaneCenterX(int detectorID)
Definition: GeomSvc.h:211
PatternDB::St23
std::set< PartTrackKey > St23
Definition: PatternDB.h:145
GeomSvc::getPlaneCenterY
double getPlaneCenterY(int detectorID)
Definition: GeomSvc.h:212
SRawEvent::getNHitsInD2
Int_t getNHitsInD2()
Definition: SRawEvent.cxx:455
GeomSvc::getUinStereoPlane
double getUinStereoPlane(int detectorID, double x, double y)
Definition: GeomSvc.h:254
KalmanDSTrk::setRawEventDebug
void setRawEventDebug(SRawEvent *event_input)
Definition: KalmanDSTrk.cxx:438
PatternDBUtil::BuildPatternDB
static int BuildPatternDB(const std::string &fin, const std::string &fout, PatternDB &db)
Definition: PatternDBUtil.cxx:67
GeomSvc::getPlaneResolution
double getPlaneResolution(int detectorID) const
Definition: GeomSvc.h:209
PropSegment::isValid
int isValid() const
isValid returns non zero if object contains vailid data
Definition: FastTracklet.cxx:210
PatternDBUtil::GetTrackletKey
static TrackletKey GetTrackletKey(const Tracklet tracklet, const PatternDB::STATION station)
Definition: PatternDBUtil.cxx:452
Tracklet::getSRecTrack
SRecTrack getSRecTrack(bool hyptest=true)
Definition: FastTracklet.cxx:1002
SRawEvent::getNHitsInD1
Int_t getNHitsInD1()
Definition: SRawEvent.cxx:444
Hit::print
void print(std::ostream &os=std::cout) const
Definition: SRawEvent.h:73
KalmanDSTrk::reduceTrackletList
int reduceTrackletList(std::list< Tracklet > &tracklets)
Definition: KalmanDSTrk.cxx:2050
SRecTrack::setKalmanStatus
void setKalmanStatus(Int_t status)
Definition: SRecEvent.h:147
SRawEvent::getNHitsInD3m
Int_t getNHitsInD3m()
Definition: SRawEvent.cxx:482
Tracklet::getMomentum
double getMomentum() const
Definition: FastTracklet.cxx:721
PHField
transient DST object for field storage and access
Definition: PHField.h:13
KalmanDSTrk::muonID_hodoAid
bool muonID_hodoAid(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1894
SRecTrack::isValid
int isValid() const
isValid returns non zero if object contains vailid data
Definition: SRecEvent.cxx:371
KalmanFitter::updateTrack
void updateTrack(KalmanTrack &_track)
Definition: KalmanFitter.cxx:123
NULL
#define NULL
Definition: Pdb.h:9
Tracklet::getExpPositionY
double getExpPositionY(double z) const
Definition: FastTracklet.cxx:657
KalmanDSTrk::muonID_comp
bool muonID_comp(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1837
SRawEvent::getPartialHitPairsInSuperDetector
std::list< SRawEvent::hit_pair > getPartialHitPairsInSuperDetector(Short_t detectorID)
Definition: SRawEvent.cxx:254
Tracklet::getCharge
int getCharge() const
Definition: FastTracklet.cxx:753
nHodoPlanes
#define nHodoPlanes
Definition: GlobalConsts.h:7
Tracklet::print
void print(std::ostream &os=std::cout)
Definition: FastTracklet.cxx:1074
KalmanDSTrk::getSagittaWindowsInSt1
void getSagittaWindowsInSt1(Tracklet &tracklet, double *pos_exp, double *window, int st1ID)
Definition: KalmanDSTrk.cxx:2115
SRawEvent::getHitsIndexInDetectors
std::list< Int_t > getHitsIndexInDetectors(std::vector< Int_t > &detectorIDs)
Definition: SRawEvent.cxx:233
Tracklet::err_y0
double err_y0
Definition: FastTracklet.h:243
TFEXIT_SUCCESS
#define TFEXIT_SUCCESS
Definition: GlobalConsts.h:17
PatternDBUtil::LoadPatternDB
static PatternDB * LoadPatternDB(const std::string &fin)
Definition: PatternDBUtil.cxx:337
PatternDB::ERR_KEY
static const TrackletKey ERR_KEY
Definition: PatternDB.h:124
KalmanTrack::printNodes
void printNodes()
Definition: KalmanTrack.cxx:851
KalmanDSTrk::processOneTracklet
SRecTrack processOneTracklet(Tracklet &tracklet)
Track fitting stuff.
Definition: KalmanDSTrk.cxx:2207
Tracklet::getExpPosErrorY
double getExpPosErrorY(double z) const
Definition: FastTracklet.cxx:662
KalmanDSTrk::buildBackPartialTracks
void buildBackPartialTracks()
Definition: KalmanDSTrk.cxx:681
KalmanDSTrk::resolveLeftRight
void resolveLeftRight(SRawEvent::hit_pair hpair, int &LR1, int &LR2)
Definition: KalmanDSTrk.cxx:1376
Tracklet::hits
std::list< SignedHit > hits
Definition: FastTracklet.h:227
PropSegment::getNHits
int getNHits() const
Definition: FastTracklet.cxx:190
PropSegment::print
void print(std::ostream &os=std::cout) const
Definition: FastTracklet.cxx:148
TrkPar::_z
double _z
Definition: KalmanUtil.h:87
SRecTrack::setPTSlope
void setPTSlope(Double_t slopeX, Double_t slopeY)
Definition: SRecEvent.h:224
KalmanDSTrk::setRawEventWorker
int setRawEventWorker(SRawEvent *event_input)
Definition: KalmanDSTrk.cxx:482
SRawEvent::getNRoadsNeg
Int_t getNRoadsNeg()
Definition: SRawEvent.h:184
KalmanDSTrk::fitTracklet
int fitTracklet(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:2008
SRawEvent::hit_pair
std::pair< Int_t, Int_t > hit_pair
Type of pair with two adjacent wires.
Definition: SRawEvent.h:169
nStations
#define nStations
Definition: GlobalConsts.h:5
Tracklet::isValid
int isValid() const
isValid returns non zero if object contains vailid data
Definition: FastTracklet.cxx:506
Tracklet::chisq
double chisq
Definition: FastTracklet.h:221
Fun4AllBase::VERBOSITY_A_LOT
Output a lot of messages.
Definition: Fun4AllBase.h:48
KalmanTrack::getNodeList
std::list< Node > & getNodeList()
Definition: KalmanTrack.h:84
Tracklet::chisq_vtx
double chisq_vtx
Definition: FastTracklet.h:224
KalmanDSTrk::Verbosity
int Verbosity() const
Definition: KalmanDSTrk.h:60
Tracklet::seg_x
PropSegment seg_x
Definition: FastTracklet.h:230
KalmanDSTrk::acceptEvent
bool acceptEvent(SRawEvent *rawEvent)
Definition: KalmanDSTrk.cxx:643
KalmanDSTrk::buildPropSegments
void buildPropSegments()
Definition: KalmanDSTrk.cxx:1943
Hit::index
Int_t index
Definition: SRawEvent.h:77
Tracklet::err_tx
double err_tx
Definition: FastTracklet.h:240
Tracklet::seg_y
PropSegment seg_y
Definition: FastTracklet.h:231
PatternDB::St3
std::set< TrackletKey > St3
Definition: PatternDB.h:144
GeomSvc::isInKMAG
bool isInKMAG(double x, double y)
Definition: GeomSvc.cxx:773
TrkPar::_state_kf
TMatrixD _state_kf
State vectors and its covariance.
Definition: KalmanUtil.h:85
KalmanDSTrk::printTimers
void printTimers()
Definition: KalmanDSTrk.cxx:2377
SRawEvent::getNRoadsPos
Int_t getNRoadsPos()
Definition: SRawEvent.h:183
Tracklet::err_invP
double err_invP
Definition: FastTracklet.h:244
TFEXIT_FAIL_GLOABL
#define TFEXIT_FAIL_GLOABL
Definition: GlobalConsts.h:24
KalmanDSTrk::acceptTracklet
bool acceptTracklet(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1684
PropSegment::getPosRef
double getPosRef(double default_val=-9999.)
Definition: FastTracklet.cxx:173
SRecTrack::setNHitsInPT
void setNHitsInPT(Int_t nHitsX, Int_t nHitsY)
Definition: SRecEvent.h:225
GeomSvc::getPlaneType
int getPlaneType(int detectorID) const
Definition: GeomSvc.h:223
GeomSvc::getPlaneScaleX
double getPlaneScaleX(int detectorID)
Definition: GeomSvc.h:204
KalmanTrack::getHitIndexList
std::list< int > & getHitIndexList()
Get the list of hits associated.
Definition: KalmanTrack.h:83
TFEXIT_FAIL_ST2_TRACKLET
#define TFEXIT_FAIL_ST2_TRACKLET
Definition: GlobalConsts.h:21
KalmanDSTrk::fitTrack
bool fitTrack(KalmanTrack &kmtrk)
Definition: KalmanDSTrk.cxx:2324
KalmanDSTrk::buildGlobalTracks
void buildGlobalTracks()
Definition: KalmanDSTrk.cxx:899
KalmanDSTrk::printAtDetectorBack
void printAtDetectorBack(int stationID, std::string outputFileName)
Definition: KalmanDSTrk.cxx:2163
GeomSvc::getPlaneScaleY
double getPlaneScaleY(int detectorID)
Definition: GeomSvc.h:205
SRawEvent::getNHitsInDetectors
Int_t getNHitsInDetectors(std::vector< Int_t > &detectorIDs)
Definition: SRawEvent.cxx:414
Tracklet::err_x0
double err_x0
Definition: FastTracklet.h:242
KalmanDSTrk::resolveSingleLeftRight
void resolveSingleLeftRight(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1253
KalmanDSTrk::KalmanDSTrk
KalmanDSTrk(const PHField *field, const TGeoManager *geom, bool enable_KF=true, int DS_level=KalmanDSTrk::NO_DS, const std::string sim_db_name="", const std::string pattern_db_name="")
Definition: KalmanDSTrk.cxx:36
Tracklet::y0
double y0
Definition: FastTracklet.h:237
KalmanDSTrk::buildTrackletsInStation
void buildTrackletsInStation(int stationID, int listID, double *pos_exp=NULL, double *window=NULL)
Tracklet finding stuff.
Definition: KalmanDSTrk.cxx:1416
KalmanTrack::getChisq
double getChisq()
Definition: KalmanTrack.h:91
nChamberPlanes
#define nChamberPlanes
Definition: GlobalConsts.h:6
SRawEvent::getNHitsInD0
Int_t getNHitsInD0()
Definition: SRawEvent.cxx:433
GeomSvc::instance
static GeomSvc * instance()
singlton instance
Definition: GeomSvc.cxx:211
Tracklet::err_ty
double err_ty
Definition: FastTracklet.h:241
PatternDB::St123
std::set< GlobTrackKey > St123
Definition: PatternDB.h:146
TFEXIT_FAIL_BACKPARTIAL
#define TFEXIT_FAIL_BACKPARTIAL
Definition: GlobalConsts.h:23
PropSegment::hodoHits
Hit hodoHits[4]
Definition: FastTracklet.h:115
TrkPar::_covar_kf
TMatrixD _covar_kf
Definition: KalmanUtil.h:86
LogInfo
#define LogInfo(message)
Definition: DbSvc.cc:14
Tracklet::ty
double ty
Definition: FastTracklet.h:235
TFEXIT_FAIL_ST3_TRACKLET
#define TFEXIT_FAIL_ST3_TRACKLET
Definition: GlobalConsts.h:22
Tracklet::invP
double invP
Definition: FastTracklet.h:238
PatternDB::St2
std::set< TrackletKey > St2
Definition: PatternDB.h:143
Tracklet::stationID
int stationID
Definition: FastTracklet.h:213
PatternDB::St1
std::set< TrackletKey > St1
Definition: PatternDB.h:142
KalmanDSTrk::chi2fit
void chi2fit(int n, double x[], double y[], double &a, double &b)
Definition: KalmanDSTrk.cxx:2409
SRecTrack::setTriggerRoad
void setTriggerRoad(Int_t roadID)
Definition: SRecEvent.h:220
SRawEvent::getNPropHitsAll
Int_t getNPropHitsAll()
Definition: SRawEvent.cxx:403
PHTimer.h
high precision timer
KalmanDSTrk::hodoMask
bool hodoMask(Tracklet &tracklet)
Definition: KalmanDSTrk.cxx:1714
GeomSvc::isInPlane
bool isInPlane(int detectorID, double x, double y)
See if a point is in a plane.
Definition: GeomSvc.cxx:752
KalmanTrack::isValid
bool isValid()
Self check to see if it is null.
Definition: KalmanTrack.cxx:96
Tracklet::getExpPositionX
double getExpPositionX(double z) const
Definition: FastTracklet.cxx:622
SRawEvent::getHitsIndexInDetector
std::list< Int_t > getHitsIndexInDetector(Short_t detectorID)
Gets.
Definition: SRawEvent.cxx:186
Tracklet::tx
double tx
Definition: FastTracklet.h:234
SRawEvent::isFPGATriggered
bool isFPGATriggered()
Definition: SRawEvent.cxx:639
PHTimer
high precision timer
Definition: PHTimer.h:24
Tracklet::merge
Tracklet merge(Tracklet &elem)
Definition: FastTracklet.cxx:875
GeomSvc::getPlanePosition
double getPlanePosition(int detectorID) const
Definition: GeomSvc.h:197
Tracklet::x0
double x0
Definition: FastTracklet.h:236
Tracklet::sortHits
void sortHits()
Definition: FastTracklet.h:142
Tracklet::getXZInfoInSt1
void getXZInfoInSt1(double &tx_st1, double &x0_st1)
Definition: FastTracklet.cxx:758
KalmanFitter::processOneTrack
int processOneTrack(KalmanTrack &_track)
Definition: KalmanFitter.cxx:51
c1
TCanvas * c1
Definition: Fun4SimTree.C:5