libTPCbase.pkg TPCbaseLinkDef.h -
[u/mrichter/AliRoot.git] / TPC / AliTPCTransformation.cxx
CommitLineData
94685752 1/*
2
3 Class AliTPCtransformation:
4 Should represent general non linear transformation. Currently tune for TPConly.
5 To be used:
6 1. Simulation-Digitization
7 2. Reconstruction - AliTPCTransform
8 3. Calibration/Alignment (KalmanFilter, Milipedde)
9 4. Set of transformation to be stored/retrieved as OCDB entry
10
11 Base functionality:
12
13 1. Double_t GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z)
14 Get correction - return the delta of coordinate coord dx or dy or dz for given volID for point at point (x,y,z)
15 All coordinates are global
16
17 2. The transformation should work only for given volIDs and detector IDs
18 Currently Bitmask is used for filtering
19
20
21
22*/
23
24#include <string.h>
25#include "TRandom.h"
26#include "TMath.h"
27#include "TBits.h"
28#include "TFormula.h"
29#include "TF1.h"
30#include "TLinearFitter.h"
31#include "TFile.h"
32#include "TObjString.h"
33
34#include "TTreeStream.h"
35#include "AliTrackPointArray.h"
36#include "AliLog.h"
37#include "AliTPCTransformation.h"
38
39ClassImp(AliTPCTransformation)
40
41
42AliTPCTransformation::GenFuncG AliTPCTransformation::fgFormulas[10000];
43TObjArray* AliTPCTransformation::fgFormulasName = new TObjArray(10000);
44
45
46void AliTPCTransformation::RegisterFormula(const char * name, GenFuncG formula){
47 //
48 // Add Formula to the list of formulas
49 //
50 Int_t last= fgFormulasName->GetEntries();
51 fgFormulasName->AddLast(new TObjString(name));
52 fgFormulas[last]=formula;
53}
54
55Int_t AliTPCTransformation::BuildBasicFormulas(){
56 //
6d438146 57
58 //
59 //build list of basic TPC formulas - corrections
94685752 60 //
61 RegisterFormula("TPCscalingRPol",(GenFuncG)(AliTPCTransformation::TPCscalingRPol));
6d438146 62 RegisterFormula("TPCscalingRIFC",(GenFuncG)(AliTPCTransformation::TPCscalingRIFC));
63 RegisterFormula("TPCscalingROFC",(GenFuncG)(AliTPCTransformation::TPCscalingROFC));
64 //
94685752 65 RegisterFormula("TPCscalingZDr",(GenFuncG)(AliTPCTransformation::TPCscalingZDr));
66 RegisterFormula("TPCscalingPhiLocal",(GenFuncG)(AliTPCTransformation::TPCscalingPhiLocal));
0e9efcbe 67 //
68 // TPC Local X and Y misalignment
69 //
70 RegisterFormula("TPClocaldLxdGX",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGX));
71 RegisterFormula("TPClocaldLxdGY",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGY));
72 RegisterFormula("TPClocaldLydGX",(GenFuncG)(AliTPCTransformation::TPClocaldLydGX));
73 RegisterFormula("TPClocaldLydGY",(GenFuncG)(AliTPCTransformation::TPClocaldLydGY));
74 //
75 // Z offset
76 //
77 RegisterFormula("TPCDeltaZ",(GenFuncG)(AliTPCTransformation::TPCDeltaZ));
78 RegisterFormula("TPCDeltaZMediumLong",(GenFuncG)(AliTPCTransformation::TPCDeltaZMediumLong));
79 RegisterFormula("TPCTiltingZ",(GenFuncG)(AliTPCTransformation::TPCTiltingZ));
94685752 80 return 0;
81}
82
83AliTPCTransformation::GenFuncG AliTPCTransformation::FindFormula(const char * name){
84 //
85 // find formula - if registered
86 //
87 if (fgFormulasName->FindObject(name)==0) return 0;
88 Int_t entries = fgFormulasName->GetEntries();
89 for (Int_t i=0;i<entries;i++){
90 if (strcmp(fgFormulasName->At(i)->GetName(), name)==0){
91 return fgFormulas[i];
92 }
93 }
94 return 0;
95}
96
97Double_t AliTPCTransformation::Eval(const char * name, const Double_t*x,const Double_t*par){
98 //
99 // Only for test purposes - very slow
100 //
101 GenFuncG fun = FindFormula(name);
102 if (!fun) return 0;
103 return fun(x,par);
104}
105
106
107AliTPCTransformation::AliTPCTransformation():
108 TNamed(),
109 fNameX(0), // x formula name
110 fNameY(0), // y formula name
111 fNameZ(0), // z formula name
112 //
113 fBitMask(0), // bitmaps - transformation only for specified volID
114 fCoordSystem(0), // coord system of output deltas
115 fParam(0), // free parameter of transformation
116 fSigma(0), // uncertainty of the parameter
6d438146 117 fSigma2Time(0), // change of the error in time - (For kalman filter)
94685752 118 fFixedParam(0), // fixed parameters of tranformation
0e9efcbe 119 fIsActive(kTRUE), // swith - On/Off
94685752 120 //
121 fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
122 fFormulaX(0), // x formula - pointer to the function
123 fFormulaY(0), // y formula - pointer to the function
124 fFormulaZ(0) // z formula - pointer to the function
125 //
126{
127 //
128 // default constructor
129 //
130}
131
132
133
6d438146 134AliTPCTransformation::AliTPCTransformation(const char *name, TBits *mask, const char *fx, const char *fy, const char *fz, Int_t coordSystem):
94685752 135 TNamed(name,name),
136 fNameX(0), // x formula name
137 fNameY(0), // y formula name
138 fNameZ(0), // z formula name
139 fBitMask(mask), // bitmaps - transformation only for specified volID
140 fCoordSystem(coordSystem), // coordinate system of output deltas
6d438146 141 fParam(0), // free parameter of transformation
142 fSigma(0),
143 fSigma2Time(0), // change of sigma in time
94685752 144 fFixedParam(0), // fixed parameters of tranformation
0e9efcbe 145 fIsActive(kTRUE), // swith - On/Off
94685752 146 //
147 fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
148 fFormulaX(0), // x formula - pointer to the function
149 fFormulaY(0), // y formula - pointer to the function
150 fFormulaZ(0) // z formula - pointer to the function
151{
152 //
153 // non default constructor
154 //
155 if (fx) fNameX= new TString(fx);
156 if (fy) fNameY= new TString(fy);
157 if (fz) fNameZ= new TString(fz);
94685752 158 Init();
159}
160
6d438146 161AliTPCTransformation::AliTPCTransformation(const AliTPCTransformation&trafo):
162 TNamed(trafo),
163 fNameX(0), // x formula name
164 fNameY(0), // y formula name
165 fNameZ(0), // z formula name
166 //
167 fBitMask(0), // bitmaps - transformation only for specified volID
168 fCoordSystem(0), // coord system of output deltas
169 fParam(trafo.fParam), // free parameter of transformation
170 fSigma(trafo.fSigma), // uncertainty of the parameter
171 fSigma2Time(trafo.fSigma2Time), // change of the error in time - (For kalman filter)
172 fFixedParam(0), // fixed parameters of tranformation
173 //
174 fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
175 fFormulaX(0), // x formula - pointer to the function
176 fFormulaY(0), // y formula - pointer to the function
177 fFormulaZ(0) // z formula - pointer to the function
178{
179 if (trafo.fNameX) fNameX = new TString(*(trafo.fNameX));
180 if (trafo.fNameY) fNameY = new TString(*(trafo.fNameY));
181 if (trafo.fNameZ) fNameZ = new TString(*(trafo.fNameZ));
182 if (trafo.fBitMask) fBitMask = new TBits(*(trafo.fBitMask));
183}
184
185AliTPCTransformation::~AliTPCTransformation(){
186 //
187 // destructor
188 //
189 delete fNameX;
190 delete fNameY;
191 delete fNameZ;
192 delete fBitMask;
193 delete fFixedParam;
194}
195void AliTPCTransformation::SetParams(Double_t param, Double_t sigma, Double_t sigma2Time, TVectorD* fixedParams){
196 //
197 // Set parameters of transformation
198 //
199 fParam = param;
200 fSigma = sigma;
201 fSigma2Time = sigma2Time;
202 if (fFixedParam) delete fFixedParam;
203 fFixedParam = new TVectorD(*fixedParams);
204 Init();
205}
94685752 206
207
208Bool_t AliTPCTransformation::Init(){
209 //
210 // associate formulas with pointer to the function
211 //
212 Bool_t isOK=kTRUE;
213 if (fNameX) {
214 fFormulaX=FindFormula(fNameX->Data());
215 if (fFormulaX==0) isOK=kFALSE;
216 }
217 if (fNameY) {
218 fFormulaY=FindFormula(fNameY->Data());
219 if (fFormulaY==0) isOK=kFALSE;
220 }
221 if (fNameZ) {
222 fFormulaZ=FindFormula(fNameZ->Data());
223 if (!fFormulaZ) isOK=kFALSE;
224 }
225 return isOK;
226}
227
0e9efcbe 228
229
94685752 230TBits * AliTPCTransformation::BitsSide(Bool_t aside){
231 //
232 TBits * bits = new TBits(72);
233 for (Int_t i=0; i<72;i++){
234 if (i%36<18 && aside) (*bits)[i]=kTRUE;
235 if (i%36<18 && (!aside)) (*bits)[i]=kFALSE;
236 if (i%36>=18 && aside) (*bits)[i]=kFALSE;
237 if (i%36>=18 && (!aside)) (*bits)[i]=kTRUE;
238 }
239 return bits;
240}
241
242TBits * AliTPCTransformation::BitsAll(){
243 //
244 //
245 //
246 TBits * bits = new TBits(72);
247 for (Int_t i=0; i<72;i++){
248 (*bits)[i]=kTRUE;
249 }
250 return bits;
251}
252
253Double_t AliTPCTransformation::GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z){
254 //
255 //
256 //
0e9efcbe 257 if (!fIsActive) return 0;
94685752 258 if (fBitMask && (!(*fBitMask)[volID])) return 0;
0e9efcbe 259 Double_t xyz[5]={x,y,z, param,volID};
94685752 260 if (fCoordSystem==0){
261 // cartezian system
262 if (coord==0 && fFormulaX) return fFormulaX(xyz,fFixedParam->GetMatrixArray());
263 if (coord==1 && fFormulaY) return fFormulaY(xyz,fFixedParam->GetMatrixArray());
0e9efcbe 264 if (coord==2 && fFormulaZ) return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
94685752 265 }
266 if (fCoordSystem==1){
267 // cylindrical system
268 if (coord==2) {
269 if (fFormulaZ==0) return 0;
270 return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
271 }
272 Double_t rrphiz[3]={0,0,0};
273 if (fFormulaX) rrphiz[0] = fFormulaX(xyz,fFixedParam->GetMatrixArray());
274 if (fFormulaY) rrphiz[1] = fFormulaY(xyz,fFixedParam->GetMatrixArray());
275 Double_t alpha = TMath::ATan2(y,x);
276 Double_t ca = TMath::Cos(alpha);
277 Double_t sa = TMath::Sin(alpha);
278 if (coord==0) return ca*rrphiz[0]-sa*rrphiz[1];
279 if (coord==1) return sa*rrphiz[0]+ca*rrphiz[1];
280 }
281 return 0;
282}
283
284Double_t AliTPCTransformation::TPCscalingRPol(Double_t *xyz, Double_t * param){
285 //
286 // Scaling and shift of TPC radius
287 // xyz[0..2] - global xyz of point
288 // xyz[3] - scale parameter
6d438146 289 // param[0] - radial scaling power
290 // param[1] - drift scaling power
291 // radius from -1(at rInner) to 1 (rOuter)
292 // driftM from -1(at 0 drift) to 1 (250 cm drift)
293
294 Double_t rInner=78.8;
295 Double_t rOuter=258.0;
296 Double_t deltaR = rOuter-rInner;
297 Double_t radius = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])-rInner)*2./deltaR;
298 Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.0;
299 Double_t delta = TMath::Power(radius,param[0])*TMath::Power(driftM,param[1]);
300 return delta*xyz[3];
94685752 301}
302
303
304Double_t AliTPCTransformation::TPCscalingZDr(Double_t *xyz, Double_t * param){
305 //
306 //
307 // Scaling and shift of TPC radius
308 // xyz[0..2] - global xyz of point
309 // xyz[3] - scale parameter
310 Double_t driftP = TMath::Power(1. - TMath::Abs(xyz[2]/250.), param[0]);
311 Double_t deltaZ = (xyz[2]>0) ? -driftP : driftP;
312 return deltaZ*xyz[3];
313}
314
315Double_t AliTPCTransformation::TPCscalingPhiLocal(Double_t *xyz, Double_t * param){
316 //
317 //
318 // Scaling if the local y -phi
319 // xyz[0..2] - global xyz of point
320 // xyz[3] - scale parameter
6d438146 321 // value = 1 for ful drift length and parameter 1
94685752 322 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
323 Double_t sector = TMath::Nint(9*alpha/TMath::Pi()-0.5);
324 Double_t localAlpha = (alpha-(sector+0.5)*TMath::Pi()/9.);
325 Double_t radius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])/250.;
6d438146 326 //
327 Double_t deltaAlpha = radius*TMath::Power(2.*9.*localAlpha/TMath::Pi(),param[0]);
94685752 328 return deltaAlpha*xyz[3];
329}
330
331
6d438146 332Double_t AliTPCTransformation::TPCscalingRIFC(Double_t *xyz, Double_t * param){
333 //
334 // inner field cage r distorion - proportinal to 1 over distance to the IFC
335 // param[0] - drift polynom order
336 // distortion at first pad row - is normalized to
337 Double_t rInner=78.8;
338 Double_t rFirst=85.2;
339 Double_t deltaR = rFirst-rInner;
340 Double_t ndistR = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])-rInner)/deltaR;
341 Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
342 Double_t value = TMath::Power(driftM,param[0])/ndistR;
343 return xyz[3]*value;
344}
345
346Double_t AliTPCTransformation::TPCscalingROFC(Double_t *xyz, Double_t * param){
347 //
348 // outer field cage r distorion - proportinal to 1 over distance to the OFC
349 // param[0] - drift polynom order
350 // driftM - from -1 to 1
351 //
352 Double_t rLast=245.8;
353 Double_t rOuter=258.0;
354 Double_t deltaR = rOuter-rLast;
355 Double_t ndistR = (rOuter-TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]))/deltaR;
356 Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
357 Double_t value = TMath::Power(driftM,param[0])/ndistR;
358 return xyz[3]*value;
359}
0e9efcbe 360
361
362//
363// TPC sector local misalignment
364//
365Double_t AliTPCTransformation:: TPClocaldLxdGX(Double_t *xyz, Double_t * /*param*/){
366 //
367 // xyz - [0..2] - position
368 // [3] - scale parameter
369 // [4] - volID
370 // return delta in global coordiante system
371 //
372 Int_t sector = TMath::Nint(xyz[4]);
373 Double_t alpha = TMath::Pi()*(sector+0.5)/9;
374 Double_t ca = TMath::Cos(alpha);
375 // Double_t sa = TMath::Sin(alpha);
376 return ca*xyz[3];
377}
378
379Double_t AliTPCTransformation::TPClocaldLxdGY(Double_t *xyz, Double_t * /*param*/){
380 //
381 // xyz - [0..2] - position
382 // [3] - scale parameter
383 // [4] - volID
384 // return delta in global coordiante system
385 //
386 Int_t sector = TMath::Nint(xyz[4]);
387 Double_t alpha = TMath::Pi()*(sector+0.5)/9;
388 //Double_t ca = TMath::Cos(alpha);
389 Double_t sa = TMath::Sin(alpha);
390 return sa*xyz[3];
391}
392
393Double_t AliTPCTransformation:: TPClocaldLydGX(Double_t *xyz, Double_t * /*param*/){
394 //
395 // xyz - [0..2] - position
396 // [3] - scale parameter
397 // [4] - volID
398 // return delta in global coordiante system
399 //
400 Int_t sector = TMath::Nint(xyz[4]);
401 Double_t alpha = TMath::Pi()*(sector+0.5)/9;
402 //Double_t ca = TMath::Cos(alpha);
403 Double_t sa = TMath::Sin(alpha);
404 return -sa*xyz[3];
405}
406
407Double_t AliTPCTransformation::TPClocaldLydGY(Double_t *xyz, Double_t * /*param*/){
408 //
409 // xyz - [0..2] - position
410 // [3] - scale parameter
411 // [4] - volID
412 // return delta in global coordiante system
413 //
414 Int_t sector = TMath::Nint(xyz[4]);
415 Double_t alpha = TMath::Pi()*(sector+0.5)/9;
416 Double_t ca = TMath::Cos(alpha);
417 //Double_t sa = TMath::Sin(alpha);
418 return ca*xyz[3];
419}
420
421
422Double_t AliTPCTransformation::TPCDeltaZ(Double_t *xyz, Double_t * /*param*/){
423 //
424 // xyz - [0..2] - position
425 // [3] - scale parameter
426 // [4] - volID
427 // return delta in global coordiante system
428 //
429 Int_t sector = TMath::Nint(xyz[4]);
430 Double_t signZ = (sector%36<18) ? 1: -1; // drift direction
431 return signZ*xyz[3]; // IROC shift
432}
433
434Double_t AliTPCTransformation::TPCDeltaZMediumLong(Double_t *xyz, Double_t * /*param*/){
435 //
436 // xyz - [0..2] - position
437 // [3] - scale parameter
438 // [4] - volID
439 // return delta in global coordinate system
440 //
441 Int_t sector = TMath::Nint(xyz[4]);
442 Double_t signZ = (sector%36<18) ? 1: -1; // drift direction
443 if (sector<36) return 0;
444 //
445 Double_t radius = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]));
446 const Double_t radiusLong = 198.1;
447 Double_t sign = (radius<radiusLong) ? 1:-1;
448 return xyz[3]*sign*signZ;
449}
450
451Double_t AliTPCTransformation::TPCTiltingZ(Double_t *xyz, Double_t * param){
452 // xyz - [0..2] - position
453 // [3] - scale parameter
454 // [4] - volID
455 // param[0] - n for cos
456 // param[1] - n for sin
457 // return delta in global coordinate system
458 Double_t radius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
459 const Double_t rFirst=85.2;
460 const Double_t rLast =245.8;
461 Double_t radiusC = (rFirst+rLast)*0.5;
462 Double_t deltaR = 2.0*(radius-radiusC)/(rLast-rFirst);
463 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
464
465 if (param[0]>0) deltaR *= TMath::Cos(param[0]*alpha);
466 if (param[1]>0) deltaR *= TMath::Sin(param[1]*alpha);
467 return deltaR*xyz[3];
468}
469