Correcting compiler warnings
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignment.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 /* $Id$ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONAlignmentV5
20 /// Alignment class fro the ALICE DiMuon spectrometer 
21 ///
22 /// MUON specific alignment class which interface to AliMillepede.   
23 /// For each track ProcessTrack calculates the local and global derivatives
24 /// at each hit and fill the corresponding local equations. Provide methods for
25 /// fixing or constraining detection elements for best results. 
26 ///
27 /// \author Bruce Becker, Javier Castillo
28 //-----------------------------------------------------------------------------
29
30 #include "AliMUONAlignment.h"
31 #include "AliMUONTrack.h"
32 #include "AliMUONRawCluster.h"
33 #include "AliMUONTrackParam.h"
34 #include "AliMUONHitForRec.h"
35 #include "AliMUONGeometryTransformer.h"
36 #include "AliMUONGeometryModuleTransformer.h"
37 #include "AliMUONGeometryDetElement.h"
38 #include "AliMUONGeometryBuilder.h"
39 #include "AliMUONConstants.h"
40 #include "AliMillepede.h"
41
42 #include "AliMpExMap.h"
43
44 #include "AliLog.h"
45
46 #include "TSystem.h"
47
48 /// \cond CLASSIMP
49 ClassImp(AliMUONAlignment)
50 /// \endcond
51
52   Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
53   Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
54   Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
55   Int_t AliMUONAlignment::fgNParCh = 3;
56   Int_t AliMUONAlignment::fgNTrkMod = 16;
57   Int_t AliMUONAlignment::fgNCh = 10;
58   Int_t AliMUONAlignment::fgNSt = 5;
59
60 AliMUONAlignment::AliMUONAlignment() 
61   : TObject(),
62     fBFieldOn(kTRUE),
63     fStartFac(16.), 
64     fResCutInitial(100.), 
65     fResCut(100.),
66     fMillepede(0),
67     fTrackParamAtHit(0),
68     fHitForRecAtHit(0),
69     fTrack(0),
70     fRecHit(0),
71     fTrackParam(0),
72     fNGlobal(fgNDetElem*fgNParCh),
73     fNLocal(4),
74     fNStdDev(3),
75     fDetElemId(0),
76     fDetElemNumber(0),
77     fPhi(0.),
78     fCosPhi(1.),
79     fSinPhi(0.),
80     fTransform(0)
81 {
82   /// Default constructor, setting define alignment parameters
83   fSigma[0] = 1.0e-1;
84   fSigma[1] = 1.0e-2;
85
86   fDoF[0] = kTRUE;  fDoF[1] = kTRUE;  fDoF[2] = kTRUE;
87   fAllowVar[0] = 0.05;  fAllowVar[1] = 0.05;  fAllowVar[2] = 0.001;
88   
89   AliInfo(Form("fAllowVar[0]: %f\t fAllowVar[1]: %f\t fPhi: %f\t fgNDetElem: %i\t fNGlobal: %i\t fNLocal: %i",fAllowVar[0],fAllowVar[1],fPhi,fgNDetElem,fNGlobal,fNLocal));
90
91   fMillepede = new AliMillepede();
92
93   Init(fNGlobal, fNLocal, fNStdDev);
94
95   ResetLocalEquation();
96   AliInfo("Parameters initialized to zero");
97
98 }
99
100 AliMUONAlignment::~AliMUONAlignment() {
101   /// Destructor
102 }
103
104 void AliMUONAlignment::Init(Int_t nGlobal,  /* number of global paramers */
105                             Int_t nLocal,   /* number of local parameters */
106                             Int_t nStdDev   /* std dev cut */ )
107 {
108   /// Initialization of AliMillepede. Fix parameters, define constraints ...
109   fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
110
111   Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
112   Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
113
114   AllowVariations(bStOnOff);
115
116   // Fix parameters or add constraints here
117   for (Int_t iSt=0; iSt<5; iSt++)
118     if (!bStOnOff[iSt]) FixStation(iSt+1);
119
120   FixHalfSpectrometer(bStOnOff,bSpecLROnOff);
121
122   ResetConstraints();
123   
124   // Define global constrains to be applied
125   // X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
126   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
127   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
128   //  AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
129
130   // Other possible way to add constrains
131   bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
132   bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
133 //   AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
134
135   bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE;
136   //  AddConstraints(bStOnOff,bVarXYT);
137   
138   // Set iterations
139   if (fStartFac>1) fMillepede->SetIterations(fStartFac);          
140 }
141
142 void AliMUONAlignment::FixStation(Int_t iSt){
143   /// Fix all detection elements of station iSt
144   Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0; 
145   Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1]; 
146   for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){    
147     FixParameter(i*fgNParCh+0, 0.0);
148     FixParameter(i*fgNParCh+1, 0.0);
149     FixParameter(i*fgNParCh+2, 0.0);
150   }
151 }
152
153 void AliMUONAlignment::FixHalfSpectrometer(Bool_t *lStOnOff, Bool_t *lSpecLROnOff){
154   /// Fix left or right detector
155   for (Int_t i = 0; i < fgNDetElem; i++){    
156     Int_t iCh=0;
157     for (iCh=1; iCh<=fgNCh; iCh++){
158       if (i<fgSNDetElemCh[iCh-1]) break;
159     }
160     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0; 
161     if (iSt){
162       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
163       if (iCh>=1 && iCh<=4){
164         if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0]){ // From track crossings
165           FixParameter(i*fgNParCh+0, 0.0);
166           FixParameter(i*fgNParCh+1, 0.0);
167           FixParameter(i*fgNParCh+2, 0.0);
168         }
169         if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1]){ // From track crossings
170           FixParameter(i*fgNParCh+0, 0.0);
171           FixParameter(i*fgNParCh+1, 0.0);
172           FixParameter(i*fgNParCh+2, 0.0);
173         }
174       }
175       if (iCh>=5 && iCh<=6){
176         if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0]){
177           FixParameter(i*fgNParCh+0, 0.0);
178           FixParameter(i*fgNParCh+1, 0.0);
179           FixParameter(i*fgNParCh+2, 0.0);
180         }
181         if (((lDetElemNumber>=0&&lDetElemNumber<=4) || 
182              (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1]){
183           FixParameter(i*fgNParCh+0, 0.0);
184           FixParameter(i*fgNParCh+1, 0.0);
185           FixParameter(i*fgNParCh+2, 0.0);
186         }
187       }
188       if (iCh>=7 && iCh<=10){
189         if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0]){
190           FixParameter(i*fgNParCh+0, 0.0);
191           FixParameter(i*fgNParCh+1, 0.0);
192           FixParameter(i*fgNParCh+2, 0.0);
193         }
194         if (((lDetElemNumber>=0&&lDetElemNumber<=6) || 
195              (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1]){
196           FixParameter(i*fgNParCh+0, 0.0);
197           FixParameter(i*fgNParCh+1, 0.0);
198           FixParameter(i*fgNParCh+2, 0.0);
199         }
200       }
201     }
202   }
203 }
204
205 void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){
206   /// Set non linear parameter flag selected stations and degrees of freedom
207   for (Int_t i = 0; i < fgNDetElem; i++){    
208     Int_t iCh=0;
209     for (iCh=1; iCh<=fgNCh; iCh++){
210       if (i<fgSNDetElemCh[iCh-1]) break;
211     }
212     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0; 
213     if (iSt){
214       if (lVarXYT[0]) { // X constraint
215         SetNonLinear(i*fgNParCh+0);
216       }
217       if (lVarXYT[1]) { // Y constraint
218         SetNonLinear(i*fgNParCh+1);
219       }
220       if (lVarXYT[2]) { // T constraint
221         SetNonLinear(i*fgNParCh+2);
222       }
223     }
224   }
225 }
226
227 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT){
228   /// Add constraint equations for selected stations and degrees of freedom 
229   for (Int_t i = 0; i < fgNDetElem; i++){    
230     Int_t iCh=0;
231     for (iCh=1; iCh<=fgNCh; iCh++){
232       if (i<fgSNDetElemCh[iCh-1]) break;
233     }
234     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0; 
235     if (iSt){
236       if (lVarXYT[0]) { // X constraint
237         fConstraintX[i*fgNParCh+0]=1.0;
238       }
239       if (lVarXYT[1]) { // Y constraint
240         fConstraintY[i*fgNParCh+1]=1.0;
241       }
242       if (lVarXYT[2]) { // T constraint
243         fConstraintP[i*fgNParCh+2]=1.0;
244       }
245     }
246   }
247   if (lVarXYT[0]) { // X constraint
248     AddConstraint(fConstraintX,0.0);
249   }
250   if (lVarXYT[1]) { // Y constraint
251     AddConstraint(fConstraintY,0.0);
252   }
253   if (lVarXYT[2]) { // T constraint
254     AddConstraint(fConstraintP,0.0);
255   }
256 }
257
258 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t *lDetTLBR, Bool_t *lSpecLROnOff){
259   /// Add constraint equations for selected stations, degrees of freedom and detector half 
260   Double_t lDetElemLocX = 0.;
261   Double_t lDetElemLocY = 0.;
262   Double_t lDetElemLocZ = 0.;
263   Double_t lDetElemGloX = 0.;
264   Double_t lDetElemGloY = 0.;
265   Double_t lDetElemGloZ = 0.;
266   Double_t lMeanY = 0.;
267   Double_t lSigmaY = 0.;
268   Double_t lMeanZ = 0.;
269   Double_t lSigmaZ = 0.;
270   Int_t lNDetElem = 0;
271   for (Int_t i = 0; i < fgNDetElem; i++){    
272     Int_t iCh=0;
273     for (iCh=1; iCh<=fgNCh; iCh++){
274       if (i<fgSNDetElemCh[iCh-1]) break;
275     }
276     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0; 
277     if (iSt){
278       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
279       Int_t lDetElemId = iCh*100+lDetElemNumber;
280       fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
281                                lDetElemGloX,lDetElemGloY,lDetElemGloZ);
282       if (iCh>=1 && iCh<=4){
283         if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings
284           lMeanY += lDetElemGloY;
285           lSigmaY += lDetElemGloY*lDetElemGloY;
286           lMeanZ += lDetElemGloZ;
287           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
288           lNDetElem++;
289         }
290         if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings
291           lMeanY += lDetElemGloY;
292           lSigmaY += lDetElemGloY*lDetElemGloY;
293           lMeanZ += lDetElemGloZ;
294           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
295           lNDetElem++;
296         }
297       }
298       if (iCh>=5 && iCh<=6){
299         if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){
300           lMeanY += lDetElemGloY;
301           lSigmaY += lDetElemGloY*lDetElemGloY;
302           lMeanZ += lDetElemGloZ;
303           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
304           lNDetElem++;
305         }
306         if (((lDetElemNumber>=0&&lDetElemNumber<=4) || 
307              (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){
308           lMeanY += lDetElemGloY;
309           lSigmaY += lDetElemGloY*lDetElemGloY;
310           lMeanZ += lDetElemGloZ;
311           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
312           lNDetElem++;
313         }
314       }
315       if (iCh>=7 && iCh<=10){
316         if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){
317           lMeanY += lDetElemGloY;
318           lSigmaY += lDetElemGloY*lDetElemGloY;
319           lMeanZ += lDetElemGloZ;
320           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
321           lNDetElem++;
322         }
323         if (((lDetElemNumber>=0&&lDetElemNumber<=6) || 
324              (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){
325           lMeanY += lDetElemGloY;
326           lSigmaY += lDetElemGloY*lDetElemGloY;
327           lMeanZ += lDetElemGloZ;
328           lSigmaZ += lDetElemGloZ*lDetElemGloZ;
329           lNDetElem++;
330         }
331       }
332     }
333   }
334   lMeanY /= lNDetElem;
335   lSigmaY /= lNDetElem;
336   lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
337   lMeanZ /= lNDetElem;
338   lSigmaZ /= lNDetElem;
339   lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
340   AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ));  
341
342   for (Int_t i = 0; i < fgNDetElem; i++){    
343     Int_t iCh=0;
344     for (iCh=1; iCh<=fgNCh; iCh++){
345       if (i<fgSNDetElemCh[iCh-1]) break;
346     }
347     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0; 
348     if (iSt){
349       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
350       Int_t lDetElemId = iCh*100+lDetElemNumber;
351       fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
352                                lDetElemGloX,lDetElemGloY,lDetElemGloZ);
353       if (lVarXYT[0]) { // X constraint
354         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
355         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
356         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
357         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
358       }
359       if (lVarXYT[1]) { // Y constraint
360         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
361         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
362         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
363         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
364       }
365       if (lVarXYT[2]) { // P constraint
366         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
367         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
368         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
369         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
370       }
371       if (lVarXYT[3]) { // X-Z shearing
372         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
373         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
374         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
375         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
376       }
377       if (lVarXYT[4]) { // Y-Z shearing
378         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
379         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
380         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
381         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
382       }
383       if (lVarXYT[5]) { // P-Z rotation
384         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
385         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
386         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
387         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
388       }
389       if (lVarXYT[6]) { // X-Y shearing
390         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
391         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
392         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
393         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
394       }
395       if (lVarXYT[7]) { // Y-Y scaling
396         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
397         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
398         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
399         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
400       }
401       if (lVarXYT[8]) { // P-Y rotation
402         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
403         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
404         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
405         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
406       }
407     }
408   }
409   if (lVarXYT[0]) { // X constraint
410     if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
411     if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
412     if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
413     if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
414   }
415   if (lVarXYT[1]) { // Y constraint
416     if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
417     if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
418     if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
419     if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
420   }
421   if (lVarXYT[2]) { // T constraint
422     if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
423     if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
424     if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
425     if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
426   }
427   if (lVarXYT[3]) { // X-Z constraint
428     if (lDetTLBR[0]) AddConstraint(fConstraintXZT,0.0); // Top half
429     if (lDetTLBR[1]) AddConstraint(fConstraintXZL,0.0); // Left half
430     if (lDetTLBR[2]) AddConstraint(fConstraintXZB,0.0); // Bottom half
431     if (lDetTLBR[3]) AddConstraint(fConstraintXZR,0.0); // Right half
432   }
433   if (lVarXYT[4]) { // Y-Z constraint
434     if (lDetTLBR[0]) AddConstraint(fConstraintYZT,0.0); // Top half
435     if (lDetTLBR[1]) AddConstraint(fConstraintYZL,0.0); // Left half
436     if (lDetTLBR[2]) AddConstraint(fConstraintYZB,0.0); // Bottom half
437     if (lDetTLBR[3]) AddConstraint(fConstraintYZR,0.0); // Right half
438   }
439   if (lVarXYT[5]) { // P-Z constraint
440     if (lDetTLBR[0]) AddConstraint(fConstraintPZT,0.0); // Top half
441     if (lDetTLBR[1]) AddConstraint(fConstraintPZL,0.0); // Left half
442     if (lDetTLBR[2]) AddConstraint(fConstraintPZB,0.0); // Bottom half
443     if (lDetTLBR[3]) AddConstraint(fConstraintPZR,0.0); // Right half
444   }
445   if (lVarXYT[6]) { // X-Y constraint
446     if (lDetTLBR[0]) AddConstraint(fConstraintXYT,0.0); // Top half
447     if (lDetTLBR[1]) AddConstraint(fConstraintXYL,0.0); // Left half
448     if (lDetTLBR[2]) AddConstraint(fConstraintXYB,0.0); // Bottom half
449     if (lDetTLBR[3]) AddConstraint(fConstraintXYR,0.0); // Right half
450   }
451   if (lVarXYT[7]) { // Y-Y constraint
452     if (lDetTLBR[0]) AddConstraint(fConstraintYYT,0.0); // Top half
453     if (lDetTLBR[1]) AddConstraint(fConstraintYYL,0.0); // Left half
454     if (lDetTLBR[2]) AddConstraint(fConstraintYYB,0.0); // Bottom half
455     if (lDetTLBR[3]) AddConstraint(fConstraintYYR,0.0); // Right half
456   }
457   if (lVarXYT[8]) { // P-Y constraint
458     if (lDetTLBR[0]) AddConstraint(fConstraintPYT,0.0); // Top half
459     if (lDetTLBR[1]) AddConstraint(fConstraintPYL,0.0); // Left half
460     if (lDetTLBR[2]) AddConstraint(fConstraintPYB,0.0); // Bottom half
461     if (lDetTLBR[3]) AddConstraint(fConstraintPYR,0.0); // Right half
462   }
463 }
464
465 void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t /*lWeight*/){
466   /// Set constrain equation for top half of spectrometer
467   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
468   if (lCh>=1 && lCh<=4){
469     if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
470       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
471     }
472   }
473   if (lCh>=5 && lCh<=6){
474     if (lDetElemNumber>=0&&lDetElemNumber<=9){
475       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
476     }
477   }
478   if (lCh>=7 && lCh<=10){
479     if (lDetElemNumber>=0&&lDetElemNumber<=13){
480       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
481     }
482   }
483 }
484
485 void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight){
486   /// Set constrain equation for left half of spectrometer
487   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
488   if (lCh>=1 && lCh<=4){
489     if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
490       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
491     }
492   }
493   if (lCh>=5 && lCh<=6){
494     if (lDetElemNumber>=5&&lDetElemNumber<=13){
495       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
496     }
497   }
498   if (lCh>=7 && lCh<=10){
499     if (lDetElemNumber>=7&&lDetElemNumber<=19){
500       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
501     }
502   }
503 }
504
505 void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t /*lWeight*/){
506   /// Set constrain equation for bottom half of spectrometer
507   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
508   if (lCh>=1 && lCh<=4){
509     if (lDetElemNumber==2 || lDetElemNumber==3){ // From track crossings
510       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
511     }
512   }
513   if (lCh>=5 && lCh<=6){
514     if ((lDetElemNumber>=9&&lDetElemNumber<=17) || 
515         (lDetElemNumber==0)){
516       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
517     }
518   }
519   if (lCh>=7 && lCh<=10){
520     if ((lDetElemNumber>=13&&lDetElemNumber<=25) || 
521         (lDetElemNumber==0)){
522       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
523     }
524   }
525 }
526
527 void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight){
528   /// Set constrain equation for right half of spectrometer
529   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
530   if (lCh>=1 && lCh<=4){
531     if (lDetElemNumber==0 || lDetElemNumber==3){ // From track crossings
532       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
533     }
534   }
535   if (lCh>=5 && lCh<=6){
536     if ((lDetElemNumber>=0&&lDetElemNumber<=4) || 
537         (lDetElemNumber>=14&&lDetElemNumber<=17)){
538       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
539     }
540   }
541   if (lCh>=7 && lCh<=10){
542     if ((lDetElemNumber>=0&&lDetElemNumber<=6) || 
543         (lDetElemNumber>=20&&lDetElemNumber<=25)){
544       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
545     }
546   }
547 }
548
549 void AliMUONAlignment::ResetConstraints(){
550   /// Reset all constraint equations
551   for (Int_t i = 0; i < fgNDetElem; i++){    
552     fConstraintX[i*fgNParCh+0]=0.0;
553     fConstraintX[i*fgNParCh+1]=0.0;
554     fConstraintX[i*fgNParCh+2]=0.0;
555     fConstraintY[i*fgNParCh+0]=0.0;
556     fConstraintY[i*fgNParCh+1]=0.0;
557     fConstraintY[i*fgNParCh+2]=0.0;
558     fConstraintP[i*fgNParCh+0]=0.0;
559     fConstraintP[i*fgNParCh+1]=0.0;
560     fConstraintP[i*fgNParCh+2]=0.0;
561     fConstraintXT[i*fgNParCh+0]=0.0;
562     fConstraintXT[i*fgNParCh+1]=0.0;
563     fConstraintXT[i*fgNParCh+2]=0.0;
564     fConstraintYT[i*fgNParCh+0]=0.0;
565     fConstraintYT[i*fgNParCh+1]=0.0;
566     fConstraintYT[i*fgNParCh+2]=0.0;
567     fConstraintPT[i*fgNParCh+0]=0.0;
568     fConstraintPT[i*fgNParCh+1]=0.0;
569     fConstraintPT[i*fgNParCh+2]=0.0;
570     fConstraintXZT[i*fgNParCh+0]=0.0;
571     fConstraintXZT[i*fgNParCh+1]=0.0;
572     fConstraintXZT[i*fgNParCh+2]=0.0;
573     fConstraintYZT[i*fgNParCh+0]=0.0;
574     fConstraintYZT[i*fgNParCh+1]=0.0;
575     fConstraintYZT[i*fgNParCh+2]=0.0;
576     fConstraintPZT[i*fgNParCh+0]=0.0;
577     fConstraintPZT[i*fgNParCh+1]=0.0;
578     fConstraintPZT[i*fgNParCh+2]=0.0;
579     fConstraintXYT[i*fgNParCh+0]=0.0;
580     fConstraintXYT[i*fgNParCh+1]=0.0;
581     fConstraintXYT[i*fgNParCh+2]=0.0;
582     fConstraintYYT[i*fgNParCh+0]=0.0;
583     fConstraintYYT[i*fgNParCh+1]=0.0;
584     fConstraintYYT[i*fgNParCh+2]=0.0;
585     fConstraintPYT[i*fgNParCh+0]=0.0;
586     fConstraintPYT[i*fgNParCh+1]=0.0;
587     fConstraintPYT[i*fgNParCh+2]=0.0;
588     fConstraintXL[i*fgNParCh+0]=0.0;
589     fConstraintXL[i*fgNParCh+1]=0.0;
590     fConstraintXL[i*fgNParCh+2]=0.0;
591     fConstraintYL[i*fgNParCh+0]=0.0;
592     fConstraintYL[i*fgNParCh+1]=0.0;
593     fConstraintYL[i*fgNParCh+2]=0.0;
594     fConstraintPL[i*fgNParCh+0]=0.0;
595     fConstraintPL[i*fgNParCh+1]=0.0;
596     fConstraintPL[i*fgNParCh+2]=0.0;
597     fConstraintXZL[i*fgNParCh+0]=0.0;
598     fConstraintXZL[i*fgNParCh+1]=0.0;
599     fConstraintXZL[i*fgNParCh+2]=0.0;
600     fConstraintYZL[i*fgNParCh+0]=0.0;
601     fConstraintYZL[i*fgNParCh+1]=0.0;
602     fConstraintYZL[i*fgNParCh+2]=0.0;
603     fConstraintPZL[i*fgNParCh+0]=0.0;
604     fConstraintPZL[i*fgNParCh+1]=0.0;
605     fConstraintPZL[i*fgNParCh+2]=0.0;
606     fConstraintXYL[i*fgNParCh+0]=0.0;
607     fConstraintXYL[i*fgNParCh+1]=0.0;
608     fConstraintXYL[i*fgNParCh+2]=0.0;
609     fConstraintYYL[i*fgNParCh+0]=0.0;
610     fConstraintYYL[i*fgNParCh+1]=0.0;
611     fConstraintYYL[i*fgNParCh+2]=0.0;
612     fConstraintPYL[i*fgNParCh+0]=0.0;
613     fConstraintPYL[i*fgNParCh+1]=0.0;
614     fConstraintPYL[i*fgNParCh+2]=0.0;
615     fConstraintXB[i*fgNParCh+0]=0.0;
616     fConstraintXB[i*fgNParCh+1]=0.0;
617     fConstraintXB[i*fgNParCh+2]=0.0;
618     fConstraintYB[i*fgNParCh+0]=0.0;
619     fConstraintYB[i*fgNParCh+1]=0.0;
620     fConstraintYB[i*fgNParCh+2]=0.0;
621     fConstraintPB[i*fgNParCh+0]=0.0;
622     fConstraintPB[i*fgNParCh+1]=0.0;
623     fConstraintPB[i*fgNParCh+2]=0.0;
624     fConstraintXZB[i*fgNParCh+0]=0.0;
625     fConstraintXZB[i*fgNParCh+1]=0.0;
626     fConstraintXZB[i*fgNParCh+2]=0.0;
627     fConstraintYZB[i*fgNParCh+0]=0.0;
628     fConstraintYZB[i*fgNParCh+1]=0.0;
629     fConstraintYZB[i*fgNParCh+2]=0.0;
630     fConstraintPZB[i*fgNParCh+0]=0.0;
631     fConstraintPZB[i*fgNParCh+1]=0.0;
632     fConstraintPZB[i*fgNParCh+2]=0.0;
633     fConstraintXYB[i*fgNParCh+0]=0.0;
634     fConstraintXYB[i*fgNParCh+1]=0.0;
635     fConstraintXYB[i*fgNParCh+2]=0.0;
636     fConstraintYYB[i*fgNParCh+0]=0.0;
637     fConstraintYYB[i*fgNParCh+1]=0.0;
638     fConstraintYYB[i*fgNParCh+2]=0.0;
639     fConstraintPYB[i*fgNParCh+0]=0.0;
640     fConstraintPYB[i*fgNParCh+1]=0.0;
641     fConstraintPYB[i*fgNParCh+2]=0.0;
642     fConstraintXR[i*fgNParCh+0]=0.0;
643     fConstraintXR[i*fgNParCh+1]=0.0;
644     fConstraintXR[i*fgNParCh+2]=0.0;
645     fConstraintYR[i*fgNParCh+0]=0.0;
646     fConstraintYR[i*fgNParCh+1]=0.0;
647     fConstraintYR[i*fgNParCh+2]=0.0;
648     fConstraintPR[i*fgNParCh+0]=0.0;
649     fConstraintPR[i*fgNParCh+1]=0.0;
650     fConstraintPR[i*fgNParCh+2]=0.0;
651     fConstraintXZR[i*fgNParCh+0]=0.0;
652     fConstraintXZR[i*fgNParCh+1]=0.0;
653     fConstraintXZR[i*fgNParCh+2]=0.0;
654     fConstraintYZR[i*fgNParCh+0]=0.0;
655     fConstraintYZR[i*fgNParCh+1]=0.0;
656     fConstraintYZR[i*fgNParCh+2]=0.0;
657     fConstraintPZR[i*fgNParCh+0]=0.0;
658     fConstraintPZR[i*fgNParCh+1]=0.0;
659     fConstraintPZR[i*fgNParCh+2]=0.0;
660     fConstraintPZR[i*fgNParCh+0]=0.0;
661     fConstraintPZR[i*fgNParCh+1]=0.0;
662     fConstraintPZR[i*fgNParCh+2]=0.0;
663     fConstraintXYR[i*fgNParCh+0]=0.0;
664     fConstraintXYR[i*fgNParCh+1]=0.0;
665     fConstraintXYR[i*fgNParCh+2]=0.0;
666     fConstraintYYR[i*fgNParCh+0]=0.0;
667     fConstraintYYR[i*fgNParCh+1]=0.0;
668     fConstraintYYR[i*fgNParCh+2]=0.0;
669     fConstraintPYR[i*fgNParCh+0]=0.0;
670     fConstraintPYR[i*fgNParCh+1]=0.0;
671     fConstraintPYR[i*fgNParCh+2]=0.0;
672   }
673 }
674
675 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
676   /// Constrain equation defined by par to value
677   fMillepede->SetGlobalConstraint(par, value);
678   AliInfo("Adding constraint");
679 }
680
681 void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
682   /// Initialize global parameters with par array
683   fMillepede->SetGlobalParameters(par);
684   AliInfo("Init Global Parameters");
685 }
686  
687 void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
688   /// Parameter iPar is encourage to vary in [-value;value]. 
689   /// If value == 0, parameter is fixed
690   fMillepede->SetParSigma(iPar, value);
691   if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
692 }
693
694 void AliMUONAlignment::ResetLocalEquation()
695 {
696   /// Reset the derivative vectors
697   for(int i=0; i<fNLocal; i++) {
698     fLocalDerivatives[i] = 0.0;
699   }
700   for(int i=0; i<fNGlobal; i++) {
701     fGlobalDerivatives[i] = 0.0;
702   }
703 }
704
705 void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) {
706   /// Set allowed variation for selected stations based on fDoF and fAllowVar
707   for (Int_t iSt=1; iSt<=5; iSt++) {
708     if (bStOnOff[iSt-1]) {
709       Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0; 
710       Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1]; 
711       for (int i=0; i<fgNParCh; i++) {
712         AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));    
713         if (fDoF[i]) {
714           for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){    
715             FixParameter(j*fgNParCh+i, fAllowVar[i]);
716           }
717         }
718       }
719     }
720   }
721 }
722
723 void AliMUONAlignment::SetNonLinear(Int_t iPar  /* set non linear flag */ ) {
724   /// Set nonlinear flag for parameter iPar
725   fMillepede->SetNonLinear(iPar);
726   AliInfo(Form("Parameter %i set to non linear", iPar));
727 }
728
729 void AliMUONAlignment::LocalEquationX() {
730   /// Define local equation for current track and hit in x coor. measurement
731   // set local derivatives
732   SetLocalDerivative(0, fCosPhi);
733   SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
734   SetLocalDerivative(2, fSinPhi);
735   SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
736
737   // set global derivatives
738   SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -fCosPhi);
739   SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fSinPhi);
740   if (fBFieldOn){
741     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
742                         -fSinPhi*(fTrackPos[0]-fDetElemPos[0]) 
743                         +fCosPhi*(fTrackPos[1]-fDetElemPos[1]));
744   }
745   else {
746     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
747                         -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*
748                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0]) 
749                         +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*
750                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
751   }
752
753   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
754 }
755
756 void AliMUONAlignment::LocalEquationY() {
757   /// Define local equation for current track and hit in y coor. measurement
758   // set local derivatives
759   SetLocalDerivative(0,-fSinPhi);
760   SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
761   SetLocalDerivative(2, fCosPhi);
762   SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
763
764   // set global derivatives
765   SetGlobalDerivative(fDetElemNumber*fgNParCh+0,  fSinPhi);
766   SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fCosPhi);
767   if (fBFieldOn){
768   SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
769                       -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
770                       -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
771   }
772   else {
773     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
774                         -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*
775                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
776                         -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*
777                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
778   }
779
780   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
781 }
782
783 void AliMUONAlignment::FillRecPointData() {
784   /// Get information of current hit
785   fClustPos[0] = fRecHit->GetNonBendingCoor();
786   fClustPos[1] = fRecHit->GetBendingCoor();
787   fClustPos[2] = fRecHit->GetZ();
788   fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
789                             fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
790 }
791
792 void AliMUONAlignment::FillTrackParamData() {
793   /// Get information of current track at current hit
794   fTrackPos[0] = fTrackParam->GetNonBendingCoor();
795   fTrackPos[1] = fTrackParam->GetBendingCoor();
796   fTrackPos[2] = fTrackParam->GetZ();
797   fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
798   fTrackSlope[1] = fTrackParam->GetBendingSlope();
799   fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
800                             fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
801 }
802
803 void AliMUONAlignment::FillDetElemData() {
804   /// Get information of current detection element
805   Double_t lDetElemLocX = 0.;
806   Double_t lDetElemLocY = 0.;
807   Double_t lDetElemLocZ = 0.;
808   fDetElemId = fRecHit->GetDetElemId();
809   fDetElemNumber = fDetElemId%100;
810   for (int iCh=0; iCh<fDetElemId/100-1; iCh++){
811     fDetElemNumber += fgNDetElemCh[iCh];
812   }
813   fTransform->Local2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
814                            fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
815 }
816
817 void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) {
818   /// Process track; Loop over hits and set local equations
819   fTrack = track;
820   // get tclones arrays.
821   fTrackParamAtHit = fTrack->GetTrackParamAtHit();
822   fHitForRecAtHit = fTrack->GetHitForRecAtHit();
823   
824   // get size of arrays
825   Int_t nTrackParam = fTrackParamAtHit->GetEntries();
826   Int_t nHitForRec = fHitForRecAtHit->GetEntries();
827   AliDebug(1,Form("Number of track param entries : %i ", nTrackParam));
828   AliDebug(1,Form("Number of hit for rec entries : %i ", nHitForRec));
829
830   for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
831     fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
832     fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
833     if (!fRecHit || !fTrackParam) continue;
834     // fill local variables for this position --> one measurement
835     FillDetElemData();
836     FillRecPointData();
837     FillTrackParamData();         
838 //     if (fDetElemId<500) continue;
839     fTrackPos0[0]      = fTrackPos[0];    
840     fTrackPos0[1]      = fTrackPos[1];    
841     fTrackPos0[2]      = fTrackPos[2];    
842     fTrackSlope0[0] = fTrackSlope[0]; 
843     fTrackSlope0[1] = fTrackSlope[1];   
844     break;
845   }
846
847   for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
848     // and get new pointers
849     fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
850     fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
851     if (!fRecHit || !fTrackParam) continue;
852     // fill local variables for this position --> one measurement
853     FillDetElemData();        
854     FillRecPointData();
855     FillTrackParamData();
856 //     if (fDetElemId<500) continue;
857     AliDebug(1,Form("cluster: %i", iHit));
858     AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
859     AliDebug(1,Form("fDetElemPos[0]: %f\t fDetElemPos[1]: %f\t fDetElemPos[2]: %f\t DetElemID: %i\t ", fDetElemPos[0],fDetElemPos[1],fDetElemPos[2], fDetElemId));
860
861     AliDebug(1,Form("Track Parameter: %i", iHit));
862     AliDebug(1,Form("x: %f\t y: %f\t z: %f\t slopex: %f\t slopey: %f", fTrackPos[0], fTrackPos[1], fTrackPos[2], fTrackSlope[0], fTrackSlope[1]));
863     AliDebug(1,Form("x0: %f\t y0: %f\t z0: %f\t slopex0: %f\t slopey0: %f", fTrackPos0[0], fTrackPos0[1], fTrackPos0[2], fTrackSlope0[0], fTrackSlope0[1]));
864     
865     fCosPhi = TMath::Cos(fPhi);
866     fSinPhi = TMath::Sin(fPhi);
867     if (fBFieldOn){
868       fMeas[0] = fTrackPos[0] - fClustPos[0];
869       fMeas[1] = fTrackPos[1] - fClustPos[1];
870     }
871     else {
872       fMeas[0] = - fClustPos[0];
873       fMeas[1] = - fClustPos[1];
874     }
875     AliDebug(1,Form("fMeas[0]: %f\t fMeas[1]: %f\t fSigma[0]: %f\t fSigma[1]: %f", fMeas[0], fMeas[1], fSigma[0], fSigma[1]));    
876     // Set local equations
877     LocalEquationX();
878     LocalEquationY();
879   }
880 }
881
882 void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
883   /// Call local fit for this tracks
884   Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
885   if (iRes && !lSingleFit) {
886     fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
887   }
888 }
889
890 void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
891   /// Call global fit; Global parameters are stored in parameters
892   fMillepede->GlobalFit(parameters,errors,pulls);
893
894   AliInfo("Done fitting global parameters!");
895   for (int iGlob=0; iGlob<fgNDetElem; iGlob++){
896     printf("%d\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+2]);
897   }
898 }
899
900 Double_t AliMUONAlignment::GetParError(Int_t iPar) {
901   /// Get error of parameter iPar
902   Double_t lErr = fMillepede->GetParError(iPar);
903   return lErr;
904 }
905
906 void AliMUONAlignment::PrintGlobalParameters() {
907   /// Print global parameters
908   fMillepede->PrintGlobalParameters();
909 }
910
911 //_________________________________________________________________________
912 TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *lMisAlignment) const
913 {
914   /// Realign given transformation by given misalignment and return the misaligned transformation
915   
916   Double_t cartMisAlig[3] = {0,0,0};
917   Double_t angMisAlig[3] = {0,0,0};
918 //   const Double_t *trans = transform.GetTranslation();
919 //   TGeoRotation *rot;
920 //   // check if the rotation we obtain is not NULL
921 //   if (transform.GetRotation()) {    
922 //     rot = transform.GetRotation();
923 //   }
924 //   else {    
925 //     rot = new TGeoRotation("rot");
926 //   }                  // default constructor.
927
928   cartMisAlig[0] = -lMisAlignment[0];
929   cartMisAlig[1] = -lMisAlignment[1];
930   angMisAlig[2] = -lMisAlignment[2]*180./TMath::Pi();
931
932   TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
933   TGeoRotation deltaRot;
934   deltaRot.RotateX(angMisAlig[0]);
935   deltaRot.RotateY(angMisAlig[1]);
936   deltaRot.RotateZ(angMisAlig[2]);
937
938   TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
939   TGeoHMatrix newTransfMat = transform * deltaTransf;
940
941   return TGeoCombiTrans(newTransfMat);
942 }
943
944 //______________________________________________________________________
945 AliMUONGeometryTransformer *
946 AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
947                             double *misAlignments, Bool_t verbose)
948                             
949 {
950   /////////////////////////////////////////////////////////////////////
951   //   Takes the internal geometry module transformers, copies them
952   // and gets the Detection Elements from them.
953   // Takes misalignment parameters and applies these
954   // to the local transform of the Detection Element
955   // Obtains the global transform by multiplying the module transformer
956   // transformation with the local transformation 
957   // Applies the global transform to a new detection element
958   // Adds the new detection element to a new module transformer
959   // Adds the new module transformer to a new geometry transformer
960   // Returns the new geometry transformer
961
962   Double_t lModuleMisAlignment[3] = {0.,0.,0.};
963   Double_t lDetElemMisAlignment[3] = {0.,0.,0.};
964   Int_t iDetElemId = 0;
965   Int_t iDetElemNumber = 0;
966
967   AliMUONGeometryTransformer *newGeometryTransformer =
968     new AliMUONGeometryTransformer(kTRUE);
969   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
970     // module transformers    
971     const AliMUONGeometryModuleTransformer *kModuleTransformer =
972       transformer->GetModuleTransformer(iMt, true);
973       
974     AliMUONGeometryModuleTransformer *newModuleTransformer =
975       new AliMUONGeometryModuleTransformer(iMt);
976     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
977     
978     TGeoCombiTrans moduleTransform =
979       TGeoCombiTrans(*kModuleTransformer->GetTransformation());
980     // New module transformation
981     TGeoCombiTrans newModuleTransform = ReAlign(moduleTransform,lModuleMisAlignment);
982     newModuleTransformer->SetTransformation(newModuleTransform);
983     
984     // Get delta transformation: 
985     // Tdelta = Tnew * Told.inverse
986     TGeoHMatrix deltaModuleTransform = 
987       AliMUONGeometryBuilder::Multiply(newModuleTransform, 
988                                        kModuleTransformer->GetTransformation()->Inverse());    
989     // Create module mis alignment matrix
990     newGeometryTransformer
991       ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
992     
993     AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
994     
995     if (verbose)
996       AliInfo(Form("%i DEs in old GeometryStore  %i",detElements->GetSize(), iMt));
997
998     for (Int_t iDe = 0; iDe < detElements->GetSize(); iDe++) {
999       // detection elements.
1000       AliMUONGeometryDetElement *detElement =
1001         (AliMUONGeometryDetElement *) detElements->GetObject(iDe);
1002       if (!detElement)
1003         AliFatal("Detection element not found.");
1004
1005       /// make a new detection element
1006       AliMUONGeometryDetElement *newDetElement =
1007         new AliMUONGeometryDetElement(detElement->GetId(),
1008                                       detElement->GetVolumePath());
1009       TString lDetElemName(detElement->GetDEName());
1010       lDetElemName.ReplaceAll("DE","");
1011       iDetElemId = lDetElemName.Atoi();
1012       iDetElemNumber = iDetElemId%100;
1013       for (int iCh=0; iCh<iDetElemId/100-1; iCh++){
1014         iDetElemNumber += fgNDetElemCh[iCh];
1015       }
1016       for (int i=0; i<fgNParCh; i++) {
1017         lDetElemMisAlignment[i] = 0.0;
1018         if (iMt<fgNTrkMod) {
1019           AliInfo(Form("iMt %i, iCh %i, iDe %i, iDeId %i, iDeNb %i, iPar %i",iMt, iDetElemId/100, iDe, iDetElemId, iDetElemNumber, iDetElemNumber*fgNParCh+i));
1020           lDetElemMisAlignment[i] =  misAlignments[iDetElemNumber*fgNParCh+i];
1021         }
1022       }
1023       // local transformation of this detection element.
1024       TGeoCombiTrans localTransform
1025         = TGeoCombiTrans(*detElement->GetLocalTransformation());
1026       TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
1027       newDetElement->SetLocalTransformation(newLocalTransform);   
1028
1029       // global transformation
1030       TGeoHMatrix newGlobalTransform =
1031         AliMUONGeometryBuilder::Multiply(newModuleTransform,
1032                                          newLocalTransform);
1033       newDetElement->SetGlobalTransformation(newGlobalTransform);
1034           
1035       // add this det element to module
1036       newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
1037                                                       newDetElement);
1038
1039       // In the Alice Alignment Framework misalignment objects store
1040       // global delta transformation
1041       // Get detection "intermediate" global transformation
1042       TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform;
1043       // Get detection element global delta transformation: 
1044       // Tdelta = Tnew * Told.inverse
1045       TGeoHMatrix  deltaGlobalTransform
1046         = AliMUONGeometryBuilder::Multiply(newGlobalTransform, 
1047                                            newOldGlobalTransform.Inverse());
1048           
1049       // Create mis alignment matrix
1050       newGeometryTransformer
1051         ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
1052     }
1053       
1054     if (verbose)
1055       AliInfo(Form("Added module transformer %i to the transformer", iMt));
1056     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1057   }
1058   return newGeometryTransformer;
1059 }
1060