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 }