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