rewrote AliMUONAlignment to - use AliMillePede2 instead of AliMillepede - use AliMill...
[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: AliMUONAlignment.cxx 51000 2011-08-08 17:58:17Z ivana $ */
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 "AliMillePede2.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 #include <TGraphErrors.h>
50
51 /// \cond CLASSIMP
52 ClassImp(AliMUONAlignment)
53 /// \endcond
54
55 //_____________________________________________________________________
56 // static variables
57 const Int_t AliMUONAlignment::fgNDetElemCh[AliMUONAlignment::fgNCh] = { 4, 4, 4, 4, 18, 18, 26, 26, 26, 26 };
58 const Int_t AliMUONAlignment::fgSNDetElemCh[AliMUONAlignment::fgNCh+1] = { 0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156 };
59
60 //_____________________________________________________________________
61 /// self initialized array, used for adding constraints
62 class Array
63 {
64
65   public:
66
67   /// contructor
68   Array( void )
69   {
70     for( Int_t i=0; i < AliMUONAlignment::fNGlobal; ++i )
71     { values[i] = 0; }
72   }
73
74   /// array
75   Double_t values[AliMUONAlignment::fNGlobal];
76
77   private:
78
79   /// Not implemented
80   Array(const Array& );
81
82   /// Not implemented
83   Array&  operator = (const Array& );
84
85 };
86
87 //_____________________________________________________________________
88 AliMUONAlignment::AliMUONAlignment()
89   : TObject(),
90     fInitialized( kFALSE ),
91     fRunNumber( 0 ),
92     fBFieldOn( kTRUE ),
93     fStartFac( 256 ),
94     fResCutInitial( 100 ),
95     fResCut( 100 ),
96     fMillepede( 0L ),
97     fCluster( 0L ),
98     fNStdDev( 3 ),
99     fDetElemNumber( 0 ),
100     fUnbias( kFALSE ),
101     fTrackRecord(),
102     fTransform( 0 ),
103     fGeoCombiTransInverse()
104 {
105   /// constructor
106   fSigma[0] = 1.5e-1;
107   fSigma[1] = 1.0e-2;
108
109   // default allowed variations
110   fAllowVar[0] = 0.5;  // x
111   fAllowVar[1] = 0.5;  // y
112   fAllowVar[2] = 0.01; // phi_z
113   fAllowVar[3] = 5;    // z
114
115   // initialize millepede
116   fMillepede = new AliMillePede2();
117
118   // initialize degrees of freedom
119   // by default all parameters are free
120   for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
121   { fGlobalParameterStatus[iPar] = kFreeParId; }
122
123   // initialize local equations
124   for(int i=0; i<fNLocal; ++i )
125   { fLocalDerivatives[i] = 0.0; }
126
127   for(int i=0; i<fNGlobal; ++i )
128   { fGlobalDerivatives[i] = 0.0; }
129
130 }
131
132 //_____________________________________________________________________
133 AliMUONAlignment::~AliMUONAlignment()
134 {
135   /// destructor
136 }
137
138 //_____________________________________________________________________
139 void AliMUONAlignment::Init( void )
140 {
141
142   /// initialize
143   /**
144   initialize millipede
145   must be called after necessary detectors have been fixed,
146   but before constrains are added and before global parameters initial value are set
147   */
148   if( fInitialized )
149   { AliFatal( "Millepede already initialized" ); }
150
151   // assign proper groupID to free parameters
152   Int_t nGlobal = 0;
153   for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
154   {
155
156     if( fGlobalParameterStatus[iPar] == kFixedParId )
157     {
158       // fixed parameters are left unchanged
159       continue;
160
161     } else if( fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId ) {
162
163       // free parameters or first element of group are assigned a new group id
164       fGlobalParameterStatus[iPar] = nGlobal++;
165       continue;
166
167     } else if( fGlobalParameterStatus[iPar] < kGroupBaseId ) {
168
169       // get detector element id from status, get chamber parameter id
170       const Int_t iDeBase( kGroupBaseId - 1 - fGlobalParameterStatus[iPar] );
171       const Int_t iParBase = iPar%fgNParCh;
172
173       // check
174       if( iDeBase < 0 || iDeBase >= iPar/fgNParCh )
175       { AliFatal( Form( "Group for parameter index %i has wrong base detector element: %i", iPar, iDeBase ) ); }
176
177       // assign identical group id to current
178       fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase*fgNParCh + iParBase];
179       AliInfo( Form( "Parameter %i grouped to detector %i (%s)", iPar, iDeBase, GetParameterMaskString( 1<<iParBase ).Data() ) );
180
181     } else AliFatal( Form( "Unrecognized parameter status for index %i: %i", iPar, fGlobalParameterStatus[iPar] ) );
182
183   }
184
185   AliInfo( Form( "Free Parameters: %i out of %i", nGlobal, fNGlobal ) );
186
187   // initialize millepede
188   fMillepede->InitMille( fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus );
189   fInitialized = kTRUE;
190
191   // some debug output
192   for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
193   {  AliInfo( Form( "fAllowVar[%i]= %f", iPar, fAllowVar[iPar] ) ); }
194
195   // set allowed variations for all parameters
196   for( Int_t iDet = 0; iDet < fgNDetElem; ++iDet )
197   {
198     for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
199     { fMillepede->SetParSigma( iDet*fgNParCh + iPar, fAllowVar[iPar] ); }
200   }
201
202   // Set iterations
203   if (fStartFac>1) fMillepede->SetIterations(fStartFac);
204
205 }
206
207 //_____________________________________________________
208 AliMillePedeRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight )
209 {
210   /// process track for alignment minimization
211   /**
212   returns the alignment records for this track.
213   They can be stored in some output for later reprocessing.
214   */
215
216   // reset track records
217   fTrackRecord.Reset();
218   if( fMillepede->GetRecord() ) fMillepede->GetRecord()->Reset();
219
220   // get number of track parameters
221   Int_t nTrackParam = track->GetTrackParamAtCluster()->GetEntries();
222
223   Bool_t first( kTRUE );
224   for( Int_t iTrackParam = 0; iTrackParam < nTrackParam; ++iTrackParam )
225   {
226
227     // get new pointers
228     AliMUONTrackParam* trackParam( (AliMUONTrackParam *) track->GetTrackParamAtCluster()->At(iTrackParam) );
229     if( !trackParam ) continue;
230
231     AliMUONVCluster* cluster = trackParam->GetClusterPtr();
232     if( !cluster ) continue;
233
234     // fill local variables for this position --> one measurement
235     FillDetElemData( cluster );
236     FillRecPointData( cluster );
237
238     // unbias and store track parameters
239     if( fUnbias && !UnbiasTrackParamData( trackParam ) ) continue;
240     FillTrackParamData( trackParam );
241
242     if( first )
243     {
244
245       // for first valid cluster, save track position as "starting" values
246       first = kFALSE;
247
248       fTrackPos0[0] = fTrackPos[0];
249       fTrackPos0[1] = fTrackPos[1];
250       fTrackPos0[2] = fTrackPos[2];
251       fTrackSlope0[0] = fTrackSlope[0];
252       fTrackSlope0[1] = fTrackSlope[1];
253
254     }
255
256     // 'inverse' (GlobalToLocal) rotation matrix
257     const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
258
259     // calculate measurements
260     if( fBFieldOn )
261     {
262
263       // use residuals (cluster - track) for measurement
264       fMeas[0] = r[0]*(fClustPos[0] - fTrackPos[0]) + r[1]*(fClustPos[1] - fTrackPos[1]);
265       fMeas[1] = r[3]*(fClustPos[0] - fTrackPos[0]) + r[4]*(fClustPos[1] - fTrackPos[1]);
266
267     } else {
268
269       // use cluster position for measurement
270       fMeas[0] = ( r[0]*fClustPos[0] + r[1]*fClustPos[1] );
271       fMeas[1] = ( r[3]*fClustPos[0] + r[4]*fClustPos[1] );
272
273     }
274
275     // Set local equations
276     LocalEquationX();
277     LocalEquationY();
278
279   }
280
281   // copy track record
282   fMillepede->SetRecordRun(fRunNumber);
283   fMillepede->SetRecordWeight(weight);
284   fTrackRecord = *fMillepede->GetRecord();
285
286   // save record data
287   if( doAlignment )
288   {
289     fMillepede->SaveRecordData();
290   }
291
292   // return record
293   return &fTrackRecord;
294
295 }
296
297 //______________________________________________________________________________
298 void AliMUONAlignment::ProcessTrack( AliMillePedeRecord* trackRecord )
299 {
300   /// process track record
301   if( !trackRecord ) return;
302
303   // make sure record storage is initialized
304   if( !fMillepede->GetRecord() ) fMillepede->InitDataRecStorage();
305
306   // copy content
307   *fMillepede->GetRecord() = *trackRecord;
308
309   // save record
310   fMillepede->SaveRecordData();
311
312   return;
313
314 }
315
316 //_____________________________________________________________________
317 void AliMUONAlignment::FixAll( UInt_t mask )
318 {
319   /// fix parameters matching mask, for all chambers
320   AliInfo( Form( "Fixing %s for all detector elements", GetParameterMaskString( mask ).Data() ) );
321
322   // fix all stations
323   for( Int_t i = 0; i < fgNDetElem; ++i )
324   {
325     if( mask & ParX )  FixParameter(i, 0);
326     if( mask & ParY )  FixParameter(i, 1);
327     if( mask & ParTZ ) FixParameter(i, 2);
328     if( mask & ParZ )  FixParameter(i, 3);
329   }
330
331 }
332
333 //_____________________________________________________________________
334 void AliMUONAlignment::FixChamber( Int_t iCh, UInt_t mask )
335 {
336   /// fix parameters matching mask, for all detector elements in a given chamber, counting from 1
337
338   // check boundaries
339   if( iCh < 1 || iCh > 10 )
340   { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
341
342   // get first and last element
343   const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
344   const Int_t iDetElemLast = fgSNDetElemCh[iCh];
345   for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
346   {
347
348     AliInfo( Form( "Fixing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
349
350     if( mask & ParX )  FixParameter(i, 0);
351     if( mask & ParY )  FixParameter(i, 1);
352     if( mask & ParTZ ) FixParameter(i, 2);
353     if( mask & ParZ )  FixParameter(i, 3);
354
355   }
356 }
357
358 //_____________________________________________________________________
359 void AliMUONAlignment::FixDetElem( Int_t iDetElemId, UInt_t mask )
360 {
361   /// fix parameters matching mask, for a given detector element, counting from 0
362   const Int_t iDet( GetDetElemNumber( iDetElemId ) );
363   if ( mask & ParX )  FixParameter(iDet, 0);
364   if ( mask & ParY )  FixParameter(iDet, 1);
365   if ( mask & ParTZ ) FixParameter(iDet, 2);
366   if ( mask & ParZ )  FixParameter(iDet, 3);
367
368 }
369
370 //_____________________________________________________________________
371 void AliMUONAlignment::FixHalfSpectrometer( const Bool_t *lChOnOff, UInt_t sidesMask, UInt_t mask )
372 {
373
374   /// Fix parameters matching mask for all detectors in selected chambers and selected sides of the spectrometer
375   for( Int_t i = 0; i < fgNDetElem; ++i )
376   {
377
378     // get chamber matching detector
379     const Int_t iCh( GetChamberId(i) );
380     if( !lChOnOff[iCh-1] ) continue;
381
382     // get detector element in chamber
383     Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
384
385     // skip detector if its side is off
386     // stations 1 and 2
387     if( iCh>=1 && iCh<=4 )
388     {
389       if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
390       if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
391       if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
392       if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
393     }
394
395     // station 3
396     if (iCh>=5 && iCh<=6)
397     {
398       if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
399       if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
400       if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
401       if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
402     }
403
404     // stations 4 and 5
405     if (iCh>=7 && iCh<=10)
406     {
407       if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
408       if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
409       if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
410       if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
411     }
412
413     // detector is accepted, fix it
414     FixDetElem( i, mask );
415
416   }
417
418 }
419
420 //______________________________________________________________________
421 void AliMUONAlignment::FixParameter( Int_t iPar )
422 {
423
424   /// fix a given parameter, counting from 0
425   if( fInitialized )
426   { AliFatal( "Millepede already initialized" ); }
427
428   fGlobalParameterStatus[iPar] = kFixedParId;
429
430 }
431
432
433 //_____________________________________________________________________
434 void AliMUONAlignment::ReleaseChamber( Int_t iCh, UInt_t mask )
435 {
436   /// release parameters matching mask, for all detector elements in a given chamber, counting from 1
437
438   // check boundaries
439   if( iCh < 1 || iCh > 10 )
440   { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
441
442   // get first and last element
443   const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
444   const Int_t iDetElemLast = fgSNDetElemCh[iCh];
445   for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
446   {
447
448     AliInfo( Form( "Releasing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
449
450     if( mask & ParX )  ReleaseParameter(i, 0);
451     if( mask & ParY )  ReleaseParameter(i, 1);
452     if( mask & ParTZ ) ReleaseParameter(i, 2);
453     if( mask & ParZ )  ReleaseParameter(i, 3);
454
455   }
456 }
457
458 //_____________________________________________________________________
459 void AliMUONAlignment::ReleaseDetElem( Int_t iDetElemId, UInt_t mask )
460 {
461   /// release parameters matching mask, for a given detector element, counting from 0
462   const Int_t iDet( GetDetElemNumber( iDetElemId ) );
463   if ( mask & ParX )  ReleaseParameter(iDet, 0);
464   if ( mask & ParY )  ReleaseParameter(iDet, 1);
465   if ( mask & ParTZ ) ReleaseParameter(iDet, 2);
466   if ( mask & ParZ )  ReleaseParameter(iDet, 3);
467
468 }
469
470 //______________________________________________________________________
471 void AliMUONAlignment::ReleaseParameter( Int_t iPar )
472 {
473
474   /// release a given parameter, counting from 0
475   if( fInitialized )
476   { AliFatal( "Millepede already initialized" ); }
477
478   fGlobalParameterStatus[iPar] = kFreeParId;
479
480 }
481
482 //_____________________________________________________________________
483 void AliMUONAlignment::GroupChamber( Int_t iCh, UInt_t mask )
484 {
485   /// group parameters matching mask for all detector elements in a given chamber, counting from 1
486   if( iCh < 1 || iCh > 10 )
487   { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
488
489   const Int_t detElemMin = 100*iCh;
490   const Int_t detElemMax = 100*iCh + fgNDetElemCh[iCh]-1;
491   GroupDetElems( detElemMin, detElemMax, mask );
492
493 }
494
495 //_____________________________________________________________________
496 void AliMUONAlignment::GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t mask )
497 {
498   /// group parameters matching mask for all detector elements between min and max
499   // check number of detector elements
500   const Int_t nDetElem = detElemMax - detElemMin + 1;
501   if( nDetElem<2 )
502   { AliFatal( Form( "Requested group of DEs %d-%d contains less than 2 DE's", detElemMin, detElemMax ) ); }
503
504   // create list
505   Int_t* detElemList = new int[nDetElem];
506   for( Int_t i = 0; i < nDetElem; ++i )
507   { detElemList[i] = detElemMin+i; }
508
509   // group
510   GroupDetElems( detElemList, nDetElem, mask );
511   delete[] detElemList;
512
513 }
514
515 //_____________________________________________________________________
516 void AliMUONAlignment::GroupDetElems( Int_t* detElemList, Int_t nDetElem, UInt_t mask )
517 {
518   /// group parameters matching mask for all detector elements in list
519   if( fInitialized )
520   { AliFatal( "Millepede already initialized" ); }
521
522   const Int_t iDeBase( GetDetElemNumber( detElemList[0] ) );
523   for( Int_t i = 0; i < nDetElem; ++i )
524   {
525     const Int_t iDeCurrent( GetDetElemNumber( detElemList[i] ) );
526     if( mask & ParX ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 0] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
527     if( mask & ParY ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 1] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
528     if( mask & ParTZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 2] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
529     if( mask & ParZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 3] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
530
531     if( i== 0 ) AliInfo( Form( "Creating new group for detector %i and variable %s", detElemList[i], GetParameterMaskString( mask ).Data() ) );
532     else AliInfo( Form( "Adding detector element %i to current group", detElemList[i] ) );
533   }
534
535 }
536
537 //______________________________________________________________________
538 void AliMUONAlignment::SetChamberNonLinear( Int_t iCh, UInt_t mask )
539 {
540   /// Set parameters matching mask as non linear, for all detector elements in a given chamber, counting from 1
541   const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
542   const Int_t iDetElemLast = fgSNDetElemCh[iCh];
543   for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
544   {
545
546       if( mask & ParX ) SetParameterNonLinear(i, 0);
547       if( mask & ParY ) SetParameterNonLinear(i, 1);
548       if( mask & ParTZ ) SetParameterNonLinear(i, 2);
549       if( mask & ParZ ) SetParameterNonLinear(i, 3);
550
551   }
552
553 }
554
555 //_____________________________________________________________________
556 void AliMUONAlignment::SetDetElemNonLinear( Int_t iDetElemId, UInt_t mask )
557 {
558   /// Set parameters matching mask as non linear, for a given detector element, counting from 0
559   const Int_t iDet( GetDetElemNumber( iDetElemId ) );
560   if ( mask & ParX )  SetParameterNonLinear(iDet, 0);
561   if ( mask & ParY )  SetParameterNonLinear(iDet, 1);
562   if ( mask & ParTZ ) SetParameterNonLinear(iDet, 2);
563   if ( mask & ParZ )  SetParameterNonLinear(iDet, 3);
564
565 }
566
567 //______________________________________________________________________
568 void AliMUONAlignment::SetParameterNonLinear( Int_t iPar )
569 {
570   /// Set nonlinear flag for parameter iPar
571   if( !fInitialized )
572   { AliFatal( "Millepede not initialized" ); }
573
574   fMillepede->SetNonLinear( iPar );
575   AliInfo( Form( "Parameter %i set to non linear", iPar ) );
576 }
577
578 //______________________________________________________________________
579 void AliMUONAlignment::AddConstraints( const Bool_t *lChOnOff, UInt_t mask )
580 {
581   /// Add constraint equations for selected chambers and degrees of freedom
582
583   Array fConstraintX;
584   Array fConstraintY;
585   Array fConstraintTZ;
586   Array fConstraintZ;
587
588   for( Int_t i = 0; i < fgNDetElem; ++i )
589   {
590
591     // get chamber matching detector
592     const Int_t iCh( GetChamberId(i) );
593     if (lChOnOff[iCh-1])
594     {
595
596       if( mask & ParX ) fConstraintX.values[i*fgNParCh+0]=1.0;
597       if( mask & ParY ) fConstraintY.values[i*fgNParCh+1]=1.0;
598       if( mask & ParTZ ) fConstraintTZ.values[i*fgNParCh+2]=1.0;
599       if( mask & ParZ ) fConstraintTZ.values[i*fgNParCh+3]=1.0;
600
601     }
602   }
603
604   if( mask & ParX ) AddConstraint(fConstraintX.values,0.0);
605   if( mask & ParY ) AddConstraint(fConstraintY.values,0.0);
606   if( mask & ParTZ ) AddConstraint(fConstraintTZ.values,0.0);
607   if( mask & ParZ ) AddConstraint(fConstraintZ.values,0.0);
608
609 }
610
611 //______________________________________________________________________
612 void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask )
613 {
614   /*
615   questions:
616   - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
617   - why is weight ignored for ConstrainT and ConstrainB
618   - why is there no constrain on z
619   */
620
621   /// Add constraint equations for selected chambers, degrees of freedom and detector half
622   Double_t lMeanY = 0.;
623   Double_t lSigmaY = 0.;
624   Double_t lMeanZ = 0.;
625   Double_t lSigmaZ = 0.;
626   Int_t lNDetElem = 0;
627
628   for( Int_t i = 0; i < fgNDetElem; ++i )
629   {
630
631     // get chamber matching detector
632     const Int_t iCh( GetChamberId(i) );
633
634     // skip detector if chamber is off
635     if( lChOnOff[iCh-1] ) continue;
636
637     // get detector element id from detector element number
638     const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
639     const Int_t lDetElemId = iCh*100+lDetElemNumber;
640
641     // skip detector if its side is off
642     // stations 1 and 2
643     if( iCh>=1 && iCh<=4 )
644     {
645       if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
646       if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
647       if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
648       if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
649     }
650
651     // station 3
652     if (iCh>=5 && iCh<=6)
653     {
654       if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
655       if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
656       if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
657       if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
658     }
659
660     // stations 4 and 5
661     if (iCh>=7 && iCh<=10)
662     {
663       if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
664       if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
665       if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
666       if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
667     }
668
669     // get global x, y and z position
670     Double_t lDetElemGloX = 0.;
671     Double_t lDetElemGloY = 0.;
672     Double_t lDetElemGloZ = 0.;
673     fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
674
675     // increment mean Y, mean Z, sigmas and number of accepted detectors
676     lMeanY += lDetElemGloY;
677     lSigmaY += lDetElemGloY*lDetElemGloY;
678     lMeanZ += lDetElemGloZ;
679     lSigmaZ += lDetElemGloZ*lDetElemGloZ;
680     lNDetElem++;
681
682   }
683
684   // calculate mean values
685   lMeanY /= lNDetElem;
686   lSigmaY /= lNDetElem;
687   lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
688   lMeanZ /= lNDetElem;
689   lSigmaZ /= lNDetElem;
690   lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
691   AliInfo( Form( "Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ ) );
692
693   // create all possible arrays
694   Array fConstraintX[4];  //Array for constraint equation X
695   Array fConstraintY[4];  //Array for constraint equation Y
696   Array fConstraintP[4];  //Array for constraint equation P
697   Array fConstraintXZ[4];  //Array for constraint equation X vs Z
698   Array fConstraintYZ[4];  //Array for constraint equation Y vs Z
699   Array fConstraintPZ[4];  //Array for constraint equation P vs Z
700
701   // do we really need these ?
702   Array fConstraintXY[4];  //Array for constraint equation X vs Y
703   Array fConstraintYY[4];  //Array for constraint equation Y vs Y
704   Array fConstraintPY[4];  //Array for constraint equation P vs Y
705
706   // fill Bool_t sides array based on masks, for convenience
707   Bool_t lDetTLBR[4];
708   lDetTLBR[0] = sidesMask & SideTop;
709   lDetTLBR[1] = sidesMask & SideLeft;
710   lDetTLBR[2] = sidesMask & SideBottom;
711   lDetTLBR[3] = sidesMask & SideRight;
712
713   for( Int_t i = 0; i < fgNDetElem; ++i )
714   {
715
716     // get chamber matching detector
717     const Int_t iCh( GetChamberId(i) );
718
719     // skip detector if chamber is off
720     if( !lChOnOff[iCh-1] ) continue;
721
722     // get detector element id from detector element number
723     const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
724     const Int_t lDetElemId = iCh*100+lDetElemNumber;
725
726     // get global x, y and z position
727     Double_t lDetElemGloX = 0.;
728     Double_t lDetElemGloY = 0.;
729     Double_t lDetElemGloZ = 0.;
730     fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
731
732     // loop over sides
733     for( Int_t iSide = 0; iSide < 4; iSide++ )
734     {
735
736       // skip if side is not selected
737       if( !lDetTLBR[iSide] ) continue;
738
739       // skip detector if it is not in the selected side
740       // stations 1 and 2
741       if( iCh>=1 && iCh<=4 )
742       {
743         if( lDetElemNumber == 0 && !(iSide == 0 || iSide == 3) ) continue; // top-right
744         if( lDetElemNumber == 1 && !(iSide == 0 || iSide == 1) ) continue; // top-left
745         if( lDetElemNumber == 2 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
746         if( lDetElemNumber == 3 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
747       }
748
749       // station 3
750       if (iCh>=5 && iCh<=6)
751       {
752         if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3) ) continue; // top-right
753         if( lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1) ) continue; // top-left
754         if( lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
755         if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
756       }
757
758       // stations 4 and 5
759       if (iCh>=7 && iCh<=10)
760       {
761         if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3) ) continue; // top-right
762         if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1) ) continue; // top-left
763         if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
764         if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
765       }
766
767       // constrain x
768       if( lVarXYT[0] ) fConstraintX[iSide].values[i*fgNParCh+0] = 1;
769
770       // constrain y
771       if( lVarXYT[1] ) fConstraintY[iSide].values[i*fgNParCh+1] = 1;
772
773       // constrain phi (rotation around z)
774       if( lVarXYT[2] ) fConstraintP[iSide].values[i*fgNParCh+2] = 1;
775
776       // x-z shearing
777       if( lVarXYT[3] ) fConstraintXZ[iSide].values[i*fgNParCh+0] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
778
779       // y-z shearing
780       if( lVarXYT[4] ) fConstraintYZ[iSide].values[i*fgNParCh+1] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
781
782       // phi-z shearing
783       if( lVarXYT[5] ) fConstraintPZ[iSide].values[i*fgNParCh+2] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
784
785       // x-y shearing
786       if( lVarXYT[6] ) fConstraintXY[iSide].values[i*fgNParCh+0] = (lDetElemGloY-lMeanY)/lSigmaY;
787
788       // y-y shearing
789       if( lVarXYT[7] ) fConstraintYY[iSide].values[i*fgNParCh+1] = (lDetElemGloY-lMeanY)/lSigmaY;
790
791       // phi-y shearing
792       if( lVarXYT[8] ) fConstraintPY[iSide].values[i*fgNParCh+2] = (lDetElemGloY-lMeanY)/lSigmaY;
793
794     }
795
796   }
797
798   // pass constraints to millepede
799   for( Int_t iSide = 0; iSide < 4; iSide++ )
800   {
801     // skip if side is not selected
802     if( !lDetTLBR[iSide] ) continue;
803
804     if( lVarXYT[0] ) AddConstraint(fConstraintX[iSide].values,0.0);
805     if( lVarXYT[1] ) AddConstraint(fConstraintY[iSide].values,0.0);
806     if( lVarXYT[2] ) AddConstraint(fConstraintP[iSide].values,0.0);
807     if( lVarXYT[3] ) AddConstraint(fConstraintXZ[iSide].values,0.0);
808     if( lVarXYT[4] ) AddConstraint(fConstraintYZ[iSide].values,0.0);
809     if( lVarXYT[5] ) AddConstraint(fConstraintPZ[iSide].values,0.0);
810     if( lVarXYT[6] ) AddConstraint(fConstraintXY[iSide].values,0.0);
811     if( lVarXYT[7] ) AddConstraint(fConstraintYY[iSide].values,0.0);
812     if( lVarXYT[8] ) AddConstraint(fConstraintPY[iSide].values,0.0);
813   }
814
815 }
816
817 //______________________________________________________________________
818 void AliMUONAlignment::InitGlobalParameters(Double_t *par)
819 {
820   /// Initialize global parameters with par array
821   if( !fInitialized )
822   { AliFatal( "Millepede is not initialized" ); }
823
824   fMillepede->SetGlobalParameters(par);
825 }
826
827 //______________________________________________________________________
828 void AliMUONAlignment::SetAllowedVariation( Int_t iPar, Double_t value )
829 {
830   /// "Encouraged" variation for degrees of freedom
831   // check initialization
832   if( fInitialized )
833   { AliFatal( "Millepede already initialized" ); }
834
835   // check initialization
836   if( !(iPar >= 0 && iPar < fgNParCh ) )
837   { AliFatal( Form( "Invalid index: %i", iPar ) ); }
838
839   fAllowVar[iPar] = value;
840 }
841
842 //______________________________________________________________________
843 void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
844 {
845
846   /// Set expected measurement resolution
847   fSigma[0] = sigmaX;
848   fSigma[1] = sigmaY;
849
850   // print
851   for( Int_t i=0; i<2; ++i )
852   { AliInfo( Form( "fSigma[%i]=%f", i, fSigma[i] ) ); }
853
854 }
855
856 //_____________________________________________________
857 void AliMUONAlignment::GlobalFit( Double_t *parameters, Double_t *errors, Double_t *pulls )
858 {
859
860   /// Call global fit; Global parameters are stored in parameters
861   fMillepede->GlobalFit( parameters, errors, pulls );
862
863   AliInfo( "Done fitting global parameters" );
864   for( int iDet=0; iDet<fgNDetElem; ++iDet )
865   {
866     AliInfo( Form( "%d\t %f\t %f\t %f\t %f",
867       iDet,
868       parameters[iDet*fgNParCh+0],parameters[iDet*fgNParCh+1],
869       parameters[iDet*fgNParCh+3],parameters[iDet*fgNParCh+2]
870       ) );
871   }
872
873 }
874
875 //_____________________________________________________
876 void AliMUONAlignment::PrintGlobalParameters() const
877 { fMillepede->PrintGlobalParameters(); }
878
879 //_____________________________________________________
880 Double_t AliMUONAlignment::GetParError(Int_t iPar) const
881 { return fMillepede->GetParError(iPar); }
882
883 //______________________________________________________________________
884 AliMUONGeometryTransformer* AliMUONAlignment::ReAlign(
885   const AliMUONGeometryTransformer * transformer,
886   const double *misAlignments, Bool_t )
887 {
888
889   /// Returns a new AliMUONGeometryTransformer with the found misalignments
890   /// applied.
891
892   // Takes the internal geometry module transformers, copies them
893   // and gets the Detection Elements from them.
894   // Takes misalignment parameters and applies these
895   // to the local transform of the Detection Element
896   // Obtains the global transform by multiplying the module transformer
897   // transformation with the local transformation
898   // Applies the global transform to a new detection element
899   // Adds the new detection element to a new module transformer
900   // Adds the new module transformer to a new geometry transformer
901   // Returns the new geometry transformer
902
903   Double_t lModuleMisAlignment[fgNParCh] = {0};
904   Double_t lDetElemMisAlignment[fgNParCh] = {0};
905   const TClonesArray* oldMisAlignArray( transformer->GetMisAlignmentData() );
906
907   AliMUONGeometryTransformer *newGeometryTransformer = new AliMUONGeometryTransformer();
908   for( Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt )
909   {
910
911     // module transformers
912     const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
913
914     AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
915     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
916
917     // get transformation
918     TGeoHMatrix deltaModuleTransform( DeltaTransform( lModuleMisAlignment ) );
919
920     // update module
921     TGeoHMatrix moduleTransform( *kModuleTransformer->GetTransformation() );
922     TGeoHMatrix newModuleTransform( AliMUONGeometryBuilder::Multiply( deltaModuleTransform, moduleTransform ) );
923     newModuleTransformer->SetTransformation(newModuleTransform);
924
925     // Get matching old alignment and update current matrix accordingly
926     if( oldMisAlignArray )
927     {
928
929       const AliAlignObjMatrix* oldAlignObj(0);
930       const Int_t moduleId( kModuleTransformer->GetModuleId() );
931       const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId );
932       for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
933       {
934
935         const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
936         if( localAlignObj && localAlignObj->GetVolUID() == volId )
937         {
938           oldAlignObj = localAlignObj;
939           break;
940         }
941
942       }
943
944       // multiply
945       if( oldAlignObj )
946       {
947
948         TGeoHMatrix oldMatrix;
949         oldAlignObj->GetMatrix( oldMatrix );
950         deltaModuleTransform.Multiply( &oldMatrix );
951
952       }
953
954     }
955
956     // Create module mis alignment matrix
957     newGeometryTransformer ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
958
959     AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
960
961     TIter next(detElements->CreateIterator());
962     AliMUONGeometryDetElement* detElement;
963     Int_t iDe(-1);
964     while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
965     {
966       ++iDe;
967       // make a new detection element
968       AliMUONGeometryDetElement *newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
969       TString lDetElemName(detElement->GetDEName());
970       lDetElemName.ReplaceAll("DE","");
971
972       // store detector element id and number
973       const Int_t iDetElemId = lDetElemName.Atoi();
974       if( !DetElemIsValid( iDetElemId ) )
975       {
976         AliInfo( Form( "Skipping invalid detector element %i", iDetElemId ) );
977         continue;
978       }
979
980       const Int_t iDetElemNumber( GetDetElemNumber( iDetElemId ) );
981
982       for( int i=0; i<fgNParCh; ++i )
983       {
984         lDetElemMisAlignment[i] = 0.0;
985         if( iMt<fgNTrkMod ) { lDetElemMisAlignment[i] =  misAlignments[iDetElemNumber*fgNParCh+i]; }
986       }
987
988       // get transformation
989       TGeoHMatrix deltaGlobalTransform( DeltaTransform( lDetElemMisAlignment ) );
990
991       // update module
992       TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
993       TGeoHMatrix newGlobalTransform( AliMUONGeometryBuilder::Multiply( deltaGlobalTransform, globalTransform ) );
994       newDetElement->SetGlobalTransformation( newGlobalTransform );
995       newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
996
997       // Get matching old alignment and update current matrix accordingly
998       if( oldMisAlignArray )
999       {
1000
1001         const AliAlignObjMatrix* oldAlignObj(0);
1002         const int detElemId( detElement->GetId() );
1003         const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
1004         for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
1005         {
1006
1007           const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
1008           if( localAlignObj && localAlignObj->GetVolUID() == volId )
1009           {
1010             oldAlignObj = localAlignObj;
1011             break;
1012           }
1013
1014         }
1015
1016         // multiply
1017         if( oldAlignObj )
1018         {
1019
1020           TGeoHMatrix oldMatrix;
1021           oldAlignObj->GetMatrix( oldMatrix );
1022           deltaGlobalTransform.Multiply( &oldMatrix );
1023
1024         }
1025
1026       }
1027
1028       // Create misalignment matrix
1029       newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
1030
1031     }
1032
1033     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1034   }
1035
1036   return newGeometryTransformer;
1037
1038 }
1039
1040 //______________________________________________________________________
1041 void AliMUONAlignment::SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY )
1042 {
1043
1044   /// Set alignment resolution to misalign objects to be stored in CDB
1045   /// if rChId is > 0 set parameters for this chamber only, counting from 1
1046   TMatrixDSym mChCorrMatrix(6);
1047   mChCorrMatrix[0][0]=chResX*chResX;
1048   mChCorrMatrix[1][1]=chResY*chResY;
1049
1050   TMatrixDSym mDECorrMatrix(6);
1051   mDECorrMatrix[0][0]=deResX*deResX;
1052   mDECorrMatrix[1][1]=deResY*deResY;
1053
1054   AliAlignObjMatrix *alignMat = 0x0;
1055
1056   for( Int_t chId = 0; chId <= 9; ++chId )
1057   {
1058
1059     // skip chamber if selection is valid, and does not match
1060     if( rChId > 0 && chId+1 != rChId ) continue;
1061
1062     TString chName1;
1063     TString chName2;
1064     if (chId<4)
1065     {
1066
1067       chName1 = Form("GM%d",chId);
1068       chName2 = Form("GM%d",chId);
1069
1070     } else {
1071
1072       chName1 = Form("GM%d",4+(chId-4)*2);
1073       chName2 = Form("GM%d",4+(chId-4)*2+1);
1074
1075     }
1076
1077     for( int i=0; i<misAlignArray->GetEntries(); ++i )
1078     {
1079
1080       alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
1081       TString volName(alignMat->GetSymName());
1082       if((volName.Contains(chName1)&&
1083         ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
1084         (volName.Length()==volName.Index(chName1)+chName1.Length())))||
1085         (volName.Contains(chName2)&&
1086         ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
1087         (volName.Length()==volName.Index(chName2)+chName2.Length()))))
1088       {
1089
1090         volName.Remove(0,volName.Last('/')+1);
1091         if (volName.Contains("GM")) alignMat->SetCorrMatrix(mChCorrMatrix);
1092         else if (volName.Contains("DE")) alignMat->SetCorrMatrix(mDECorrMatrix);
1093
1094       }
1095
1096     }
1097
1098   }
1099
1100 }
1101
1102
1103 //_____________________________________________________
1104 void AliMUONAlignment::FillDetElemData( AliMUONVCluster* cluster )
1105 {
1106
1107   /// Get information of current detection element
1108   // get detector element number from Alice ID
1109   const Int_t detElemId = cluster->GetDetElemId();
1110   fDetElemNumber = GetDetElemNumber( detElemId );
1111
1112   // get detector element
1113   const AliMUONGeometryDetElement* detElement = fTransform->GetDetElement( detElemId );
1114
1115   /*
1116   get the global transformation matrix and store its inverse, in order to manually perform
1117   the global to Local transformations needed to calculate the derivatives
1118   */
1119   fGeoCombiTransInverse = detElement->GetGlobalTransformation()->Inverse();
1120
1121 }
1122
1123 //______________________________________________________________________
1124 void AliMUONAlignment::FillRecPointData( AliMUONVCluster* cluster )
1125 {
1126
1127   /// Get information of current cluster
1128   fClustPos[0] = cluster->GetX();
1129   fClustPos[1] = cluster->GetY();
1130   fClustPos[2] = cluster->GetZ();
1131
1132 }
1133
1134 //______________________________________________________________________
1135 void AliMUONAlignment::FillTrackParamData( AliMUONTrackParam* trackParam )
1136 {
1137
1138   /// Get information of current track at current cluster
1139   fTrackPos[0] = trackParam->GetNonBendingCoor();
1140   fTrackPos[1] = trackParam->GetBendingCoor();
1141   fTrackPos[2] = trackParam->GetZ();
1142   fTrackSlope[0] = trackParam->GetNonBendingSlope();
1143   fTrackSlope[1] = trackParam->GetBendingSlope();
1144
1145 }
1146
1147 //______________________________________________________________________
1148 Bool_t AliMUONAlignment::UnbiasTrackParamData( AliMUONTrackParam* trackParam ) const
1149 {
1150
1151   /**
1152   calculate unbiased track parameters at a given detector, that is,
1153   after taking out the contribution of the detector's cluster from the track
1154   */
1155
1156   // check track parameters
1157   if( !trackParam ) return kFALSE;
1158
1159   // Remove cluster contibution from smoothed track param
1160   TMatrixD smoothParameters( trackParam->GetSmoothParameters() );
1161   TMatrixD smoothCovariances( trackParam->GetSmoothCovariances() );
1162
1163   AliMUONVCluster* cluster = trackParam->GetClusterPtr();
1164   // p' = p + K(m - H*p)
1165   // K  = C H (-V + H C H^t)^-1
1166   // C' = C - K H C
1167   // where p,C are smoothed param,cov at cluster position
1168   // H is the matrix converting the state vector to measurement
1169   // V is the measurement cov.matrix
1170   // m is the measurement vector
1171   static TMatrixD H(2,5);
1172   H.Zero();
1173   H(0,0)=H(1,2) = 1.;
1174
1175   // (-Vk+H_k C^n_k H_k^T)^-1
1176   static TMatrixD df(2,2);
1177   df.Zero();
1178   df(0,0) = smoothCovariances(0,0) - cluster->GetErrX2();
1179   df(1,1) = smoothCovariances(2,2) - cluster->GetErrY2();
1180   df(0,1) = smoothCovariances(0,2);
1181   df(1,0) = smoothCovariances(2,0);
1182
1183   if (df.Determinant() != 0) df.Invert();
1184   else {
1185     AliInfo( "Determinant = 0\n" );
1186     return kFALSE;
1187   }
1188
1189   // gain matrix
1190   TMatrixD kTmp( smoothCovariances, TMatrixD::kMultTranspose, H );
1191   TMatrixD K(kTmp, TMatrixD::kMult, df);
1192
1193   TMatrixD dfc(2,1);
1194   dfc.Zero();
1195   dfc(0,0) = cluster->GetX()-smoothParameters(0,0);
1196   dfc(1,0) = cluster->GetY()-smoothParameters(2,0);
1197   TMatrixD tmp0(K,TMatrixD::kMult, dfc);
1198   smoothParameters += tmp0;
1199
1200   TMatrixD tmp1(K,   TMatrixD::kMult, H);
1201   TMatrixD tmp2(tmp1,TMatrixD::kMult, smoothCovariances);
1202   smoothCovariances -= tmp2;
1203
1204   // update track parameters
1205   trackParam->SetParameters( smoothParameters );
1206   trackParam->SetCovariances( smoothCovariances );
1207
1208   return kTRUE;
1209
1210 }
1211
1212 //______________________________________________________________________
1213 void AliMUONAlignment::LocalEquationX( void )
1214 {
1215   /// local equation along X
1216
1217   // 'inverse' (GlobalToLocal) rotation matrix
1218   const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
1219
1220   // local derivatives
1221   SetLocalDerivative( 0, r[0] );
1222   SetLocalDerivative( 1, r[0]*(fTrackPos[2] - fTrackPos0[2]) );
1223   SetLocalDerivative( 2, r[1] );
1224   SetLocalDerivative( 3, r[1]*(fTrackPos[2] - fTrackPos0[2]) );
1225
1226   // global derivatives
1227   /*
1228   alignment parameters are
1229   0: delta_x
1230   1: delta_y
1231   2: delta_phiz
1232   3: delta_z
1233   */
1234
1235   SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[0] );
1236   SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[1] );
1237
1238   if( fBFieldOn )
1239   {
1240
1241     // use local position for derivatives vs 'delta_phi_z'
1242     SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*fTrackPos[0] + r[0]*fTrackPos[1] );
1243
1244     // use local slopes for derivatives vs 'delta_z'
1245     SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope[0] + r[1]*fTrackSlope[1] );
1246
1247   } else {
1248
1249     // local copy of extrapolated track positions
1250     const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1251     const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
1252
1253     // use properly extrapolated position for derivatives vs 'delta_phi_z'
1254     SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*trackPosX + r[0]*trackPosY );
1255
1256     // use slopes at origin for derivatives vs 'delta_z'
1257     SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope0[0] + r[1]*fTrackSlope0[1] );
1258
1259   }
1260
1261   // store local equation
1262   fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0] );
1263
1264 }
1265
1266 //______________________________________________________________________
1267 void AliMUONAlignment::LocalEquationY( void )
1268 {
1269   /// local equation along Y
1270
1271   // 'inverse' (GlobalToLocal) rotation matrix
1272   const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
1273
1274   // store local derivatives
1275   SetLocalDerivative( 0, r[3] );
1276   SetLocalDerivative( 1, r[3]*(fTrackPos[2] - fTrackPos0[2] ) );
1277   SetLocalDerivative( 2, r[4] );
1278   SetLocalDerivative( 3, r[4]*(fTrackPos[2] - fTrackPos0[2] ) );
1279
1280   // set global derivatives
1281   SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[3]);
1282   SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[4]);
1283
1284   if( fBFieldOn )
1285   {
1286
1287     // use local position for derivatives vs 'delta_phi'
1288     SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*fTrackPos[0] + r[3]*fTrackPos[1]);
1289
1290     // use local slopes for derivatives vs 'delta_z'
1291     SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope[0]+r[4]*fTrackSlope[1] );
1292
1293   } else {
1294
1295     // local copy of extrapolated track positions
1296     const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1297     const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
1298
1299     // use properly extrapolated position for derivatives vs 'delta_phi'
1300     SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*trackPosX + r[3]*trackPosY );
1301
1302     // use slopes at origin for derivatives vs 'delta_z'
1303     SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope0[0]+r[4]*fTrackSlope0[1] );
1304
1305   }
1306
1307   // store local equation
1308   fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1] );
1309
1310 }
1311
1312 //_________________________________________________________________________
1313 TGeoCombiTrans AliMUONAlignment::DeltaTransform( const double *lMisAlignment) const
1314 {
1315   /// Get Delta Transformation, based on alignment parameters
1316
1317   // translation
1318   const TGeoTranslation deltaTrans( lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
1319
1320   // rotation
1321   TGeoRotation deltaRot;
1322   deltaRot.RotateZ(lMisAlignment[2]*180./TMath::Pi());
1323
1324   // combined rotation and translation.
1325   return TGeoCombiTrans(deltaTrans,deltaRot);
1326
1327 }
1328
1329 //______________________________________________________________________
1330 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value)
1331 {
1332   /// Constrain equation defined by par to value
1333   if( !fInitialized )
1334   { AliFatal( "Millepede is not initialized" ); }
1335
1336   fMillepede->SetGlobalConstraint(par, value);
1337 }
1338
1339 //______________________________________________________________________
1340 Bool_t AliMUONAlignment::DetElemIsValid( Int_t iDetElemId ) const
1341 {
1342   /// return true if given detector element is valid (and belongs to muon tracker)
1343   const Int_t iCh = iDetElemId/100;
1344   const Int_t iDet = iDetElemId%100;
1345   return ( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] );
1346 }
1347
1348 //______________________________________________________________________
1349 Int_t AliMUONAlignment::GetDetElemNumber( Int_t iDetElemId ) const
1350 {
1351   /// get det element number from ID
1352   // get chamber and element number in chamber
1353   const Int_t iCh = iDetElemId/100;
1354   const Int_t iDet = iDetElemId%100;
1355
1356   // make sure detector index is valid
1357   if( !( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] ) )
1358   { AliFatal( Form( "Invalid detector element id: %i", iDetElemId ) ); }
1359
1360   // add number of detectors up to this chamber
1361   return iDet + fgSNDetElemCh[iCh-1];
1362
1363 }
1364
1365 //______________________________________________________________________
1366 Int_t AliMUONAlignment::GetChamberId( Int_t iDetElemNumber ) const
1367 {
1368   /// get chamber (counting from 1) matching a given detector element id
1369   Int_t iCh( 0 );
1370   for( iCh=0; iCh<fgNCh; iCh++ )
1371   { if( iDetElemNumber < fgSNDetElemCh[iCh] ) break; }
1372
1373   return iCh;
1374 }
1375
1376 //______________________________________________________________________
1377 TString AliMUONAlignment::GetParameterMaskString( UInt_t mask ) const
1378 {
1379   TString out;
1380   if( mask & ParX ) out += "X";
1381   if( mask & ParY ) out += "Y";
1382   if( mask & ParZ ) out += "Z";
1383   if( mask & ParTZ ) out += "T";
1384   return out;
1385 }
1386
1387 //______________________________________________________________________
1388 TString AliMUONAlignment::GetSidesMaskString( UInt_t mask ) const
1389 {
1390   TString out;
1391   if( mask & SideTop ) out += "T";
1392   if( mask & SideLeft ) out += "L";
1393   if( mask & SideBottom ) out += "B";
1394   if( mask & SideRight ) out += "R";
1395   return out;
1396 }