]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSAlignmentModule.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update info about cvs installation using cvs account
[u/mrichter/AliRoot.git] / ITS / AliITSAlignmentModule.cxx
CommitLineData
e8189707 1/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
3/* $Id$ */
4/* $Author$ */
5/* $Date$ */
6/* $Name$ */
7/* $Header$ */
8/*
9 $Log$
10 Revision 1.1.2.3 2000/06/04 16:35:09 nilsen
11 One more try to fix log comments.
12
13 Revision 1.1.2.2 2000/03/04 23:39:36 nilsen
14 Fixed the logs???
15 */
16/* Revision 1.1.2.1 2000/03/02 20:12:23 nilsen */
17/* A new class usefull for ITS detector Alignment Studdies */
18/* */
19/* $Revision$ */
20
21// Standard C & C++ libraries
22#include <TObject.h>
23
24// Standard Root Libraries
25#include <TParticle.h>
26
27// ITS libraries
28#include "AliITSgeom.h"
29#include "AliITSAlignmentTrack.h"
30#include "AliITSAlignmentModule.h"
31
32ClassImp(AliITSAlignmentModule)
33
34//______________________________________________________________________
35AliITSAlignmentModule::AliITSAlignmentModule(){
36
37 findex = -1;
38 flay = -1;
39 flad = -1;
40 fdet = -1;
41 fx0[0] = 0.0;
42 fx0[1] = 0.0;
43 fx0[2] = 0.0;
44 fangles[0] = 0.0;
45 fangles[1] = 0.0;
46 fangles[2] = 0.0;
47 fM[0][0] = 0.0;
48 fM[0][1] = 0.0;
49 fM[0][2] = 0.0;
50 fM[1][0] = 0.0;
51 fM[1][1] = 0.0;
52 fM[1][2] = 0.0;
53 fM[2][0] = 0.0;
54 fM[2][1] = 0.0;
55 fM[2][2] = 0.0;
56 ftrksM =0;
57 fChi2 = 0.0;
58}
59//______________________________________________________________________
60AliITSAlignmentModule::AliITSAlignmentModule(Int_t index,AliITSgeom *gm,
61 Int_t ntrk,AliITSAlignmentTrack *trk){
62 Float_t x,y,z,n;
63 Int_t i,j;
64
65 findex = index;
66 gm->GetModuleId(index,flay,flad,fdet);
67 gm->GetRotMatrix(index,(Double_t *) &(fM[0][0]));
68 gm->GetTrans(flay,flad,fdet,x,y,z);
69 fx0[0] = (Double_t )x;
70 fx0[1] = (Double_t )y;
71 fx0[2] = (Double_t )z;
72 gm->GetAngles(flay,flad,fdet,x,y,z);
73 fangles[0] = (Double_t )x;
74 fangles[1] = (Double_t )y;
75 fangles[2] = (Double_t )z;
76 ftrksM = new TObjArray();
77 fChi2 = 0.0;
78
79 for(i=0;i<ntrk;i++){
80 n = 0;
81 for(j=0;j<trk[i].GetNumberOfClustersSl();j++){
82 if(trk[i].GetIndex(j,index)==findex){
83 ftrksM->AddAtFree((TObject *) &trk[i]);
84 break; // break out of the j loop
85 } // end if
86 } // end for j
87 } // end for i
88}
89//______________________________________________________________________
90Double_t AliITSAlignmentModule::ComputeChi2(){
91 Float_t n;
92 Int_t i,j,ntrk,ntr;
93
94 ntrk = ftrksM->GetEntriesFast();
95 fChi2 = 0.0;
96 for(i=0;i<ntrk;i++){
97 n = 0;
98 ntr =((AliITSAlignmentTrack *)(ftrksM->At(i)))->GetNumberOfClustersSl();
99 for(j=0;j<ntr;j++){
100 fChi2 += (Double_t) ((AliITSAlignmentTrack *)(ftrksM->At(i)))->
101 GetChi2();
102 n++;
103 } // end for j
104 if(n==0){
105 fChi2 = -1.0;
106 }else{
107 fChi2 /= (Double_t) n;
108 } // end if n==0
109 } // end for i
110 return fChi2;
111}
112//______________________________________________________________________
113void AliITSAlignmentModule::lnsrch(Int_t npar,Double_t *xold,Double_t fold,
114 Double_t *g,Double_t *p,Double_t *x,
115 Double_t &f, Double_t stpmax,Int_t &check){
116 Double_t ALF = 1.0e-4, TOLX = 1.0E-7;
117
118 Int_t i;
119 Double_t a,alam,alam2=0.0,alamin,b,disc,f2=0.0,rhs1,rhs2,slope,sum,temp,
120 test,tmplam;
121
122 check = 0;
123 for(sum=0.0,i=0;i<npar;i++) sum += p[i]*p[i];
124 sum = TMath::Sqrt(sum);
125 if(sum>stpmax) for(i=0;i<npar;i++) p[i] *= stpmax/sum;
126 for(slope=0.0,i=0;i<npar;i++) slope += g[i]*p[i];
127 if(slope >=0.0) printf("Error: round off problem in lnsrch.\n");
128 test = 0.0;
129 for(i=0;i<npar;i++){
130 temp = TMath::Abs(p[i])/TMath::Max(TMath::Abs(xold[i]),1.0);
131 if(temp > test) test = temp;
132 } // end for i
133 alamin = TOLX/test;
134 alam = 1.0;
135 for(;;){
136 for(i=0;i<npar;i++) x[i] = xold[i] + alam*p[i];
137 f = Chi2(x);
138 if(alam < alamin){
139 for(i=0;i<npar;i++) x[i] = xold[i];
140 check = 1;
141 return;
142 }else if(f <= fold+ALF*alam*slope) return;
143 else{
144 if(alam == 1.0) tmplam = -slope/(2.0*(f-fold-slope));
145 else{
146 rhs1 = f-fold-alam*slope;
147 rhs2 = f2-fold-alam2*slope;
148 a = (rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
149 b = (-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam*alam);
150 if(a==0.0) tmplam = -slope/(2.0*b);
151 else{
152 disc = b*b - 3.0*a*slope;
153 if(disc<0.0) tmplam = 0.5*alam;
154 else if(b<=0.0) tmplam = (-b+TMath::Sqrt(disc))/(3.0*a);
155 else tmplam = -slope/(b+TMath::Sqrt(disc));
156 } // end if a == 0.0
157 if(tmplam > 0.5*alam) tmplam = 0.5*alam;
158 } // end if alam == 1.0
159 } // end if alam < alamin
160 alam2 = alam;
161 f2 = f;
162 alam = TMath::Max(tmplam,0.1*alam);
163 } // end for ever loop
164
165}
166//______________________________________________________________________
167void AliITSAlignmentModule::MRVMminimization(Int_t npar,Double_t *p,
168 Double_t &fret,Double_t gtol,
169 Int_t &iter){
170 Int_t ITMAX = 200;
171 Double_t EPS = 3.0e-8, TOLX = 4.0*EPS, STPMX = 100.0;
172
173 Int_t check,i,its,j;
174 Double_t den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
175 Double_t *dg,*g,*hdg,**hessin,*pnew,*xi;
176
177 // allocate space.
178 dg = new Double_t[npar];
179 g = new Double_t[npar];
180 hdg = new Double_t[npar];
181 hessin = new Double_t * [npar];
182 for(i=0;i<npar;i++) hessin[i] = new Double_t[npar];
183 pnew = new Double_t[npar];
184 xi = new Double_t[npar];
185
186 // init function values
187 fp = Chi2(p);
188 dChi2(p,g);
189
190 for(i=0;i<npar;i++){
191 for(j=0;j<npar;j++) hessin[i][j] = 0.0;
192 hessin[i][i] = 1.0;
193 xi[i] = -g[i];
194 sum += p[i]*p[i];
195 } // end for i
196 stpmax = STPMX*TMath::Max(TMath::Sqrt(sum),(Double_t)(npar+1));
197
198 for(its=0;its<ITMAX;its++){
199 iter = its;
200 lnsrch(npar,p,fp,g,xi,pnew,fret,stpmax,check);
201 fp = fret;
202 for(i=0;i<npar;i++){
203 xi[i] = pnew[i] - p[i];
204 p[i] = pnew[i];
205 } // end for i
206 test = 0.0;
207 for(i=0;i<npar;i++){
208 temp = TMath::Abs(xi[i])/TMath::Max(TMath::Abs(p[i]),1.0);
209 if(temp > test) test = temp;
210 } // end for i
211 if(test<TOLX) break;
212 for(i=0;i<npar;i++) dg[i] = g[i];
213 dChi2(p,g);
214 test = 0.0;
215 den = TMath::Max(fret,1.0);
216 for(i=0;i<npar;i++){
217 temp = TMath::Abs(g[i])*TMath::Max(TMath::Abs(p[i]),1.0)/den;
218 if(temp>test) test=temp;
219 } // end for i
220 if(test<gtol) break;
221 for(i=0;i<npar;i++) dg[i] = g[i] - dg[i];
222 for(i=0;i<npar;i++){
223 hdg[i] = 0.0;
224 for(j=0;j<npar;j++) hdg[i] += hessin[i][j]*dg[j];
225 } // end for i
226 fac = fae = sumdg = sumxi = 0.0;
227 for(i=0;i<npar;i++){
228 fac += dg[i]*xi[i];
229 fae += dg[i]*hdg[i];
230 sumdg += TMath::Sqrt(dg[i]);
231 sumxi += TMath::Sqrt(xi[i]);
232 } // end for i
233 if(fac>TMath::Sqrt(EPS*sumdg*sumxi)){
234 fac = 1.0/fac;
235 fad = 1.0/fae;
236 for(i=0;i<npar;i++) dg[i] = fac*xi[i]-fad*hdg[i];
237 for(i=0;i<npar;i++) for(j=i;j<npar;j++){
238 hessin[i][j] += fac*xi[i]*xi[j] - fad*hdg[i]*hdg[j] +
239 fae*dg[i]*dg[j];
240 hessin[j][i] = hessin[i][j];
241 } // end for i,j
242 } // end if fac>...
243 for(i=0;i<npar;i++){
244 xi[i] = 0.0;
245 for(j=0;j<npar;j++) xi[i] -= hessin[i][j]*g[j];
246 } // end for i
247 } // end for its
248 if(its==ITMAX) printf("Error: too many iterations\n");
249 delete[] dg;
250 delete[] g;
251 delete[] hdg;
252 delete[] pnew;
253 delete[] xi;
254 for(i=0;i<npar;i++) delete[] hessin[i];
255 delete[] hessin;
256}
257//______________________________________________________________________
258void AliITSAlignmentModule::SetByAngles(Double_t *th){
259 Double_t sx,cx,sy,cy,sz,cz;
260
261 sx = TMath::Sin(th[0]); cx = TMath::Cos(th[0]);
262 sy = TMath::Sin(th[1]); cy = TMath::Cos(th[1]);
263 sz = TMath::Sin(th[2]); cz = TMath::Cos(th[2]);
264 for(Int_t i=0;i<3;i++) fangles[i] = th[i];
265 fM[0][0] = cz*cy;
266 fM[0][1] = -cz*sy*sx - sz*cx;
267 fM[0][2] = -cz*sy*cx + sz*sx;
268 fM[1][0] = sz*cy;
269 fM[1][1] = -sz*sy*sx + cz*cx;
270 fM[1][2] = -sz*sy*cx - cz*sx;
271 fM[2][0] = sy;
272 fM[2][1] = cy*sx;
273 fM[2][2] = cy*cx;
274}
275//______________________________________________________________________
276void AliITSAlignmentModule::dfMdthx(Double_t dfMx[3][3]){
277 Double_t sx,cx,sy,cy,sz,cz;
278
279 sx = TMath::Sin(fangles[0]); cx = TMath::Cos(fangles[0]);
280 sy = TMath::Sin(fangles[1]); cy = TMath::Cos(fangles[1]);
281 sz = TMath::Sin(fangles[2]); cz = TMath::Cos(fangles[2]);
282 dfMx[0][0] = 0.0;
283 dfMx[0][1] = -cz*sy*cx + sz*sx;
284 dfMx[0][2] = cz*sy*sx + sz*cx;
285 dfMx[1][0] = 0.0;
286 dfMx[1][1] = -sz*sy*cx - cz*sx;
287 dfMx[1][2] = sz*sy*sx - cz*cx;
288 dfMx[2][0] = 0.0;
289 dfMx[2][1] = cy*cx;
290 dfMx[2][2] = -cy*sx;
291}
292//______________________________________________________________________
293void AliITSAlignmentModule::dfMdthy(Double_t dfMx[3][3]){
294 Double_t sx,cx,sy,cy,sz,cz;
295
296 sx = TMath::Sin(fangles[0]); cx = TMath::Cos(fangles[0]);
297 sy = TMath::Sin(fangles[1]); cy = TMath::Cos(fangles[1]);
298 sz = TMath::Sin(fangles[2]); cz = TMath::Cos(fangles[2]);
299 dfMx[0][0] = -cz*sy;
300 dfMx[0][1] = -cz*cy*sx;
301 dfMx[0][2] = -cz*cy*cx;
302 dfMx[1][0] = -sz*sy;
303 dfMx[1][1] = -sz*cy*sx;
304 dfMx[1][2] = -sz*cy*cx;
305 dfMx[2][0] = cy;
306 dfMx[2][1] = -sy*sx;
307 dfMx[2][2] = -sy*cx;
308}
309//______________________________________________________________________
310void AliITSAlignmentModule::dfMdthz(Double_t dfMx[3][3]){
311 Double_t sx,cx,sy,cy,sz,cz;
312
313 sx = TMath::Sin(fangles[0]); cx = TMath::Cos(fangles[0]);
314 sy = TMath::Sin(fangles[1]); cy = TMath::Cos(fangles[1]);
315 sz = TMath::Sin(fangles[2]); cz = TMath::Cos(fangles[2]);
316 dfMx[0][0] = -sz*cy;
317 dfMx[0][1] = sz*sy*sx - cz*cx;
318 dfMx[0][2] = sz*sy*cx + cz*sx;
319 dfMx[1][0] = cz*cy;
320 dfMx[1][1] = -cz*sy*sx - sz*cx;
321 dfMx[1][2] = -cz*sy*cx + sz*sx;
322 dfMx[2][0] = 0.0;
323 dfMx[2][1] = 0.0;
324 dfMx[2][2] = 0.0;
325}
326//______________________________________________________________________
327void AliITSAlignmentModule::LtoG(Double_t xl[],Double_t xg[]){
328
329 Int_t i;
330 for(i=0;i<3;i++) xg[i] = fx0[i];
331 for(i=0;i<3;i++)for(Int_t j=0;j<3;j++) xg[i] += fM[j][i]*xl[j];
332}
333//______________________________________________________________________
334void AliITSAlignmentModule::GtoL(Double_t xg[],Double_t xl[]){
335 Int_t i;
336 for(i=0;i<3;i++) xl[i] = 0.0;
337 for(i=0;i<3;i++)for(Int_t j=0;j<3;j++)xl[i]+=fM[i][j]*(xg[j]-fx0[j]);
338}
339//______________________________________________________________________
340Double_t AliITSAlignmentModule::Chi2(Double_t p[]){
341 Int_t i,j,k,l,n=0,indx;
342 Double_t chi=0.0,xo[3],xi[3],xg[3],Ex[3][3];
343 AliITSAlignmentTrack *tr;
344
345 for(i=0;i<3;i++) fx0[i] = p[i];
346 SetByAngles(&(p[3]));
347
348 for(i=0;i<ftrksM->GetEntriesFast();i++) {
349 tr = (AliITSAlignmentTrack *)(ftrksM->At(i));
350 for(j=0;j<tr->GetNumberOfClustersSl();j++){
351 tr->GetIndex(j,indx);
352 if(indx==findex){
353 n++;
354 tr->GetPointG(j,(Double_t *)xi);
355 tr->func(xi,xo);
356 tr->GetPointL(j,(Double_t *)xi);
357 tr->GetErrorG(j,(Double_t **)Ex);
358 LtoG(xi,xg);
359 for(k=0;k<3;k++)for(l=0;l<3;l++)
360 chi += (xg[k] - xo[k])*Ex[k][l]*(xg[l]-xo[l]);
361 } // end if indx==findex
362 } // end for j
363 } // end for i
364 if(n<7) return chi;
365 return chi/(Double_t)(n-6);
366}
367//______________________________________________________________________
368void AliITSAlignmentModule::dChi2(Double_t p[],Double_t dChi2[]){
369 Int_t i,j,k,l,m,n=0,indx;
370 Double_t chi[6]={0.0,0.0,0.0,0.0,0.0,0.0},xo[3],xi[3],xg[3],Ex[3][3];
371 Double_t dxdp[3][6],fMx[3][3],fMy[3][3],fMz[3][3];
372 AliITSAlignmentTrack *tr;
373
374 for(i=0;i<3;i++) fx0[i] = p[i];
375 SetByAngles(&(p[3]));
376 dfMdthx(fMx);
377 dfMdthy(fMy);
378 dfMdthz(fMz);
379 for(i=0;i<3;i++)for(j=0;j<6;j++) dxdp[i][j] = 0.0;
380 dxdp[0][0] = 1.0; // dx/dx
381 dxdp[1][1] = 1.0; // dy/dy
382 dxdp[2][2] = 1.0; // dz/dz
383
384 for(i=0;i<ftrksM->GetEntriesFast();i++) {
385 tr = (AliITSAlignmentTrack *)(ftrksM->At(i));
386 for(j=0;j<tr->GetNumberOfClustersSl();j++){
387 tr->GetIndex(j,indx);
388 if(indx==findex){
389 n++;
390 tr->GetPointG(j,(Double_t *)xi);
391 tr->func(xi,xo);
392 tr->GetPointL(j,(Double_t *)xi);
393 tr->GetErrorG(j,(Double_t **)Ex);
394 LtoG(xi,xg);
395 for(m=0;m<3;m++) for(k=0;k<3;k++){
396 dxdp[m][3] += fMx[m][k]*xi[k];
397 dxdp[m][4] += fMy[m][k]*xi[k];
398 dxdp[m][5] += fMz[m][k]*xi[k];
399 } // end for m
400 for(m=0;m<6;m++){
401 for(k=0;k<3;k++)for(l=0;l<3;l++)
402 chi[m] += (xg[k] - xo[k])*Ex[k][l]*dxdp[l][m];
403 } // end for m
404 } // end if indx==findex
405 } // end for j
406 } // end for i
407 if(n<7) return ;
408 for(m=0;m<6;m++) chi[m] /= (Double_t)(n-6);
409 return;
410}
411//______________________________________________________________________
412void AliITSAlignmentModule::AlignModule(){
413 static Int_t npar=6;
414 Int_t iter,i;
415 //Double_t p[npar],fret,gtol=1.0E-5;
416 Double_t p[6],fret,gtol=1.0E-5;
417
418 for(i=0;i<3;i++) {p[i] = fx0[i]; p[i+3] = fangles[i];}
419 MRVMminimization(npar,(Double_t *)p,fret,gtol,iter);
420 for(i=0;i<3;i++) {fx0[i] = p[i]; fangles[i] = p[i+3];}
421 printf("AlignModule #%d:Xt=(%e,%e,%e) cm angles=(%e,%e,%e)rad,"
422 " Chi2=%e loops=%d\n",findex,
423 fx0[0],fx0[1],fx0[2],fangles[0],fangles[1],fangles[2],fret,iter);
424}
425//______________________________________________________________________
426void AliITSAlignmentModule::Streamer(TBuffer &R__b){
427 // Stream an object of class AliITSAlignmentModule.
428
429 if (R__b.IsReading()) {
430 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
431 TObject::Streamer(R__b);
432 R__b >> findex;
433 R__b >> flay;
434 R__b >> flad;
435 R__b >> fdet;
436 R__b >> ftrksM;
437 R__b >> fChi2;
438// R__b.ReadStaticArray(fx0);
439// R__b.ReadStaticArray((double*)fM);
440// R__b.ReadStaticArray(fangles);
441 } else {
442 R__b.WriteVersion(AliITSAlignmentModule::IsA());
443 TObject::Streamer(R__b);
444 R__b << findex;
445 R__b << flay;
446 R__b << flad;
447 R__b << fdet;
448 R__b << ftrksM;
449 R__b << fChi2;
450// R__b.WriteArray(fx0, 3);
451// R__b.WriteArray((double*)fM, 9);
452// R__b.WriteArray(fangles, 3);
453 }
454}