New version of alignment framework.
[u/mrichter/AliRoot.git] / STEER / AliTrackResidualsFast.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 //-----------------------------------------------------------------
17 //   Implementation of the derived class for track residuals
18 //   based on linear chi2 minimization (in approximation of
19 //   small alignment angles and translations)
20 //
21 //-----------------------------------------------------------------
22
23 #include <TMinuit.h>
24 #include <TGeoMatrix.h>
25
26 #include "AliAlignObj.h"
27 #include "AliTrackPointArray.h"
28 #include "AliTrackResidualsFast.h"
29
30 ClassImp(AliTrackResidualsFast)
31
32 //______________________________________________________________________________
33 AliTrackResidualsFast::AliTrackResidualsFast():AliTrackResiduals()
34 {
35   // Default constructor
36   for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
37   fSumR = 0;
38 }
39
40 //______________________________________________________________________________
41 AliTrackResidualsFast::AliTrackResidualsFast(Int_t ntracks):
42   AliTrackResiduals(ntracks)
43 {
44   // Constructor
45   for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
46   fSumR = 0;
47 }
48  
49 //______________________________________________________________________________
50 AliTrackResidualsFast::AliTrackResidualsFast(const AliTrackResidualsFast &res):
51   AliTrackResiduals(res)
52 {
53   // Copy constructor
54   for (Int_t i = 0; i < 16; i++) fSum[i] = res.fSum[i];
55   fSumR = res.fSumR;
56 }
57
58 //______________________________________________________________________________
59 AliTrackResidualsFast &AliTrackResidualsFast::operator= (const AliTrackResidualsFast& res)
60 {
61   // Assignment operator
62  ((AliTrackResiduals *)this)->operator=(res);
63  for (Int_t i = 0; i < 16; i++) fSum[i] = res.fSum[i];
64  fSumR = res.fSumR;
65
66  return *this;
67 }
68
69 //______________________________________________________________________________
70 Bool_t AliTrackResidualsFast::Minimize()
71 {
72   // Implementation of fast linear Chi2
73   // based minimization of track residuals sum
74   for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
75   fSumR = 0;
76
77   AliTrackPoint p1,p2;
78
79   for (Int_t itrack = 0; itrack < fLast; itrack++) {
80     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
81     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
82       fVolArray[itrack]->GetPoint(p1,ipoint);
83       fTrackArray[itrack]->GetPoint(p2,ipoint);
84       AddPoints(p1,p2);
85     }
86   }
87
88   return Update();
89
90   // debug info
91 //   Float_t chi2 = 0;
92 //   for (Int_t itrack = 0; itrack < fLast; itrack++) {
93 //     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
94 //     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
95 //       fVolArray[itrack]->GetPoint(p1,ipoint);
96 //       fAlignObj->Transform(p1);
97 //       fTrackArray[itrack]->GetPoint(p2,ipoint);
98 //       Float_t residual = p2.GetResidual(p1,kFALSE);
99 //       chi2 += residual;
100 //     }
101 //   }
102 //   printf("Final chi2 = %f\n",chi2);
103 }
104
105 //______________________________________________________________________________
106 void AliTrackResidualsFast::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
107 {
108   // Update the sums used for
109   // the linear chi2 minimization
110   Float_t xyz[3],xyzp[3];
111   p.GetXYZ(xyz); pprime.GetXYZ(xyzp);
112   Double_t weight = 1;
113   weight *= weight;
114   fSum[0] += weight;
115   fSum[1] += xyz[0]*weight;
116   fSum[2] += xyz[1]*weight;
117   fSum[3] += xyz[2]*weight;
118   fSum[4] += (xyz[0]*xyz[0]+xyz[1]*xyz[1])*weight;
119   fSum[5] += (xyz[0]*xyz[0]+xyz[2]*xyz[2])*weight;
120   fSum[6] += (xyz[1]*xyz[1]+xyz[2]*xyz[2])*weight;
121   fSum[7] -= xyz[0]*xyz[1]*weight;
122   fSum[8] -= xyz[0]*xyz[2]*weight;
123   fSum[9] -= xyz[1]*xyz[2]*weight;
124   fSum[10] += (xyzp[0]-xyz[0])*weight;
125   fSum[11] += (xyzp[1]-xyz[1])*weight;
126   fSum[12] += (xyzp[2]-xyz[2])*weight;
127   fSum[13] += (xyzp[2]*xyz[1]-xyzp[1]*xyz[2])*weight;
128   fSum[14] += (xyzp[0]*xyz[2]-xyzp[2]*xyz[0])*weight;
129   fSum[15] += (xyzp[1]*xyz[0]-xyzp[0]*xyz[1])*weight;
130
131   fSumR += ((xyzp[0]-xyz[0])*(xyzp[0]-xyz[0])+
132             (xyzp[1]-xyz[1])*(xyzp[1]-xyz[1])+
133             (xyzp[2]-xyz[2])*(xyzp[2]-xyz[2]))*weight;
134
135   fNdf += 3;
136 }
137
138 //______________________________________________________________________________
139 Bool_t AliTrackResidualsFast::Update()
140 {
141   // Find the alignment parameters
142   // by using the already accumulated
143   // sums
144   TMatrixDSym     smatrix(6);
145   TMatrixD        sums(1,6);
146
147   smatrix(0,0) = fSum[0]; smatrix(1,1) = fSum[0]; smatrix(2,2)=fSum[0];
148   smatrix(3,3) = fSum[6]; smatrix(4,4) = fSum[5]; smatrix(5,5)=fSum[4];
149
150   smatrix(0,1) = 0; smatrix(0,2) = 0;
151   smatrix(0,3) = 0; smatrix(0,4) = fSum[3]; smatrix(0,5) = -fSum[2];
152   smatrix(1,2) = 0;
153   smatrix(1,3) = -fSum[3]; smatrix(1,4) = 0; smatrix(1,5) = fSum[1];
154   smatrix(2,3) = fSum[2]; smatrix(2,4) = -fSum[1]; smatrix(2,5) = 0;
155
156   smatrix(3,4) = fSum[7]; smatrix(3,5) = fSum[8];
157   smatrix(4,5) = fSum[9];
158
159   sums(0,0)    = fSum[10]; sums(0,1) = fSum[11]; sums(0,2) = fSum[12];
160   sums(0,3)    = fSum[13]; sums(0,4) = fSum[14]; sums(0,5) = fSum[15];
161
162   smatrix.Invert();
163   if (!smatrix.IsValid()) return kFALSE;
164
165   TMatrixD res = sums*smatrix;
166   fAlignObj->SetPars(res(0,0),res(0,1),res(0,2),
167                      TMath::RadToDeg()*res(0,3),
168                      TMath::RadToDeg()*res(0,4),
169                      TMath::RadToDeg()*res(0,5));
170   TMatrixD  tmp = res*sums.T();
171   fChi2 = fSumR - tmp(0,0);
172   fNdf -= 6;
173
174   return kTRUE;
175 }