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