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