StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFtpcTrackingParams.cc
1 // $Id: StFtpcTrackingParams.cc,v 1.39 2013/02/18 16:30:42 fisyak Exp $
2 // $Log: StFtpcTrackingParams.cc,v $
3 // Revision 1.39 2013/02/18 16:30:42 fisyak
4 // gufld => agufld
5 //
6 // Revision 1.38 2012/11/07 23:30:18 fisyak
7 // Supress warnings
8 //
9 // Revision 1.37 2009/08/25 19:41:19 fine
10 // fix the compilation issues under SL5_64_bits gcc 4.3.2
11 //
12 // Revision 1.36 2007/12/13 10:35:21 jcs
13 // standardize logger messages
14 //
15 // Revision 1.35 2007/12/12 12:55:18 jcs
16 // Markus Oldenburg replaced assert() with a return code which can be tested in StFtpcTrackMaker
17 // replaced 'return 1' with 'return kStOK'
18 // replaced 'return 0' with 'return KStErr'
19 //
20 // Revision 1.34 2007/11/13 10:30:40 jcs
21 // add code (commented out) which enables testing rotation values without changing the database
22 //
23 // Revision 1.33 2007/05/08 10:47:33 jcs
24 // replace StMagUtilities with StarMagField as requested by Yuri
25 //
26 // Revision 1.32 2007/02/06 11:42:17 jcs
27 // move unessential output messages from INFO to DEBUG
28 //
29 // Revision 1.31 2007/01/15 08:23:02 jcs
30 // replace printf, cout and gMesMgr with Logger commands
31 //
32 // Revision 1.30 2006/09/25 14:03:42 jcs
33 // move the reconstruction parameters maxDcaVertex,minNumTracks from code to CodeParams
34 //
35 // Revision 1.29 2004/06/04 11:05:26 jcs
36 // replaced StarDb/ftpc/fdepars/fdepar with StarDb/ftpc/ftpcdEdxPars
37 //
38 // Revision 1.28 2004/04/29 03:34:49 perev
39 // minor cleanup
40 //
41 // Revision 1.27 2004/04/06 19:03:02 oldi
42 // Code cleanup.
43 //
44 // Revision 1.26 2004/04/05 06:38:46 oldi
45 // Reported problem fixed (delete -> delete[]).
46 //
47 // Revision 1.25 2004/01/28 01:41:32 jeromel
48 // *** empty log message ***
49 //
50 // Revision 1.24 2003/10/07 14:12:54 jcs
51 // use gufld to determine magnetic field factor (this is where we started from a long,long time ago)
52 //
53 // Revision 1.23 2003/10/02 00:10:37 perev
54 // Zeroing of members added and bug in ResetMagField fixed
55 //
56 // Revision 1.22 2003/09/30 00:10:29 oldi
57 // mMagField set to 0.
58 //
59 // Revision 1.21 2003/09/26 06:08:57 oldi
60 // Check if the magentic field was reversed 'by hand' with a chain option.
61 // If yes, multiply the scaleFactor of the field with -1.
62 //
63 // Revision 1.20 2003/09/11 21:31:30 jeromel
64 // removed inline as it would leave a few undefined reference
65 //
66 // Revision 1.19 2003/09/02 17:58:17 perev
67 // gcc 3.2 updates + WarnOff
68 //
69 // Revision 1.18 2003/06/27 13:11:25 putschke
70 // *** empty log message ***
71 //
72 // Revision 1.17 2003/05/23 21:10:44 oldi
73 // Cosmetics.
74 //
75 // Revision 1.16 2003/05/21 09:47:35 putschke
76 // Include rotation around y-axis for FTPC east and west
77 //
78 // Revision 1.15 2003/05/20 18:35:04 oldi
79 // Cuts for vertex estimation introduced (globDca < 1 cm, multiplicity >= 200).
80 //
81 // Revision 1.14 2003/01/29 17:24:37 oldi
82 // Message added which will be printed when StMagUtilities is initialized.
83 //
84 // Revision 1.13 2003/01/21 10:04:13 jcs
85 // initialize variables to eliminate compiler warnings for NODEBUG=yes
86 //
87 // Revision 1.12 2003/01/16 18:04:34 oldi
88 // Bugs eliminated. Now it compiles on Solaris again.
89 // Split residuals for global and primary fit.
90 //
91 // Revision 1.11 2002/11/21 15:46:25 oldi
92 // Enabled rotation for FTPC west. If there is an observed shift of the vertex
93 // position in y-direction (in FTPC west), just fill this offset into the Db.
94 // Up to now this offset is set to 0., i.e. only FTPC east is rotated (because
95 // the offset is at 0.3427 cm).
96 //
97 // Revision 1.10 2002/11/19 12:45:09 oldi
98 // A new database entry (installationPointY[east/west]) was introduced. Now
99 // the rotation of FTPC east is done around the correct axis, which isn't
100 // measured but comes from the drawings. The measurements used before were true
101 // measurements but had nothing to do with the rotation axis, unfortunately.
102 // Anyway, the difference is rather small since a typical cluster is rotated
103 // by less than 0.1mm.
104 // Some code cleanup done.
105 //
106 // Revision 1.9 2002/11/06 13:47:34 oldi
107 // All current database values hardcoded (for stand alone usage).
108 // Code clean ups.
109 //
110 // Revision 1.8 2002/10/31 13:42:26 oldi
111 // Everything read from database now.
112 //
113 // Revision 1.7 2002/10/11 15:45:43 oldi
114 // Get FTPC geometry and dimensions from database.
115 // No field fit activated: Returns momentum = 0 but fits a helix.
116 // Bug in TrackMaker fixed (events with z_vertex > outer_ftpc_radius were cut).
117 // QA histograms corrected (0 was supressed).
118 // Code cleanup (several lines of code changed due to *params -> Instance()).
119 // cout -> gMessMgr.
120 //
121 // Revision 1.6 2002/10/03 10:34:06 oldi
122 // Usage of gufld removed.
123 // Magnetic field is read by StMagUtilities, now.
124 //
125 // Revision 1.5 2002/08/02 11:19:32 oldi
126 // MaxDCA is set to 100 cm now. Therefore 'every' track is fitted with the
127 // primary vertex in addition to the global fit, which is perfpormed anyway.
128 // The cut, which was set to 2.5 cm before, has to be applied during the
129 // analysis, now.
130 // The y-offest of FTPC east changed slightly due to the 'new' t0 of 2.96 mus.
131 //
132 // Revision 1.4 2002/06/07 06:00:39 oldi
133 // New value for rotation angle of FTPC east after temperature offset was corrected.
134 //
135 
136 //----------Author: Markus D. Oldenburg
137 //----------Last Modified: 25.04.2002
138 //----------Copyright: &copy MDO Production 2002
139 
140 #include "StFtpcTrackingParams.hh"
141 #include "SystemOfUnits.h"
142 #include "StMessMgr.h"
143 
144 #ifndef ST_NO_NAMESPACES
145 using namespace units;
146 #endif
147 
148 #include "TMath.h"
149 #include <Stiostream.h>
150 
152 // StFtpcTrackingParams //
153 // //
154 // This class provides the necessary parameters for FTPC tracking. //
156 
157 ClassImp(StFtpcTrackingParams)
158 
159 class St_ftpcdEdxPars;
160 
161 // Initialization of instance
162 StFtpcTrackingParams* StFtpcTrackingParams::mInstance = 0;
163 
164 
165 // FTPC geometry
166 Double_t StFtpcTrackingParams::InnerRadius() { return mInnerRadius; }
167 Double_t StFtpcTrackingParams::OuterRadius() { return mOuterRadius; }
168  Int_t StFtpcTrackingParams::NumberOfPadRows() { return mNumberOfPadRows; }
169  Int_t StFtpcTrackingParams::NumberOfPadRowsPerSide() { return mNumberOfPadRowsPerSide; }
170 Double_t StFtpcTrackingParams::PadRowPosZ(Int_t row) { return mPadRowPosZ[row]; }
171 
172 // Vertex position
173 Double_t StFtpcTrackingParams::MaxVertexPosZWarning() { return mMaxVertexPosZWarning; }
174 Double_t StFtpcTrackingParams::MaxVertexPosZError() { return mMaxVertexPosZError; }
175 
176 // Vertex reconstruction
177  Int_t StFtpcTrackingParams::HistoBins() { return mHistoBins; }
178 Double_t StFtpcTrackingParams::HistoMin() { return mHistoMin; }
179 Double_t StFtpcTrackingParams::HistoMax() { return mHistoMax; }
180 Double_t StFtpcTrackingParams::MaxDcaVertex() { return mMaxDcaVertex; }
181  Int_t StFtpcTrackingParams::MinNumTracks() { return mMinNumTracks; }
182 
183 // Tracker
184 Int_t StFtpcTrackingParams::RowSegments() { return mRowSegments; }
185 Int_t StFtpcTrackingParams::PhiSegments() { return mPhiSegments; }
186 Int_t StFtpcTrackingParams::EtaSegments() { return mEtaSegments; }
187 
188 // Tracking
189  Bool_t StFtpcTrackingParams::Laser(Int_t tracking_method) { return mLaser[tracking_method]; }
190  Bool_t StFtpcTrackingParams::VertexConstraint(Int_t tracking_method) { return mVertexConstraint[tracking_method]; }
191  Int_t StFtpcTrackingParams::MaxTrackletLength(Int_t tracking_method) { return mMaxTrackletLength[tracking_method]; }
192  Int_t StFtpcTrackingParams::MinTrackLength(Int_t tracking_method) { return mMinTrackLength[tracking_method]; }
193  Int_t StFtpcTrackingParams::RowScopeTracklet(Int_t tracking_method) { return mRowScopeTracklet[tracking_method]; }
194  Int_t StFtpcTrackingParams::RowScopeTrack(Int_t tracking_method) { return mRowScopeTrack[tracking_method]; }
195  Int_t StFtpcTrackingParams::PhiScope(Int_t tracking_method) { return mPhiScope[tracking_method]; }
196  Int_t StFtpcTrackingParams::EtaScope(Int_t tracking_method) { return mEtaScope[tracking_method]; }
197 Double_t StFtpcTrackingParams::MaxDca(Int_t tracking_method) { return mMaxDca[tracking_method]; }
198 
199 // Tracklets
200 Double_t StFtpcTrackingParams::MaxAngleTracklet(Int_t tracking_method) { return mMaxAngleTracklet[tracking_method]; }
201 
202 // Tracks
203 Double_t StFtpcTrackingParams::MaxAngleTrack(Int_t tracking_method) { return mMaxAngleTrack[tracking_method]; }
204 Double_t StFtpcTrackingParams::MaxCircleDist(Int_t tracking_method) { return mMaxCircleDist[tracking_method]; }
205 Double_t StFtpcTrackingParams::MaxLengthDist(Int_t tracking_method) { return mMaxLengthDist[tracking_method]; }
206 
207 // Split tracks
208 Double_t StFtpcTrackingParams::MaxDist() { return mMaxDist; }
209 Double_t StFtpcTrackingParams::MinPointRatio() { return mMinPointRatio; }
210 Double_t StFtpcTrackingParams::MaxPointRatio() { return mMaxPointRatio; }
211 
212 // dE/dx
213 Int_t StFtpcTrackingParams::DebugLevel() { return mDebugLevel; }
214 Int_t StFtpcTrackingParams::IdMethod() { return mIdMethod; }
215 Int_t StFtpcTrackingParams::NoAngle() { return mNoAngle; }
216 Int_t StFtpcTrackingParams::MaxHit() { return mMaxHit; }
217 Int_t StFtpcTrackingParams::MinHit() { return mMinHit; }
218 Int_t StFtpcTrackingParams::MaxTrack() { return mMaxTrack; }
219 
220 Double_t StFtpcTrackingParams::PadLength() { return mPadLength; }
221 Double_t StFtpcTrackingParams::FracTrunc() { return mFracTrunc; }
222 Double_t StFtpcTrackingParams::Aip() { return mAip; }
223 Double_t StFtpcTrackingParams::ALargeNumber() { return mALargeNumber; }
224 
225 // transformation due to rotated and displaced TPC
226  StMatrixD StFtpcTrackingParams::TpcToGlobalRotation() { return mTpcToGlobalRotation; }
227  StMatrixD StFtpcTrackingParams::GlobalToTpcRotation() { return mGlobalToTpcRotation; }
228 StThreeVectorD StFtpcTrackingParams::TpcPositionInGlobal() { return mTpcPositionInGlobal; }
229 
230 StMatrixD StFtpcTrackingParams::FtpcRotation(Int_t i) { return *mFtpcRotation[i]; }
231 StMatrixD StFtpcTrackingParams::FtpcRotationInverse(Int_t i) { return *mFtpcRotationInverse[i]; }
232 StMatrixD StFtpcTrackingParams::FtpcRotationX(Int_t i) { return *mFtpcRotationX[i]; }
233 StMatrixD StFtpcTrackingParams::FtpcRotationXInverse(Int_t i) { return *mFtpcRotationXInverse[i]; }
234 StMatrixD StFtpcTrackingParams::FtpcRotationY(Int_t i) { return *mFtpcRotationY[i]; }
235 StMatrixD StFtpcTrackingParams::FtpcRotationYInverse(Int_t i) { return *mFtpcRotationYInverse[i]; }
236 
237  Double_t StFtpcTrackingParams::InstallationPointX(Int_t i) { return (i>=0 && i<=1) ? mInstallationPointX[i] : 0.; }
238  Double_t StFtpcTrackingParams::InstallationPointY(Int_t i) { return (i>=0 && i<=1) ? mInstallationPointY[i] : 0.; }
239  Double_t StFtpcTrackingParams::InstallationPointZ(Int_t i) { return (i>=0 && i<=1) ? mInstallationPointZ[i] : 0.; }
240  Double_t StFtpcTrackingParams::ObservedVertexOffsetY(Int_t i) { return (i>=0 && i<=1) ? mObservedVertexOffsetY[i] : 0.; }
241  Double_t StFtpcTrackingParams::ObservedVertexOffsetX(Int_t i) { return (i>=0 && i<=1) ? mObservedVertexOffsetX[i] : 0.; }
242 
243 
244 StFtpcTrackingParams* StFtpcTrackingParams::Instance(Bool_t debug,
245  St_ftpcTrackingPars *trackPars,
246  St_ftpcdEdxPars *dEdxPars,
247  St_ftpcDimensions *dimensions,
248  St_ftpcPadrowZ *padrow_z) {
249  // makes new instance or returns old one if it exists already
250 
251  if (!mInstance) {
252  mInstance = new StFtpcTrackingParams(trackPars, dEdxPars, dimensions, padrow_z);
253  }
254 
255  return mInstance;
256 }
257 
258 
259 StFtpcTrackingParams* StFtpcTrackingParams::Instance(Bool_t debug, St_ftpcCoordTrans *ftpcCoordTrans) {
260  // updates magnetic field, if necessary
261 
262  mInstance->InitCoordTransformation(ftpcCoordTrans->GetTable()); // Has to be invoked here, because it could change from run to run.
263  mInstance->InitSpaceTransformation(); // Has to be invoked here, since gStTpcDb isn't set before.
264 
265  if (debug) {
266  mInstance->PrintParams();
267  }
268 
269  return mInstance;
270 }
271 
272 
273 StFtpcTrackingParams* StFtpcTrackingParams::Instance() {
274  // return instance only
275 
276  return mInstance;
277 }
278 
279 
280 StFtpcTrackingParams::StFtpcTrackingParams(St_ftpcTrackingPars *trackPars,
281  St_ftpcdEdxPars *dEdxPars,
282  St_ftpcDimensions *dimensions,
283  St_ftpcPadrowZ *zrow)
284  : mReturnCode(0), mTpcToGlobalRotation(3, 3, 1), mGlobalToTpcRotation(3, 3, 1)
285 {
286  // zero everythig what possible
287  memset(&mStart,0,&mEnd-&mStart+1);
288 
289  // default constructor
290 
291  mFtpcRotation[0] = new StMatrixD(3, 3, 1);
292  mFtpcRotation[1] = new StMatrixD(3, 3, 1);
293  mFtpcRotationInverse[0] = new StMatrixD(3, 3, 1);
294  mFtpcRotationInverse[1] = new StMatrixD(3, 3, 1);
295  mFtpcRotationX[0] = new StMatrixD(3, 3, 1);
296  mFtpcRotationX[1] = new StMatrixD(3, 3, 1);
297  mFtpcRotationXInverse[0] = new StMatrixD(3, 3, 1);
298  mFtpcRotationXInverse[1] = new StMatrixD(3, 3, 1);
299  mFtpcRotationY[0] = new StMatrixD(3, 3, 1);
300  mFtpcRotationY[1] = new StMatrixD(3, 3, 1);
301  mFtpcRotationYInverse[0] = new StMatrixD(3, 3, 1);
302  mFtpcRotationYInverse[1] = new StMatrixD(3, 3, 1);
303 
304 
305  InitTrackingParams(trackPars->GetTable());
306  InitdEdx(dEdxPars->GetTable());
307  InitDimensions(dimensions->GetTable());
308  InitPadRows(zrow->GetTable());
309  InitCoordTransformation();
310 }
311 
312 
313 StFtpcTrackingParams::StFtpcTrackingParams()
314  : mReturnCode(0), mTpcToGlobalRotation(3, 3, 1), mGlobalToTpcRotation(3, 3, 1)
315 {
316  // zero everythig what is possible
317  memset(&mStart,0,&mEnd-&mStart+1);
318 
319  // Initialization with hardcoded values.
320 
321  mFtpcRotation[0] = new StMatrixD(3, 3, 1);
322  mFtpcRotation[1] = new StMatrixD(3, 3, 1);
323  mFtpcRotationInverse[0] = new StMatrixD(3, 3, 1);
324  mFtpcRotationInverse[1] = new StMatrixD(3, 3, 1);
325  mFtpcRotationX[0] = new StMatrixD(3, 3, 1);
326  mFtpcRotationX[1] = new StMatrixD(3, 3, 1);
327  mFtpcRotationXInverse[0] = new StMatrixD(3, 3, 1);
328  mFtpcRotationXInverse[1] = new StMatrixD(3, 3, 1);
329  mFtpcRotationY[0] = new StMatrixD(3, 3, 1);
330  mFtpcRotationY[1] = new StMatrixD(3, 3, 1);
331  mFtpcRotationYInverse[0] = new StMatrixD(3, 3, 1);
332  mFtpcRotationYInverse[1] = new StMatrixD(3, 3, 1);
333 
334  // FTPC geometry
335  mInnerRadius = 7.73 * centimeter;
336  mOuterRadius = 30.05 * centimeter;
337 
338  mNumberOfPadRows = 20;
339  mNumberOfPadRowsPerSide = 10;
340 
341  mPadRowPosZ = new Double_t[mNumberOfPadRows];
342 
343  mPadRowPosZ[0] = 162.75 * centimeter;
344  mPadRowPosZ[1] = 171.25 * centimeter;
345  mPadRowPosZ[2] = 184.05 * centimeter;
346  mPadRowPosZ[3] = 192.55 * centimeter;
347  mPadRowPosZ[4] = 205.35 * centimeter;
348  mPadRowPosZ[5] = 213.85 * centimeter;
349  mPadRowPosZ[6] = 226.65 * centimeter;
350  mPadRowPosZ[7] = 235.15 * centimeter;
351  mPadRowPosZ[8] = 247.95 * centimeter;
352  mPadRowPosZ[9] = 256.45 * centimeter;
353  mPadRowPosZ[10] = -162.75 * centimeter;
354  mPadRowPosZ[11] = -171.25 * centimeter;
355  mPadRowPosZ[12] = -184.05 * centimeter;
356  mPadRowPosZ[13] = -192.55 * centimeter;
357  mPadRowPosZ[14] = -205.35 * centimeter;
358  mPadRowPosZ[15] = -213.85 * centimeter;
359  mPadRowPosZ[16] = -226.65 * centimeter;
360  mPadRowPosZ[17] = -235.15 * centimeter;
361  mPadRowPosZ[18] = -247.95 * centimeter;
362  mPadRowPosZ[19] = -256.45 * centimeter;
363 
364  // Vertex position
365  mMaxVertexPosZWarning = 50 * centimeter;
366  mMaxVertexPosZError = 100 * centimeter;
367 
368  // Vertex reconstruction
369  mHistoBins = 300;
370  mHistoMin = -75. * centimeter;
371  mHistoMax = 75. * centimeter;
372  mMaxDcaVertex = 100. * centimeter;
373  mMinNumTracks = 1; // must be >0 !
374 
375  // Tracker
376  mRowSegments = 20;
377  mPhiSegments = 100;
378  mEtaSegments = 200;
379 
380  // Tracking
381  mLaser[0] = 0;
382  mLaser[1] = 0;
383  mLaser[2] = 1;
384  mLaser[3] = 1;
385 
386  mVertexConstraint[0] = 1;
387  mVertexConstraint[1] = 0;
388  mVertexConstraint[2] = 1;
389  mVertexConstraint[3] = 0;
390 
391  mMaxTrackletLength[0] = 3;
392  mMaxTrackletLength[1] = 3;
393  mMaxTrackletLength[2] = 10;
394  mMaxTrackletLength[3] = 10;
395 
396  mMinTrackLength[0] = 5;
397  mMinTrackLength[1] = 5;
398  mMinTrackLength[2] = 5;
399  mMinTrackLength[3] = 5;
400 
401  mRowScopeTracklet[0] = 2;
402  mRowScopeTracklet[1] = 2;
403  mRowScopeTracklet[2] = 2;
404  mRowScopeTracklet[3] = 3;
405 
406  mRowScopeTrack[0] = 3;
407  mRowScopeTrack[1] = 3;
408  mRowScopeTrack[2] = 3;
409  mRowScopeTrack[3] = 3;
410 
411  mPhiScope[0] = 1;
412  mPhiScope[1] = 1;
413  mPhiScope[2] = 1;
414  mPhiScope[3] = 2;
415 
416  mEtaScope[0] = 1;
417  mEtaScope[1] = 3;
418  mEtaScope[2] = 3;
419  mEtaScope[3] = 15;
420 
421  mMaxDca[0] = 100. * centimeter;
422  mMaxDca[1] = 100. * centimeter;
423  mMaxDca[2] = 100. * centimeter;
424  mMaxDca[3] = 100. * centimeter;
425 
426  // Tracklets
427  mMaxAngleTracklet[0] = 0.015 * radian;
428  mMaxAngleTracklet[1] = 0.03 * radian;
429  mMaxAngleTracklet[2] = 0.03 * radian;
430  mMaxAngleTracklet[3] = 0.05 * radian;
431 
432  // Tracks
433  mMaxAngleTrack[0] = 0.03 * radian;
434  mMaxAngleTrack[1] = 0.08 * radian;
435  mMaxAngleTrack[2] = 0.007 * radian;
436  mMaxAngleTrack[3] = 0.007 * radian;
437 
438  mMaxCircleDist[0] = 0.05 * 1./centimeter;
439  mMaxCircleDist[1] = 0.05 * 1./centimeter;
440  mMaxCircleDist[2] = 0.03 * 1./centimeter;
441  mMaxCircleDist[3] = 0.03 * 1./centimeter;
442 
443  mMaxLengthDist[0] = 30. * centimeter;
444  mMaxLengthDist[1] = 70. * centimeter;
445  mMaxLengthDist[2] = 30. * centimeter;
446  mMaxLengthDist[3] = 30. * centimeter;
447 
448  // Split tracks
449  mMaxDist = 0.11 * centimeter;
450  mMinPointRatio = 0.5;
451  mMaxPointRatio = 0.5;
452 
453  // dE/dx
454  mDebugLevel = 10;
455  mIdMethod = 0;
456  mNoAngle = 0;
457  mMaxHit = 10;
458  mMinHit = 4;
459  mMaxTrack = 100000;
460 
461  mPadLength = 0.02;
462  mFracTrunc = 0.8;
463  mAip = 2.6e-08;
464  mALargeNumber = 1.0e+10;
465 
466  // transformation due to rotated and displaced TPC
467  mTpcPositionInGlobal.setX(-0.256 * centimeter);
468  mTpcPositionInGlobal.setY(-0.082 * centimeter);
469  mTpcPositionInGlobal.setZ(-0.192 * centimeter);
470 
471  Double_t phi = 0.0 * radian; //large uncertainty, so set to 0
472  Double_t theta = -0.000381 * radian;
473  Double_t psi = -0.000156 * radian;
474 
475  mGlobalToTpcRotation(1, 1) = TMath::Cos(theta) * TMath::Cos(phi);
476  mGlobalToTpcRotation(1, 2) = TMath::Cos(theta) * TMath::Sin(phi);
477  mGlobalToTpcRotation(1, 3) = - TMath::Sin(theta);
478  mGlobalToTpcRotation(2, 1) = TMath::Sin(psi) * TMath::Sin(theta) * TMath::Cos(phi) - TMath::Cos(psi) * TMath::Sin(phi);
479  mGlobalToTpcRotation(2, 2) = TMath::Sin(psi) * TMath::Sin(theta) * TMath::Sin(phi) + TMath::Cos(psi) * TMath::Cos(phi);
480  mGlobalToTpcRotation(2, 3) = TMath::Cos(theta) * TMath::Sin(psi);
481  mGlobalToTpcRotation(3, 1) = TMath::Cos(psi) * TMath::Sin(theta) * TMath::Cos(phi) + TMath::Sin(psi) * TMath::Sin(phi);
482  mGlobalToTpcRotation(3, 2) = TMath::Cos(psi) * TMath::Sin(theta) * TMath::Sin(phi) - TMath::Sin(psi) * TMath::Cos(phi);
483  mGlobalToTpcRotation(3, 3) = TMath::Cos(theta) * TMath::Cos(psi);
484 
485  size_t ierr;
486  mTpcToGlobalRotation = mGlobalToTpcRotation.inverse(ierr);
487 
488  if (ierr!=0) {
489  LOG_WARN << "Cant invert rotation matrix!" << endm;
490  LOG_WARN << "Global to TPC rotation matrix:" << mGlobalToTpcRotation << endm;
491  LOG_WARN << "TPC to global rotation matrix:" << mTpcToGlobalRotation << endm;
492  }
493 
494  // internal FTPC rotation [has do be done before local -> global]
495  mInstallationPointX[0] = 0.0 * centimeter;
496  mInstallationPointX[1] = 0.0 * centimeter;
497  mInstallationPointY[0] = -22.25 * centimeter;
498  mInstallationPointY[1] = -22.25 * centimeter;
499  mInstallationPointZ[0] = -234.325 * centimeter;
500  mInstallationPointZ[1] = 234.325 * centimeter;
501 
502  mObservedVertexOffsetY[0] = 0.1 * centimeter;
503  mObservedVertexOffsetY[1] = -0.034 * centimeter;
504 
505  for (Int_t i = 0; i <= 1; i++) { // east and west rotation
506 
507  // define rotation angle
508  Double_t zShift, phi0, phi1, alpha;
509  Double_t pq = TMath::Sqrt(TMath::Power(mObservedVertexOffsetY[i], 2.)
510  - 2.*mObservedVertexOffsetY[i]*mInstallationPointY[i]
511  + TMath::Power(mInstallationPointZ[i], 2.)); // p-q-formula
512 
513  // Initialize variable to avoid compiler warning
514  alpha = 0;
515 
516  if (i == 0) { // east
517  zShift = (mInstallationPointZ[i] + pq) * centimeter; // take correct solution of p-q-formula
518  phi0 = TMath::ATan((mInstallationPointY[i] - mObservedVertexOffsetY[i]) / mInstallationPointZ[i]) * radian;
519  phi1 = TMath::ATan(mInstallationPointY[i] / (mInstallationPointZ[i] - zShift)) * radian;
520  alpha = phi0 - phi1;
521  }
522 
523  else if (i == 1) { //west
524  zShift = (mInstallationPointZ[i] - pq) * centimeter; // take correct solution of p-q-formula
525  phi0 = TMath::ATan((mObservedVertexOffsetY[i] - mInstallationPointY[i]) / mInstallationPointZ[i]) * radian;
526  phi1 = TMath::ATan(mInstallationPointY[i] / (zShift - mInstallationPointZ[i]));
527  alpha = phi1 - phi0;
528  }
529 
530  // define rotation axis
531  // simplify to rotation about x-axis because of very small y-z-offset
532  Double_t rx = 1.0 * centimeter;
533  Double_t ry = 0.0 * centimeter;
534  Double_t rz = 0.0 * centimeter;
535 
536  // take the normal vector as rotation vector
537  Double_t norm_r = TMath::Sqrt(rx*rx + ry*ry + rz*rz);
538  rx = rx/norm_r;
539  ry = ry/norm_r;
540  rz = rz/norm_r;
541 
542  // rotation maxtrix : rotation about r(rx, ry, rz) with angle alpha
543  // before that the coordinate system has to be transformed to x_installation,
544  // y_installation, z_installation as new coordinate system origin
545  (*mFtpcRotationY[i])(1, 1) = rx * rx * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
546  (*mFtpcRotationY[i])(1, 2) = ry * rx * (1 - TMath::Cos(alpha)) - rz * TMath::Sin(alpha);
547  (*mFtpcRotationY[i])(1, 3) = rz * rx * (1 - TMath::Cos(alpha)) + ry * TMath::Sin(alpha);
548  (*mFtpcRotationY[i])(2, 1) = rx * ry * (1 - TMath::Cos(alpha)) + rz * TMath::Sin(alpha);
549  (*mFtpcRotationY[i])(2, 2) = ry * ry * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
550  (*mFtpcRotationY[i])(2, 3) = rz * ry * (1 - TMath::Cos(alpha)) - rx * TMath::Sin(alpha);
551  (*mFtpcRotationY[i])(3, 1) = rx * ry * (1 - TMath::Cos(alpha)) - ry * TMath::Sin(alpha);
552  (*mFtpcRotationY[i])(3, 2) = ry * rz * (1 - TMath::Cos(alpha)) + rx * TMath::Sin(alpha);
553  (*mFtpcRotationY[i])(3, 3) = rz * rz * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
554 
555  (*mFtpcRotationYInverse[i]) = (*mFtpcRotationY[i]).inverse(ierr);
556 
557  if (ierr!=0) {
558  LOG_WARN << "Can't invert FTPC ";
559  if (i == 0) { LOG_WARN << " east rotation matrix Y!" << endm; }
560  else { LOG_WARN << " west rotation matrix Y!" << endm; }
561  LOG_WARN << "FTPC rotation matrix Y:" << (*mFtpcRotationY[i]) << endm;
562  LOG_WARN << "Inverse FTPC rotation matrix Y:" << (*mFtpcRotationYInverse[i]) << endm;
563  }
564 
565  // define rotation axis
566 
567  mObservedVertexOffsetX[0] = 0. * centimeter;
568  mObservedVertexOffsetX[1] = -0.08 * centimeter;
569 
570  rx = 0.0 * centimeter;
571  ry = 1.0 * centimeter;
572  rz = 0.0 * centimeter;
573 
574  norm_r = TMath::Sqrt(rx*rx + ry*ry + rz*rz);
575 
576  rx = rx/norm_r;
577  ry = ry/norm_r;
578  rz = rz/norm_r;
579 
580  // calculate beta analog to alpha
581 
582  pq = TMath::Sqrt(TMath::Power(mObservedVertexOffsetX[i], 2.)
583  - 2.*mObservedVertexOffsetX[i]*mInstallationPointX[i]
584  + TMath::Power(mInstallationPointZ[i], 2.)); // p-q-formula
585 
586  Double_t beta=0;
587 
588  if (i == 0) { // east
589  zShift = (mInstallationPointZ[i] + pq) * centimeter; // take correct solution of p-q-formula
590  phi0 = TMath::ATan((mInstallationPointX[i] - mObservedVertexOffsetX[i]) / mInstallationPointZ[i]) * radian;
591  phi1 = TMath::ATan(mInstallationPointX[i] / (mInstallationPointZ[i] - zShift)) * radian;
592  beta = phi0 - phi1;
593  }
594 
595  else if (i == 1) { //west
596  zShift = (mInstallationPointZ[i] - pq) * centimeter; // take correct solution of p-q-formula
597  phi0 = TMath::ATan((mObservedVertexOffsetX[i] - mInstallationPointX[i]) / mInstallationPointZ[i]) * radian;
598  phi1 = TMath::ATan(mInstallationPointX[i] / (zShift - mInstallationPointZ[i]));
599  beta = phi1 - phi0;
600  }
601 
602  (*mFtpcRotationX[i])(1, 1) = rx * rx * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
603  (*mFtpcRotationX[i])(1, 2) = ry * rx * (1 - TMath::Cos(beta)) - rz * TMath::Sin(beta);
604  (*mFtpcRotationX[i])(1, 3) = rz * rx * (1 - TMath::Cos(beta)) + ry * TMath::Sin(beta);
605  (*mFtpcRotationX[i])(2, 1) = rx * ry * (1 - TMath::Cos(beta)) + rz * TMath::Sin(beta);
606  (*mFtpcRotationX[i])(2, 2) = ry * ry * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
607  (*mFtpcRotationX[i])(2, 3) = rz * ry * (1 - TMath::Cos(beta)) - rx * TMath::Sin(beta);
608  (*mFtpcRotationX[i])(3, 1) = rx * ry * (1 - TMath::Cos(beta)) - ry * TMath::Sin(beta);
609  (*mFtpcRotationX[i])(3, 2) = ry * rz * (1 - TMath::Cos(beta)) + rx * TMath::Sin(beta);
610  (*mFtpcRotationX[i])(3, 3) = rz * rz * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
611 
612  (*mFtpcRotationXInverse[i]) = (*mFtpcRotationX[i]).inverse(ierr);
613 
614  if (ierr!=0) {
615  LOG_WARN << "Can't invert FTPC ";
616  if (i == 0) {LOG_WARN << " east rotation matrix X !" << endm;}
617  else {LOG_WARN << " west rotation matrix X !" << endm;}
618  LOG_WARN << "FTPC rotation matrix X:" << (*mFtpcRotationX[i]) << endm;
619  LOG_WARN << "Inverse FTPC rotation matrix X:" << (*mFtpcRotationXInverse[i]) << endm;
620  }
621 
622  // combine Y and X rotation to rotation matrix
623 
624  (*mFtpcRotation[i]) = (*mFtpcRotationY[i])*(*mFtpcRotationX[i]);
625  (*mFtpcRotationInverse[i])=(*mFtpcRotation[i]).inverse(ierr);
626 
627 
628  if (ierr!=0) {
629  LOG_WARN << "Can't invert FTPC ";
630  if (i == 0) {LOG_WARN << " east rotation matrix!" << endm;}
631  else {LOG_WARN << " west rotation matrix!" << endm;}
632  LOG_WARN << "FTPC rotation matrix:" << (*mFtpcRotation[i]) << endm;
633  LOG_WARN << "Inverse FTPC rotation matrix:" << (*mFtpcRotationInverse[i]) << endm;
634  }
635  }
636 }
637 
638 StFtpcTrackingParams::~StFtpcTrackingParams() {
639  // delete created pointers
640 
641  for (Int_t i = 0; i < 1; i++) {
642  delete mFtpcRotation[i];
643  delete mFtpcRotationInverse[i];
644  delete mFtpcRotationX[i];
645  delete mFtpcRotationXInverse[i];
646  delete mFtpcRotationY[i];
647  delete mFtpcRotationYInverse[i];
648  }
649 
650  delete[] mPadRowPosZ;
651 
652  mInstance = 0;
653 }
654 
655 
656 Int_t StFtpcTrackingParams::InitTrackingParams(ftpcTrackingPars_st *trackParsTable) {
657  // Sets tracking parameters.
658  if (trackParsTable) {
659 
660  // Vertex position
661  mMaxVertexPosZWarning = trackParsTable->maxVertexPosZWarning * centimeter;
662  mMaxVertexPosZError = trackParsTable->maxVertexPosZError * centimeter;
663 
664  // Vertex reconstruction
665  mHistoBins = trackParsTable->histoBins;
666  mHistoMin = trackParsTable->histoMin * centimeter;
667  mHistoMax = trackParsTable->histoMax * centimeter;
668  mMaxDcaVertex = trackParsTable->maxDcaVertex * centimeter;
669  mMinNumTracks = trackParsTable->minNumTracks; // must be >0 !
670 
671 
672  // Tracker
673  mRowSegments = trackParsTable->rowSegments;
674  mPhiSegments = trackParsTable->phiSegments;
675  mEtaSegments = trackParsTable->etaSegments;
676 
677  // Tracking
678  // the 4 indizes represent: 0: main vertex tracking
679  // 1: non vertex tracking
680  // 2: no field tracking
681  // 3: laser tracking
682 
683  for (Int_t i = 0; i < 4; i++) {
684 
685  mLaser[i] = trackParsTable->laser[i] ? kTRUE : kFALSE;
686  mVertexConstraint[i] = trackParsTable->vertexConstraint[i] ? kTRUE : kFALSE;
687 
688  mMaxTrackletLength[i] = trackParsTable->maxTrackletLength[i];
689  mMinTrackLength[i] = trackParsTable->minTrackLength[i];
690  mRowScopeTracklet[i] = trackParsTable->rowScopeTracklet[i];
691 
692  mRowScopeTrack[i] = trackParsTable->rowScopeTrack[i];
693  mPhiScope[i] = trackParsTable->phiScope[i];
694  mEtaScope[i] = trackParsTable->etaScope[i];
695 
696  mMaxDca[i] = trackParsTable->maxDca[i] * centimeter;
697 
698  // Tracklets
699  mMaxAngleTracklet[i] = trackParsTable->maxAngleTracklet[i] * radian;
700 
701  // Tracks
702  mMaxAngleTrack[i] = trackParsTable->maxAngleTrack[i] * radian;
703  mMaxCircleDist[i] = trackParsTable->maxCircleDist[i] * 1./centimeter;
704  mMaxLengthDist[i] = trackParsTable->maxLengthDist[i] * centimeter;
705  }
706 
707  // Split tracks
708  mMaxDist = trackParsTable->maxDist * centimeter;
709  mMinPointRatio = trackParsTable->minPointRatio;
710  mMaxPointRatio = trackParsTable->maxPointRatio;
711 
712  return kStOK;
713  }
714 
715  else {
716  LOG_ERROR << "No data in table class St_TrackingPars." << endm;
717  mReturnCode += 1;
718  return kStErr;
719  }
720 }
721 
722 
723 Int_t StFtpcTrackingParams::InitdEdx(ftpcdEdxPars_st *dEdxParsTable) {
724  // Sets dEdx parameters
725 
726  if (dEdxParsTable) {
727  mDebugLevel = dEdxParsTable[0].debug_level;
728  mIdMethod = dEdxParsTable[0].id_method;
729  mNoAngle = dEdxParsTable[0].no_angle;
730  mMaxHit = dEdxParsTable[0].max_hit;
731  mMinHit = dEdxParsTable[0].min_hit;
732  mMaxTrack = dEdxParsTable[0].max_track;
733  mPadLength = dEdxParsTable[0].pad_length/100.; // from cm in um/keV (microsecond/keV = 100)
734  mFracTrunc = dEdxParsTable[0].frac_trun;
735  mAip = dEdxParsTable[0].a_i_p * 1.0e-9; // in GeV
736  mALargeNumber = dEdxParsTable[0].a_large_number;
737 
738  return kStOK;
739  }
740 
741  else {
742  LOG_ERROR << "No data in table class St_ftpcdEdxPars." << endm;
743  mReturnCode += 2;
744 
745  return kStErr;
746  }
747 }
748 
749 
750 Int_t StFtpcTrackingParams::InitDimensions(ftpcDimensions_st* dimensionsTable) {
751  // Sets FTPC geometry.
752 
753  if (dimensionsTable) {
754  mInnerRadius = dimensionsTable->innerRadiusSensitiveVolume * centimeter;
755  mOuterRadius = dimensionsTable->outerRadiusSensitiveVolume * centimeter;
756  mNumberOfPadRows = dimensionsTable->totalNumberOfPadrows;
757  mNumberOfPadRowsPerSide = dimensionsTable->numberOfPadrowsPerSide;
758  mInstallationPointX[0] = dimensionsTable->installationPointX[0] * centimeter;
759  mInstallationPointX[1] = dimensionsTable->installationPointX[1] * centimeter;
760  mInstallationPointY[0] = dimensionsTable->installationPointY[0] * centimeter;
761  mInstallationPointY[1] = dimensionsTable->installationPointY[1] * centimeter;
762  mInstallationPointZ[0] = dimensionsTable->installationPointZ[0] * centimeter;
763  mInstallationPointZ[1] = dimensionsTable->installationPointZ[1] * centimeter;
764 
765  return kStOK;
766  }
767 
768  else {
769  LOG_ERROR << "No data in table class St_ftpcDimensions." << endm;
770  mReturnCode += 4;
771 
772  return kStErr;
773  }
774 }
775 
776 
777 Int_t StFtpcTrackingParams::InitPadRows(ftpcPadrowZ_st* padrowzTable) {
778  // Set z-positione of pad rows.
779 
780  if (padrowzTable) {
781  mPadRowPosZ = new Double_t[mNumberOfPadRows];
782 
783  for (Int_t i = 0; i < mNumberOfPadRows; i++) {
784  mPadRowPosZ[i] = ((Float_t *)padrowzTable->z)[i];
785  }
786 
787  return kStOK;
788  }
789 
790  else {
791  LOG_ERROR << "No data in table class St_ftpcPadrowZ." << endm;
792  mReturnCode += 8;
793 
794  return kStErr;
795  }
796 }
797 
798 
799 Int_t StFtpcTrackingParams::InitCoordTransformation() {
800  // Fill dummies.
801 
802  mObservedVertexOffsetY[0] = -9999.;
803  mObservedVertexOffsetY[1] = -9999.;
804  mObservedVertexOffsetX[0] = -9999.;
805  mObservedVertexOffsetX[1] = -9999.;
806 
807  return kStOK;
808 }
809 
810 
811 Int_t StFtpcTrackingParams::InitCoordTransformation(ftpcCoordTrans_st* ftpcCoordTrans) {
812  // Set rotation values.
813 
814  if (ftpcCoordTrans) {
815 
816  Double_t oldY[2] = {mObservedVertexOffsetY[0] * centimeter, mObservedVertexOffsetY[1] * centimeter};
817  Double_t oldX[2] = {mObservedVertexOffsetX[0] * centimeter, mObservedVertexOffsetX[1] * centimeter};
818 
819  mObservedVertexOffsetY[0] = ftpcCoordTrans->observedVertexOffsetY[0] * centimeter;
820  mObservedVertexOffsetY[1] = ftpcCoordTrans->observedVertexOffsetY[1] * centimeter;
821 
822  mObservedVertexOffsetX[0] = ftpcCoordTrans->observedVertexOffsetX[0] * centimeter;
823  mObservedVertexOffsetX[1] = ftpcCoordTrans->observedVertexOffsetX[1] * centimeter;
824 
825 // For testing without changing the database, set rotation values here
826 // mObservedVertexOffsetY[0] = 0.0;
827 // mObservedVertexOffsetY[1] = 0.0;
828 // mObservedVertexOffsetX[0] = 0.0;
829 // mObservedVertexOffsetX[1] = 0.0;
830 
831  if ((mObservedVertexOffsetY[0] != oldY[0] || mObservedVertexOffsetY[1] != oldY[1]) &&
832  (oldY[0] < -9990. || oldY[1] < -9990.)) {
833 
834  LOG_DEBUG << "Observed vertex offset in y direction has changed. Changed from "
835  << oldY[0] << " to " << mObservedVertexOffsetY[0] << " (east) and from "
836  << oldY[1] << " to " << mObservedVertexOffsetY[1] << " (west)." << endm;
837  }
838 
839 
840  if ((mObservedVertexOffsetX[0] != oldX[0] || mObservedVertexOffsetX[1] != oldX[1]) &&
841  (oldX[0] < -9990. || oldX[1] < -9990.)) {
842 
843  LOG_DEBUG << "Observed vertex offset in x direction has changed. Changed from "
844  << oldX[0] << " to " << mObservedVertexOffsetX[0] << " (east) and from "
845  << oldX[1] << " to " << mObservedVertexOffsetX[1] << " (west)." << endm;
846  }
847 
848  return kStOK;
849  }
850 
851  else {
852  LOG_ERROR << "No data in table class St_ftpcCoordTrans." << endm;
853  mReturnCode += 16;
854 
855  return kStErr;
856  }
857 }
858 
859 
860 Int_t StFtpcTrackingParams::InitSpaceTransformation() {
861 
862  // transformation due to rotated and displaced TPC
863  // lines commented out go in if TPC database is available
864 
865  if (gStTpcDb->GlobalPosition()) {
866 
867  Double_t phi = 0.0 * radian; //large uncertainty, so set to 0
868  Double_t theta = gStTpcDb->GlobalPosition()->TpcRotationAroundGlobalAxisY() * radian;
869  Double_t psi = gStTpcDb->GlobalPosition()->TpcRotationAroundGlobalAxisX() * radian;
870 
871  mGlobalToTpcRotation(1, 1) = TMath::Cos(theta) * TMath::Cos(phi);
872  mGlobalToTpcRotation(1, 2) = TMath::Cos(theta) * TMath::Sin(phi);
873  mGlobalToTpcRotation(1, 3) = - TMath::Sin(theta);
874  mGlobalToTpcRotation(2, 1) = TMath::Sin(psi) * TMath::Sin(theta) * TMath::Cos(phi) - TMath::Cos(psi) * TMath::Sin(phi);
875  mGlobalToTpcRotation(2, 2) = TMath::Sin(psi) * TMath::Sin(theta) * TMath::Sin(phi) + TMath::Cos(psi) * TMath::Cos(phi);
876  mGlobalToTpcRotation(2, 3) = TMath::Cos(theta) * TMath::Sin(psi);
877  mGlobalToTpcRotation(3, 1) = TMath::Cos(psi) * TMath::Sin(theta) * TMath::Cos(phi) + TMath::Sin(psi) * TMath::Sin(phi);
878  mGlobalToTpcRotation(3, 2) = TMath::Cos(psi) * TMath::Sin(theta) * TMath::Sin(phi) - TMath::Sin(psi) * TMath::Cos(phi);
879  mGlobalToTpcRotation(3, 3) = TMath::Cos(theta) * TMath::Cos(psi);
880 
881  size_t ierr;
882  mTpcToGlobalRotation = mGlobalToTpcRotation.inverse(ierr);
883 
884  if (ierr!=0) {
885  LOG_WARN << "Cant invert rotation matrix!" << endm;
886  LOG_WARN << "Global to TPC rotation matrix:" << mGlobalToTpcRotation << endm;
887  LOG_WARN << "TPC to global rotation matrix:" << mTpcToGlobalRotation << endm;
888  }
889 
890  mTpcPositionInGlobal.setX(gStTpcDb->GlobalPosition()->TpcCenterPositionX() * centimeter);
891  mTpcPositionInGlobal.setY(gStTpcDb->GlobalPosition()->TpcCenterPositionY() * centimeter);
892  mTpcPositionInGlobal.setZ(gStTpcDb->GlobalPosition()->TpcCenterPositionZ() * centimeter);
893  }
894 
895  else {
896  LOG_ERROR << "StTpcDb IS INCOMPLETE! Cannot contstruct Coordinate transformation." << endm;
897  mReturnCode += 32;
898  return kStErr;
899  }
900 
901  // internal FTPC rotation [has do be done before local -> global]
902  for (Int_t i = 0; i <= 1; i++) { // east and west rotation
903 
904  // define rotation angle
905  Double_t zShift, phi0, phi1, alpha;
906  Double_t pq = TMath::Sqrt(TMath::Power(mObservedVertexOffsetY[i], 2.)
907  - 2.*mObservedVertexOffsetY[i]*mInstallationPointY[i]
908  + TMath::Power(mInstallationPointZ[i], 2.)); // p-q-formula
909 
910  // Initialize variable to avoid compiler warning
911  alpha = 0;
912 
913  if (i == 0) { // east
914  zShift = (mInstallationPointZ[i] + pq) * centimeter; // take correct solution of p-q-formula
915  phi0 = TMath::ATan((mInstallationPointY[i] - mObservedVertexOffsetY[i]) / mInstallationPointZ[i]) * radian;
916  phi1 = TMath::ATan(mInstallationPointY[i] / (mInstallationPointZ[i] - zShift)) * radian;
917  alpha = phi0 - phi1;
918  }
919 
920  else if (i == 1) { //west
921  zShift = (mInstallationPointZ[i] - pq) * centimeter; // take correct solution of p-q-formula
922  phi0 = TMath::ATan((mObservedVertexOffsetY[i] - mInstallationPointY[i]) / mInstallationPointZ[i]) * radian;
923  phi1 = TMath::ATan(mInstallationPointY[i] / (zShift - mInstallationPointZ[i]));
924  alpha = phi1 - phi0;
925  }
926 
927  // define rotation axis
928  // simplify to rotation about x-axis because of very small y-z-offset
929  Double_t rx = 1.0 * centimeter;
930  Double_t ry = 0.0 * centimeter;
931  Double_t rz = 0.0 * centimeter;
932 
933  // take the normal vector as rotation vector
934  Double_t norm_r = TMath::Sqrt(rx*rx + ry*ry + rz*rz);
935  rx = rx/norm_r;
936  ry = ry/norm_r;
937  rz = rz/norm_r;
938 
939  // rotation maxtrix : rotation about r(rx, ry, rz) with angle alpha
940  // before that the coordinate system has to be transformed to x_installation,
941  // y_installation, z_installation as new coordinate system origin
942  (*mFtpcRotationY[i])(1, 1) = rx * rx * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
943  (*mFtpcRotationY[i])(1, 2) = ry * rx * (1 - TMath::Cos(alpha)) - rz * TMath::Sin(alpha);
944  (*mFtpcRotationY[i])(1, 3) = rz * rx * (1 - TMath::Cos(alpha)) + ry * TMath::Sin(alpha);
945  (*mFtpcRotationY[i])(2, 1) = rx * ry * (1 - TMath::Cos(alpha)) + rz * TMath::Sin(alpha);
946  (*mFtpcRotationY[i])(2, 2) = ry * ry * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
947  (*mFtpcRotationY[i])(2, 3) = rz * ry * (1 - TMath::Cos(alpha)) - rx * TMath::Sin(alpha);
948  (*mFtpcRotationY[i])(3, 1) = rx * ry * (1 - TMath::Cos(alpha)) - ry * TMath::Sin(alpha);
949  (*mFtpcRotationY[i])(3, 2) = ry * rz * (1 - TMath::Cos(alpha)) + rx * TMath::Sin(alpha);
950  (*mFtpcRotationY[i])(3, 3) = rz * rz * (1 - TMath::Cos(alpha)) + TMath::Cos(alpha);
951 
952  size_t ierr;
953 
954  (*mFtpcRotationYInverse[i]) = (*mFtpcRotationY[i]).inverse(ierr);
955 
956  if (ierr!=0) {
957  LOG_WARN << "Can't invert FTPC ";
958  if (i == 0) {LOG_WARN << " east rotation matrix! Y" << endm;}
959  else {LOG_WARN << " west rotation matrix! Y" << endm;}
960  LOG_WARN << "FTPC rotation matrix Y:" << (*mFtpcRotationY[i]) << endm;
961  LOG_WARN << "Inverse FTPC rotation matrix Y:" << (*mFtpcRotationYInverse[i]) << endm;
962  }
963 
964  // define rotation axis
965 
966  rx = 0.0 * centimeter;
967  ry = 1.0 * centimeter;
968  rz = 0.0 * centimeter;
969 
970  norm_r = TMath::Sqrt(rx*rx + ry*ry + rz*rz);
971 
972  rx = rx/norm_r;
973  ry = ry/norm_r;
974  rz = rz/norm_r;
975 
976  // calculate beta analog to alpha
977 
978  pq = TMath::Sqrt(TMath::Power(mObservedVertexOffsetX[i], 2.)
979  - 2.*mObservedVertexOffsetX[i]*mInstallationPointX[i]
980  + TMath::Power(mInstallationPointZ[i], 2.)); // p-q-formula
981 
982  Double_t beta=0;
983 
984  if (i == 0) { // east
985  zShift = (mInstallationPointZ[i] + pq) * centimeter; // take correct solution of p-q-formula
986  phi0 = TMath::ATan((mInstallationPointX[i] - mObservedVertexOffsetX[i]) / mInstallationPointZ[i]) * radian;
987  phi1 = TMath::ATan(mInstallationPointX[i] / (mInstallationPointZ[i] - zShift)) * radian;
988  beta = phi0 - phi1;
989  }
990 
991  else if (i == 1) { //west
992  zShift = (mInstallationPointZ[i] - pq) * centimeter; // take correct solution of p-q-formula
993  phi0 = TMath::ATan((mObservedVertexOffsetX[i] - mInstallationPointX[i]) / mInstallationPointZ[i]) * radian;
994  phi1 = TMath::ATan(mInstallationPointX[i] / (zShift - mInstallationPointZ[i]));
995  beta = phi1 - phi0;
996  }
997 
998  (*mFtpcRotationX[i])(1, 1) = rx * rx * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
999  (*mFtpcRotationX[i])(1, 2) = ry * rx * (1 - TMath::Cos(beta)) - rz * TMath::Sin(beta);
1000  (*mFtpcRotationX[i])(1, 3) = rz * rx * (1 - TMath::Cos(beta)) + ry * TMath::Sin(beta);
1001  (*mFtpcRotationX[i])(2, 1) = rx * ry * (1 - TMath::Cos(beta)) + rz * TMath::Sin(beta);
1002  (*mFtpcRotationX[i])(2, 2) = ry * ry * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
1003  (*mFtpcRotationX[i])(2, 3) = rz * ry * (1 - TMath::Cos(beta)) - rx * TMath::Sin(beta);
1004  (*mFtpcRotationX[i])(3, 1) = rx * ry * (1 - TMath::Cos(beta)) - ry * TMath::Sin(beta);
1005  (*mFtpcRotationX[i])(3, 2) = ry * rz * (1 - TMath::Cos(beta)) + rx * TMath::Sin(beta);
1006  (*mFtpcRotationX[i])(3, 3) = rz * rz * (1 - TMath::Cos(beta)) + TMath::Cos(beta);
1007 
1008  (*mFtpcRotationXInverse[i]) = (*mFtpcRotationX[i]).inverse(ierr);
1009 
1010  if (ierr!=0) {
1011  LOG_WARN << "Can't invert FTPC ";
1012  if (i == 0) {LOG_WARN << " east rotation matrix X !" << endm;}
1013  else {LOG_WARN << " west rotation matrix X !" << endm;}
1014  LOG_WARN << "FTPC rotation matrix X:" << (*mFtpcRotationX[i]) << endm;
1015  LOG_WARN << "Inverse FTPC rotation matrix X:" << (*mFtpcRotationXInverse[i]) << endm;
1016  }
1017 
1018  // combine Y and X rotation to rotation matrix
1019 
1020  (*mFtpcRotation[i]) = (*mFtpcRotationY[i])*(*mFtpcRotationX[i]);
1021  (*mFtpcRotationInverse[i])=(*mFtpcRotation[i]).inverse(ierr);
1022 
1023 
1024  if (ierr!=0) {
1025  LOG_WARN << "Can't invert FTPC ";
1026  if (i == 0) {LOG_WARN << " east rotation matrix!" << endm;}
1027  else {LOG_WARN << " west rotation matrix!" << endm;}
1028  LOG_WARN << "FTPC rotation matrix:" << (*mFtpcRotation[i]) << endm;
1029  LOG_WARN << "Inverse FTPC rotation matrix:" << (*mFtpcRotationInverse[i]) << endm;
1030  }
1031 
1032  }
1033 
1034  return kStOK;
1035 }
1036 
1037 
1038 void StFtpcTrackingParams::PrintParams() {
1039  // prints params to screen
1040 
1041  LOG_INFO << endm;
1042 
1043  LOG_INFO << "Used parameters for FTPC tracking" << endm;
1044  LOG_INFO << "---------------------------------" << endm;
1045 
1046  LOG_INFO << endm;
1047  LOG_INFO << "FTPC geometry" << endm;
1048  LOG_INFO << "Inner radius (cm)...........: " << mInnerRadius << endm;
1049  LOG_INFO << "Outer radius (cm)...........: " << mOuterRadius << endm;
1050  LOG_INFO << "Total number of padows......: " << mNumberOfPadRows << endm;
1051  LOG_INFO << "Number of padows per FTPC...: " << mNumberOfPadRowsPerSide << endm;
1052 
1053  for (Int_t i = 0; i < NumberOfPadRows(); i++) {
1054  LOG_DEBUG << "z-position of padrow ";
1055  if (i<10) {LOG_DEBUG << "z-position of padrow "<< i << " (cm): " << PadRowPosZ(i) << endm;}
1056  else {LOG_DEBUG << "z-position of padrow "<< i << " (cm): " << PadRowPosZ(i) << endm;}
1057  }
1058 
1059  LOG_INFO << endm;
1060  LOG_INFO << "Vertex position" << endm;
1061  LOG_INFO << "Max. vertex z-position to do tracking w/o warning (cm): " << mMaxVertexPosZWarning << endm;
1062  LOG_INFO << "Max. vertex z-position to do tracking at all (cm).....: " << mMaxVertexPosZError << endm;
1063 
1064  LOG_INFO << endm;
1065  LOG_INFO << "Vertex reconstruction (with FTPC hits)" << endm;
1066  LOG_INFO << "# of histogram bins.............: " << mHistoBins << endm;
1067  LOG_INFO << "lower boundary of histogram (cm): " << mHistoMin << endm;
1068  LOG_INFO << "upper boundary of histogram (cm): " << mHistoMax << endm;
1069  LOG_INFO << "Dca cut (cm)....................: " << mMaxDcaVertex << endm;
1070  LOG_INFO << "min. # of tracks required.......: " << mMinNumTracks << endm;
1071 
1072  LOG_INFO << endm;
1073  LOG_INFO << "Tracker settings" << endm;
1074  LOG_INFO << "# of row segments: " << mRowSegments << endm;
1075  LOG_INFO << "# of phi segments: " << mPhiSegments << endm;
1076  LOG_INFO << "# of eta segments: " << mEtaSegments << endm;
1077 
1078  LOG_INFO << endm;
1079  LOG_INFO << "Settings for tracking" << endm;
1080  LOG_INFO << "Tracking method main vtx non vtx no fld laser" << endm;
1081  LOG_INFO << "Laser tracking switch: " << (Int_t)mLaser[0] << " "
1082  << (Int_t)mLaser[1] << " " << (Int_t)mLaser[2] << " "
1083  << (Int_t)mLaser[3] << endm;
1084  LOG_INFO << "Vertex constraint....: " << (Int_t)mVertexConstraint[0] << " "
1085  << (Int_t)mVertexConstraint[1] << " " << (Int_t)mVertexConstraint[2] << " "
1086  << (Int_t)mVertexConstraint[3] << endm;
1087  LOG_INFO << "Max. tracklet length.: " << mMaxTrackletLength[0] << " "
1088  << mMaxTrackletLength[1] << " " << mMaxTrackletLength[2] << " "
1089  << mMaxTrackletLength[3] << endm;
1090  LOG_INFO << "Min. track length....: " << mMinTrackLength[0] << " "
1091  << mMinTrackLength[1] << " " << mMinTrackLength[2] << " "
1092  << mMinTrackLength[3] << endm;
1093  LOG_INFO << "Tracklet row scope...: " << mRowScopeTracklet[0] << " "
1094  << mRowScopeTracklet[1] << " " << mRowScopeTracklet[2] << " "
1095  << mRowScopeTracklet[3] << endm;
1096  LOG_INFO << "Track row scope......: " << mRowScopeTrack[0] << " "
1097  << mRowScopeTrack[1] << " " << mRowScopeTrack[2] << " "
1098  << mRowScopeTrack[3] << endm;
1099  LOG_INFO << "Phi scope............: " << mPhiScope[0] << " "
1100  << mPhiScope[1] << " " << mPhiScope[2] << " "
1101  << mPhiScope[3] << endm;
1102  LOG_INFO << "Eta scope............: " << mEtaScope[0] << " "
1103  << mEtaScope[1] << " " << mEtaScope[2] << " "
1104  << mEtaScope[3] << endm;
1105  LOG_INFO << "Max. DCA (cm)........: " << mMaxDca[0] << " "
1106  << mMaxDca[1] << " " << mMaxDca[2] << " "
1107  << mMaxDca[3] << endm;
1108 
1109  LOG_INFO << endm;
1110  LOG_INFO << "Cuts for tracking" << endm;
1111  LOG_INFO << "Tracking method main vtx non vtx no fld laser" << endm;
1112  LOG_INFO << "Max. tracklet angle (rad).......: " << mMaxAngleTracklet[0] << " "
1113  << mMaxAngleTracklet[1] << " " << mMaxAngleTracklet[2] << " "
1114  << mMaxAngleTracklet[3] << endm;
1115  LOG_INFO << "Max. track angle (rad)..........: " << mMaxAngleTrack[0] << " "
1116  << mMaxAngleTrack[1] << " " << mMaxAngleTrack[2] << " "
1117  << mMaxAngleTrack[3] << endm;
1118  LOG_INFO << "Max. dist. in circle fit (1/cm).: " << mMaxCircleDist[0] << " "
1119  << mMaxCircleDist[1] << " " << mMaxCircleDist[2] << " "
1120  << mMaxCircleDist[3] << endm;
1121  LOG_INFO << "Max. dist. in length fit (cm)...: " << mMaxLengthDist[0] << " "
1122  << mMaxLengthDist[1] << " " << mMaxLengthDist[2] << " "
1123  << mMaxLengthDist[3] << endm;
1124 
1125  LOG_INFO << endm;
1126  LOG_INFO << "Settings for split track merging" << endm;
1127  LOG_INFO << "Max. distance (cm): " << mMaxDist << endm;
1128  LOG_INFO << "Min. point ratio..: " << mMinPointRatio << endm;
1129  LOG_INFO << "Max. point ratio..: " << mMaxPointRatio << endm;
1130 
1131  LOG_DEBUG << endm;
1132  LOG_DEBUG << "dE/dx" << endm;
1133  LOG_DEBUG << "Debugging level..............: " << mDebugLevel << endm;
1134  LOG_DEBUG << "Method Id....................: " << mIdMethod << endm;
1135  LOG_DEBUG << "Switch for dip/cross angles..: " << mNoAngle << endm;
1136  LOG_DEBUG << "Max. allowable hits per track: " << mMaxHit << endm;
1137  LOG_DEBUG << "Min. no. of hits required....: " << mMinHit << endm;
1138  LOG_DEBUG << "Max. no. of tracks to be used: " << mMaxTrack << endm;
1139  LOG_DEBUG << "Length of pad (us/keV).......: " << mPadLength << endm;
1140  LOG_DEBUG << "Fraction for trunc. mean.....: " << mFracTrunc << endm;
1141  LOG_DEBUG << "a i p (GeV)..................: " << mAip << endm;
1142  LOG_DEBUG << "A large number...............: " << mALargeNumber << endm;
1143 
1144  LOG_INFO << endm;
1145  LOG_INFO << "FTPC to global transformation" << endm;
1146  LOG_INFO << "Installation point x, y, z (east, cm).: " << InstallationPointX(0) << ", "<<InstallationPointY(0) << ", " << InstallationPointZ(0) << endm;
1147  LOG_INFO << "Installation point x, y, z (west, cm).: " << InstallationPointX(1) << ", "<< InstallationPointY(1) << ", " << InstallationPointZ(1) << endm;
1148  LOG_INFO << "Observed vertex offset y (east, cm): " << ObservedVertexOffsetY(0) << endm;
1149  LOG_INFO << "Observed vertex offset y (west, cm): " << ObservedVertexOffsetY(1) << endm;
1150  LOG_INFO << "Observed vertex offset x (east, cm): " << ObservedVertexOffsetX(0) << endm;
1151  LOG_INFO << "Observed vertex offset x (west, cm): " << ObservedVertexOffsetX(1) << endm;
1152  LOG_DEBUG << "FTPC east to global rotation Y: " << FtpcRotationY(0) << endm;
1153  LOG_DEBUG << "Global to FTPC east rotation Y: " << FtpcRotationYInverse(0) << endm;
1154  LOG_DEBUG << "FTPC west to global rotation Y: " << FtpcRotationY(1) << endm;
1155  LOG_DEBUG << "Global to FTPC west rotation Y: " << FtpcRotationYInverse(1) << endm;
1156  LOG_DEBUG << "FTPC east to global rotation X: " << FtpcRotationX(0) << endm;
1157  LOG_DEBUG << "Global to FTPC east rotation X: " << FtpcRotationXInverse(0) << endm;
1158  LOG_DEBUG << "FTPC west to global rotation X: " << FtpcRotationX(1) << endm;
1159  LOG_DEBUG << "Global to FTPC west rotation X: " << FtpcRotationXInverse(1) << endm;
1160  LOG_DEBUG << "FTPC east to global rotation: " << FtpcRotation(0) << endm;
1161  LOG_DEBUG << "Global to FTPC east rotation: " << FtpcRotationInverse(0) << endm;
1162  LOG_DEBUG << "FTPC west to global rotation: " << FtpcRotation(1) << endm;
1163  LOG_DEBUG << "Global to FTPC west rotation: " << FtpcRotationInverse(1) << endm;
1164 
1165  LOG_DEBUG << "TPC to global transformation" << endm;
1166  LOG_INFO << "Position of TPC (cm)..: " << TpcPositionInGlobal() <<endm;
1167  LOG_DEBUG << "TPC to global rotation: " << TpcToGlobalRotation() << endm;
1168  LOG_DEBUG << "Global to TPC rotation: " << GlobalToTpcRotation()<< endm;
1169 
1170  return;
1171 }
StFtpcTrackingParams::TpcToGlobalRotation
StMatrixD TpcToGlobalRotation()
transformation due to rotated and displaced TPC
Definition: StFtpcTrackingParams.cc:226
kStOK
Definition: Stypes.h:40
StFtpcTrackingParams::MaxAngleTracklet
Double_t MaxAngleTracklet(Int_t tracking_method)
Tracklets.
Definition: StFtpcTrackingParams.cc:200
kStErr
Definition: Stypes.h:44
Track
C++ STL includes.
Definition: AgUStep.h:47
Generated by   doxygen 1.8.5