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