New convention for the solenoid field
[u/mrichter/AliRoot.git] / STEER / AliTrackResidualsFast.cxx
CommitLineData
46ae650f 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
30ClassImp(AliTrackResidualsFast)
31
32//______________________________________________________________________________
33AliTrackResidualsFast::AliTrackResidualsFast():AliTrackResiduals()
34{
35 // Default constructor
36 for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
37 fSumR = 0;
38}
39
40//______________________________________________________________________________
41AliTrackResidualsFast::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//______________________________________________________________________________
50AliTrackResidualsFast::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//______________________________________________________________________________
59AliTrackResidualsFast &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//______________________________________________________________________________
70Bool_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//______________________________________________________________________________
106void 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//______________________________________________________________________________
139Bool_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}