next50 trigger mask in AliHLTGlobalEsdConverterComponent
[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 AliMUONAlignment
20 /// Alignment class for 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 cluster and fill the corresponding local equations. Provide methods
25 /// for 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 "AliMUONTrackParam.h"
33 #include "AliMUONVCluster.h"
34 #include "AliMUONGeometryTransformer.h"
35 #include "AliMUONGeometryModuleTransformer.h"
36 #include "AliMUONGeometryDetElement.h"
37 #include "AliMUONGeometryBuilder.h"
38 #include "AliMillepede.h"
39
40 #include "AliMpExMap.h"
41 #include "AliMpExMapIterator.h"
42
43 #include "AliAlignObjMatrix.h"
44 #include "AliLog.h"
45
46 #include "TMath.h"
47 #include "TMatrixDSym.h"
48 #include "TClonesArray.h"
49
50 /// \cond CLASSIMP
51 ClassImp(AliMUONAlignment)
52 /// \endcond
53
54 //_____________________________________________________________________
55 // static variables
56 Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
57 Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
58 Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
59 Int_t AliMUONAlignment::fgNParCh = 4;
60 Int_t AliMUONAlignment::fgNTrkMod = 16;
61 Int_t AliMUONAlignment::fgNCh = 10;
62 Int_t AliMUONAlignment::fgNSt = 5;
63
64 //_____________________________________________________________________
65 AliMUONAlignment::AliMUONAlignment()
66   : TObject(),
67     fBFieldOn(kTRUE),
68     fStartFac(256.),
69     fResCutInitial(100.),
70     fResCut(100.),
71     fMillepede(0),
72     fTrack(0),
73     fCluster(0),
74     fTrackParam(0),
75     fNGlobal(fgNDetElem*fgNParCh),
76     fNLocal(4),
77     fNStdDev(3),
78     fDetElemId(0),
79     fDetElemNumber(0),
80     fPhi(0.),
81     fCosPhi(1.),
82     fSinPhi(0.),
83     fTrackRecord(),
84     fTransform(0)
85 {
86   /// Default constructor       
87   fSigma[0] = 1.5e-1;
88   fSigma[1] = 1.0e-2;
89
90   AliInfo(Form("fSigma[0]: %f\t fSigma[1]: %f",fSigma[0],fSigma[1]));
91
92   fDoF[0] = kTRUE;  fDoF[1] = kTRUE;  fDoF[2] = kTRUE;  fDoF[3] = kTRUE;
93   fAllowVar[0] = 0.05;  fAllowVar[1] = 0.05;  fAllowVar[2] = 0.001;  fAllowVar[3] = 0.5;
94
95   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));
96
97   fMillepede = new AliMillepede();
98
99   Init(fNGlobal, fNLocal, fNStdDev);
100
101   ResetLocalEquation();
102   AliInfo("Parameters initialized to zero");
103
104 }
105
106 //_____________________________________________________________________
107 AliMUONAlignment::AliMUONAlignment(TRootIOCtor* /*dummy*/)
108 : TObject(),
109 fBFieldOn(kFALSE),
110 fStartFac(0.),
111 fResCutInitial(0.),
112 fResCut(0.),
113 fMillepede(0),
114 fTrack(0),
115 fCluster(0),
116 fTrackParam(0),
117 fNGlobal(0),
118 fNLocal(0),
119 fNStdDev(0),
120 fDetElemId(0),
121 fDetElemNumber(0),
122 fPhi(0.),
123 fCosPhi(0.),
124 fSinPhi(0.),
125 fTrackRecord(),
126 fTransform(0)
127 {
128   /// Root IO constructor
129   ResetConstraints();
130   for (Int_t iCh=0; iCh<10; iCh++) { 
131     fChOnOff[iCh] = kFALSE; 
132   }
133   fSpecLROnOff[0] = kFALSE; fSpecLROnOff[1] = kFALSE;
134   for (Int_t iDoF=0; iDoF<4; iDoF++) {
135     fDoF[iDoF] = kFALSE;
136     fAllowVar[iDoF] = kFALSE;
137   }
138   for (Int_t i=0; i<3; i++) {
139     fClustPos[i] = 0;
140     fClustPosLoc[i] = 0;
141     fTrackPos0[i] = 0;  
142     fTrackPos[i] = 0;   
143     fTrackPosLoc[i] = 0;
144         fDetElemPos[i] = 0;
145   }
146   fTrackSlope0[0]=0;    fTrackSlope0[1]=0;
147   fTrackSlope[0]=0;     fTrackSlope[1]=0;
148   fMeas[0]=0;   fMeas[1]=0;
149   fSigma[0]=0;  fSigma[1]=0;
150   for (Int_t iLPar=0; iLPar<4; iLPar++) {
151     fLocalDerivatives[iLPar]=0; 
152   }
153   for (Int_t iGPar=0; iGPar<624; iGPar++) {
154     fGlobalDerivatives[iGPar]=0;        
155   }
156 }
157
158 //_____________________________________________________________________
159 AliMUONAlignment::~AliMUONAlignment()
160 {}
161
162 //_____________________________________________________________________
163 void AliMUONAlignment::Init(
164   Int_t nGlobal,  /* number of global paramers */
165   Int_t nLocal,   /* number of local parameters */
166   Int_t nStdDev   /* std dev cut */ )
167 {
168   /// Initialization of AliMillepede. Fix parameters, define constraints ...
169   fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
170
171   //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
172   //  Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
173   //  Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
174
175   //  AllowVariations(bChOnOff);
176
177   // Fix parameters or add constraints here
178   //   for (Int_t iSt=0; iSt<5; iSt++)
179   //   { if (!bStOnOff[iSt]) FixStation(iSt+1); }
180
181   //   for (Int_t iCh=0; iCh<10; iCh++)
182   //   { if (!bChOnOff[iCh]) FixChamber(iCh+1); }
183
184 //   FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
185
186   ResetConstraints();
187
188   // Define global constrains to be applied
189   // X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
190   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
191   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
192   //  AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
193
194   // Other possible way to add constrains
195   bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
196   bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
197   //  AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
198
199   bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE;
200   //  AddConstraints(bStOnOff,bVarXYT);
201
202   // Set iterations
203   if (fStartFac>1) fMillepede->SetIterations(fStartFac);
204 }
205
206 //_____________________________________________________________________
207 void AliMUONAlignment::FixStation(Int_t iSt)
208 {
209   /// Fix all detection elements of station iSt
210   Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
211   Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
212   for (Int_t i = iDetElemFirst; i < iDetElemLast; i++)
213   {
214     FixParameter(i*fgNParCh+0, 0.0);
215     FixParameter(i*fgNParCh+1, 0.0);
216     FixParameter(i*fgNParCh+2, 0.0);
217     FixParameter(i*fgNParCh+3, 0.0);
218   }
219
220 }
221
222 //_____________________________________________________________________
223 void AliMUONAlignment::FixChamber(Int_t iCh)
224 {
225   /// Fix all detection elements of chamber iCh
226   Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0;
227   Int_t iDetElemLast = fgSNDetElemCh[iCh-1];
228   for (Int_t i = iDetElemFirst; i < iDetElemLast; i++)
229   {
230     FixParameter(i*fgNParCh+0, 0.0);
231     FixParameter(i*fgNParCh+1, 0.0);
232     FixParameter(i*fgNParCh+2, 0.0);
233     FixParameter(i*fgNParCh+3, 0.0);
234   }
235 }
236
237 //_____________________________________________________________________
238 void AliMUONAlignment::FixDetElem(Int_t iDetElemId, TString sVarXYT)
239 {
240
241   /// Fix a given detection element
242   Int_t iDetElemNumber = iDetElemId%100;
243   for (int iCh=0; iCh<iDetElemId/100-1; iCh++)
244   {
245     iDetElemNumber += fgNDetElemCh[iCh];
246   }
247
248   if (sVarXYT.Contains("X"))
249   {
250     // X constraint
251     FixParameter(iDetElemNumber*fgNParCh+0, 0.0);
252   }
253
254   if (sVarXYT.Contains("Y"))
255   {
256     // Y constraint
257     FixParameter(iDetElemNumber*fgNParCh+1, 0.0);
258   }
259
260   if (sVarXYT.Contains("T"))
261   {
262     // T constraint
263     FixParameter(iDetElemNumber*fgNParCh+2, 0.0);
264   }
265
266   if (sVarXYT.Contains("Z"))
267   {
268     // T constraint
269     FixParameter(iDetElemNumber*fgNParCh+3, 0.0);
270   }
271
272 }
273
274 //_____________________________________________________________________
275 void AliMUONAlignment::FixHalfSpectrometer(const Bool_t *lChOnOff, const Bool_t *lSpecLROnOff)
276 {
277
278   /// Fix left or right detector
279   for (Int_t i = 0; i < fgNDetElem; i++)
280   {
281
282     Int_t iCh=0;
283     for (iCh=1; iCh<=fgNCh; iCh++)
284     { if (i<fgSNDetElemCh[iCh-1]) break; }
285
286     if (lChOnOff[iCh-1])
287     {
288
289       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
290       if (iCh>=1 && iCh<=4)
291       {
292
293         if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0])
294         {
295           // From track crossings
296           FixParameter(i*fgNParCh+0, 0.0);
297           FixParameter(i*fgNParCh+1, 0.0);
298           FixParameter(i*fgNParCh+2, 0.0);
299           FixParameter(i*fgNParCh+3, 0.0);
300         }
301
302         if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1])
303         {
304           // From track crossings
305           FixParameter(i*fgNParCh+0, 0.0);
306           FixParameter(i*fgNParCh+1, 0.0);
307           FixParameter(i*fgNParCh+2, 0.0);
308           FixParameter(i*fgNParCh+3, 0.0);
309         }
310
311       }
312
313       if (iCh>=5 && iCh<=6)
314       {
315
316         if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0])
317         {
318           FixParameter(i*fgNParCh+0, 0.0);
319           FixParameter(i*fgNParCh+1, 0.0);
320           FixParameter(i*fgNParCh+2, 0.0);
321           FixParameter(i*fgNParCh+3, 0.0);
322         }
323
324         if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
325           (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1])
326         {
327
328           FixParameter(i*fgNParCh+0, 0.0);
329           FixParameter(i*fgNParCh+1, 0.0);
330           FixParameter(i*fgNParCh+2, 0.0);
331           FixParameter(i*fgNParCh+3, 0.0);
332         }
333
334       }
335
336       if (iCh>=7 && iCh<=10)
337       {
338
339         if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0])
340         {
341           FixParameter(i*fgNParCh+0, 0.0);
342           FixParameter(i*fgNParCh+1, 0.0);
343           FixParameter(i*fgNParCh+2, 0.0);
344           FixParameter(i*fgNParCh+3, 0.0);
345         }
346
347         if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
348           (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1])
349         {
350           FixParameter(i*fgNParCh+0, 0.0);
351           FixParameter(i*fgNParCh+1, 0.0);
352           FixParameter(i*fgNParCh+2, 0.0);
353           FixParameter(i*fgNParCh+3, 0.0);
354         }
355
356       }
357
358     }
359
360   }
361
362 }
363
364 //______________________________________________________________________
365 void AliMUONAlignment::SetNonLinear(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
366 {
367
368   /// Set non linear parameter flag selected chambers and degrees of freedom
369   for (Int_t i = 0; i < fgNDetElem; i++)
370   {
371
372     Int_t iCh=0;
373     for (iCh=1; iCh<=fgNCh; iCh++)
374     { if (i<fgSNDetElemCh[iCh-1]) break; }
375
376     if (lChOnOff[iCh-1])
377     {
378
379       if (lVarXYT[0])
380       {
381         // X constraint
382         SetNonLinear(i*fgNParCh+0);
383       }
384
385       if (lVarXYT[1])
386       {
387         // Y constraint
388         SetNonLinear(i*fgNParCh+1);
389       }
390
391       if (lVarXYT[2])
392       {
393         // T constraint
394         SetNonLinear(i*fgNParCh+2);
395       }
396
397       if (lVarXYT[3])
398       {
399         // Z constraint
400         SetNonLinear(i*fgNParCh+3);
401       }
402
403     }
404
405   }
406
407 }
408
409 //______________________________________________________________________
410 void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
411 {
412
413   /// Add constraint equations for selected chambers and degrees of freedom
414   for (Int_t i = 0; i < fgNDetElem; i++)
415   {
416
417     Int_t iCh=0;
418     for (iCh=1; iCh<=fgNCh; iCh++)
419     {
420       if (i<fgSNDetElemCh[iCh-1]) break;
421     }
422
423     if (lChOnOff[iCh-1])
424     {
425       if (lVarXYT[0])
426       {
427         // X constraint
428         fConstraintX[i*fgNParCh+0]=1.0;
429       }
430
431       if (lVarXYT[1])
432       {
433         // Y constraint
434         fConstraintY[i*fgNParCh+1]=1.0;
435       }
436
437       if (lVarXYT[2])
438       {
439         // T constraint
440         fConstraintP[i*fgNParCh+2]=1.0;
441       }
442       //       if (lVarXYT[3]) { // Z constraint
443       //        fConstraintP[i*fgNParCh+3]=1.0;
444       //       }
445
446     }
447   }
448
449   if (lVarXYT[0])
450   {
451     // X constraint
452     AddConstraint(fConstraintX,0.0);
453   }
454
455   if (lVarXYT[1])
456   {
457     // Y constraint
458     AddConstraint(fConstraintY,0.0);
459   }
460
461   if (lVarXYT[2])
462   {
463     // T constraint
464     AddConstraint(fConstraintP,0.0);
465   }
466
467 //   if (lVarXYT[3]) { // Z constraint
468 //     AddConstraint(fConstraintP,0.0);
469 //   }
470 }
471
472 //______________________________________________________________________
473 void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, const Bool_t *lDetTLBR, const Bool_t *lSpecLROnOff)
474 {
475   /// Add constraint equations for selected chambers, degrees of freedom and detector half
476   Double_t lDetElemLocX = 0.;
477   Double_t lDetElemLocY = 0.;
478   Double_t lDetElemLocZ = 0.;
479   Double_t lDetElemGloX = 0.;
480   Double_t lDetElemGloY = 0.;
481   Double_t lDetElemGloZ = 0.;
482   Double_t lMeanY = 0.;
483   Double_t lSigmaY = 0.;
484   Double_t lMeanZ = 0.;
485   Double_t lSigmaZ = 0.;
486   Int_t lNDetElem = 0;
487   for (Int_t i = 0; i < fgNDetElem; i++)
488   {
489
490     Int_t iCh=0;
491     for (iCh=1; iCh<=fgNCh; iCh++){
492       if (i<fgSNDetElemCh[iCh-1]) break;
493     }
494     if (lChOnOff[iCh-1]){
495       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
496       Int_t lDetElemId = iCh*100+lDetElemNumber;
497       fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
498              lDetElemGloX,lDetElemGloY,lDetElemGloZ);
499       if (iCh>=1 && iCh<=4){
500   if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings
501     lMeanY += lDetElemGloY;
502     lSigmaY += lDetElemGloY*lDetElemGloY;
503     lMeanZ += lDetElemGloZ;
504     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
505     lNDetElem++;
506   }
507   if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings
508     lMeanY += lDetElemGloY;
509     lSigmaY += lDetElemGloY*lDetElemGloY;
510     lMeanZ += lDetElemGloZ;
511     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
512     lNDetElem++;
513   }
514       }
515       if (iCh>=5 && iCh<=6){
516   if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){
517     lMeanY += lDetElemGloY;
518     lSigmaY += lDetElemGloY*lDetElemGloY;
519     lMeanZ += lDetElemGloZ;
520     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
521     lNDetElem++;
522   }
523   if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
524        (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){
525     lMeanY += lDetElemGloY;
526     lSigmaY += lDetElemGloY*lDetElemGloY;
527     lMeanZ += lDetElemGloZ;
528     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
529     lNDetElem++;
530   }
531       }
532       if (iCh>=7 && iCh<=10){
533   if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){
534     lMeanY += lDetElemGloY;
535     lSigmaY += lDetElemGloY*lDetElemGloY;
536     lMeanZ += lDetElemGloZ;
537     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
538     lNDetElem++;
539   }
540   if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
541        (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){
542     lMeanY += lDetElemGloY;
543     lSigmaY += lDetElemGloY*lDetElemGloY;
544     lMeanZ += lDetElemGloZ;
545     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
546     lNDetElem++;
547   }
548       }
549     }
550   }
551   if (lNDetElem) {
552     lMeanY /= lNDetElem;
553     lSigmaY /= lNDetElem;
554     lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
555     lMeanZ /= lNDetElem;
556     lSigmaZ /= lNDetElem;
557     lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
558      AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ));
559   } else {
560     AliError("No detection elements to constrain!!!");
561     return;    
562   }
563                 
564   for (Int_t i = 0; i < fgNDetElem; i++){
565     Int_t iCh=0;
566     for (iCh=1; iCh<=fgNCh; iCh++){
567       if (i<fgSNDetElemCh[iCh-1]) break;
568     }
569     if (lChOnOff[iCh-1]){
570       Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
571       Int_t lDetElemId = iCh*100+lDetElemNumber;
572       fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
573              lDetElemGloX,lDetElemGloY,lDetElemGloZ);
574       if (lVarXYT[0]) { // X constraint
575   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
576   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
577   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
578   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
579       }
580       if (lVarXYT[1]) { // Y constraint
581   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
582   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
583   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
584   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
585       }
586       if (lVarXYT[2]) { // P constraint
587   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
588   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
589   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
590   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
591       }
592       if (lVarXYT[3]) { // X-Z shearing
593   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
594   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
595   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
596   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
597       }
598       if (lVarXYT[4]) { // Y-Z shearing
599   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
600   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
601   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
602   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
603       }
604       if (lVarXYT[5]) { // P-Z rotation
605   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
606   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
607   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
608   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
609       }
610       if (lVarXYT[6]) { // X-Y shearing
611   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
612   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
613   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
614   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
615       }
616       if (lVarXYT[7]) { // Y-Y scaling
617   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
618   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
619   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
620   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
621       }
622       if (lVarXYT[8]) { // P-Y rotation
623   if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
624   if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
625   if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
626   if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
627       }
628     }
629   }
630   if (lVarXYT[0]) { // X constraint
631     if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
632     if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
633     if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
634     if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
635   }
636   if (lVarXYT[1]) { // Y constraint
637     if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
638     if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
639     if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
640     if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
641   }
642   if (lVarXYT[2]) { // T constraint
643     if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
644     if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
645     if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
646     if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
647   }
648   if (lVarXYT[3]) { // X-Z constraint
649     if (lDetTLBR[0]) AddConstraint(fConstraintXZT,0.0); // Top half
650     if (lDetTLBR[1]) AddConstraint(fConstraintXZL,0.0); // Left half
651     if (lDetTLBR[2]) AddConstraint(fConstraintXZB,0.0); // Bottom half
652     if (lDetTLBR[3]) AddConstraint(fConstraintXZR,0.0); // Right half
653   }
654   if (lVarXYT[4]) { // Y-Z constraint
655     if (lDetTLBR[0]) AddConstraint(fConstraintYZT,0.0); // Top half
656     if (lDetTLBR[1]) AddConstraint(fConstraintYZL,0.0); // Left half
657     if (lDetTLBR[2]) AddConstraint(fConstraintYZB,0.0); // Bottom half
658     if (lDetTLBR[3]) AddConstraint(fConstraintYZR,0.0); // Right half
659   }
660   if (lVarXYT[5]) { // P-Z constraint
661     if (lDetTLBR[0]) AddConstraint(fConstraintPZT,0.0); // Top half
662     if (lDetTLBR[1]) AddConstraint(fConstraintPZL,0.0); // Left half
663     if (lDetTLBR[2]) AddConstraint(fConstraintPZB,0.0); // Bottom half
664     if (lDetTLBR[3]) AddConstraint(fConstraintPZR,0.0); // Right half
665   }
666   if (lVarXYT[6]) { // X-Y constraint
667     if (lDetTLBR[0]) AddConstraint(fConstraintXYT,0.0); // Top half
668     if (lDetTLBR[1]) AddConstraint(fConstraintXYL,0.0); // Left half
669     if (lDetTLBR[2]) AddConstraint(fConstraintXYB,0.0); // Bottom half
670     if (lDetTLBR[3]) AddConstraint(fConstraintXYR,0.0); // Right half
671   }
672   if (lVarXYT[7]) { // Y-Y constraint
673     if (lDetTLBR[0]) AddConstraint(fConstraintYYT,0.0); // Top half
674     if (lDetTLBR[1]) AddConstraint(fConstraintYYL,0.0); // Left half
675     if (lDetTLBR[2]) AddConstraint(fConstraintYYB,0.0); // Bottom half
676     if (lDetTLBR[3]) AddConstraint(fConstraintYYR,0.0); // Right half
677   }
678   if (lVarXYT[8]) { // P-Y constraint
679     if (lDetTLBR[0]) AddConstraint(fConstraintPYT,0.0); // Top half
680     if (lDetTLBR[1]) AddConstraint(fConstraintPYL,0.0); // Left half
681     if (lDetTLBR[2]) AddConstraint(fConstraintPYB,0.0); // Bottom half
682     if (lDetTLBR[3]) AddConstraint(fConstraintPYR,0.0); // Right half
683   }
684 }
685
686 void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t /*lWeight*/) const{
687   /// Set constrain equation for top half of spectrometer
688   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
689   if (lCh>=1 && lCh<=4){
690     if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
691       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
692     }
693   }
694   if (lCh>=5 && lCh<=6){
695     if (lDetElemNumber>=0&&lDetElemNumber<=9){
696       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
697     }
698   }
699   if (lCh>=7 && lCh<=10){
700     if (lDetElemNumber>=0&&lDetElemNumber<=13){
701       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
702     }
703   }
704 }
705
706 //______________________________________________________________________
707 void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight) const{
708   /// Set constrain equation for left half of spectrometer
709   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
710   if (lCh>=1 && lCh<=4){
711     if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
712       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
713     }
714   }
715   if (lCh>=5 && lCh<=6){
716     if (lDetElemNumber>=5&&lDetElemNumber<=13){
717       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
718     }
719   }
720   if (lCh>=7 && lCh<=10){
721     if (lDetElemNumber>=7&&lDetElemNumber<=19){
722       lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
723     }
724   }
725 }
726
727 //______________________________________________________________________
728 void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t /*lWeight*/) const{
729   /// Set constrain equation for bottom half of spectrometer
730   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
731   if (lCh>=1 && lCh<=4){
732     if (lDetElemNumber==2 || lDetElemNumber==3){ // From track crossings
733       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
734     }
735   }
736   if (lCh>=5 && lCh<=6){
737     if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
738   (lDetElemNumber==0)){
739       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
740     }
741   }
742   if (lCh>=7 && lCh<=10){
743     if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
744   (lDetElemNumber==0)){
745       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
746     }
747   }
748 }
749
750 //______________________________________________________________________
751 void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight) const{
752   /// Set constrain equation for right half of spectrometer
753   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
754   if (lCh>=1 && lCh<=4){
755     if (lDetElemNumber==0 || lDetElemNumber==3){ // From track crossings
756       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
757     }
758   }
759   if (lCh>=5 && lCh<=6){
760     if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
761   (lDetElemNumber>=14&&lDetElemNumber<=17)){
762       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
763     }
764   }
765   if (lCh>=7 && lCh<=10){
766     if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
767   (lDetElemNumber>=20&&lDetElemNumber<=25)){
768       lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
769     }
770   }
771 }
772
773 //______________________________________________________________________
774 void AliMUONAlignment::ResetConstraints(){
775   /// Reset all constraint equations
776   for (Int_t i = 0; i < fgNDetElem; i++){
777         fConstraintX3[i*fgNParCh+0]=0.0; 
778         fConstraintY3[i*fgNParCh+0]=0.0; 
779         fConstraintX4[i*fgNParCh+0]=0.0; 
780         fConstraintY4[i*fgNParCh+0]=0.0; 
781         fConstraintP4[i*fgNParCh+0]=0.0; 
782         fConstraintX5[i*fgNParCh+0]=0.0; 
783         fConstraintY5[i*fgNParCh+0]=0.0; 
784         fConstraintX[i*fgNParCh+0]=0.0;
785     fConstraintX[i*fgNParCh+1]=0.0;
786     fConstraintX[i*fgNParCh+2]=0.0;
787     fConstraintY[i*fgNParCh+0]=0.0;
788     fConstraintY[i*fgNParCh+1]=0.0;
789     fConstraintY[i*fgNParCh+2]=0.0;
790     fConstraintP[i*fgNParCh+0]=0.0;
791     fConstraintP[i*fgNParCh+1]=0.0;
792     fConstraintP[i*fgNParCh+2]=0.0;
793     fConstraintXT[i*fgNParCh+0]=0.0;
794     fConstraintXT[i*fgNParCh+1]=0.0;
795     fConstraintXT[i*fgNParCh+2]=0.0;
796     fConstraintYT[i*fgNParCh+0]=0.0;
797     fConstraintYT[i*fgNParCh+1]=0.0;
798     fConstraintYT[i*fgNParCh+2]=0.0;
799     fConstraintPT[i*fgNParCh+0]=0.0;
800     fConstraintPT[i*fgNParCh+1]=0.0;
801     fConstraintPT[i*fgNParCh+2]=0.0;
802     fConstraintXZT[i*fgNParCh+0]=0.0;
803     fConstraintXZT[i*fgNParCh+1]=0.0;
804     fConstraintXZT[i*fgNParCh+2]=0.0;
805     fConstraintYZT[i*fgNParCh+0]=0.0;
806     fConstraintYZT[i*fgNParCh+1]=0.0;
807     fConstraintYZT[i*fgNParCh+2]=0.0;
808     fConstraintPZT[i*fgNParCh+0]=0.0;
809     fConstraintPZT[i*fgNParCh+1]=0.0;
810     fConstraintPZT[i*fgNParCh+2]=0.0;
811     fConstraintXYT[i*fgNParCh+0]=0.0;
812     fConstraintXYT[i*fgNParCh+1]=0.0;
813     fConstraintXYT[i*fgNParCh+2]=0.0;
814     fConstraintYYT[i*fgNParCh+0]=0.0;
815     fConstraintYYT[i*fgNParCh+1]=0.0;
816     fConstraintYYT[i*fgNParCh+2]=0.0;
817     fConstraintPYT[i*fgNParCh+0]=0.0;
818     fConstraintPYT[i*fgNParCh+1]=0.0;
819     fConstraintPYT[i*fgNParCh+2]=0.0;
820     fConstraintXL[i*fgNParCh+0]=0.0;
821     fConstraintXL[i*fgNParCh+1]=0.0;
822     fConstraintXL[i*fgNParCh+2]=0.0;
823     fConstraintYL[i*fgNParCh+0]=0.0;
824     fConstraintYL[i*fgNParCh+1]=0.0;
825     fConstraintYL[i*fgNParCh+2]=0.0;
826     fConstraintPL[i*fgNParCh+0]=0.0;
827     fConstraintPL[i*fgNParCh+1]=0.0;
828     fConstraintPL[i*fgNParCh+2]=0.0;
829     fConstraintXZL[i*fgNParCh+0]=0.0;
830     fConstraintXZL[i*fgNParCh+1]=0.0;
831     fConstraintXZL[i*fgNParCh+2]=0.0;
832     fConstraintYZL[i*fgNParCh+0]=0.0;
833     fConstraintYZL[i*fgNParCh+1]=0.0;
834     fConstraintYZL[i*fgNParCh+2]=0.0;
835     fConstraintPZL[i*fgNParCh+0]=0.0;
836     fConstraintPZL[i*fgNParCh+1]=0.0;
837     fConstraintPZL[i*fgNParCh+2]=0.0;
838     fConstraintXYL[i*fgNParCh+0]=0.0;
839     fConstraintXYL[i*fgNParCh+1]=0.0;
840     fConstraintXYL[i*fgNParCh+2]=0.0;
841     fConstraintYYL[i*fgNParCh+0]=0.0;
842     fConstraintYYL[i*fgNParCh+1]=0.0;
843     fConstraintYYL[i*fgNParCh+2]=0.0;
844     fConstraintPYL[i*fgNParCh+0]=0.0;
845     fConstraintPYL[i*fgNParCh+1]=0.0;
846     fConstraintPYL[i*fgNParCh+2]=0.0;
847     fConstraintXB[i*fgNParCh+0]=0.0;
848     fConstraintXB[i*fgNParCh+1]=0.0;
849     fConstraintXB[i*fgNParCh+2]=0.0;
850     fConstraintYB[i*fgNParCh+0]=0.0;
851     fConstraintYB[i*fgNParCh+1]=0.0;
852     fConstraintYB[i*fgNParCh+2]=0.0;
853     fConstraintPB[i*fgNParCh+0]=0.0;
854     fConstraintPB[i*fgNParCh+1]=0.0;
855     fConstraintPB[i*fgNParCh+2]=0.0;
856     fConstraintXZB[i*fgNParCh+0]=0.0;
857     fConstraintXZB[i*fgNParCh+1]=0.0;
858     fConstraintXZB[i*fgNParCh+2]=0.0;
859     fConstraintYZB[i*fgNParCh+0]=0.0;
860     fConstraintYZB[i*fgNParCh+1]=0.0;
861     fConstraintYZB[i*fgNParCh+2]=0.0;
862     fConstraintPZB[i*fgNParCh+0]=0.0;
863     fConstraintPZB[i*fgNParCh+1]=0.0;
864     fConstraintPZB[i*fgNParCh+2]=0.0;
865     fConstraintXYB[i*fgNParCh+0]=0.0;
866     fConstraintXYB[i*fgNParCh+1]=0.0;
867     fConstraintXYB[i*fgNParCh+2]=0.0;
868     fConstraintYYB[i*fgNParCh+0]=0.0;
869     fConstraintYYB[i*fgNParCh+1]=0.0;
870     fConstraintYYB[i*fgNParCh+2]=0.0;
871     fConstraintPYB[i*fgNParCh+0]=0.0;
872     fConstraintPYB[i*fgNParCh+1]=0.0;
873     fConstraintPYB[i*fgNParCh+2]=0.0;
874     fConstraintXR[i*fgNParCh+0]=0.0;
875     fConstraintXR[i*fgNParCh+1]=0.0;
876     fConstraintXR[i*fgNParCh+2]=0.0;
877     fConstraintYR[i*fgNParCh+0]=0.0;
878     fConstraintYR[i*fgNParCh+1]=0.0;
879     fConstraintYR[i*fgNParCh+2]=0.0;
880     fConstraintPR[i*fgNParCh+0]=0.0;
881     fConstraintPR[i*fgNParCh+1]=0.0;
882     fConstraintPR[i*fgNParCh+2]=0.0;
883     fConstraintXZR[i*fgNParCh+0]=0.0;
884     fConstraintXZR[i*fgNParCh+1]=0.0;
885     fConstraintXZR[i*fgNParCh+2]=0.0;
886     fConstraintYZR[i*fgNParCh+0]=0.0;
887     fConstraintYZR[i*fgNParCh+1]=0.0;
888     fConstraintYZR[i*fgNParCh+2]=0.0;
889     fConstraintPZR[i*fgNParCh+0]=0.0;
890     fConstraintPZR[i*fgNParCh+1]=0.0;
891     fConstraintPZR[i*fgNParCh+2]=0.0;
892     fConstraintPZR[i*fgNParCh+0]=0.0;
893     fConstraintPZR[i*fgNParCh+1]=0.0;
894     fConstraintPZR[i*fgNParCh+2]=0.0;
895     fConstraintXYR[i*fgNParCh+0]=0.0;
896     fConstraintXYR[i*fgNParCh+1]=0.0;
897     fConstraintXYR[i*fgNParCh+2]=0.0;
898     fConstraintYYR[i*fgNParCh+0]=0.0;
899     fConstraintYYR[i*fgNParCh+1]=0.0;
900     fConstraintYYR[i*fgNParCh+2]=0.0;
901     fConstraintPYR[i*fgNParCh+0]=0.0;
902     fConstraintPYR[i*fgNParCh+1]=0.0;
903     fConstraintPYR[i*fgNParCh+2]=0.0;
904   }
905 }
906
907 //______________________________________________________________________
908 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
909   /// Constrain equation defined by par to value
910   fMillepede->SetGlobalConstraint(par, value);
911   AliInfo("Adding constraint");
912 }
913
914 //______________________________________________________________________
915 void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
916   /// Initialize global parameters with par array
917   fMillepede->SetGlobalParameters(par);
918   AliInfo("Init Global Parameters");
919 }
920
921 //______________________________________________________________________
922 void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
923   /// Parameter iPar is encourage to vary in [-value;value].
924   /// If value == 0, parameter is fixed
925   fMillepede->SetParSigma(iPar, value);
926   if (TMath::Abs(value)<1e-4) AliInfo(Form("Parameter %i Fixed", iPar));
927 }
928
929 //______________________________________________________________________
930 void AliMUONAlignment::ResetLocalEquation()
931 {
932   /// Reset the derivative vectors
933   for(int i=0; i<fNLocal; i++) {
934     fLocalDerivatives[i] = 0.0;
935   }
936   for(int i=0; i<fNGlobal; i++) {
937     fGlobalDerivatives[i] = 0.0;
938   }
939 }
940
941 //______________________________________________________________________
942 void AliMUONAlignment::AllowVariations(const Bool_t *bChOnOff)
943 {
944
945   /// Set allowed variation for selected chambers based on fDoF and fAllowVar
946   for (Int_t iCh=1; iCh<=10; iCh++)
947   {
948     if (bChOnOff[iCh-1])
949     {
950
951       Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0;
952       Int_t iDetElemLast = fgSNDetElemCh[iCh-1];
953       for (int i=0; i<fgNParCh; i++)
954       {
955         AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
956         if (fDoF[i])
957         {
958           for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
959             FixParameter(j*fgNParCh+i, fAllowVar[i]);
960           }
961         }
962
963       }
964
965     }
966
967   }
968
969 }
970
971 //______________________________________________________________________
972 void AliMUONAlignment::SetNonLinear(Int_t iPar  /* set non linear flag */ )
973 {
974   /// Set nonlinear flag for parameter iPar
975   fMillepede->SetNonLinear(iPar);
976   AliInfo(Form("Parameter %i set to non linear", iPar));
977 }
978
979 //______________________________________________________________________
980 void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
981 {
982
983   /// Set expected measurement resolution
984   fSigma[0] = sigmaX;   fSigma[1] = sigmaY;
985   AliInfo(Form("Using fSigma[0]=%f and fSigma[1]=%f",fSigma[0],fSigma[1]));
986 }
987
988
989 //______________________________________________________________________
990 void AliMUONAlignment::LocalEquationX( Bool_t doAlignment )
991 {
992
993
994   // local cluster record
995   AliMUONAlignmentClusterRecord clusterRecord;
996
997   // store detector and measurement
998   clusterRecord.SetDetElemId( fDetElemId );
999   clusterRecord.SetDetElemNumber( fDetElemNumber );
1000   clusterRecord.SetMeas( fMeas[0] );
1001   clusterRecord.SetSigma( fSigma[0] );
1002
1003   // store local derivatives
1004   clusterRecord.SetLocalDerivative( 0, fCosPhi );
1005   clusterRecord.SetLocalDerivative( 1, fCosPhi*(fTrackPos[2] - fTrackPos0[2]) );
1006   clusterRecord.SetLocalDerivative( 2, fSinPhi );
1007   clusterRecord.SetLocalDerivative( 3, fSinPhi*(fTrackPos[2] - fTrackPos0[2]) );
1008
1009   // store global derivatives
1010   clusterRecord.SetGlobalDerivative( 0, -fCosPhi );
1011   clusterRecord.SetGlobalDerivative( 1, -fSinPhi );
1012
1013   if (fBFieldOn)
1014   {
1015
1016     clusterRecord.SetGlobalDerivative(
1017       2,
1018       -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
1019       +fCosPhi*(fTrackPos[1]-fDetElemPos[1]) );
1020
1021   } else {
1022
1023     clusterRecord.SetGlobalDerivative(
1024       2,
1025       -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
1026       +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]) );
1027
1028   }
1029
1030   clusterRecord.SetGlobalDerivative( 3, fCosPhi*fTrackSlope0[0]+fSinPhi*fTrackSlope0[1] );
1031
1032   // append to trackRecord
1033   fTrackRecord.AddClusterRecord( clusterRecord );
1034
1035   // store local equation
1036   if( doAlignment ) LocalEquation( clusterRecord );
1037
1038 }
1039
1040 //______________________________________________________________________
1041 void AliMUONAlignment::LocalEquationY(Bool_t doAlignment )
1042 {
1043
1044   // local cluster record
1045   AliMUONAlignmentClusterRecord clusterRecord;
1046
1047   // store detector and measurement
1048   clusterRecord.SetDetElemId( fDetElemId );
1049   clusterRecord.SetDetElemNumber( fDetElemNumber );
1050   clusterRecord.SetMeas( fMeas[1] );
1051   clusterRecord.SetSigma( fSigma[1] );
1052
1053   // store local derivatives
1054   clusterRecord.SetLocalDerivative( 0, -fSinPhi );
1055   clusterRecord.SetLocalDerivative( 1, -fSinPhi*(fTrackPos[2] - fTrackPos0[2] ) );
1056   clusterRecord.SetLocalDerivative( 2, fCosPhi );
1057   clusterRecord.SetLocalDerivative( 3, fCosPhi*(fTrackPos[2] - fTrackPos0[2] ) );
1058
1059   // set global derivatives
1060   clusterRecord.SetGlobalDerivative( 0,  fSinPhi);
1061   clusterRecord.SetGlobalDerivative( 1, -fCosPhi);
1062
1063   if (fBFieldOn)
1064   {
1065
1066     clusterRecord.SetGlobalDerivative(
1067        2,
1068       -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
1069       -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
1070
1071   } else {
1072
1073     clusterRecord.SetGlobalDerivative(
1074        2,
1075       -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
1076       -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
1077   }
1078
1079   clusterRecord.SetGlobalDerivative( 3, -fSinPhi*fTrackSlope0[0]+fCosPhi*fTrackSlope0[1]);
1080
1081   // append to trackRecord
1082   fTrackRecord.AddClusterRecord( clusterRecord );
1083
1084   // store local equation
1085   if( doAlignment ) LocalEquation( clusterRecord );
1086
1087 }
1088
1089 //_____________________________________________________
1090 void AliMUONAlignment::LocalEquation( const AliMUONAlignmentClusterRecord& clusterRecord )
1091 {
1092
1093   // copy to local variables
1094   for( Int_t index = 0; index < 4; ++index )
1095   {
1096     SetLocalDerivative( index, clusterRecord.GetLocalDerivative( index ) );
1097     SetGlobalDerivative( clusterRecord.GetDetElemNumber()*fgNParCh + index, clusterRecord.GetGlobalDerivative( index ) );
1098   }
1099
1100   // pass equation parameters to millepede
1101   fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, clusterRecord.GetMeas(), clusterRecord.GetSigma() );
1102
1103 }
1104
1105 //______________________________________________________________________
1106 void AliMUONAlignment::FillRecPointData()
1107 {
1108
1109   /// Get information of current cluster
1110   fClustPos[0] = fCluster->GetX();
1111   fClustPos[1] = fCluster->GetY();
1112   fClustPos[2] = fCluster->GetZ();
1113   fTransform->Global2Local(
1114     fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
1115     fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
1116
1117 }
1118
1119 //______________________________________________________________________
1120 void AliMUONAlignment::FillTrackParamData()
1121 {
1122
1123   /// Get information of current track at current cluster
1124   fTrackPos[0] = fTrackParam->GetNonBendingCoor();
1125   fTrackPos[1] = fTrackParam->GetBendingCoor();
1126   fTrackPos[2] = fTrackParam->GetZ();
1127   fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
1128   fTrackSlope[1] = fTrackParam->GetBendingSlope();
1129   fTransform->Global2Local(
1130     fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
1131     fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
1132
1133 }
1134
1135 //_____________________________________________________
1136 void AliMUONAlignment::FillDetElemData()
1137 {
1138
1139   /// Get information of current detection element
1140   Double_t lDetElemLocX = 0.;
1141   Double_t lDetElemLocY = 0.;
1142   Double_t lDetElemLocZ = 0.;
1143   fDetElemId = fCluster->GetDetElemId();
1144   fDetElemNumber = fDetElemId%100;
1145   for (int iCh=0; iCh<fDetElemId/100-1; iCh++)
1146   { fDetElemNumber += fgNDetElemCh[iCh]; }
1147
1148   fTransform->Local2Global(
1149     fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
1150     fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
1151
1152 }
1153
1154 //_____________________________________________________
1155 AliMUONAlignmentTrackRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment )
1156 {
1157   /// Process track and set Local Equations
1158   // store current track in running member.
1159   fTrack = track;
1160
1161   // clear track record
1162   fTrackRecord.Clear();
1163
1164   // get number of tracks
1165   Int_t nTrackParam = fTrack->GetTrackParamAtCluster()->GetEntries();
1166   AliDebug(1,Form("Number of track param entries : %i ", nTrackParam));
1167
1168   Bool_t first( kTRUE );
1169   for(Int_t iCluster=0; iCluster<nTrackParam; iCluster++)
1170   {
1171
1172     // and get new pointers
1173     fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtCluster()->At(iCluster);
1174     if ( ! fTrackParam ) continue;
1175     fCluster = fTrackParam->GetClusterPtr();
1176     if ( ! fCluster ) continue;
1177     // if (fDetElemId<500) continue;
1178
1179     // fill local variables for this position --> one measurement
1180     FillDetElemData();
1181     FillRecPointData();
1182     FillTrackParamData();
1183
1184     if( first )
1185     {
1186
1187       // for first valid cluster, save track position as "starting" values
1188       first = kFALSE;
1189
1190       fTrackPos0[0] = fTrackPos[0];
1191       fTrackPos0[1] = fTrackPos[1];
1192       fTrackPos0[2] = fTrackPos[2];
1193       fTrackSlope0[0] = fTrackSlope[0];
1194       fTrackSlope0[1] = fTrackSlope[1];
1195
1196     }
1197
1198     // calculate measurements
1199     fCosPhi = TMath::Cos(fPhi);
1200     fSinPhi = TMath::Sin(fPhi);
1201     if( fBFieldOn )
1202     {
1203
1204       fMeas[0] = fTrackPos[0] - fClustPos[0];
1205       fMeas[1] = fTrackPos[1] - fClustPos[1];
1206
1207     } else {
1208
1209       fMeas[0] = - fClustPos[0];
1210       fMeas[1] = - fClustPos[1];
1211
1212     }
1213
1214     // soùe debugging
1215     AliDebug(1,Form("cluster: %i", iCluster));
1216     AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
1217     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));
1218
1219     AliDebug(1,Form("Track Parameter: %i", iCluster));
1220     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]));
1221     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]));
1222
1223     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]));
1224
1225     // Set local equations
1226     LocalEquationX( doAlignment );
1227     LocalEquationY( doAlignment );
1228
1229   }
1230
1231   return &fTrackRecord;
1232
1233 }
1234
1235 //______________________________________________________________________________
1236 void AliMUONAlignment::ProcessTrack( AliMUONAlignmentTrackRecord* track, Bool_t doAlignment )
1237 {
1238   /// Process track (from record) and set Local Equations
1239   if( !( track && doAlignment ) ) return;
1240
1241   // loop over clusters
1242   for( Int_t index = 0; index < track->GetNRecords(); ++index )
1243   { if( AliMUONAlignmentClusterRecord* record = track->GetRecord( index ) ) LocalEquation( *record ); }
1244
1245   return;
1246
1247 }
1248
1249 //_____________________________________________________
1250 void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit)
1251 {
1252
1253   /// Call local fit for this tracks
1254   Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
1255   if (iRes && !lSingleFit)
1256   { fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1); }
1257
1258 }
1259
1260 //_____________________________________________________
1261 void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
1262 {
1263
1264   /// Call global fit; Global parameters are stored in parameters
1265   fMillepede->GlobalFit(parameters,errors,pulls);
1266
1267   AliInfo("Done fitting global parameters!");
1268   for (int iGlob=0; iGlob<fgNDetElem; iGlob++)
1269   { printf("%d\t %f\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+3],parameters[iGlob*fgNParCh+2]); }
1270
1271
1272 }
1273
1274 //_____________________________________________________
1275 Double_t AliMUONAlignment::GetParError(Int_t iPar)
1276 {
1277   /// Get error of parameter iPar
1278   Double_t lErr = fMillepede->GetParError(iPar);
1279   return lErr;
1280 }
1281
1282 //_____________________________________________________
1283 void AliMUONAlignment::PrintGlobalParameters()
1284 {
1285   /// Print global parameters
1286   fMillepede->PrintGlobalParameters();
1287 }
1288
1289 //_________________________________________________________________________
1290 TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, const double *lMisAlignment) const
1291 {
1292   /// Realign given transformation by given misalignment and return the misaligned transformation
1293
1294   Double_t cartMisAlig[3] = {0,0,0};
1295   Double_t angMisAlig[3] = {0,0,0};
1296 //   const Double_t *trans = transform.GetTranslation();
1297 //   TGeoRotation *rot;
1298 //   // check if the rotation we obtain is not NULL
1299 //   if (transform.GetRotation()) {
1300 //     rot = transform.GetRotation();
1301 //   }
1302 //   else {
1303 //     rot = new TGeoRotation("rot");
1304 //   }                  // default constructor.
1305
1306   cartMisAlig[0] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0])*lMisAlignment[0];
1307   cartMisAlig[1] = -TMath::Sign(1.0,transform.GetRotationMatrix()[4])*lMisAlignment[1];
1308   cartMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[8])*lMisAlignment[3];
1309   angMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0]*transform.GetRotationMatrix()[4])*lMisAlignment[2]*180./TMath::Pi();
1310
1311   TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
1312   TGeoRotation deltaRot;
1313   deltaRot.RotateX(angMisAlig[0]);
1314   deltaRot.RotateY(angMisAlig[1]);
1315   deltaRot.RotateZ(angMisAlig[2]);
1316
1317   TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
1318   TGeoHMatrix newTransfMat = transform * deltaTransf;
1319
1320   return TGeoCombiTrans(newTransfMat);
1321 }
1322
1323 //______________________________________________________________________
1324 AliMUONGeometryTransformer *
1325 AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
1326           const double *misAlignments, Bool_t verbose)
1327
1328 {
1329   /// Returns a new AliMUONGeometryTransformer with the found misalignments
1330   /// applied.
1331
1332   // Takes the internal geometry module transformers, copies them
1333   // and gets the Detection Elements from them.
1334   // Takes misalignment parameters and applies these
1335   // to the local transform of the Detection Element
1336   // Obtains the global transform by multiplying the module transformer
1337   // transformation with the local transformation
1338   // Applies the global transform to a new detection element
1339   // Adds the new detection element to a new module transformer
1340   // Adds the new module transformer to a new geometry transformer
1341   // Returns the new geometry transformer
1342
1343   Double_t lModuleMisAlignment[4] = {0.,0.,0.,0.};
1344   Double_t lDetElemMisAlignment[4] = {0.,0.,0.,0.};
1345   Int_t iDetElemId = 0;
1346   Int_t iDetElemNumber = 0;
1347
1348   AliMUONGeometryTransformer *newGeometryTransformer =
1349     new AliMUONGeometryTransformer();
1350   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
1351     // module transformers
1352     const AliMUONGeometryModuleTransformer *kModuleTransformer =
1353       transformer->GetModuleTransformer(iMt, true);
1354
1355     AliMUONGeometryModuleTransformer *newModuleTransformer =
1356       new AliMUONGeometryModuleTransformer(iMt);
1357     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1358
1359     TGeoCombiTrans moduleTransform =
1360       TGeoCombiTrans(*kModuleTransformer->GetTransformation());
1361     // New module transformation
1362     TGeoCombiTrans newModuleTransform = ReAlign(moduleTransform,lModuleMisAlignment);
1363     newModuleTransformer->SetTransformation(newModuleTransform);
1364
1365     // Get delta transformation:
1366     // Tdelta = Tnew * Told.inverse
1367     TGeoHMatrix deltaModuleTransform =
1368       AliMUONGeometryBuilder::Multiply(newModuleTransform,
1369                kModuleTransformer->GetTransformation()->Inverse());
1370     // Create module mis alignment matrix
1371     newGeometryTransformer
1372       ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
1373
1374     AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
1375
1376     if (verbose)
1377       AliInfo(Form("%i DEs in old GeometryStore  %i",detElements->GetSize(), iMt));
1378
1379     TIter next(detElements->CreateIterator());
1380     AliMUONGeometryDetElement* detElement;
1381     Int_t iDe(-1);
1382     while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
1383     {
1384       ++iDe;
1385       // make a new detection element
1386       AliMUONGeometryDetElement *newDetElement =
1387   new AliMUONGeometryDetElement(detElement->GetId(),
1388               detElement->GetVolumePath());
1389       TString lDetElemName(detElement->GetDEName());
1390       lDetElemName.ReplaceAll("DE","");
1391       iDetElemId = lDetElemName.Atoi();
1392       iDetElemNumber = iDetElemId%100;
1393       for (int iCh=0; iCh<iDetElemId/100-1; iCh++){
1394   iDetElemNumber += fgNDetElemCh[iCh];
1395       }
1396       for (int i=0; i<fgNParCh; i++) {
1397   lDetElemMisAlignment[i] = 0.0;
1398   if (iMt<fgNTrkMod) {
1399     AliInfo(Form("iMt %i, iCh %i, iDe %i, iDeId %i, iDeNb %i, iPar %i",iMt, iDetElemId/100, iDe, iDetElemId, iDetElemNumber, iDetElemNumber*fgNParCh+i));
1400     lDetElemMisAlignment[i] =  misAlignments[iDetElemNumber*fgNParCh+i];
1401   }
1402       }
1403       // local transformation of this detection element.
1404       TGeoCombiTrans localTransform
1405   = TGeoCombiTrans(*detElement->GetLocalTransformation());
1406       TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
1407       newDetElement->SetLocalTransformation(newLocalTransform);
1408
1409       // global transformation
1410       TGeoHMatrix newGlobalTransform =
1411   AliMUONGeometryBuilder::Multiply(newModuleTransform,
1412            newLocalTransform);
1413       newDetElement->SetGlobalTransformation(newGlobalTransform);
1414
1415       // add this det element to module
1416       newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
1417                   newDetElement);
1418
1419       // In the Alice Alignment Framework misalignment objects store
1420       // global delta transformation
1421       // Get detection "intermediate" global transformation
1422       TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform;
1423       // Get detection element global delta transformation:
1424       // Tdelta = Tnew * Told.inverse
1425       TGeoHMatrix  deltaGlobalTransform
1426   = AliMUONGeometryBuilder::Multiply(newGlobalTransform,
1427              newOldGlobalTransform.Inverse());
1428
1429       // Create mis alignment matrix
1430       newGeometryTransformer
1431   ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
1432     }
1433
1434     if (verbose)
1435       AliInfo(Form("Added module transformer %i to the transformer", iMt));
1436     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1437   }
1438   return newGeometryTransformer;
1439 }
1440
1441 //______________________________________________________________________
1442 void AliMUONAlignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY)
1443 {
1444
1445   /// Set alignment resolution to misalign objects to be stored in CDB
1446   Int_t chIdMin = (rChId<0)? 0 : rChId;
1447   Int_t chIdMax = (rChId<0)? 9 : rChId;
1448   Double_t chResX = rChResX;
1449   Double_t chResY = rChResY;
1450   Double_t deResX = rDeResX;
1451   Double_t deResY = rDeResY;
1452
1453   TMatrixDSym mChCorrMatrix(6);
1454   mChCorrMatrix[0][0]=chResX*chResX;
1455   mChCorrMatrix[1][1]=chResY*chResY;
1456   //  mChCorrMatrix.Print();
1457
1458   TMatrixDSym mDECorrMatrix(6);
1459   mDECorrMatrix[0][0]=deResX*deResX;
1460   mDECorrMatrix[1][1]=deResY*deResY;
1461   //  mDECorrMatrix.Print();
1462
1463   AliAlignObjMatrix *alignMat = 0x0;
1464
1465   for(Int_t chId=chIdMin; chId<=chIdMax; chId++) {
1466     TString chName1;
1467     TString chName2;
1468     if (chId<4){
1469       chName1 = Form("GM%d",chId);
1470       chName2 = Form("GM%d",chId);
1471     } else {
1472       chName1 = Form("GM%d",4+(chId-4)*2);
1473       chName2 = Form("GM%d",4+(chId-4)*2+1);
1474     }
1475
1476     for (int i=0; i<misAlignArray->GetEntries(); i++) {
1477       alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
1478       TString volName(alignMat->GetSymName());
1479       if((volName.Contains(chName1)&&
1480     ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
1481      (volName.Length()==volName.Index(chName1)+chName1.Length())))||
1482    (volName.Contains(chName2)&&
1483     ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
1484      (volName.Length()==volName.Index(chName2)+chName2.Length())))){
1485   volName.Remove(0,volName.Last('/')+1);
1486   if (volName.Contains("GM")) {
1487     //  alignMat->Print("NULL");
1488     alignMat->SetCorrMatrix(mChCorrMatrix);
1489   } else if (volName.Contains("DE")) {
1490     //  alignMat->Print("NULL");
1491     alignMat->SetCorrMatrix(mDECorrMatrix);
1492   }
1493       }
1494     }
1495   }
1496 }