MUON/AliMUONAlignment.cxx

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