MUON/AliMUONAlignment.cxx

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