1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Class to make a internal alignemnt of TPC chambers //
21 // Different linear tranformation investigated
22 // 12 parameters - arbitrary transformation
23 // 9 parameters - scaling fixed to 1
26 //// Only the 12 parameter version works so far...
31 #include "TLinearFitter.h"
32 #include "AliTPCcalibAlign.h"
33 #include "AliExternalTrackParam.h"
34 #include "AliTPCTracklet.h"
41 ClassImp(AliTPCcalibAlign)
43 AliTPCcalibAlign::AliTPCcalibAlign()
44 : fDphiHistArray(72*72),
45 fDthetaHistArray(72*72),
48 fFitterArray12(72*72),
55 for (Int_t i=0;i<72*72;++i) {
60 AliTPCcalibAlign::~AliTPCcalibAlign() {
66 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
68 AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
70 // TObjArray trackletsL=
71 // AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
73 // TObjArray trackletsQ=
74 // AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
76 // TObjArray trackletsR=
77 // AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
80 // trackletsL.SetOwner();
81 // trackletsQ.SetOwner();
82 // trackletsR.SetOwner();
83 if (tracklets.GetEntries()==2) {
84 AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
85 AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
86 if (t1->GetSector()>t2->GetSector()) {
87 AliTPCTracklet* tmp=t1;
91 AliExternalTrackParam *common1=0,*common2=0;
92 if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
93 ProcessTracklets(*common1,*common2,t1->GetSector(),t2->GetSector());
100 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
101 const AliExternalTrackParam &tp2,
104 if (s2-s1==36) {//only inner-outer
105 if (!fDphiHistArray[s1*72+s2]) {
108 name<<"hist_phi_"<<s1<<"_"<<s2;
109 title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
110 fDphiHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.01,0.01); // +/- 10 mrad
111 ((TH1D*)fDphiHistArray[s1*72+s2])->SetDirectory(0);
113 if (!fDthetaHistArray[s1*72+s2]) {
116 name<<"hist_theta_"<<s1<<"_"<<s2;
117 title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
118 fDthetaHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.01,0.01); // +/- 10 mrad
119 ((TH1D*)fDthetaHistArray[s1*72+s2])->SetDirectory(0);
121 if (!fDyHistArray[s1*72+s2]) {
124 name<<"hist_y_"<<s1<<"_"<<s2;
125 title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
126 fDyHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.3,0.3); // +/- 3 mm
127 ((TH1D*)fDyHistArray[s1*72+s2])->SetDirectory(0);
129 if (!fDzHistArray[s1*72+s2]) {
132 name<<"hist_z_"<<s1<<"_"<<s2;
133 title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
134 fDzHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.3,0.3); // +/- 3 mm
135 ((TH1D*)fDzHistArray[s1*72+s2])->SetDirectory(0);
137 static_cast<TH1D*>(fDphiHistArray[s1*72+s2])
138 ->Fill(TMath::ASin(tp1.GetSnp())
139 -TMath::ASin(tp2.GetSnp()));
140 static_cast<TH1D*>(fDthetaHistArray[s1*72+s2])
141 ->Fill(TMath::ATan(tp1.GetTgl())
142 -TMath::ATan(tp2.GetTgl()));
143 static_cast<TH1D*>(fDyHistArray[s1*72+s2])
146 static_cast<TH1D*>(fDzHistArray[s1*72+s2])
155 // Process function to fill fitters
157 Double_t t1[5],t2[5];
158 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
159 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3]/*, &dzdx2=t2[4]*/;
163 Double_t snp1=tp1.GetSnp();
164 dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
165 Double_t tgl1=tp1.GetTgl();
166 // dz/dx = 1/(cos(theta)*cos(phi))
167 dzdx1=1./TMath::Sqrt((1.+tgl1*tgl1)*(1.-snp1*snp1));
171 Double_t snp2=tp2.GetSnp();
172 dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
173 Double_t tgl2=tp2.GetTgl();
174 dzdx1=1./TMath::Sqrt((1.+tgl2*tgl2)*(1.-snp2*snp2));
176 Process12(t1,t2,GetOrMakeFitter12(s1,s2));
177 Process9(t1,t2,GetOrMakeFitter12(s1,s2));
178 Process6(t1,t2,GetOrMakeFitter12(s1,s2));
182 void AliTPCcalibAlign::Process12(Double_t *t1,
184 TLinearFitter *fitter) {
185 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
186 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
187 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
188 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
189 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
191 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
192 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
193 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
194 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
195 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
206 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
207 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
208 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
209 for (Int_t i=0; i<12;i++) p[i]=0.;
214 p[0+1] = y1*dydx2; // a01
215 p[0+2] = z1*dydx2; // a02
216 p[9+0] = dydx2; // a03
218 fitter->AddPoint(p,value,sy);
220 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
221 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
222 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
223 for (Int_t i=0; i<12;i++) p[i]=0.;
228 p[0+1] = y1*dzdx2; // a01
229 p[0+2] = z1*dzdx2; // a02
230 p[9+0] = dzdx2; // a03
232 fitter->AddPoint(p,value,sz);
234 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
235 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
236 for (Int_t i=0; i<12;i++) p[i]=0.;
238 p[3+1] = dydx1; // a11
239 p[3+2] = dzdx1; // a12
240 p[0+0] = -dydx2; // a00
241 p[0+1] = -dydx1*dydx2; // a01
242 p[0+2] = -dzdx1*dydx2; // a02
244 fitter->AddPoint(p,value,sdydx);
246 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
247 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
248 for (Int_t i=0; i<12;i++) p[i]=0.;
250 p[6+1] = dydx1; // a21
251 p[6+2] = dzdx1; // a22
252 p[0+0] = -dzdx2; // a00
253 p[0+1] = -dydx1*dzdx2; // a01
254 p[0+2] = -dzdx1*dzdx2; // a02
256 fitter->AddPoint(p,value,sdzdx);
259 void AliTPCcalibAlign::Process9(Double_t *t1,
261 TLinearFitter *fitter) {
262 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
263 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
264 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
265 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
266 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
268 // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
269 // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
270 // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
271 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
272 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
283 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
284 // y2 = a10*x1 + a11*y1 + a12*z1 + a13
285 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a03)*dydx2
286 for (Int_t i=0; i<12;i++) p[i]=0.;
295 fitter->AddPoint(p,value,sy);
297 // x2 = a00*x1 + a01*y1 + a02*z1 + a03
298 // z2 = a20*x1 + a21*y1 + a22*z1 + a23
299 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a03)*dzdx2;
300 for (Int_t i=0; i<12;i++) p[i]=0.;
309 fitter->AddPoint(p,value,sz);
311 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
312 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
313 for (Int_t i=0; i<12;i++) p[i]=0.;
318 p[0+1] = -dydx1*dydx2;
319 p[0+2] = -dzdx1*dydx2;
321 fitter->AddPoint(p,value,sdydx);
323 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
324 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
325 for (Int_t i=0; i<12;i++) p[i]=0.;
330 p[0+1] = -dydx1*dzdx2;
331 p[0+2] = -dzdx1*dzdx2;
333 fitter->AddPoint(p,value,sdzdx);
336 void AliTPCcalibAlign::Process6(Double_t *t1,
338 TLinearFitter *fitter) {
339 // x2 = 1 *x1 +-a01*y1 + 0 +a03
340 // y2 = a01*x1 + 1 *y1 + 0 +a13
341 // z2 = a20*x1 + a21*y1 + 1 *z1 +a23
342 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
343 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
345 // 1 -a01 0 a03 x -p[0] x p[3]
346 // a10 1 0 a13 ==> p[0] x x p[4]
347 // a20 a21 1 a23 p[1] p[2] x p[5]
348 Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
349 Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
360 // x2 = 1 *x1 +-a01*y1 + 0 +a03
361 // y2 = a01*x1 + 1 *y1 + 0 +a13
362 // y2' = a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
363 for (Int_t i=0; i<12;i++) p[i]=0.;
372 fitter->AddPoint(p,value,sy);
374 // x2 = 1 *x1 +-a01*y1 + 0 + a03
375 // z2 = a20*x1 + a21*y1 + 1 *z1 + a23
376 // z2' = a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 +a03)*dzdx2;
377 for (Int_t i=0; i<12;i++) p[i]=0.;
386 fitter->AddPoint(p,value,sz);
388 // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
389 // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
390 for (Int_t i=0; i<12;i++) p[i]=0.;
395 p[0+1] = -dydx1*dydx2;
396 p[0+2] = -dzdx1*dydx2;
398 fitter->AddPoint(p,value,sdydx);
400 // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
401 // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
402 for (Int_t i=0; i<12;i++) p[i]=0.;
407 p[0+1] = -dydx1*dzdx2;
408 p[0+2] = -dzdx1*dzdx2;
410 fitter->AddPoint(p,value,sdzdx);
413 void AliTPCcalibAlign::Eval() {
415 for (Int_t s1=0;s1<72;++s1)
416 for (Int_t s2=0;s2<72;++s2)
417 if ((f=GetFitter12(s1,s2))&&fPoints[72*s1+s2]>12) {
418 // cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
420 cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
427 chi212 = align->GetChisquare()/(4.*entries);
431 align->GetParameters(par);
432 align->GetCovarianceMatrix(mat);
436 for (Int_t i=0; i<12;i++){
437 palign12(i)= par(i+1);
438 for (Int_t j=0; j<12;j++){
439 pcovar12(i,j) = mat(i+1,j+1);
440 pcovar12(i,j) *= chi212;
444 for (Int_t i=0; i<12;i++){
445 psigma12(i) = TMath::Sqrt(pcovar12(i,i));
446 palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
447 for (Int_t j=0; j<12;j++){
448 pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
454 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
455 if (!GetFitter12(s1,s2))
460 GetFitter12(s1,s2)->GetParameters(p);
483 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
484 if (!GetFitter9(s1,s2))
488 GetFitter9(s1,s2)->GetParameters(p);
510 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
511 if (!GetFitter6(s1,s2))
516 GetFitter6(s1,s2)->GetParameters(p);