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