StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiKalmanTrackNode.cxx
1 //StiKalmanTrack.cxx
2 /*
3  * $Id: StiKalmanTrackNode.cxx,v 2.182 2018/11/10 00:22:03 perev Exp $
4  *
5  * /author Claude Pruneau
6  *
7  * $Log: StiKalmanTrackNode.cxx,v $
8  * Revision 2.182 2018/11/10 00:22:03 perev
9  * 1. implementation of StiKalmanTrackNode::initialize(const double dirg[3])
10  * 2. Replace Step to Path in THelixTrack call
11  *
12  * Revision 2.181 2018/06/29 21:46:27 smirnovd
13  * Revert iTPC-related changes committed on 2018-06-20 through 2018-06-28
14  *
15  * Revert "NoDead option added"
16  * Revert "Fill mag field more carefully"
17  * Revert "Assert commented out"
18  * Revert "Merging with TPC group code"
19  * Revert "Remove too strong assert"
20  * Revert "Restore removed by mistake line"
21  * Revert "Remove not used anymore file"
22  * Revert "iTPCheckIn"
23  *
24  * Revision 2.178 2018/04/10 11:38:34 smirnovd
25  * Replace thrown exceptions with runtime asserts
26  *
27  * Revision 2.177 2018/04/10 11:32:09 smirnovd
28  * Minor corrections across multiple files
29  *
30  * - Remove ClassImp macro
31  * - Change white space
32  * - Correct windows newlines to unix
33  * - Remove unused debugging
34  * - Correct StTpcRTSHitMaker header guard
35  * - Remove unused preprocessor directives in StiCA
36  * - Minor changes in status and debug print out
37  * - Remove using std namespace from StiKalmanTrackFinder
38  * - Remove includes for unused headers
39  *
40  * Revision 2.176 2018/04/10 11:31:24 smirnovd
41  * Remove dead code
42  *
43  * Revision 2.175 2018/01/16 22:46:09 smirnovd
44  * Remove inline attribute to match the declaration
45  *
46  * Let compiler decide whether to inline or not
47  *
48  * Revision 2.174 2016/11/07 23:58:03 perev
49  * More accurate tracking when in refit track sometimes missed the vollume.
50  * It is related to bug #3243
51  *
52  * Revision 2.173 2016/07/08 16:11:32 perev
53  * Method print enhancened
54  *
55  * Revision 2.172 2016/06/29 18:37:39 perev
56  * 1. Removed debug codes non actual now.
57  *
58  * Revision 2.169.2.1 2016/06/02 16:45:42 smirnovd
59  * Squashed changes on MAIN branch after StiCA_2016 was brached off
60  *
61  * commit 0b534582b5bf40a64870088f6864387a7941a9be
62  * Author: perev <perev>
63  * Date: Tue May 31 17:11:46 2016 +0000
64  *
65  * Coverity
66  *
67  * commit cbfeeef5e8f9a6e24ddd7329ff5770086e535493
68  * Author: perev <perev>
69  * Date: Tue Apr 19 01:58:39 2016 +0000
70  *
71  * Assignment out of array boundary removed(J.Lauret)
72  *
73  * commit a49f5f23dc613c1ee8ab61c543e713f776d3c7fe
74  * Author: perev <perev>
75  * Date: Tue Apr 19 01:37:22 2016 +0000
76  *
77  * WarnOff
78  *
79  * commit 48ca225cc052db66cd8a3934f15c46345c9862c6
80  * Author: perev <perev>
81  * Date: Fri Apr 15 20:47:42 2016 +0000
82  *
83  * Warnoff
84  *
85  * commit b1b0f73cef0f5675bd84106241067329e0221079
86  * Author: perev <perev>
87  * Date: Fri Apr 15 20:13:06 2016 +0000
88  *
89  * Warnoff
90  *
91  * commit 393adde57febc06a90d054f71e621e8efd082e10
92  * Author: perev <perev>
93  * Date: Wed Apr 13 23:08:44 2016 +0000
94  *
95  * -opt2 proble solved. Array A[1] removed
96  *
97  * commit 1c105bdc0cbde40ccec63fdbf40e79dfb3e7f0e0
98  * Author: perev <perev>
99  * Date: Mon Mar 28 00:17:55 2016 +0000
100  *
101  * 1st hit must be not used at all
102  *
103  * commit 1eca42192ef93788d149625ecebc8390f8b0bc3a
104  * Author: perev <perev>
105  * Date: Mon Mar 28 00:15:53 2016 +0000
106  *
107  * Add max number of tracks assigned to one hit
108  *
109  * commit b349ba99342bc38eaa82f3d2a8d25aa29ba73c29
110  * Author: genevb <genevb>
111  * Date: Thu Feb 25 23:04:50 2016 +0000
112  *
113  * kSsdId => kSstId
114  *
115  * commit a06d8162931b223b4a405ea5714e703b1cad14e3
116  * Author: perev <perev>
117  * Date: Mon Dec 28 23:50:27 2015 +0000
118  *
119  * Remove assert temporary
120  *
121  * commit f8646d17ed86b9be5b5fa940691f9871346a5ee2
122  * Author: perev <perev>
123  * Date: Mon Dec 21 19:41:31 2015 +0000
124  *
125  * bug #3166 assert vertex closer to 0,0 <9 removed
126  *
127  * commit 48a6813db30f593a90a79beb688c27d0e8946bfa
128  * Author: perev <perev>
129  * Date: Sat Dec 19 03:40:50 2015 +0000
130  *
131  * assert rxy<4 ==> <9 temporary
132  *
133  * commit d49576f25ba887ba4ff82c3bf1ffcc760c8da6b2
134  * Author: perev <perev>
135  * Date: Fri Dec 18 03:50:06 2015 +0000
136  *
137  * *** empty log message ***
138  *
139  * commit 23e9c0447bd41151e45728a6f4dd3cc554be1cfb
140  * Author: perev <perev>
141  * Date: Thu Dec 3 19:12:24 2015 +0000
142  *
143  * Remove redundant GTrack error: mFlag: is Negative
144  *
145  * Revision 2.171 2016/04/15 20:13:06 perev
146  * Warnoff
147  *
148  * Revision 2.170 2016/04/13 23:08:44 perev
149  * -opt2 proble solved. Array A[1] removed
150  *
151  * Revision 2.169 2015/07/29 01:28:05 smirnovd
152  * Added std:: to resolve ambiguity for isnan in g++ (4.8.2)
153  *
154  * Also switched to standard c++ header <cmath>
155  *
156  * Revision 2.168 2015/04/09 22:54:40 perev
157  * new method evaluateChi2Info added to recalculate Xi2 for testing only
158  *
159  * Revision 2.167 2015/02/25 20:10:20 perev
160  * In StiKalmanTrackNode::propagateError() recov(1) is called when length is bigger kBigLen
161  *
162  * Revision 2.166 2015/02/25 02:35:26 perev
163  * Some outdated asserts and debug lines removed
164  * In propagateError recov is colled
165  * When length is big recov(1) is called, otherwice recov(0).
166  * recov(1) is slower and negative error matrix happens with big length,
167  * about 100.
168  *
169  * Revision 2.165 2015/02/23 19:54:06 perev
170  * Bug #3048 fixed. Very starnge mistype. instead of if(something) it was if(!something)
171  * How it was happened I have no idea (VP)
172  *
173  * Revision 2.164 2015/02/21 04:48:03 perev
174  * All asserts with sign() replaced for zign() for speedup
175  * Some outdated asserts removed
176  * Some outdated debug prints removed as well
177  *
178  * Revision 2.163 2015/02/17 01:42:34 perev
179  * Fix bug #3034. It is the conisedence of -1 (return local()) and kEnded.
180  * Now if (local()==-1) volume is missed, return kFailed
181  *
182  * Revision 2.162 2015/02/16 22:12:46 perev
183  * Roll back all changes related to nudge problem wall x in particular
184  *
185  * Revision 2.153 2015/01/15 20:05:27 perev
186  * Simplified check of simplified locate()
187  *
188  * Revision 2.152 2015/01/15 19:23:26 perev
189  * Method locate() simplified. Redundunt info removed.
190  * For instance, which part of detector track missed. This is not used anyway.
191  * Some debug added/removed
192  *
193  * Revision 2.151 2014/11/10 21:48:03 perev
194  * Zero field accounting using isZeroH(0 methot
195  *
196  * Revision 2.150 2014/11/03 20:53:08 perev
197  * For the zero field defined minimum non zero field eqaal 1GeV radius 1km
198  * Such notation was before, but because zero field is not used too often
199  * it was disappeared. Now fixed again
200  *
201  * Revision 2.149 2014/10/30 15:03:54 jeromel
202  * Reverted to Oct 2nd
203  *
204  * Revision 2.142 2014/09/30 15:44:51 perev
205  * Added StELoss class to keep ELoss info
206  *
207  * Revision 2.141 2014/09/18 18:45:00 perev
208  * Debug++
209  * Using new cylCross method
210  * More checks for out of region
211  *
212  * Revision 2.140 2014/09/05 21:55:29 perev
213  * bug #2903 fixed. x0,x0p,x0Gas initialised now tp 1e11 instead of -1
214  * Many asserts adde temporary. Some of them time consuming
215  *
216  * Revision 2.139 2014/08/27 01:33:59 perev
217  * ::print bug fixed (printed local coordinates instead of global ones)
218  * added print of rxy and direction of track, outside +ve, inside -ve
219  *
220  * Revision 2.138 2014/08/22 16:25:20 perev
221  * Fix old bug double counting of density
222  * Fix ELoss bug. dEdX(density) ==> dEdX(density,material)
223  * Save calculated ELoss in StiNode for technical analisys
224  *
225  * Revision 2.137 2014/06/03 16:48:38 genevb
226  * Reduce visibility of inactive materials
227  *
228  * Revision 2.136 2013/04/10 22:09:01 fisyak
229  * Roll back to version 04/04/2013
230  *
231  * Revision 2.134 2013/01/14 22:19:48 fisyak
232  * Set Bz = 0 for laser tracks
233  *
234  * Revision 2.133 2011/11/21 17:05:26 fisyak
235  * Correct no. of possible point for CA case
236  *
237  * Revision 2.132 2010/09/07 18:37:31 fisyak
238  * Restore Sti logic before TPCCATracker
239  *
240  * Revision 2.131 2010/09/06 18:20:48 fisyak
241  * Add TPCCATracker
242  *
243  * Revision 1.4 2010/08/23 21:56:28 ikulakov
244  * Fix - alinghment bag.
245  *
246  * Revision 1.3 2010/08/04 19:46:43 ikulakov
247  * Fix - edge (don't cut 2 cm from detector)
248  *
249  * Revision 1.2 2010/07/29 16:19:11 fisyak
250  * GSI CA tracker
251  *
252  * Revision 2.130 2010/07/14 18:45:10 perev
253  * MoreComments
254  *
255  * Revision 2.129 2010/04/03 04:02:29 perev
256  * Account field=0
257  *
258  * Revision 2.128 2010/04/01 20:25:34 perev
259  * Remove test to small error _cPP
260  *
261  * Revision 2.127 2010/04/01 00:27:33 perev
262  * Zero fied = 2e-6
263  *
264  * Revision 2.126 2009/11/05 17:37:52 fine
265  * remove the compilation warnings
266  *
267  * Revision 2.125 2009/10/18 22:48:58 perev
268  * remove STAR LOG in print()
269  *
270  * Revision 2.124 2009/08/19 21:19:37 perev
271  * getTime() is a const now
272  *
273  * Revision 2.123 2009/04/23 02:39:03 perev
274  * GetTime defence sin <1 added
275  *
276  * Revision 2.122 2009/04/01 19:20:17 perev
277  * Replace asserts to error condition
278  *
279  * Revision 2.121 2009/03/16 13:50:15 fisyak
280  * Move out all Sti Chairs into StDetectorDb
281  *
282  * Revision 2.120 2008/12/26 15:18:00 fisyak
283  * Enlarge fitting volume from 200 => 250 cm
284  *
285  * Revision 2.119 2008/06/11 22:04:37 fisyak
286  * Add dead material
287  *
288  * Revision 2.118 2008/06/09 20:12:09 perev
289  * BigStepBack
290  *
291  * Revision 2.115 2008/04/03 20:03:36 fisyak
292  * Straighten out DB access via chairs
293  *
294  * Revision 2.114 2008/03/25 18:02:53 perev
295  * remove field field from everythere
296  *
297  * Revision 2.113 2007/09/10 21:26:52 perev
298  * getPt non positive bug fix. introduces 3 month ago
299  *
300  * Revision 2.112 2007/08/30 19:13:27 fine
301  * replace the repmaining cout with LOG_DEBUG
302  *
303  * Revision 2.111 2007/08/16 20:21:24 fine
304  * replace printf with logger
305  *
306  * Revision 2.110 2007/07/12 00:21:00 perev
307  * Normal radius instead of layer one
308  *
309  * Revision 2.109 2007/06/25 19:31:52 perev
310  * Init of _sinCA and _cosCA non zeros now
311  *
312  * Revision 2.108 2007/06/07 20:13:42 perev
313  * BugFix in getPt()
314  *
315  * Revision 2.107 2007/06/06 04:03:03 perev
316  * getTime() cleanup
317  *
318  * Revision 2.106 2007/04/30 19:53:16 fisyak
319  * Correct time of flight calculation, add time of flight corrrection for Laser
320  *
321  * Revision 2.105 2007/03/21 17:47:24 fisyak
322  * Time of Flight
323  *
324  * Revision 2.104 2006/12/24 02:16:36 perev
325  * _inf=0 added in copy constructor
326  *
327  * Revision 2.103 2006/12/18 01:17:41 perev
328  * Info block added and filled for pulls
329  *
330  * Revision 2.102 2006/10/16 20:29:35 fisyak
331  * Clean up useless classes
332  *
333  * Revision 2.101 2006/10/09 15:47:07 fisyak
334  * take out Central represantation, remove StiDedxCalculator
335  *
336  * Revision 2.100 2006/05/31 03:58:06 fisyak
337  * Add Victor's dca track parameters, clean up
338  *
339  * Revision 2.99 2006/04/15 23:12:10 perev
340  * Supress printout
341  *
342  * Revision 2.98 2006/04/07 18:01:56 perev
343  * Back to the latest Sti
344  *
345  * Revision 2.96 2006/02/16 20:44:50 perev
346  * assert changed back
347  *
348  * Revision 2.95 2006/02/15 19:07:18 perev
349  * assrt in nudge cos < 1
350  *
351  * Revision 2.94 2006/02/14 18:10:41 perev
352  * getGlobalHitErrors added+CleanUp
353  *
354  * Revision 2.93 2005/12/31 01:37:12 perev
355  * Primary node perpendicular to track
356  *
357  * Revision 2.92 2005/12/20 00:41:21 perev
358  * unassigned sind fixed(thanxYF)
359  *
360  * Revision 2.91 2005/12/18 23:41:46 perev
361  * Dependency from StiKalmanTrackNode removed
362  *
363  * Revision 2.90 2005/12/08 22:05:45 perev
364  * nudge assert replaced by print. But very strangeStiKalmanTrackNode.cxx
365  *
366  * Revision 2.89 2005/12/07 22:29:27 perev
367  * Big Cleanup
368  *
369  * Revision 2.88 2005/08/09 14:55:34 perev
370  * extend()/reduce() of node
371  *
372  * Revision 2.87 2005/08/04 03:52:54 perev
373  * Cleanup
374  *
375  * Revision 2.86 2005/07/20 17:24:25 perev
376  * Nudge actions in evaluateChi2 added
377  *
378  * Revision 2.85 2005/06/14 22:22:46 perev
379  * Replace assert to error return
380  *
381  * Revision 2.84 2005/06/03 19:57:04 perev
382  * Bug fix, violation of array size
383  *
384  * Revision 2.83 2005/06/02 17:27:41 perev
385  * More weak assert in nudge()
386  *
387  * Revision 2.82 2005/05/31 16:47:56 perev
388  * technical reorganization
389  *
390  * Revision 2.81 2005/05/12 18:10:04 perev
391  * dL/dCurv more accurate
392  *
393  * Revision 2.80 2005/05/04 19:33:00 perev
394  * Supress assert
395  *
396  * Revision 2.79 2005/04/30 20:45:18 perev
397  * Less strong test for assert in propagateError
398  *
399  * Revision 2.78 2005/04/25 20:20:25 fisyak
400  * replace assert by print out
401  *
402  * Revision 2.77 2005/04/12 14:35:39 fisyak
403  * Add print out for dE/dx
404  *
405  * Revision 2.76 2005/04/11 22:48:30 perev
406  * assert removed
407  *
408  * Revision 2.75 2005/04/11 17:33:55 perev
409  * Wrong sorting accounted, check for accuracy inctreased
410  *
411  * Revision 2.74 2005/04/11 14:32:18 fisyak
412  * Use gdrelx from GEANT for dE/dx calculation with accouning density effect
413  *
414  * Revision 2.73 2005/03/30 21:01:43 perev
415  * asserts replaced to prints
416  *
417  * Revision 2.72 2005/03/28 05:52:40 perev
418  * Reorganization of node container
419  *
420  * Revision 2.71 2005/03/24 19:28:35 perev
421  * Switch off DerivTest
422  *
423  * Revision 2.70 2005/03/24 18:05:07 perev
424  * Derivatives and their test fixed to eta==Psi model
425  *
426  * Revision 2.69 2005/03/19 00:20:33 perev
427  * Assert for zero determinant ==> print
428  *
429  * Revision 2.68 2005/03/18 17:35:38 perev
430  * some asserts removed
431  *
432  * Revision 2.67 2005/03/18 17:13:07 perev
433  * assert in rotate fix
434  *
435  * Revision 2.66 2005/03/17 06:24:52 perev
436  * A lot of changes. _eta now is Psi
437  *
438  * Revision 2.65 2005/02/25 17:05:41 perev
439  * Scaling for errors added
440  *
441  * Revision 2.64 2005/02/19 20:23:37 perev
442  * Cleanup
443  *
444  * Revision 2.63 2005/02/18 19:02:55 fisyak
445  * Add debug print out for extendToVertex
446  *
447  * Revision 2.62 2005/02/17 23:19:02 perev
448  * NormalRefangle + Error reseting
449  *
450  * Revision 2.61 2005/02/17 19:58:06 fisyak
451  * Add debug print out flags
452  *
453  * Revision 2.60 2005/02/16 17:47:16 perev
454  * assert in nudge 1==>5
455  *
456  * Revision 2.59 2005/02/07 18:33:42 fisyak
457  * Add VMC dead material
458  *
459  * Revision 2.58 2005/01/20 16:51:32 perev
460  * Remove redundant print
461  *
462  * Revision 2.57 2005/01/17 01:31:25 perev
463  * New parameter model
464  *
465  * Revision 2.56 2005/01/06 00:59:41 perev
466  * Initial errors tuned
467  *
468  * Revision 2.55 2005/01/04 01:37:47 perev
469  * minor bug fix
470  *
471  * Revision 2.54 2004/12/23 18:15:46 perev
472  * Cut for -ve cosCA added
473  *
474  * Revision 2.53 2004/12/14 17:10:17 perev
475  * Propagate for 0 not called
476  *
477  * Revision 2.52 2004/12/13 22:52:23 perev
478  * Off testError
479  *
480  * Revision 2.51 2004/12/13 20:01:38 perev
481  * old version of testError temporary activated
482  *
483  * Revision 2.50 2004/12/12 01:34:24 perev
484  * More smart testError, partial error reset
485  *
486  * Revision 2.49 2004/12/11 22:17:49 pruneau
487  * new eloss calculation
488  *
489  * Revision 2.48 2004/12/11 04:31:36 perev
490  * set of bus fixed
491  *
492  * Revision 2.47 2004/12/10 15:51:44 fisyak
493  * Remove fudge factor from eloss calculation, add more debug printout and tests, reorder calculation of cov. matrix for low triangular form
494  *
495  * Revision 2.46 2004/12/08 16:56:16 fisyak
496  * Fix sign in dE/dx; move from upper to lower triangular matrix convention (StEvent) for px,py,pz
497  *
498  * Revision 2.45 2004/12/05 00:39:07 fisyak
499  * Add test suit for matrix manipulation debugging under overall CPPFLAGS=-DSti_DEBUG
500  *
501  * Revision 2.44 2004/12/01 14:04:57 pruneau
502  * z propagation fix
503  *
504  * Revision 2.43 2004/11/24 17:59:26 fisyak
505  * Set ionization potential for Ar in eloss calculateion instead 5
506  *
507  * Revision 2.42 2004/11/22 19:43:06 pruneau
508  * commented out offending cout statement
509  *
510  * Revision 2.41 2004/11/22 19:23:20 pruneau
511  * minor changes
512  *
513  * Revision 2.40 2004/11/10 21:46:02 pruneau
514  * added extrapolation function; minor change to updateNode function
515  *
516  * Revision 2.39 2004/11/08 15:32:54 pruneau
517  * 3 sets of modifications
518  * (1) Changed the StiPlacement class to hold keys to both the radial and angle placement. Propagated the use
519  * of those keys in StiSvt StiTpc StiSsd and all relevant Sti classes.
520  * (2) Changed the StiKalmanTrackFinder::find(StiTrack*) function's algorithm for the navigation of the
521  * detector volumes. The new code uses an iterator to visit all relevant volumes. The code is now more robust and compact
522  * as well as much easier to read and maintain.
523  * (3) Changed the chi2 calculation in StiKalmanTrack::getChi2 and propagated the effects of this change
524  * in both StiTrackingPlots and StiStEventFiller classes.
525  *
526  * Revision 2.38 2004/10/27 03:25:49 perev
527  * Version V3V
528  *
529  * Revision 2.37 2004/10/26 21:53:23 pruneau
530  * No truncation but bad hits dropped
531  *
532  * Revision 2.36 2004/10/26 06:45:37 perev
533  * version V2V
534  *
535  * Revision 2.35 2004/10/25 14:15:56 pruneau
536  * various changes to improve track quality.
537  *
538  * Revision 2.34 2004/03/24 22:01:07 pruneau
539  * Removed calls to center representation and replaced by normal representation
540  *
541  * Revision 2.33 2004/03/17 21:01:53 andrewar
542  * Trapping for negative track error (^2) values _cYY and _cZZ. This should
543  * be a temporary fix until the root of the problem is found. Problem seems
544  * localized to trackNodes without hits.
545  * Also trapping for asin(x), x>1 in ::length; point to point cord length
546  * on the helix is greater than twice radius of curvature. This should also be
547  * resovled.
548  *
549  * Revision 2.32 2004/01/30 21:40:21 pruneau
550  * some clean up of the infinite checks
551  *
552  * Revision 2.31 2003/09/02 17:59:41 perev
553  * gcc 3.2 updates + WarnOff
554  *
555  * Revision 2.30 2003/08/13 21:04:21 pruneau
556  * transfered relevant tracking pars to detector builders
557  *
558  * Revision 2.29 2003/08/02 08:23:10 pruneau
559  * best performance so far
560  *
561  * Revision 2.28 2003/07/30 19:18:58 pruneau
562  * sigh
563  *
564  * Revision 2.26 2003/07/15 13:56:19 andrewar
565  * Revert to previous version to remove bug.
566  *
567  * Revision 2.24 2003/05/22 18:42:33 andrewar
568  * Changed max eloss correction from 1% to 10%.
569  *
570  * Revision 2.23 2003/05/09 22:07:57 pruneau
571  * Added protection to avoid 90deg tracks and ill defined eloss
572  *
573  * Revision 2.22 2003/05/09 14:57:20 pruneau
574  * Synching
575  *
576  * Revision 2.21 2003/05/08 18:49:09 pruneau
577  * fudge=1
578  *
579  * Revision 2.20 2003/05/07 03:01:39 pruneau
580  * *** empty log message ***
581  *
582  * Revision 2.19 2003/05/03 14:37:22 pruneau
583  * *** empty log message ***
584  *
585  * Revision 2.18 2003/05/01 20:46:47 pruneau
586  * changed error parametrization
587  *
588  * Revision 2.17 2003/04/22 21:20:17 pruneau
589  * Added hit filter
590  * Tuning og finder pars
591  * Tuning of KalmanTrackNode
592  *
593  * Revision 2.16 2003/04/17 22:49:36 andrewar
594  * Fixed getPhase function to conform to StHelixModel convention.
595  *
596  * Revision 2.15 2003/03/31 17:18:56 pruneau
597  * various
598  *
599  * Revision 2.14 2003/03/13 21:21:27 pruneau
600  * getPhase() fixed. MUST inclde -helicity()*pi/2
601  *
602  * Revision 2.13 2003/03/13 18:59:13 pruneau
603  * various updates
604  *
605  * Revision 2.12 2003/03/12 17:57:31 pruneau
606  * Elss calc updated.
607  *
608  * Revision 2.11 2003/03/04 21:31:05 pruneau
609  * Added getX0() and getGasX0() conveninence methods.
610  *
611  * Revision 2.10 2003/03/04 18:41:27 pruneau
612  * Fixed StiHit to use global coordinates as well as locals.
613  * Fixed Logic Bug in StiKalmanTrackFinder
614  *
615  * Revision 2.9 2003/03/04 15:25:48 andrewar
616  * Added several functions for radlength calculation.
617  *
618  */
619 
620 #include <cassert>
621 #include <Stiostream.h>
622 #include <stdexcept>
623 #include <math.h>
624 #include <stdio.h>
625 #include <assert.h>
626 using namespace std;
627 #include "TCernLib.h"
628 #include "StiHit.h"
629 #include "StiDetector.h"
630 #include "StiPlacement.h"
631 #include "StiMaterial.h"
632 #include "StiShape.h"
633 #include "StiPlanarShape.h"
634 #include "StiCylindricalShape.h"
635 #include "StiKalmanTrackNode.h"
636 #include "StiElossCalculator.h"
637 #include "StDetectorDbMaker/StiTrackingParameters.h"
638 #include "StDetectorDbMaker/StiKalmanTrackFinderParameters.h"
639 #include "StDetectorDbMaker/StiHitErrorCalculator.h"
640 #include "StiTrackNodeHelper.h"
641 #include "StiFactory.h"
642 #include "StiUtilities/StiDebug.h"
643 
644 #include "TString.h"
645 #if ROOT_VERSION_CODE < 331013
646 #include "TCL.h"
647 #else
648 #include "TCernLib.h"
649 #endif
650 #include "THelixTrack.h"
651 #include "StThreeVector.hh"
652 #include "StThreeVectorF.hh"
653 #include "StarMagField.h"
654 #include "TMath.h"
655 #include "StMessMgr.h"
656 
657 #define PrP(A) { LOG_INFO << "\t" << (#A) << " = \t" << ( A ) }
658 #define PrPP(A,B) {LOG_INFO << "=== StiKalmanTrackNode::" << (#A); PrP((B)); LOG_INFO << endm;}
659 // Local Track Model
660 //
661 // x[0] = y coordinate
662 // x[1] = z position along beam axis
663 // x[2] = (Psi)
664 // x[3] = C (local) curvature of the track
665 // x[4] = tan(l)
666 
667 static const double kMaxEta = 1.25; // 72 degrees for laser tracks
668 static const double kMaxSinEta = sin(kMaxEta);
669 static const double kMaxCur = 0.2;
670 static const double kFarFromBeam = 10.;
671 static const Double_t kMaxZ = 250;
672 static const Double_t kMaxR = 250;
673 StiNodeStat StiKalmanTrackNode::mgP;
674 
675 
676 static const int idx33[3][3] = {{0,1,3},{1,2,4},{3,4,5}};
677 static const int idx55[5][5] =
678  {{0,1,3,6,10},{1,2,4,7,11},{3,4,5, 8,12},{6,7, 8, 9,13},{10,11,12,13,14}};
679 static const int idx55tpt[5][5] =
680  {{0,1,2,3, 4},{1,5,6,7, 8},{2,6,9,10,11},{3,7,10,12,13},{ 4, 8,11,13,14}};
681 
682 static const int idx66[6][6] =
683  {{ 0, 1, 3, 6,10,15},{ 1, 2, 4, 7,11,16},{ 3, 4, 5, 8,12,17}
684  ,{ 6, 7, 8, 9,13,18},{10,11,12,13,14,19},{15,16,17,18,19,20}};
685 
686 bool StiKalmanTrackNode::useCalculatedHitError = true;
687 TString StiKalmanTrackNode::comment("Legend: \tE - extapolation\tM Multiple scattering\tV at Vertex\tB at beam\tR at Radius\tU Updated\n");
688 TString StiKalmanTrackNode::commentdEdx("");
689 //debug vars
690 //#define STI_ERROR_TEST
691 //#define STI_DERIV_TEST
692 #ifdef STI_DERIV_TEST
693 int StiKalmanTrackNode::fDerivTestOn=0;
694 #endif
695 #ifndef STI_DERIV_TEST
696 int StiKalmanTrackNode::fDerivTestOn=-10;
697 #endif
698 
699 double StiKalmanTrackNode::fDerivTest[kNPars][kNPars];
700 int gCurrShape=0;
701 
702 /* bit mask for debug printout
703  0 => 1 - covariance and propagate matrices
704  1 => 2 - hit associated with the node
705  2 => 4 - test matrix manipulation
706  3 => 8 - test locate
707  */
708 int StiKalmanTrackNode::_debug = 0;
709 int StiKalmanTrackNode::_laser = 0;
710 
711 //______________________________________________________________________________
712 void StiKalmanTrackNode::reset()
713 {
714 static int myCount=0;
715  StiTrackNode::reset();
716  memset(_beg,0,_end-_beg+1);
717  _ext=0; _inf=0;
718  mId = ++myCount;
719  mHz = 999;
720  StiDebug::Break(mId);
721 }
722 //______________________________________________________________________________
723 void StiKalmanTrackNode::unset()
724 {
725  reduce();
726  if (_inf) BFactory::Free(_inf); _inf=0;
727 }
728 //______________________________________________________________________________
729 void StiKalmanTrackNode::resetError(double fak)
730 {
731 static const double DY=0.3,DZ=0.3,DEta=0.03,DPti=1.,DTan=0.05;
732 
733  if (!fak) {
734  mFE.reset();
735  mFE._cYY=DY*DY;
736  mFE._cZZ=DZ*DZ;
737  mFE._cEE=DEta*DEta;
738  mFE._cPP=DPti*DPti;
739  mFE._cTT=DTan*DTan;
740  } else {
741  for (int i=0;i<kNErrs;i++) mFE.G()[i] *=fak;
742  }
743  mPE() = mFE;
744 }
745 //_____________________________________________________________
750 //______________________________________________________________________________
751 void StiKalmanTrackNode::setState(const StiKalmanTrackNode * n)
752 {
753  _state = n->_state;
754  _alpha = n->_alpha;
755  mFP = n->mFP;
756  mFE = n->mFE;
757  mFP.hz()=0;
758  nullCount = n->nullCount;
759  contiguousHitCount = n->contiguousHitCount;
760  contiguousNullCount = n->contiguousNullCount;
761  setChi2(1e62);
762 }
763 
772 //______________________________________________________________________________
773 void StiKalmanTrackNode::get(double& alpha,
774  double& xRef,
775  double x[kNPars],
776  double e[kNErrs],
777  double& chi2)
778 {
779  alpha = _alpha;
780  xRef = getRefPosition();
781  memcpy(x,mFP.P,kNPars*sizeof(mFP.x()));
782  memcpy(e,mFE.G(),sizeof(mFE));
783  chi2 = getChi2();
784 }
785 
786 //______________________________________________________________________________
794 //______________________________________________________________________________
795 double StiKalmanTrackNode::getPt() const
796 {
797  assert(!std::isnan(mFP.ptin()));
798  return (fabs(mFP.ptin())<1e-3) ? 1e3: 1./fabs(mFP.ptin());
799 }
800 //______________________________________________________________________________
801 void StiKalmanTrackNode::propagateCurv(const StiKalmanTrackNode *parent)
802 {
803  mFP.ptin()=parent->mFP.ptin();
804  mFP.curv()=getHz()*mFP.ptin();
805 }
806 //______________________________________________________________________________
813 double StiKalmanTrackNode::getHz() const
814 {
815 
816 static const double EC = 2.99792458e-4,ZEROHZ = 2e-6;
817  if (fabs(mHz)<999) return mHz;
818  if (! _laser) {
819  double h[3];
820  StarMagField::Instance()->BField(&(getGlobalPoint().x()),h);
821  h[2] = EC*h[2];
822  if (fabs(h[2]) < ZEROHZ) h[2]=ZEROHZ;
823  mHz = h[2];
824  } else mHz = ZEROHZ;
825  assert(mHz);
826  return mHz;
827 }
828 //______________________________________________________________________________
858 //______________________________________________________________________________
859 void StiKalmanTrackNode::getMomentum(double p[3], double e[6]) const
860 {
861 // keep in mind that _eta == CA
862 // keep in mind that pt == SomeCoef/rho
863 enum {jX=0,jY,jZ,jE,jP,jT};
864 
865  double pt = getPt();
866  p[0] = pt*mFP._cosCA;
867  p[1] = pt*mFP._sinCA;
868  p[2] = pt*mFP.tanl();
869 
870 // if e==0, error calculation is not needed, then return
871  if (!e) return;
872 
873  double F[3][kNPars]; memset(F,0,sizeof(F));
874  double dPtdPi = -pt*pt; if (mFP.ptin()<0) dPtdPi=-dPtdPi;
875  F[jX][jE] = -pt *mFP._sinCA;
876  F[jX][jP] = dPtdPi*mFP._cosCA;
877  F[jX][jT] = 0;
878 
879  F[jY][jE] = pt *mFP._cosCA;
880  F[jY][jP] = dPtdPi*mFP._sinCA;
881  F[jY][jT] = 0;
882 
883  F[jZ][jE] = 0;
884  F[jZ][jP] = dPtdPi*mFP.tanl();
885  F[jZ][jT] = pt;
886 
887 
888  TCL::trasat(F[0],mFE.G(),e,3,kNPars);
889 }
890 //______________________________________________________________________________
909 //______________________________________________________________________________
910 void StiKalmanTrackNode::getGlobalRadial(double x[6],double e[15])
911 {
912  enum {jRad=0,jPhi,jZ,jTan,jPsi,jCur, kX=0,kY,kZ,kE,kC,kT};
913  double alpha,xRef,chi2;
914  double xx[kNPars],ee[kNErrs];
915 
916  get(alpha,xRef,xx,ee,chi2);
917 
918  x[jRad] = sqrt(pow(xx[kX],2)+pow(xx[kY],2));
919  x[jPhi] = atan2(xx[kY],xx[kX]) + alpha;
920  x[jZ ] = xx[kZ];
921  x[jTan] = xx[kT];
922  x[jPsi] = xx[kE] + alpha;
923  x[jCur] = xx[kC];
924  if (!e) return;
925 
926  double F[kNErrs][kNErrs]; memset(F,0,sizeof(F));
927  F[jPhi][kX] = -1e5;
928  F[jPhi][kY] = 1e5;
929  if (fabs(xx[kY])>1e-5) F[jPhi][kX] = -1./(xx[kY]);
930  if (fabs(xx[kX])>1e-5) F[jPhi][kY] = 1./(xx[kX]);
931  F[jZ][kZ] = 1.;
932  F[jTan][kT] = 1;
933  F[jPsi][kE] = 1;
934  F[jCur][kC] = 1;
935  memset(e,0,sizeof(*e)*15);
936  for (int k1=0;k1<kNPars;k1++) {
937  for (int k2=0;k2<kNPars;k2++) {
938  double cc = mFE.G()[idx66[k1][k2]];
939  for (int j1=jPhi;j1<= 5;j1++){
940  for (int j2=jPhi;j2<=j1;j2++){
941  e[idx55[j1-1][j2-1]]+= cc*F[j1][k1]*F[j2][k2];
942  }}}}
943 
944 }
945 //______________________________________________________________________________
978 //______________________________________________________________________________
979 void StiKalmanTrackNode::getGlobalTpt(float x[6],float e[15])
980 {
981  enum {jRad=0,jPhi,jZ,jTan,jPsi,jCur,jPt=jCur};
982 static const double DEG = 180./M_PI;
983 static double fak[6] = {1,0,1,1,DEG,0};
984 
985  double xx[6],ee[15];
986  getGlobalRadial(xx,ee);
987  double pt = getPt();
988  fak[jPhi] = DEG*xx[jRad];
989  fak[jPt] = (double(getCharge())/pt)/xx[jCur];
990 
991  for (int i=0;i<6;i++) {x[i] = (float)(fak[i]*xx[i]);}
992  if (!e) return;
993 
994  for (int j1=jPhi;j1<= 5;j1++){
995  for (int j2=jPhi;j2<=j1;j2++){
996  e[idx55tpt[j1-1][j2-1]] = (float)fak[j1]*fak[j2]*ee[idx55[j1-1][j2-1]];
997  }}
998 
999 }
1000 //______________________________________________________________________________
1001 double StiKalmanTrackNode::getPhase() const
1002 {
1006  return getPsi()-getHelicity()*M_PI/2;
1007 
1008 }
1009 //______________________________________________________________________________
1010 double StiKalmanTrackNode::getPsi() const
1011 {
1012  return nice(mFP.eta()+_alpha);
1013 }
1014 
1015 //______________________________________________________________________________
1029 
1030 //______________________________________________________________________________
1031 void StiKalmanTrackNode::getGlobalMomentum(double p[3], double e[6]) const
1032 {
1033  // first get p & e in the local ref frame
1034  enum {jXX=0,jXY,jYY};
1035 
1036  getMomentum(p,e);
1037  // now rotate the p & e in the global ref frame
1038  // for the time being, assume an azimuthal rotation
1039  // by alpha is sufficient.
1040  // transformation matrix - needs to be set
1041  double px=p[0];
1042  double py=p[1];
1043  double cosAlpha = cos(_alpha);
1044  double sinAlpha = sin(_alpha);
1045  p[0] = cosAlpha*px - sinAlpha*py;
1046  p[1] = sinAlpha*px + cosAlpha*py;
1047  if (e==0) return;
1048 
1049  // original error matrix
1050 
1051  double cXX = e[jXX];
1052  double cXY = e[jXY];
1053  double cYY = e[jYY];
1054  double cc = cosAlpha*cosAlpha;
1055  double ss = sinAlpha*sinAlpha;
1056  double cs = cosAlpha*sinAlpha;
1057  e[jXX] = cc*cXX - 2.*cs*cXY + ss*cYY;
1058  e[jYY] = ss*cXX + 2.*cs*cXY + cc*cYY;
1059  e[jXY] = cs*cXX + (cc-ss)*cXY - cs*cYY;
1060 }
1061 
1062 
1063 //______________________________________________________________________________
1078 //______________________________________________________________________________
1079 int StiKalmanTrackNode::propagate(StiKalmanTrackNode *pNode,
1080  const StiDetector * tDet,int dir)
1081 {
1082 static int nCall=0; nCall++;
1083 StiDebug::Break(nCall);
1084  int position = 0;
1085  setState(pNode);
1086  setDetector(tDet);
1087  if (mFP._cosCA <-1e-5) return -1;
1088  if (debug()) ResetComment(::Form("%40s ",tDet->getName().c_str()));
1089 
1090  StiPlacement * place = tDet->getPlacement();
1091  double nNormalRadius = place->getNormalRadius();
1092 
1093  StiShape * sh = tDet->getShape();
1094  int shapeCode = sh->getShapeCode();
1095  double endVal,dAlpha;
1096  switch (shapeCode) {
1097 
1098  case kPlanar: endVal = nNormalRadius;
1099  { //flat volume
1100  dAlpha = place->getNormalRefAngle();
1101  dAlpha = nice(dAlpha - _alpha);
1102  // bail out if the rotation fails...
1103  position = rotate(dAlpha);
1104  if (position) return -10;
1105  }
1106  break;
1107  case kDisk:
1108  case kCylindrical: endVal = nNormalRadius;
1109  {
1110  double out[2][3];
1111  double rxy = nNormalRadius;
1112  double rxy2P = mFP.rxy2();
1113  int outside = (rxy2P>rxy*rxy);
1114  int nSol = StiTrackNode::cylCross(mFP.P,&mFP._cosCA,mFP.curv(),rxy,dir,out);
1115  if (!nSol) return -11;
1116  double *ou = out[0];
1117  if (nSol==2) {
1118  int kaze = outside + 2*dir;
1119  switch (kaze) {
1120  case 0: break;
1121  case 1: ou = out[1]; break;
1122  case 2: ou = out[1]; break;
1123  case 3: return -99;
1124  } }
1125 
1126  dAlpha = atan2(ou[1],ou[0]);
1127  position = rotate(dAlpha);
1128  if (position) return -11;
1129  }
1130  break;
1131  default: assert(0);
1132  }
1133 
1134  position = propagate(endVal,shapeCode,dir);
1135 
1136  if (position) return position;
1137  assert(mFP.x() > 0.);
1138  if (mFP[0]*mFP._cosCA+mFP[1]*mFP._sinCA<0) return kEnded;
1139  propagateError();
1140  if (debug() & 8) { PrintpT("E");}
1141 
1142  // Multiple scattering
1143  if (StiKalmanTrackFinderParameters::instance()->mcsCalculated() && getHz())
1144  propagateMCS(pNode,tDet);
1145  if (debug() & 8) { PrintpT("M");}
1146  return position;
1147 }
1148 
1149 //______________________________________________________________________________
1160 bool StiKalmanTrackNode::propagate(const StiKalmanTrackNode *parentNode, StiHit * vertex,int dir)
1161 {
1162 static int nCall=0; nCall++;
1163 StiDebug::Break(nCall);
1164 
1165  setState(parentNode);
1166  TCircle tc(&mFP.x(),&mFP._cosCA,mFP.curv());
1167  double xy[2]; xy[0]=vertex->x(),xy[1]=vertex->y();
1168  double s = tc.Path(xy); tc.Move(s);
1169  double ang = atan2(tc.Dir()[1],tc.Dir()[0]);
1170  vertex->rotate(ang);
1171  rotate(ang);
1172  if (debug()) ResetComment(::Form("Vtx:%8.3f %8.3f %8.3f",vertex->x(),vertex->y(),vertex->z()));
1173  if (propagate(vertex->x(),kPlanar,dir)) return false; // track does not reach vertex "plane"
1174  propagateError();
1175  setHit(vertex);
1176  setDetector(0);
1177  return true;
1178 }
1179 
1180 //______________________________________________________________________________
1183 bool StiKalmanTrackNode::propagateToBeam(const StiKalmanTrackNode *parentNode,int dir)
1184 {
1185  setState(parentNode);
1186  if (debug()) {
1187  if (parentNode->getDetector())
1188  ResetComment(::Form("%40s ",parentNode->getDetector()->getName().c_str()));
1189  else ResetComment("Unknown Detector");
1190  }
1191  if (propagate(0., kPlanar,dir)) return false; // track does not reach vertex "plane"
1192 
1193  propagateError();
1194  if (mFE.zign()<0) return false;
1195  if (debug() & 8) { PrintpT("B"); PrintStep();}
1196  setHit(0);
1197  setDetector(0);
1198  return true;
1199 }
1200 
1201 //______________________________________________________________________________
1204 int StiKalmanTrackNode::propagateToRadius(StiKalmanTrackNode *pNode, double radius,int dir)
1205 {
1206  int position = 0;
1207  setState(pNode);
1208  if (debug()) ResetComment(::Form("%40s ",pNode->getDetector()->getName().c_str()));
1209  position = propagate(radius,kCylindrical,dir);
1210  if (position<0) return position;
1211  propagateError();
1212  if (debug() & 8) { PrintpT("R"); PrintStep();}
1213  _detector = 0;
1214  return position;
1215 }
1216 
1217 
1218 //______________________________________________________________________________
1227 //______________________________________________________________________________
1228 int StiKalmanTrackNode::propagate(double xk, int option,int dir)
1229 {
1230 static int nCall=0; nCall++;
1231 StiDebug::Break(nCall);
1232 
1233  _state = kTNProBeg;
1234 // numeDeriv(xk,1,option,dir);
1235  mgP.x1=mFP.x(); mgP.y1=mFP.y(); mgP.cosCA1 =mFP._cosCA; mgP.sinCA1 =mFP._sinCA;
1236  double rho = mFP.curv();
1237  mgP.x2 = xk;
1238 
1239  mgP.dx=mgP.x2-mgP.x1;
1240  double test = (dir)? mgP.dx:-mgP.dx;
1241 // if track is coming back stop tracking
1242 //VP if (test<0) return -3; //Unfortunatelly correct order not garanteed
1243 
1244  double dsin = mFP.curv()*mgP.dx;
1245  mgP.sinCA2=mgP.sinCA1 + dsin;
1246 // Orientation is bad. Fit is non reliable
1247  if (fabs(mgP.sinCA2)>kMaxSinEta) return -4;
1248  mgP.cosCA2 = ::sqrt((1.-mgP.sinCA2)*(1.+mgP.sinCA2));
1249 // Check what sign of cosCA2 must be
1250  test = (2*dir-1)*mgP.dx*mgP.cosCA1;
1251  if (test<0) mgP.cosCA2 = -mgP.cosCA2;
1252 
1253  int nIt = (mgP.cosCA2 <0)? 2:1;
1254  int ians = 0;
1255  StiNodePars save = mFP;
1256  for (int iIt=0; iIt<nIt; iIt++) {//try 2 cases, +ve and -ve cosCA
1257  ians = -1;
1258  mFP = save;
1259  mgP.cosCA2 = (!iIt)? fabs(mgP.cosCA2):-fabs(mgP.cosCA2);
1260  mgP.sumSin = mgP.sinCA1+mgP.sinCA2;
1261  mgP.sumCos = mgP.cosCA1+mgP.cosCA2;
1262  if (fabs(mgP.sumCos)<1e-6) continue;
1263  mgP.dy = mgP.dx*(mgP.sumSin/mgP.sumCos);
1264  mgP.y2 = mgP.y1+mgP.dy;
1265 
1266 
1267  mgP.dl0 = mgP.cosCA1*mgP.dx+mgP.sinCA1*mgP.dy;
1268  double sind = mgP.dl0*rho;
1269 
1270  if (fabs(dsin) < 0.02 ) { //tiny angle
1271  mgP.dl = mgP.dl0*(1.+sind*sind/6);
1272 
1273  } else {
1274  double cosd = mgP.cosCA2*mgP.cosCA1+mgP.sinCA2*mgP.sinCA1;
1275  mgP.dl = atan2(sind,cosd)/rho;
1276  }
1277  if (mgP.y2*mgP.y2+mgP.x2*mgP.x2>kMaxR*kMaxR) return -5;
1278  mFP.z() += mgP.dl*mFP.tanl();
1279  if (fabs(mFP.z()) > kMaxZ) return -6;
1280  mFP.y() = mgP.y2;
1281  mFP.eta() = nice(mFP.eta()+rho*mgP.dl); /*VP*/
1282  mFP.x() = mgP.x2;
1283  mFP._sinCA = mgP.sinCA2;
1284  mFP._cosCA = mgP.cosCA2;
1285  ians = locate();
1286  if (!ians) break;
1287  }
1288  if (ians) return kFailed;
1289  if (fabs(mFP.eta())>kMaxEta) return kFailed;
1290  if (mFP.x()> kFarFromBeam) {
1291  if (mFP.x()*mgP.cosCA2+mFP.y()*mgP.sinCA2<=0) return kEnded;
1292  }
1293  mFP.hz() = getHz();
1294  mFP.curv() = mFP.hz()*mFP.ptin();
1295 
1296  mPP() = mFP;
1297  _state = kTNProEnd;
1298  return ians;
1299 }
1300 
1301 //______________________________________________________________________________
1302 int StiKalmanTrackNode::nudge(StiHit *hitp)
1303 {
1304  StiHit *hit = hitp;
1305  if (!hit) hit = getHit();
1306  double deltaX = 0;
1307  if (hit) { deltaX = hit->x()-mFP.x();}
1308  else { if (_detector) deltaX = _detector->getPlacement()->getNormalRadius()-mFP.x();}
1309  if(fabs(deltaX)>5) return -1;
1310  if (fabs(deltaX) <1.e-3) return 0;
1311  double deltaS = mFP.curv()*(deltaX);
1312  double sCA2 = mFP._sinCA + deltaS;
1313  if (fabs(sCA2)>0.99) return -2;
1314  double cCA2,deltaY,deltaL,sind;
1315  if (fabs(deltaS) < 1e-3 && fabs(mFP.eta())<1) { //Small angle approx
1316  cCA2= mFP._cosCA - mFP._sinCA/mFP._cosCA*deltaS;
1317  if (cCA2> 1) cCA2= 1;
1318  if (cCA2<-1) cCA2=-1;
1319  deltaY = deltaX*(mFP._sinCA+sCA2)/(mFP._cosCA+cCA2);
1320  deltaL = deltaX*mFP._cosCA + deltaY*mFP._sinCA;
1321  sind = deltaL*mFP.curv();
1322  deltaL = deltaL*(1.+sind*sind/6);
1323  } else {
1324  cCA2= sqrt((1.-sCA2)*(1.+sCA2));
1325  if (mFP._cosCA <0) cCA2 = -cCA2;
1326  deltaY = deltaX*(mFP._sinCA+sCA2)/(mFP._cosCA+cCA2);
1327  deltaL = deltaX*mFP._cosCA + deltaY*mFP._sinCA;
1328  sind = deltaL*mFP.curv();
1329  deltaL = asin(sind)/mFP.curv();
1330  }
1331  double deltaZ = mFP.tanl()*(deltaL);
1332  mFP._sinCA = mgP.sinCA2 = sCA2;
1333  mFP._cosCA = mgP.cosCA2 = cCA2;
1334  mgP.sumSin = mgP.sinCA1+mgP.sinCA2;
1335  mgP.sumCos = mgP.cosCA1+mgP.cosCA2;
1336  mFP.x() += deltaX;
1337  mFP.y() += deltaY;
1338  mFP.z() += deltaZ;
1339  mFP.eta() += deltaL*mFP.curv();
1340  mgP.dx += deltaX;
1341  mgP.dy += deltaY;
1342  mgP.dl0 += deltaL;
1343  mgP.dl += deltaL;
1344 
1345 
1346 // assert(fabs(mFP._sinCA) < 1.);
1347  if (fabs(mFP._sinCA)>=1) {
1348  LOG_DEBUG << Form("StiKalmanTrackNode::nudge WRONG WRONG WRONG sinCA=%g",mFP._sinCA)
1349  << endm;
1350  mFP.print();
1351  return -13;
1352  }
1353  assert(fabs(mFP._cosCA) <= 1.);
1354  mPP() = mFP;
1355  return 0;
1356 }
1357 //______________________________________________________________________________
1360 void StiKalmanTrackNode::propagateMtx()
1361 {
1362 // fYE == dY/dEta
1363  double fYE= mgP.dx*(1.+mgP.cosCA1*mgP.cosCA2+mgP.sinCA1*mgP.sinCA2)/(mgP.sumCos*mgP.cosCA2);
1364 // fEC == dEta/dRho
1365  double fEC = mgP.dx/mgP.cosCA2;
1366 // fYC == dY/dRho
1367  double fYC=(mgP.dy*mgP.sinCA2+mgP.dx*mgP.cosCA2)/mgP.sumCos*fEC;
1368 // fZC == dZ/dRho
1369  double dang = mgP.dl*mFP.curv();
1370  double C2LDX = mgP.dl*mgP.dl*(
1371  0.5*mgP.sinCA2*pow((1+pow(dang/2,2)*sinX(dang/2)),2) +
1372  mgP.cosCA2*dang*sinX(dang));
1373 
1374  double fZC = mFP.tanl()*C2LDX/mgP.cosCA2;
1375 // fZE == dZ/dEta
1376  double dLdEta = mgP.dy/mgP.cosCA2;
1377  double fZE = mFP.tanl()*dLdEta;
1378 
1379 // fZT == dZ/dTanL;
1380  double fZT= mgP.dl;
1381  double hz = getHz(); fEC *=hz; fYC*=hz; fZC*=hz;
1382 
1383  double ca =1, sa=0;
1384  if (mMtx().A[0][0]) { ca = mMtx().A[0][0]+1.;sa = mMtx().A[0][1];}
1385  mMtx().reset();
1386 // X related derivatives
1387  mMtx().A[0][0] = -1;
1388  mMtx().A[1][0] = -mgP.sinCA2/mgP.cosCA2;
1389  mMtx().A[2][0] = -mFP.tanl() /mgP.cosCA2;
1390  mMtx().A[3][0] = -mFP.curv() /mgP.cosCA2;
1391 
1392  mMtx().A[1][3]=fYE; mMtx().A[1][4]=fYC; mMtx().A[2][3]=fZE;
1393  mMtx().A[2][4]=fZC; mMtx().A[2][5]=fZT; mMtx().A[3][4]=fEC;
1394  if (sa) {
1395  double fYX = mMtx().A[1][0];
1396  mMtx().A[1][0] = fYX*ca-sa;
1397  mMtx().A[1][1] = fYX*sa+ca-1;
1398  }
1399 }
1400 
1401 
1402 
1403 //______________________________________________________________________________
1406 void StiKalmanTrackNode::propagateError()
1407 {
1408  static int nCall=0; nCall++;
1409  StiDebug::Break(nCall);
1410  propagateMtx();
1411  errPropag6(mFE.G(),mMtx().A,kNPars);
1412  int force = (fabs(mgP.dl) > StiNodeErrs::kBigLen);
1413  mFE.recov(force);
1414  mFE._cXX = mFE._cYX= mFE._cZX = mFE._cEX = mFE._cPX = mFE._cTX = 0;
1415 // now set hiterrors
1416  if (_hit) setHitErrors();
1417 
1418 // set state node is ready
1419 
1420  mPE() = mFE;
1421  _state = kTNReady;
1422 }
1423 
1424 //______________________________________________________________________________
1440 //delta(mgP.dx,dy,dz) = here - there
1441 double StiKalmanTrackNode::pathLToNode(const StiKalmanTrackNode * const oNode)
1442 {
1443  const StThreeVector<double> delta =
1444  getGlobalPoint() - oNode->getGlobalPoint();
1445  double R = getCurvature();
1446  // s = 2c * asin( t/(2c)); t=::sqrt(mgP.dx^2+dy^2+dz^2)
1447  return length(delta, R);
1448 }
1449 
1450 //______________________________________________________________________________
1451 double StiKalmanTrackNode::length(const StThreeVector<double>& delta, double curv)
1452 {
1453 
1454  double m = delta.perp();
1455  double as = 0.5*m*fabs(curv);
1456  if (as >= 1) {
1457  if (as > 1.1)
1458  cout << "StiKalmanTrackNode::length m = " << m << " curv = " << curv << " as = " << as << " illegal. reset to 1" << endl;
1459  as = 0.99;
1460  }
1461  double lxy=0;
1462  if (fabs(as)<0.01) { lxy = m*(1.+as*as/24);}
1463  else { lxy = 2.*asin(as)/fabs(curv);}
1464  return sqrt(lxy*lxy+delta.z()*delta.z());
1465 }
1466 
1467 //______________________________________________________________________________
1481 double StiKalmanTrackNode::evaluateChi2(const StiHit * hit)
1482 {
1483  double r00, r01,r11;
1484  //If required, recalculate the errors of the detector hits.
1485  //Do not attempt this calculation for the main vertex.
1486  double dsin =mFP.curv()*(hit->x()-mFP.x());
1487  if (fabs(mFP._sinCA+dsin)>0.99 ) return 1e41;
1488  if (fabs(mFP.eta()) >kMaxEta) return 1e41;
1489  if (fabs(mFP.curv()) >kMaxCur) return 1e41;
1490  if (mHrr.hYY>1000*mFE._cYY
1491  && mHrr.hZZ>1000*mFE._cZZ) return 1e41;
1492  setHitErrors(hit);
1493  r00=mHrr.hYY+mFE._cYY;
1494  r01=mHrr.hZY+mFE._cZY;
1495  r11=mHrr.hZZ+mFE._cZZ;
1496 
1497  _det=r00*r11 - r01*r01;
1498  if (_det<r00*r11*1.e-5) {
1499  LOG_DEBUG << Form("StiKalmanTrackNode::evalChi2 *** zero determinant %g",_det)<< endm;
1500  return 1e60;
1501  }
1502  double tmp=r00; r00=r11; r11=tmp; r01=-r01;
1503  double deltaX = hit->x()-mFP.x();
1504  double deltaL = deltaX/mFP._cosCA;
1505  double deltaY = mFP._sinCA *deltaL;
1506  double deltaZ = mFP.tanl() *deltaL;
1507  double dyt=(mFP.y()-hit->y()) + deltaY;
1508  double dzt=(mFP.z()-hit->z()) + deltaZ;
1509  double cc= (dyt*r00*dyt + 2*r01*dyt*dzt + dzt*r11*dzt)/_det;
1510  if (debug() & 8) {comment += Form(" chi2 = %6.2f",cc);}
1511  return cc;
1512 }
1513 //______________________________________________________________________________
1514 int StiKalmanTrackNode::isEnded() const
1515 {
1516 
1517  if(fabs(mFP.eta() )<=kMaxEta) return 0;
1518  return 1;
1519 }
1520 //______________________________________________________________________________
1521 int StiKalmanTrackNode::isDca() const
1522 {
1523  return (fabs(mFP.x())<=0);
1524 }
1525 
1526 //______________________________________________________________________________
1535 void StiKalmanTrackNode::propagateMCS(StiKalmanTrackNode * previousNode, const StiDetector * tDet)
1536 {
1537 static const int keepElossBug = StiDebug::iFlag("keepElossBug");
1538 
1539 static int nCall=0; nCall++;
1540  propagateCurv(previousNode);
1541  double pt = getPt();
1542  if (pt>=1e2) {
1543  mPP() = mFP; mPE() = mFE;
1544  return;
1545  }
1546  double relRadThickness;
1547  // Half path length in previous node
1548  double pL1,pL2,pL3,d1,d2,d3,dxEloss=0;
1549  pL1=previousNode->pathlength()/2.;
1550  // Half path length in this node
1551  pL3=pathlength()/2.;
1552  // Gap path length
1553  pL2= pathLToNode(previousNode);
1554  if (pL1<0) pL1=0;
1555  if (pL2<0) pL2=0;
1556  if (pL3<0) pL3=0;
1557  double x0p =1e11,x0Gas=1e11,x0=1e11;
1558 
1559  if (pt > 0.350 && TMath::Abs(getHz()) < 1e-3) pt = 0.350;
1560  double p2=(1.+mFP.tanl()*mFP.tanl())*pt*pt;
1561  double m=StiKalmanTrackFinderParameters::instance()->massHypothesis();
1562  double m2=m*m;
1563  double e2=p2+m2;
1564  double beta2=p2/e2;
1565 
1566 if (keepElossBug) { //Old Eloss bug prezerved
1567  d1 = previousNode->getDensity();
1568  x0p = previousNode->getX0();
1569  d3 = tDet->getMaterial()->getDensity();
1570  x0 = tDet->getMaterial()->getX0();
1571 
1572 
1573  if (pL2> (pL1+pL3))
1574  {
1575  pL2=pL2-pL1-pL3;
1576  if (mgP.dx>0)
1577  {
1578  x0Gas = tDet->getGas()->getX0();
1579  d2 = tDet->getGas()->getDensity();
1580  }
1581  else
1582  {
1583  x0Gas = previousNode->getGasX0();
1584  d2 = previousNode->getGasDensity();
1585  }
1586  relRadThickness = 0.;
1587  dxEloss = 0;
1588  if (x0p>0.)
1589  {
1590  relRadThickness += pL1/x0p;
1591  dxEloss += d1*pL1;
1592  }
1593  if (x0Gas>0.)
1594  {
1595  relRadThickness += pL2/x0Gas;
1596  dxEloss += d2*pL2;
1597  }
1598  if (x0>0.)
1599  {
1600  relRadThickness += pL3/x0;
1601  dxEloss += d3*pL3;
1602  }
1603  }
1604  else
1605  {
1606  relRadThickness = 0.;
1607  dxEloss = 0;
1608  if (x0p>0.)
1609  {
1610  relRadThickness += pL1/x0p;
1611  dxEloss += d1*pL1;
1612  }
1613  if (x0>0.)
1614  {
1615  relRadThickness += pL3/x0;
1616  dxEloss += d3*pL3;
1617  }
1618  }
1619  //cout << " m2:"<<m2<<" p2:"<<p2<<" beta2:"<<beta2;
1620  double theta2=mcs2(relRadThickness,beta2,p2);
1621  //cout << " theta2:"<<theta2;
1622  double pti = mFP.ptin(), tanl = mFP.tanl();
1623 
1624  double cos2Li = (1.+ tanl*tanl); // 1/cos(lamda)**2
1625 
1626  mFE._cEE += cos2Li *theta2;
1627  mFE._cPP += tanl*tanl*pti*pti *theta2;
1628  mFE._cTP += pti*tanl*cos2Li *theta2;
1629  mFE._cTT += cos2Li*cos2Li *theta2;
1630 
1631  double dE=0;
1632  double sign = (mgP.dx>0)? 1:-1;
1633 
1634 // const static double I2Ar = (15.8*18) * (15.8*18) * 1e-18; // GeV**2
1635  StiElossCalculator * calculator = tDet->getMaterial()->getElossCalculator();
1636 assert(calculator);
1637  double eloss = calculator->calculate(1.,m, beta2);
1638  dE = sign*dxEloss*eloss;
1639 
1640 // save detLoss and gasLoss for investigation only
1641 // setELoss(2*sign*d3*eloss,sign*d2*eloss);
1642 
1643  if (TMath::Abs(dE)>0)
1644  {
1645  if (debug()) {
1646  commentdEdx = Form("%6.3g cm(%5.2f) %6.3g keV %6.3f GeV",mgP.dx,100*relRadThickness,1e6*dE,TMath::Sqrt(e2)-m);
1647  }
1648  double correction =1. + ::sqrt(e2)*dE/p2;
1649  if (correction>1.1) correction = 1.1;
1650  else if (correction<0.9) correction = 0.9;
1651  mFP.curv() = mFP.curv()*correction;
1652  mFP.ptin() = mFP.ptin()*correction;
1653  }
1654  mPP() = mFP; mPE() = mFE;
1655 
1656 } else { // ELoss bug fixed =====================================================================
1657 
1658  const StiDetector *preDet = previousNode->getDetector();
1659  const StiMaterial *preMat = preDet->getMaterial();
1660  const StiElossCalculator *preLos = preMat->getElossCalculator();
1661  d1 =(preLos) ? preLos->calculate(1.,m, beta2):0;
1662  x0p = preMat->getX0();
1663 
1664  const StiMaterial *curMat = tDet->getMaterial();
1665  const StiElossCalculator *curLos = curMat->getElossCalculator();
1666  d3 =(curLos) ? curLos->calculate(1.,m, beta2):0;
1667  x0 = curMat->getX0();
1668  double sign = (mgP.dx>0)? 1:-1;
1669 
1670 // Gas of detector is placed UNDER of it
1671  const StiMaterial *gasMat = (sign>0)? tDet->getGas() : preDet->getGas();
1672  x0Gas = gasMat->getX0();
1673  const StiElossCalculator *gasLos = gasMat->getElossCalculator();
1674  d2 =(gasLos) ? gasLos->calculate(1.,m, beta2):0;
1675 
1676  pL2=pL2-pL1-pL3; if (pL2<0) pL2=0;
1677  relRadThickness = pL1/x0p+pL2/x0Gas+pL3/x0;
1678  dxEloss = d1*pL1+ d2*pL2 + d3*pL3;
1679 
1680  //cout << " m2:"<<m2<<" p2:"<<p2<<" beta2:"<<beta2;
1681  double theta2=mcs2(relRadThickness,beta2,p2);
1682  //cout << " theta2:"<<theta2;
1683  double pti = mFP.ptin(), tanl = mFP.tanl();
1684 
1685  double cos2Li = (1.+ tanl*tanl); // 1/cos(lamda)**2
1686 
1687  mFE._cEE += cos2Li *theta2;
1688  mFE._cPP += tanl*tanl*pti*pti *theta2;
1689  mFE._cTP += pti*tanl*cos2Li *theta2;
1690  mFE._cTT += cos2Li*cos2Li *theta2;
1691 assert(mFE._cEE>0);
1692 assert(mFE._cPP>0);
1693 assert(mFE._cTT>0);
1694 
1695 assert(mFE.zign()>0);
1696 
1697  double dE = sign*dxEloss;
1698 // save detLoss and gasLoss for investigation only
1699 
1700  mELoss[0].mELoss = 2*d3*pL3;
1701  mELoss[0].mLen = 2*pL3;
1702  mELoss[0].mDens = curMat->getDensity();
1703  mELoss[0].mX0 = x0;
1704  mELoss[1].mELoss = 2*d2*pL2;
1705  mELoss[1].mLen = 2*pL2;
1706  mELoss[1].mDens = gasMat->getDensity();
1707  mELoss[1].mX0 = x0Gas;
1708 
1709  if (fabs(dE)>0)
1710  {
1711  double correction =1. + ::sqrt(e2)*dE/p2;
1712  if (correction>1.1) correction = 1.1;
1713  else if (correction<0.9) correction = 0.9;
1714  mFP.curv() = mFP.curv()*correction;
1715  mFP.ptin() = mFP.ptin()*correction;
1716  }
1717  mPP() = mFP; mPE() = mFE;
1718 }
1719 
1720 
1721 }
1722 
1723 //______________________________________________________________________________
1743 //______________________________________________________________________________
1744 int StiKalmanTrackNode::updateNode()
1745 {
1746 static int nCall=0; nCall++;
1747  _state = kTNFitBeg;
1748 assert(mFE.zign()>0);
1749  assert(mFE._cXX<1e-8);
1750  double r00,r01,r11;
1751  r00 = mHrr.hYY + mFE._cYY;
1752  r01 = mHrr.hZY + mFE._cZY;
1753  r11 = mHrr.hZZ + mFE._cZZ;
1754  _det=r00*r11 - r01*r01;
1755  if (!finite(_det) || _det<(r00*r11)*1.e-5) {
1756  LOG_DEBUG << Form("StiKalmanTrackNode::updateNode *** zero determinant %g",_det)
1757  << endm;
1758  return -11;
1759  }
1760  // inverse matrix
1761  double tmp=r00; r00=r11/_det; r11=tmp/_det; r01=-r01/_det;
1762  // update error matrix
1763  double k00=mFE._cYY*r00+mFE._cZY*r01, k01=mFE._cYY*r01+mFE._cZY*r11;
1764  double k10=mFE._cZY*r00+mFE._cZZ*r01, k11=mFE._cZY*r01+mFE._cZZ*r11;
1765  double k20=mFE._cEY*r00+mFE._cEZ*r01, k21=mFE._cEY*r01+mFE._cEZ*r11;
1766  double k30=mFE._cPY*r00+mFE._cPZ*r01, k31=mFE._cPY*r01+mFE._cPZ*r11;
1767  double k40=mFE._cTY*r00+mFE._cTZ*r01, k41=mFE._cTY*r01+mFE._cTZ*r11;
1768  double dyt = getHit()->y() - mFP.y();
1769  double dzt = getHit()->z() - mFP.z();
1770  double dp3 = k30*dyt + k31*dzt;
1771  double dp2 = k20*dyt + k21*dzt;
1772  double dp4 = k40*dyt + k41*dzt;
1773  double eta = nice(mFP.eta() + dp2);
1774  if (fabs(eta)>kMaxEta) return -14;
1775  double pti = mFP.ptin() + dp3;
1776  if (fabs(pti) < 1e-3) pti=1e-3;
1777  double cur = pti*getHz();
1778  if (fabs(cur)>kMaxCur) return -16;
1779  assert(finite(cur));
1780  double tanl = mFP.tanl() + dp4;
1781  // update Kalman state
1782  double p0 = mFP.y() + k00*dyt + k01*dzt;
1783 //VP mFP.y() += k00*dy + k01*dz;
1784  if (fabs(p0)>kMaxR)
1785  {
1786  LOG_DEBUG << "updateNode()[1] -W- _y:"<<mFP.y()<<" _z:"<<mFP.z()<<endm;
1787  return -12;
1788  }
1789  double p1 = mFP.z() + k10*dyt + k11*dzt;
1790  if (fabs(p1)>kMaxZ)
1791  {
1792  LOG_DEBUG << "updateNode()[2] -W- _y:"<<mFP.y()<<" _z:"<<mFP.z()<<endm;
1793  return -13;
1794  }
1795  //mFP.tanl() += k40*dyt + k41*dzt;
1796  double sinCA = sin(eta);
1797  // The following test introduces a track propagation error but happens
1798  // only when the track should be aborted so we don't care...
1799  mFP.y() = p0;
1800  mFP.z() = p1;
1801  mFP.hz() = getHz();
1802  if (mFP.isZeroH()) { mFP.ready(); pti = mFP.ptin(); cur = mFP.curv();}
1803  mFP.eta() = eta;
1804  mFP.ptin() = pti;
1805  mFP.curv() = cur;
1806  mFP.tanl() = tanl;
1807  mFP._sinCA = sinCA;
1808  mFP._cosCA = ::sqrt((1.-mFP._sinCA)*(1.+mFP._sinCA));
1809  assert(!_detector || mFP.x() > 0.);
1810 
1811 // update error matrix
1812  double c00=mFE._cYY;
1813  double c10=mFE._cZY, c11=mFE._cZZ;
1814  double c20=mFE._cEY, c21=mFE._cEZ;//, c22=mFE._cEE;
1815  double c30=mFE._cPY, c31=mFE._cPZ;//, c32=mFE._cPE, c33=mFE._cPP;
1816  double c40=mFE._cTY, c41=mFE._cTZ;//, c42=mFE._cTE, c43=mFE._cTP, c44=mFE._cTT;
1817  mFE._cYY-=k00*c00+k01*c10;
1818  mFE._cZY-=k10*c00+k11*c10;mFE._cZZ-=k10*c10+k11*c11;
1819  mFE._cEY-=k20*c00+k21*c10;mFE._cEZ-=k20*c10+k21*c11;mFE._cEE-=k20*c20+k21*c21;
1820  mFE._cPY-=k30*c00+k31*c10;mFE._cPZ-=k30*c10+k31*c11;mFE._cPE-=k30*c20+k31*c21;mFE._cPP-=k30*c30+k31*c31;
1821  mFE._cTY-=k40*c00+k41*c10;mFE._cTZ-=k40*c10+k41*c11;mFE._cTE-=k40*c20+k41*c21;mFE._cTP-=k40*c30+k41*c31;mFE._cTT-=k40*c40+k41*c41;
1822  mFE.recov(1);
1823 
1824  _state = kTNFitEnd;
1825  return 0;
1826 }
1827 
1828 //______________________________________________________________________________
1840 int StiKalmanTrackNode::rotate (double alpha)
1841 {
1842  mMtx().A[0][0]=0;
1843  if (fabs(alpha)<1.e-6) return 0;
1844  _state = kTNRotBeg;
1845  _alpha += alpha;
1846  _alpha = nice(_alpha);
1847  //cout << " new _alpha:"<< 180.*_alpha/3.1415927<<endl;
1848 
1849  double xt1=mFP.x();
1850  double yt1=mFP.y();
1851  mgP.sinCA1 = mFP._sinCA;
1852  mgP.cosCA1 = mFP._cosCA;
1853  double ca = cos(alpha);
1854  double sa = sin(alpha);
1855  mFP.x() = xt1*ca + yt1*sa;
1856  mFP.y()= -xt1*sa + yt1*ca;
1857  mFP._cosCA = mgP.cosCA1*ca+mgP.sinCA1*sa;
1858  mFP._sinCA = -mgP.cosCA1*sa+mgP.sinCA1*ca;
1859  double nor = 0.5*(mFP._sinCA*mFP._sinCA+mFP._cosCA*mFP._cosCA +1);
1860  mFP._cosCA /= nor;
1861  mFP._sinCA /= nor;
1862 
1863  mFP.eta()= nice(mFP.eta()-alpha); /*VP*/
1864  mFP._sinCA = sin(mFP.eta());
1865  mFP._cosCA = cos(mFP.eta());
1866 //cout << " mFP._sinCA:"<<mFP._sinCA<<endl;
1867  assert(fabs(mFP._sinCA)<=1.);
1868  assert(fabs(mFP._cosCA)<=1.);
1869  memset(mMtx().A,0,sizeof(mMtx()));
1870  mMtx().A[0][0]= ca-1;
1871  mMtx().A[0][1]= sa;
1872  mMtx().A[1][0]=-sa;
1873  mMtx().A[1][1]= ca-1;
1874 
1875  _state = kTNRotEnd;
1876  mPP() = mFP;
1877  return 0;
1878 }
1879 
1880 
1881 //_____________________________________________________________________________
1884 ostream& operator<<(ostream& os, const StiKalmanTrackNode& n)
1885 {
1886  const StiDetector *det = n.getDetector();
1887  if (det) os <<"Det:"<<n.getDetector()->getName();
1888  else os << "Det:UNknown";
1889  os << " a:" << 180*n._alpha/M_PI<<" degs"
1890  << "\tx:" << n.mFP.x()
1891  << " p0:" << n.mFP.y()
1892  << " p1:" << n.mFP.z()
1893  << " p2:" << n.mFP.eta()
1894  << " p3:" << n.mFP.ptin()
1895  << " p4:" << n.mFP.tanl()
1896  << " c00:" <<n.mFE._cYY<< " c11:"<<n.mFE._cZZ
1897  << " pT:" << n.getPt() << endl;
1898  if (n.debug() & 2) {
1899  StiHit * hit = n.getHit();
1900  if (hit) os << "\thas hit with chi2 = " << n.getChi2()
1901  << " n:"<<n.hitCount
1902  << " null:"<<n.nullCount
1903  << endl<<"\t hit:"<<*hit;
1904  }
1905  else os << endl;
1906  return os;
1907 }
1908 
1909 //______________________________________________________________________________
1910 double StiKalmanTrackNode::getWindowY()
1911 {
1912  const StiDetector * detector = getDetector();
1913  const StiTrackingParameters * parameters = detector->getTrackingParameters();
1914  double searchWindowScale = parameters->getSearchWindowScale();
1915  double minSearchWindow = parameters->getMinSearchWindow();
1916  double maxSearchWindow = parameters->getMaxSearchWindow();
1917 
1918  const StiHitErrorCalculator * calc = detector->getHitErrorCalculator();
1919  double myEyy,myEzz;
1920  calc->calculateError(&mFP,myEyy,myEzz);
1921  double window = searchWindowScale*::sqrt(mFE._cYY+myEyy);
1922  if (window<minSearchWindow) window = minSearchWindow;
1923  else if (window>maxSearchWindow) window = maxSearchWindow;
1924  return window;
1925 }
1926 
1927 //_____________________________________________________________________________
1928 double StiKalmanTrackNode::getWindowZ()
1929 {
1930  const StiDetector * detector = getDetector();
1931  const StiTrackingParameters * parameters = detector->getTrackingParameters();
1932  double searchWindowScale = parameters->getSearchWindowScale();
1933  double minSearchWindow = parameters->getMinSearchWindow();
1934  double maxSearchWindow = parameters->getMaxSearchWindow();
1935 
1936  const StiHitErrorCalculator * calc = detector->getHitErrorCalculator();
1937  double myEyy,myEzz;
1938  calc->calculateError(&mFP,myEyy,myEzz);
1939 
1940  double window = searchWindowScale*::sqrt(mFE._cZZ+myEzz);
1941  if (window<minSearchWindow) window = minSearchWindow;
1942  else if (window>maxSearchWindow) window = maxSearchWindow;
1943  return window;
1944 }
1945 
1946 //______________________________________________________________________________
1947 StThreeVector<double> StiKalmanTrackNode::getHelixCenter() const
1948 {
1949  assert(mFP.curv());
1950  double xt0 = mFP.x()-mFP._sinCA/mFP.curv(); /*VP*/
1951  double yt0 = mFP.y()+mFP._cosCA/(mFP.curv());
1952  double zt0 = mFP.z()+mFP.tanl()*asin(mFP._sinCA)/mFP.curv();
1953  double cosAlpha = cos(_alpha);
1954  double sinAlpha = sin(_alpha);
1955  return (StThreeVector<double>(cosAlpha*xt0-sinAlpha*yt0,sinAlpha*xt0+cosAlpha*yt0,zt0));
1956 }
1957 #if 1
1958 //______________________________________________________________________________
1959 int StiKalmanTrackNode::locate()
1960 {
1961  double yOff, zOff,ang;
1962  //fast way out for projections going out of fiducial volume
1963  const StiDetector *tDet = getDetector();
1964  if (!tDet) return 0;
1965  const StiPlacement *place = tDet->getPlacement();
1966  const StiShape *sh = tDet->getShape();
1967 
1968  if (fabs(mFP.z())>kMaxZ || mFP.rxy()> kMaxR) return -1;
1969 
1970 
1971  //YF edge is tolerance when we consider that detector is hit. // edge = 0; //VP the meaning of edge is not clear
1972  Int_t shapeCode = sh->getShapeCode();
1973  switch (shapeCode) {
1974  case kDisk:
1975  case kCylindrical: // cylinder
1976  break;
1977  case kSector: // cylinder sector
1978  ang = atan2(mFP.y(),mFP.x());
1979  yOff = nice(ang +_alpha - place->getLayerAngle());
1980  if (fabs(yOff)>sh->getOpeningAngle()/2) return -1;
1981  break;
1982  case kPlanar:
1983  yOff = mFP.y() - place->getNormalYoffset();
1984  if (fabs(yOff)> sh->getHalfWidth()) return -1;
1985  break;
1986  default: assert(0 && "Wrong Shape code");
1987  }
1988  zOff = mFP.z() - place->getZcenter();
1989  if (fabs(zOff)>sh->getHalfDepth()) return -1;
1990  return 0;
1991  }
1992 #endif //1
1993 
1994 //______________________________________________________________________________
1995 void StiKalmanTrackNode::initialize(const double dirg[3])
1996 {
1997 // dir - direction normalised dir[0]**2 +dir[1]**2 = 1
1998 
1999  double cosp = _detector->getCos();
2000  double sinp = _detector->getSin();
2001  double dirl[3] = {0,0,dirg[2]};
2002  dirl[0] = cosp*dirg[0] + sinp*dirg[1];
2003  dirl[1] =-sinp*dirg[0] + cosp*dirg[1];
2004  if (fabs(dirl[0])>=1.) {dirl[0] = 1; dirl[1]=0;}
2005  if (fabs(dirl[1])>=1.) {dirl[1] = 1; dirl[0]=0;}
2006  mFP._cosCA = dirl[0];
2007  mFP._sinCA = dirl[1];
2008  mFP.P[StiNodePars::kPhi] = atan2(dirl[1],dirl[0]);
2009  mFP.P[StiNodePars::kTan] = dirl[2];
2010 
2011 // assert(_hit->x_g()*dirg[0]+_hit->y_g()*dirg[1]>0);
2012 // assert(mFP.x()*mFP._cosCA + mFP.y()*mFP._sinCA>0);
2013 //
2014 //
2015  mPP() = mFP;
2016 }
2017 //______________________________________________________________________________
2018 void StiKalmanTrackNode::initialize(StiHit *h)
2019 {
2020  reset();
2021  setHit(h);
2022  _detector = h->detector();
2023  _alpha = _detector->getPlacement()->getNormalRefAngle();
2024  mFP._sinCA = 0.6;
2025  mFP._cosCA = 0.8;
2026  mFP.x() = h->x();
2027  mFP.y() = h->y();
2028  mFP.z() = h->z();
2029  mFP.hz() = getHz();
2030  resetError();
2031  mPP() = mFP;
2032  setHitErrors();
2033  _state = kTNInit;
2034  setChi2(0.1);
2035 }
2036 //______________________________________________________________________________
2037 void StiKalmanTrackNode::initialize(StiDetector *d)
2038 {
2039  reset();
2040  _detector = d;
2041  _alpha = _detector->getPlacement()->getNormalRefAngle();
2042  _state = kTNInit;
2043  setChi2(1e10);
2044 }
2045 //______________________________________________________________________________
2046 StiKalmanTrackNode::StiKalmanTrackNode(const StiKalmanTrackNode &n)
2047 {
2048  _ext=0; _inf=0 ; *this = n;
2049 }
2050 //______________________________________________________________________________
2051 const StiKalmanTrackNode& StiKalmanTrackNode::operator=(const StiKalmanTrackNode &n)
2052 {
2053  StiTrackNode::operator=(n);
2054  memcpy(_beg,n._beg,_end-_beg+1);
2055  if (n._ext) { extend();*_ext = *n._ext;}
2056  else { if(_ext) _ext->reset(); }
2057  if (n._inf) { extinf();*_inf = *n._inf;}
2058  else { if(_inf) {BFactory::Free(_inf); _inf=0;}}
2059  return *this;
2060 }
2061 //______________________________________________________________________________
2062 void StiKalmanTrackNode::setHitErrors(const StiHit *hit)
2063 {
2064  if (!hit) hit = _hit;
2065  assert(hit);
2066  StiTrackNodeHelper::getHitErrors(hit,&mFP,&mHrr);
2067 }
2068 //______________________________________________________________________________
2069 StiHitErrs StiKalmanTrackNode::getGlobalHitErrs(const StiHit *hit) const
2070 {
2071  StiHitErrs hr;
2072  StiTrackNodeHelper::getHitErrors(hit,&mFP,&hr);
2073  hr.rotate(_alpha);
2074  return hr;
2075 }
2076 #if 1
2077 //______________________________________________________________________________
2078 int StiKalmanTrackNode::testError(double *emx, int begend)
2079 {
2080 // Test and correct error matrix. Output : number of fixes
2081 // DO NOT IMPROVE weird if() here. This accounts NaN
2082 
2083 
2084  static int nCall=0; nCall++;
2085  static const double dia[6] = { 1000.,1000., 1000.,1000.,1000,1000.};
2086  static double emxBeg[kNErrs];
2087 //return 0;
2088  if (!begend) { memcpy(emxBeg,emx,sizeof(emxBeg));}
2089  int ians=0,j1,j2,jj;
2090  for (j1=0; j1<5;j1++){
2091  jj = idx55[j1][j1];
2092  if (!(emx[jj]>0)) {;
2093  LOG_DEBUG << Form("<StiKalmanTrackNode::testError> Negative diag %g[%d][%d]",emx[jj],j1,j1)
2094  << endm;
2095  continue;}
2096  if (emx[jj]<=10*dia[j1] ) continue;
2097  assert(finite(emx[jj]));
2098  LOG_DEBUG << Form("<StiKalmanTrackNode::testError> Huge diag %g[%d][%d]",emx[jj],j1,j1)
2099  << endm;
2100  continue;
2101 // ians++; emx[jj]=dia[j1];
2102 // for (j2=0; j2<5;j2++){if (j1!=j2) emx[idx55[j1][j2]]=0;}
2103  }
2104  for (j1=0; j1< 5;j1++){
2105  for (j2=0; j2<j1;j2++){
2106  jj = idx55[j1][j2];
2107  assert(finite(emx[jj]));
2108  double cormax = emx[idx55[j1][j1]]*emx[idx55[j2][j2]];
2109  if (emx[jj]*emx[jj]<cormax) continue;
2110  cormax= sqrt(cormax);
2111 // ians++;emx[jj]= (emx[jj]<0) ? -cormax:cormax;
2112  }}
2113  return ians;
2114 }
2115 #endif//0
2116 //______________________________________________________________________________
2117 void StiKalmanTrackNode::numeDeriv(double val,int kind,int shape,int dir)
2118 {
2119  double maxStep[kNPars]={0.1,0.1,0.1,0.01,0.001,0.01};
2120  if (fDerivTestOn<0) return;
2121  gCurrShape = shape;
2122  fDerivTestOn=-1;
2123  double save[20];
2124  StiKalmanTrackNode myNode;
2125  double *pars = &myNode.mFP.x();
2126  int state=0;
2127  saveStatics(save);
2128  if (fabs(mFP.curv())> 0.02) goto FAIL;
2129  int ipar;
2130  for (ipar=1;ipar<kNPars;ipar++)
2131  {
2132  for (int is=-1;is<=1;is+=2) {
2133  myNode = *this;
2134  backStatics(save);
2135  double step = 0.1*sqrt((mFE.G())[idx66[ipar][ipar]]);
2136  if (step>maxStep[ipar]) step = maxStep[ipar];
2137 // if (step>0.1*fabs(pars[ipar])) step = 0.1*pars[ipar];
2138 // if (fabs(step)<1.e-7) step = 1.e-7;
2139  pars[ipar] +=step*is;
2140 // Update sinCA & cosCA
2141  myNode.mFP._sinCA = sin(myNode.mFP.eta());
2142  if (fabs(myNode.mFP._sinCA) > 0.9) goto FAIL;
2143  myNode.mFP._cosCA = cos(myNode.mFP.eta());
2144 
2145  switch (kind) {
2146  case 1: //propagate
2147  state = myNode.propagate(val,shape,dir); break;
2148  case 2: //rotate
2149  state = myNode.rotate(val);break;
2150  default: assert(0);
2151  }
2152  if(state ) goto FAIL;
2153 
2154  for (int jpar=1;jpar<kNPars;jpar++) {
2155  if (is<0) {
2156  fDerivTest[jpar][ipar]= pars[jpar];
2157  } else {
2158  double tmp = fDerivTest[jpar][ipar];
2159  fDerivTest[jpar][ipar] = (pars[jpar]-tmp)/(2.*step);
2160  if (ipar==jpar) fDerivTest[jpar][ipar]-=1.;
2161  }
2162  }
2163  }
2164  }
2165  fDerivTestOn=1;backStatics(save);return;
2166 FAIL:
2167  fDerivTestOn=0;backStatics(save);return;
2168 }
2169 //______________________________________________________________________________
2170 int StiKalmanTrackNode::testDeriv(double *der)
2171 {
2172  if (fDerivTestOn!=1) return 0;
2173  double *num = fDerivTest[0];
2174  int nerr = 0;
2175  for (int i=1;i<kNErrs;i++) {
2176  int ipar = i/kNPars;
2177  int jpar = i%kNPars;
2178  if (ipar==jpar) continue;
2179  if (ipar==0) continue;
2180  if (jpar==0) continue;
2181  double dif = fabs(num[i]-der[i]);
2182  if (fabs(dif) <= 1.e-5) continue;
2183  if (fabs(dif) <= 0.2*0.5*fabs(num[i]+der[i])) continue;
2184  nerr++;
2185  LOG_DEBUG << Form("***Wrong deriv [%d][%d] = %g %g %g",ipar,jpar,num[i],der[i],dif)
2186  << endm;
2187  }
2188  fDerivTestOn=0;
2189  return nerr;
2190 }
2191 //______________________________________________________________________________
2192 void StiKalmanTrackNode::saveStatics(double *sav)
2193 {
2194  sav[ 0]=mgP.x1;
2195  sav[ 1]=mgP.x2;
2196  sav[ 2]=mgP.y1;
2197  sav[ 3]=mgP.y2;
2198  sav[ 5]=mgP.dx;
2199  sav[ 6]=mgP.cosCA1;
2200  sav[ 7]=mgP.sinCA1;
2201  sav[ 8]=mgP.cosCA2;
2202  sav[ 9]=mgP.sinCA2;
2203  sav[10]=mgP.sumSin;
2204  sav[11]=mgP.sumCos;
2205  sav[14]=mgP.dl;
2206  sav[15]=mgP.dl0;
2207  sav[16]=mgP.dy;
2208 }
2209 //______________________________________________________________________________
2210 void StiKalmanTrackNode::backStatics(double *sav)
2211 {
2212  mgP.x1= sav[ 0];
2213  mgP.x2= sav[ 1];
2214  mgP.y1= sav[ 2];
2215  mgP.y2= sav[ 3];
2216  mgP.dx= sav[ 5];
2217  mgP.cosCA1= sav[ 6];
2218  mgP.sinCA1= sav[ 7];
2219  mgP.cosCA2= sav[ 8];
2220  mgP.sinCA2= sav[ 9];
2221  mgP.sumSin= sav[10];
2222  mgP.sumCos= sav[11];
2223  mgP.dl= sav[14];
2224  mgP.dl0= sav[15];
2225  mgP.dy= sav[16];
2226 }
2227 //________________________________________________________________________________
2228 void StiKalmanTrackNode::PrintpT(const Char_t *opt) const {
2229  // opt = "E" extapolation
2230  // "M" Multiple scattering
2231  // "V" at Vertex
2232  // "B" at beam
2233  // "R" at Radius
2234  // "U" Updated
2235  // mFP fit parameters
2236  // mFE fit errors
2237  // _ext->mPP
2238  // _ext->mPE
2239  // _ext->mMtx
2240  Double_t dpTOverpT = 100*TMath::Sqrt(mFE._cPP/(mFP.ptin()*mFP.ptin()));
2241  if (dpTOverpT > 9999.9) dpTOverpT = 9999.9;
2242  comment += ::Form(" %s pT %8.3f+-%6.1f sy %6.4f",opt,getPt(),dpTOverpT,TMath::Sqrt(mFE._cYY));
2243 }
2244 //________________________________________________________________________________
2245 void StiKalmanTrackNode::PrintStep() {
2246  LOG_INFO << comment.Data() << "\t" << commentdEdx.Data() << endm;
2247  ResetComment();
2248 }
2249 //________________________________________________________________________________
2250 int StiKalmanTrackNode::print(const char *opt) const
2251 {
2252 static const char *txt = "xyzeptchrXYZED";
2253 //locals xyz, e=Psi,p=1/pt, c=curvature, h=mag field, r=rxy
2254 //global XYZ, E=Psi D= direction +=outside
2255 
2256 static const char *hhh = "uvwUVW";
2257 static const char *HHH = "xyzXYZ";
2258  if (!opt || !opt[0]) opt = "2xh";
2259  StiHit *hit = getHit();
2260  TString ts;
2261  if (!isValid()) ts+="*";
2262  if (hit) {ts+=(getChi2()>1e3)? "h":"H";}
2263  if (hit && strchr(opt,'H')) {
2264  printf("%p(%s=%p)",(void*)this,ts.Data(),hit);
2265  } else {
2266  printf("%p(%s)",(void*)this,ts.Data());
2267  }
2268 
2269  if (strchr(opt,'2')) printf("\t%s=%g","ch2",getChi2());
2270  double val=-9999;
2271  for (int i=0;txt[i];i++) {
2272  if (!strchr(opt,txt[i])) continue;
2273  switch(i) {
2274  case 0:;case 1:;case 2:; case 3:;case 4:;case 5:;case 6:;case 7:;
2275  val = mFP[i]; break;
2276  case 8: val = mFP.rxy(); break;
2277  case 9: val = x_g(); break;
2278  case 10: val = y_g(); break;
2279  case 11: val = z_g(); break;
2280  case 12: val = getPsi(); break;
2281  case 13: val = mFP._cosCA*mFP[0]+mFP._sinCA*mFP[1];
2282 
2283  }
2284  printf("\t%c=%g",txt[i],val);
2285  }//end for i
2286 
2287  for (int i=0;hit && hhh[i];i++) {
2288  if (!strchr(opt,hhh[i])) continue;
2289  switch(i) {
2290  case 0:val = hit->x(); break;
2291  case 1:val = hit->y(); break;
2292  case 2:val = hit->z(); break;
2293  case 3:val = hit->x_g(); break;
2294  case 4:val = hit->y_g(); break;
2295  case 5:val = hit->z_g(); break;
2296  }
2297  printf("\th%c=%g",HHH[i],val);
2298  }
2299  if (isValid()) printf(" valid");
2300  if (getDetector()) {printf(" %s",getDetector()->getName().c_str());}
2301  printf("\n");
2302  return 1;
2303 }
2304 //________________________________________________________________________________
2305 StiNodeExt *StiKalmanTrackNode::nodeExtInstance()
2306 {
2307 static StiFactory<StiNodeExt,StiNodeExt> *extFactory=0;
2308  if (!extFactory) {
2309  extFactory = StiFactory<StiNodeExt,StiNodeExt>::myInstance();
2310  extFactory->setMaxIncrementCount(400000);
2311  extFactory->setFastDelete();
2312  }
2313  return extFactory->getInstance();
2314 }
2315 //________________________________________________________________________________
2316 StiNodeInf *StiKalmanTrackNode::nodeInfInstance()
2317 {
2318 static StiFactory<StiNodeInf,StiNodeInf> *infFactory=0;
2319  if (!infFactory) {
2320  infFactory = StiFactory<StiNodeInf,StiNodeInf>::myInstance();
2321  infFactory->setMaxIncrementCount(40000);
2322  infFactory->setFastDelete();
2323  }
2324  return infFactory->getInstance();
2325 }
2326 //________________________________________________________________________________
2327 void StiKalmanTrackNode::extend()
2328 {
2329  if(_ext) return;
2330  _ext = nodeExtInstance();
2331 }
2332 //________________________________________________________________________________
2333 void StiKalmanTrackNode::extinf()
2334 {
2335  if(_inf) return;
2336  _inf = nodeInfInstance();
2337 }
2338 //________________________________________________________________________________
2339 void StiKalmanTrackNode::saveInfo(int kase)
2340 {
2341  if (kase){};
2342  extinf();
2343  _inf->mPP = mPP();
2344  _inf->mPE = mPE();
2345  _inf->mHrr = mHrr;
2346 }
2347 
2348 //________________________________________________________________________________
2349 void StiKalmanTrackNode::reduce()
2350 {
2351  if(!_ext) return;
2352  BFactory::Free(_ext); _ext=0;
2353 }
2354 
2355 
2356 //________________________________________________________________________________
2357 StThreeVector<double> StiKalmanTrackNode::getPoint() const
2358 {
2359  return StThreeVector<double>(mFP.x(),mFP.y(),mFP.z());
2360 }
2361 
2362 //________________________________________________________________________________
2363 StThreeVector<double> StiKalmanTrackNode::getGlobalPoint() const
2364 {
2365  double ca = cos(_alpha),sa=sin(_alpha);
2366  return StThreeVector<double>(ca*mFP.x()-sa*mFP.y(), sa*mFP.x()+ca*mFP.y(), mFP.z());
2367 }
2368 
2369 //________________________________________________________________________________
2370 double StiKalmanTrackNode::x_g() const
2371 {
2372  return cos(_alpha)*mFP.x()-sin(_alpha)*mFP.y();
2373 }
2374 
2375 //________________________________________________________________________________
2376 double StiKalmanTrackNode::y_g() const
2377 {
2378  return sin(_alpha)*mFP.x()+cos(_alpha)*mFP.y();
2379 }
2380 
2381 //________________________________________________________________________________
2382 double StiKalmanTrackNode::z_g() const
2383 {
2384  return mFP.z();
2385 }
2386 
2387 //________________________________________________________________________________
2388 void StiKalmanTrackNode::setUntouched()
2389 {
2390  mUnTouch.set(mPP(),mPE());
2391 }
2392 //________________________________________________________________________________
2393 double StiKalmanTrackNode::getTime() const
2394 {
2395  static const double smax = 1e3;
2396  double time = 0;
2397  if (! _laser) {
2398  double d = sqrt(mFP.x()*mFP.x()+mFP.y()*mFP.y());
2399  double sn = fabs(mFP._cosCA*mFP.y() - mFP._sinCA*mFP.x())/d;
2400  if (sn> 0.99) sn = 0.99;
2401  if (sn<0.2) {
2402  d *= (1.+sn*sn/6);
2403  } else {
2404  d *= asin(sn)/sn;
2405  }
2406  d *= sqrt(1.+mFP.tanl()*mFP.tanl());
2407  double beta = 1;
2408  double pt = fabs(mFP.ptin());
2409  if (pt>0.1) {
2410  pt = 1./pt;
2411  double p2=(1.+mFP.tanl()*mFP.tanl())*pt*pt;
2412  double m=StiKalmanTrackFinderParameters::instance()->massHypothesis();
2413  double m2=m*m;
2414  double e2=p2+m2;
2415  double beta2=p2/e2;
2416  beta = TMath::Sqrt(beta2);
2417  }
2418  time = d/(TMath::Ccgs()*beta*1e-6); // mksec
2419  } else {
2420  if (TMath::Abs(mFP.z()) > 20.0) {
2421 static const double Radius = 197.;
2422 static const int nsurf = 6;
2423 static const double surf[6] = {-Radius*Radius, 0, 0, 0, 1, 1};
2424  double dir[3] = {mFP._cosCA,mFP._sinCA,mFP.tanl()};
2425  THelixTrack tc(mFP.P,dir,mFP.curv());
2426  double s = tc.Path(smax, surf, nsurf,0,0,1);
2427  if (TMath::Abs(s) < smax)
2428  time = TMath::Abs(s)/(TMath::Ccgs()*1e-6); // mksec
2429  }
2430  }
2431  return time;
2432 }
2433 //______________________________________________________________________________
2435 double StiKalmanTrackNode::evaluateChi2Info(const StiHit * hit) const
2436 {
2437  if (!_inf) return 1e41;
2438  const StiNodePars &myFP = _inf->mPP;
2439  const StiNodeErrs &myFE = _inf->mPE;
2440  StiHitErrs myHrr;
2441  StiTrackNodeHelper::getHitErrors(hit,&myFP,&myHrr);
2442 
2443  double r00, r01,r11,det;
2444  //If required, recalculate the errors of the detector hits.
2445  //Do not attempt this calculation for the main vertex.
2446  double dsin =myFP.curv()*(hit->x()-myFP.x());
2447  if (fabs(myFP._sinCA+dsin)>0.99) return 1e41;
2448  if (fabs(myFP.eta()) >kMaxEta) return 1e41;
2449  if (fabs(myFP.curv()) >kMaxCur) return 1e41;
2450  if (myHrr.hYY>1000*myFE._cYY
2451  && myHrr.hZZ>1000*myFE._cZZ) return 1e41;
2452  r00=myHrr.hYY+myFE._cYY;
2453  r01=myHrr.hZY+myFE._cZY;
2454  r11=myHrr.hZZ+myFE._cZZ;
2455 
2456  det=r00*r11 - r01*r01;
2457  if (det<r00*r11*1.e-5) {
2458  LOG_DEBUG << Form("StiKalmanTrackNode::evalChi2 *** zero determinant %g",_det)<< endm;
2459  return 1e60;
2460  }
2461  double tmp=r00; r00=r11; r11=tmp; r01=-r01;
2462  double deltaX = hit->x()-myFP.x();
2463  double deltaL = deltaX/myFP._cosCA;
2464  double deltaY = myFP._sinCA *deltaL;
2465  double deltaZ = myFP.tanl() *deltaL;
2466  double dyt=(myFP.y()-hit->y()) + deltaY;
2467  double dzt=(myFP.z()-hit->z()) + deltaZ;
2468  double cc= (dyt*r00*dyt + 2*r01*dyt*dzt + dzt*r11*dzt)/det;
2469  return cc;
2470 }
2471 
2472 //________________________________________________________________________________
2476 double StiKalmanTrackNode::StiKalmanTrackNode::pathlength() const
2477 {
2478  const StiDetector * det = getDetector();
2479  if (!det) return 0.;
2480  double c = mFP._cosCA;
2481  if (det->getShape()->getShapeCode()!=kPlanar) {
2482  double CA = mFP.eta()-atan2(mFP.y(),mFP.x());
2483  c = cos(CA);
2484  }
2485  double thickness = det->getShape()->getThickness();
2486  return (thickness*::sqrt(1.+mFP.tanl()*mFP.tanl())) / fabs(c);
2487 }
StiKalmanTrackNode::setState
void setState(const StiKalmanTrackNode *node)
Sets the Kalman state of this node equal to that of the given node.
Definition: StiKalmanTrackNode.cxx:751
StiKalmanTrackNode::propagate
int propagate(StiKalmanTrackNode *p, const StiDetector *tDet, int dir)
Propagates a track encapsulated by the given node &quot;p&quot; to the given detector &quot;tDet&quot;.
Definition: StiKalmanTrackNode.cxx:1079
StiKalmanTrackNode::getX0
double getX0() const
Definition: StiKalmanTrackNode.h:363
StiKalmanTrackNode::propagateToRadius
int propagateToRadius(StiKalmanTrackNode *pNode, double radius, int dir)
Definition: StiKalmanTrackNode.cxx:1204
StiKalmanTrackNode::getGlobalRadial
void getGlobalRadial(double x[6], double e[15])
Extract state information from this node in Radial representation.
Definition: StiKalmanTrackNode.cxx:910
StiKalmanTrackNode::propagateToBeam
bool propagateToBeam(const StiKalmanTrackNode *p, int dir)
Definition: StiKalmanTrackNode.cxx:1183
StiHit
Definition: StiHit.h:51
StiHit::x
const Float_t & x() const
Return the local x, y, z values.
Definition: StiHit.h:67
StiKalmanTrackNode::reset
void reset()
Resets the node to a &quot;null&quot; un-used state.
Definition: StiKalmanTrackNode.cxx:712
StiHit::rotate
void rotate(double angle)
Definition: StiHit.cxx:46
StiHit::detector
const StiDetector * detector() const
Definition: StiHit.h:96
StiKalmanTrackNode::pathLToNode
double pathLToNode(const StiKalmanTrackNode *const oNode)
Definition: StiKalmanTrackNode.cxx:1441
StiKalmanTrackNode::rotate
int rotate(double alpha)
Definition: StiKalmanTrackNode.cxx:1840
StiHit::x_g
Float_t x_g() const
Return the global x, y, z values.
Definition: StiHit.h:73
StiHitErrorCalculator::calculateError
virtual void calculateError(Double_t _z, Double_t _eta, Double_t _tanl, Double_t &ecross, Double_t &edip, Double_t fudgeFactor=1) const
coeff[6] = 0:intrinsicY 1: driftY 2: crossY 3:intrinsicZ 4: driftZ 5: crossZ
Definition: StiChairs.cxx:7
StiMaterial::getElossCalculator
StiElossCalculator * getElossCalculator() const
Get Eloss calculator.
Definition: StiMaterial.cxx:102
StiKalmanTrackNode::getGasX0
double getGasX0() const
Definition: StiKalmanTrackNode.h:373
StiKalmanTrackNode::print
int print(const char *opt) const
rotation angle of local coordinates wrt global coordinates
Definition: StiKalmanTrackNode.cxx:2250
StiKalmanTrackNode::getHelixCenter
StThreeVector< double > getHelixCenter() const
Return center of helix circle in global coordinates.
Definition: StiKalmanTrackNode.cxx:1947
StiKalmanTrackNode::getPt
double getPt() const
Calculates and returns the transverse momentum of the track at this node.
Definition: StiKalmanTrackNode.cxx:795
StiElossCalculator::calculate
double calculate(double charge2, double m, double beta2) const
Definition: StiElossCalculator.cxx:56
BFactory::Free
static void Free(void *obj)
Free an object for reuse.
Definition: Factory.h:73
TCL::trasat
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:540
StiNodePars::_cosCA
double _cosCA
sine and cosine of cross angle
Definition: StiNodePars.h:51
StiKalmanTrackNode::resetError
void resetError(double fak=0)
Resets errors for refit.
Definition: StiKalmanTrackNode.cxx:729
StiFactory::getInstance
Abstract * getInstance()
Get a pointer to instance of objects served by this factory.
Definition: StiFactory.h:114
StiKalmanTrackNode::evaluateChi2
double evaluateChi2(const StiHit *hit)
Definition: StiKalmanTrackNode.cxx:1481
StiKalmanTrackNode::getGlobalTpt
void getGlobalTpt(float x[6], float e[15])
Extract state information from this node in TPT representation.
Definition: StiKalmanTrackNode.cxx:979
StiMaterial::getX0
double getX0() const
Get the radiation length in centimeter.
Definition: StiMaterial.h:37
TDataSet::Path
virtual TString Path() const
return the full path of this data set
Definition: TDataSet.cxx:626
StiKalmanTrackNode::propagateMCS
void propagateMCS(StiKalmanTrackNode *previousNode, const StiDetector *tDet)
Definition: StiKalmanTrackNode.cxx:1535
StiKalmanTrackNode::get
void get(double &alpha, double &xRef, double x[kNPars], double cc[kNErrs], double &chi2)
Extract state information from this node.
Definition: StiKalmanTrackNode.cxx:773
StiKalmanTrackNode::mFE
StiNodeErrs mFE
covariance matrix of the track parameters
Definition: StiKalmanTrackNode.h:284
StiKalmanTrackNode::evaluateChi2Info
double evaluateChi2Info(const StiHit *hit) const
Definition: StiKalmanTrackNode.cxx:2435
StiMaterial::getDensity
double getDensity() const
Get the material density in grams/cubic centimeter.
Definition: StiMaterial.h:35
StiKalmanTrackNode::initialize
void initialize(StiHit *h)
Initialize this node with the given hit information.
Definition: StiKalmanTrackNode.cxx:2018
Named::getName
const string & getName() const
Get the name of the object.
Definition: Named.cxx:22
Generated by   doxygen 1.8.5