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 "AliMUONGeometryStore.h"
  39 #include "AliMUONGeometryBuilder.h"
  40 #include "AliLog.h"
  41 #include "TSystem.h"
  42 #include "AliMUONConstants.h"
  43
  44 #include "AliMillepede.h"
  45
  46 ClassImp(AliMUONAlignment)
  47   Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
  48   Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
  49   Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
  50   Int_t AliMUONAlignment::fgNParCh = 3;
  51   Int_t AliMUONAlignment::fgNCh = 10;
  52   Int_t AliMUONAlignment::fgNSt = 5;
  53
  54 AliMUONAlignment::AliMUONAlignment() : TObject()
  55 {
  56   /// Default constructor, setting define alignment parameters
  57   fSigma[0] = 1.0e-1;
  58   fSigma[1] = 1.0e-2;
  59   fPhi = 0;
  60   fNGlobal = fgNDetElem*fgNParCh;
  61   fNLocal = 4;
  62   fNStdDev = 3;
  63
  64   fDoF[0] = kTRUE;  fDoF[1] = kTRUE;  fDoF[2] = kTRUE;
  65   fAllowVar[0] = 0.05;  fAllowVar[1] = 0.05;  fAllowVar[2] = 0.001;
  66   fStartFac = 16.;
  67   fResCut = 100.;
  68   fResCutInitial = 100.;
  69
  70   fBFieldOn = kTRUE;
  71 //   fBFieldOn = kFALSE;
  72
  73   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));
  74
  75   fMillepede = new AliMillepede();
  76
  77   Init(fNGlobal, fNLocal, fNStdDev);
  78
  79   ResetLocalEquation();
  80   AliInfo("Parameters initialized to zero");
  81
  82 }
  83
  84 AliMUONAlignment::~AliMUONAlignment() {
  85   /// Destructor
  86 }
  87
  88 void AliMUONAlignment::Init(Int_t nGlobal,  /* number of global paramers */
  89                             Int_t nLocal,   /* number of local parameters */
  90                             Int_t nStdDev   /* std dev cut */ )
  91 {
  92   /// Initialization of AliMillepede. Fix parameters, define constraints ...
  93   fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
  94
  95   Bool_t bStOnOff[5] = {kFALSE,kFALSE,kTRUE,kTRUE,kTRUE};
  96
  97   AllowVariations(bStOnOff);
  98
  99   // Fix parameters or add constraints here
 100   for (Int_t iSt=0; iSt<5; iSt++)
 101     if (!bStOnOff[iSt]) FixStation(iSt+1);
 102
 103   Bool_t bVarXYT[3] = {kFALSE,kTRUE,kFALSE};
 104   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
 105   ResetConstraints();
 106   AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
 107   bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
 108   bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
 109   AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
 110   bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
 111   AddConstraints(bStOnOff,bVarXYT);
 112
 113   // Set iterations
 114   if (fStartFac>1) fMillepede->SetIterations(fStartFac);
 115 }
 116
 117 void AliMUONAlignment::FixStation(Int_t iSt){
 118   /// Fix all detection elements of station iSt
 119   Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
 120   Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
 121   for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){
 122     FixParameter(i*fgNParCh+0, 0.0);
 123     FixParameter(i*fgNParCh+1, 0.0);
 124     FixParameter(i*fgNParCh+2, 0.0);
 125   }
 126 }
 127
 128 void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){
 129   /// Set non linear parameter flag selected stations and degrees of freedom
 130   for (Int_t i = 0; i < fgNDetElem; i++){
 131     Int_t iCh=0;
 132     for (iCh=1; iCh<=fgNCh; iCh++){
 133       if (i<fgSNDetElemCh[iCh-1]) break;
 134     }
 135     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
 136     if (iSt){
 137       if (lVarXYT[0]) { // X constraint
 138         SetNonLinear(i*fgNParCh+0);
 139       }
 140       if (lVarXYT[1]) { // Y constraint
 141         SetNonLinear(i*fgNParCh+1);
 142       }
 143       if (lVarXYT[2]) { // T constraint
 144         SetNonLinear(i*fgNParCh+2);
 145       }
 146     }
 147   }
 148 }
 149
 150 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT){
 151   /// Add constraint equations for selected stations and degrees of freedom
 152   for (Int_t i = 0; i < fgNDetElem; i++){
 153     Int_t iCh=0;
 154     for (iCh=1; iCh<=fgNCh; iCh++){
 155       if (i<fgSNDetElemCh[iCh-1]) break;
 156     }
 157     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
 158     if (iSt){
 159       if (lVarXYT[0]) { // X constraint
 160         fConstraintX[i*fgNParCh+0]=1.0;
 161       }
 162       if (lVarXYT[1]) { // Y constraint
 163         fConstraintY[i*fgNParCh+1]=1.0;
 164       }
 165       if (lVarXYT[2]) { // T constraint
 166         fConstraintP[i*fgNParCh+2]=1.0;
 167       }
 168     }
 169   }
 170   if (lVarXYT[0]) { // X constraint
 171     AddConstraint(fConstraintX,0.0);
 172   }
 173   if (lVarXYT[1]) { // Y constraint
 174     AddConstraint(fConstraintY,0.0);
 175   }
 176   if (lVarXYT[2]) { // T constraint
 177     AddConstraint(fConstraintP,0.0);
 178   }
 179 }
 180
 181 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t *lDetTLBR){
 182   /// Add constraint equations for selected stations, degrees of freedom detector half
 183   for (Int_t i = 0; i < fgNDetElem; i++){
 184     Int_t iCh=0;
 185     for (iCh=1; iCh<=fgNCh; iCh++){
 186       if (i<fgSNDetElemCh[iCh-1]) break;
 187     }
 188     Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
 189     if (iSt){
 190       if (lVarXYT[0]) { // X constraint
 191         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
 192         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
 193         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
 194         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
 195       }
 196       if (lVarXYT[1]) { // X constraint
 197         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
 198         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
 199         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
 200         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
 201       }
 202       if (lVarXYT[2]) { // X constraint
 203         if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
 204         if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
 205         if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
 206         if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
 207       }
 208     }
 209   }
 210   if (lVarXYT[0]) { // X constraint
 211     if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
 212     if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
 213     if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
 214     if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
 215   }
 216   if (lVarXYT[1]) { // Y constraint
 217     if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
 218     if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
 219     if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
 220     if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
 221   }
 222   if (lVarXYT[2]) { // T constraint
 223     if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
 224     if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
 225     if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
 226     if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
 227   }
 228 }
 229
 230 void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar){
 231   /// Set constrain equation for top half of spectrometer
 232   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
 233   if (lCh>=1 && lCh<=4){
 234     if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
 235       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
 236     }
 237   }
 238   if (lCh>=5 && lCh<=6){
 239     if (lDetElemNumber>=0&&lDetElemNumber<=9){
 240       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
 241     }
 242   }
 243   if (lCh>=7 && lCh<=10){
 244     if (lDetElemNumber>=0&&lDetElemNumber<=13){
 245       lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
 246     }
 247   }
 248 }
 249
 250 void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar){
 251   /// Set constrain equation for left half of spectrometer
 252   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
 253   if (lCh>=1 && lCh<=4){
 254     if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
 255       lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
 256     }
 257   }
 258   if (lCh>=5 && lCh<=6){
 259     if (lDetElemNumber>=5&&lDetElemNumber<=13){
 260       lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
 261     }
 262   }
 263   if (lCh>=7 && lCh<=10){
 264     if (lDetElemNumber>=7&&lDetElemNumber<=19){
 265       lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
 266     }
 267   }
 268 }
 269
 270 void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar){
 271   /// Set constrain equation for bottom half of spectrometer
 272   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
 273   if (lCh>=1 && lCh<=4){
 274     if (lDetElemNumber==2 && lDetElemNumber==3){ // From track crossings
 275       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
 276     }
 277   }
 278   if (lCh>=5 && lCh<=6){
 279     if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
 280         (lDetElemNumber==0)){
 281       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
 282     }
 283   }
 284   if (lCh>=7 && lCh<=10){
 285     if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
 286         (lDetElemNumber==0)){
 287       lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
 288     }
 289   }
 290 }
 291
 292 void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar){
 293   /// Set constrain equation for right half of spectrometer
 294   Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
 295   if (lCh>=1 && lCh<=4){
 296     if (lDetElemNumber==0 && lDetElemNumber==3){ // From track crossings
 297       lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
 298     }
 299   }
 300   if (lCh>=5 && lCh<=6){
 301     if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
 302         (lDetElemNumber>=14&&lDetElemNumber<=17)){
 303       lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
 304     }
 305   }
 306   if (lCh>=7 && lCh<=10){
 307     if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
 308         (lDetElemNumber>=20&&lDetElemNumber<=25)){
 309       lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
 310     }
 311   }
 312 }
 313
 314 void AliMUONAlignment::ResetConstraints(){
 315   /// Reset all constraint equations
 316   for (Int_t i = 0; i < fgNDetElem; i++){
 317       fConstraintX[i*fgNParCh+0]=0.0;
 318       fConstraintX[i*fgNParCh+1]=0.0;
 319       fConstraintX[i*fgNParCh+2]=0.0;
 320       fConstraintY[i*fgNParCh+0]=0.0;
 321       fConstraintY[i*fgNParCh+1]=0.0;
 322       fConstraintY[i*fgNParCh+2]=0.0;
 323       fConstraintP[i*fgNParCh+0]=0.0;
 324       fConstraintP[i*fgNParCh+1]=0.0;
 325       fConstraintP[i*fgNParCh+2]=0.0;
 326       fConstraintXT[i*fgNParCh+0]=0.0;
 327       fConstraintXT[i*fgNParCh+1]=0.0;
 328       fConstraintXT[i*fgNParCh+2]=0.0;
 329       fConstraintYT[i*fgNParCh+0]=0.0;
 330       fConstraintYT[i*fgNParCh+1]=0.0;
 331       fConstraintYT[i*fgNParCh+2]=0.0;
 332       fConstraintPT[i*fgNParCh+0]=0.0;
 333       fConstraintPT[i*fgNParCh+1]=0.0;
 334       fConstraintPT[i*fgNParCh+2]=0.0;
 335       fConstraintXL[i*fgNParCh+0]=0.0;
 336       fConstraintXL[i*fgNParCh+1]=0.0;
 337       fConstraintXL[i*fgNParCh+2]=0.0;
 338       fConstraintYL[i*fgNParCh+0]=0.0;
 339       fConstraintYL[i*fgNParCh+1]=0.0;
 340       fConstraintYL[i*fgNParCh+2]=0.0;
 341       fConstraintPL[i*fgNParCh+0]=0.0;
 342       fConstraintPL[i*fgNParCh+1]=0.0;
 343       fConstraintPL[i*fgNParCh+2]=0.0;
 344       fConstraintXB[i*fgNParCh+0]=0.0;
 345       fConstraintXB[i*fgNParCh+1]=0.0;
 346       fConstraintXB[i*fgNParCh+2]=0.0;
 347       fConstraintYB[i*fgNParCh+0]=0.0;
 348       fConstraintYB[i*fgNParCh+1]=0.0;
 349       fConstraintYB[i*fgNParCh+2]=0.0;
 350       fConstraintPB[i*fgNParCh+0]=0.0;
 351       fConstraintPB[i*fgNParCh+1]=0.0;
 352       fConstraintPB[i*fgNParCh+2]=0.0;
 353       fConstraintXR[i*fgNParCh+0]=0.0;
 354       fConstraintXR[i*fgNParCh+1]=0.0;
 355       fConstraintXR[i*fgNParCh+2]=0.0;
 356       fConstraintYR[i*fgNParCh+0]=0.0;
 357       fConstraintYR[i*fgNParCh+1]=0.0;
 358       fConstraintYR[i*fgNParCh+2]=0.0;
 359       fConstraintPR[i*fgNParCh+0]=0.0;
 360       fConstraintPR[i*fgNParCh+1]=0.0;
 361       fConstraintPR[i*fgNParCh+2]=0.0;
 362   }
 363 }
 364
 365 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
 366   /// Constrain equation defined by par to value
 367   fMillepede->SetGlobalConstraint(par, value);
 368   AliInfo("Adding constraint");
 369 }
 370
 371 void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
 372   /// Initialize global parameters with par array
 373   fMillepede->SetGlobalParameters(par);
 374   AliInfo("Init Global Parameters");
 375 }
 376
 377 void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
 378   /// Parameter iPar is encourage to vary in [-value;value].
 379   /// If value == 0, parameter is fixed
 380   fMillepede->SetParSigma(iPar, value);
 381   if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
 382 }
 383
 384 void AliMUONAlignment::ResetLocalEquation()
 385 {
 386   /// Reset the derivative vectors
 387   for(int i=0; i<fNLocal; i++) {
 388     fLocalDerivatives[i] = 0.0;
 389   }
 390   for(int i=0; i<fNGlobal; i++) {
 391     fGlobalDerivatives[i] = 0.0;
 392   }
 393 }
 394
 395 void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) {
 396   /// Set allowed variation for selected stations based on fDoF and fAllowVar
 397   for (Int_t iSt=1; iSt<=5; iSt++) {
 398     if (bStOnOff[iSt-1]) {
 399       Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
 400       Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
 401       for (int i=0; i<fgNParCh; i++) {
 402         AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
 403         if (fDoF[i]) {
 404           for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
 405             FixParameter(j*fgNParCh+i, fAllowVar[i]);
 406           }
 407         }
 408       }
 409     }
 410   }
 411 }
 412
 413 void AliMUONAlignment::SetNonLinear(Int_t iPar  /* set non linear flag */ ) {
 414   /// Set nonlinear flag for parameter iPar
 415   fMillepede->SetNonLinear(iPar);
 416   AliInfo(Form("Parameter %i set to non linear", iPar));
 417 }
 418
 419 void AliMUONAlignment::LocalEquationX() {
 420   /// Define local equation for current track and hit in x coor. measurement
 421   // set local derivatives
 422   SetLocalDerivative(0, fCosPhi);
 423   SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
 424   SetLocalDerivative(2, fSinPhi);
 425   SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
 426
 427   // set global derivatives
 428   SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -1.);
 429   SetGlobalDerivative(fDetElemNumber*fgNParCh+1,  0.);
 430   if (fBFieldOn){
 431     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
 432                         -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
 433                         +fCosPhi*(fTrackPos[1]-fDetElemPos[1]));
 434   }
 435   else {
 436     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
 437                         -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*
 438                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
 439                         +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*
 440                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
 441   }
 442
 443   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
 444 }
 445
 446 void AliMUONAlignment::LocalEquationY() {
 447   /// Define local equation for current track and hit in y coor. measurement
 448   // set local derivatives
 449   SetLocalDerivative(0,-fSinPhi);
 450   SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
 451   SetLocalDerivative(2, fCosPhi);
 452   SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
 453
 454   // set global derivatives
 455   SetGlobalDerivative(fDetElemNumber*fgNParCh+0,  0.);
 456   SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -1.);
 457   if (fBFieldOn){
 458   SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
 459                       -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
 460                       -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
 461   }
 462   else {
 463     SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
 464                         -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*
 465                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
 466                         -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*
 467                                   (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
 468   }
 469
 470   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
 471 }
 472
 473 void AliMUONAlignment::FillRecPointData() {
 474   /// Get information of current hit
 475   fClustPos[0] = fRecHit->GetNonBendingCoor();
 476   fClustPos[1] = fRecHit->GetBendingCoor();
 477   fClustPos[2] = fRecHit->GetZ();
 478   fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
 479                             fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
 480 }
 481
 482 void AliMUONAlignment::FillTrackParamData() {
 483   /// Get information of current track at current hit
 484   fTrackPos[0] = fTrackParam->GetNonBendingCoor();
 485   fTrackPos[1] = fTrackParam->GetBendingCoor();
 486   fTrackPos[2] = fTrackParam->GetZ();
 487   fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
 488   fTrackSlope[1] = fTrackParam->GetBendingSlope();
 489   fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
 490                             fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
 491 }
 492
 493 void AliMUONAlignment::FillDetElemData() {
 494   /// Get information of current detection element
 495   Double_t lDetElemLocX = 0.;
 496   Double_t lDetElemLocY = 0.;
 497   Double_t lDetElemLocZ = 0.;
 498   fDetElemId = fRecHit->GetDetElemId();
 499   fDetElemNumber = fDetElemId%100;
 500   for (int iCh=0; iCh<fDetElemId/100-1; iCh++){
 501     fDetElemNumber += fgNDetElemCh[iCh];
 502   }
 503   fTransform->Local2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
 504                            fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
 505   if (fDetElemId/100 < 5){
 506     fSigma[0] = 3.0e-2;
 507     fSigma[1] = 3.0e-2;
 508   }
 509   else {
 510     fSigma[0] = 1.0e-1;
 511     fSigma[1] = 1.0e-2;
 512   }
 513 }
 514
 515 void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) {
 516   /// Process track; Loop over hits and set local equations
 517   fTrack = track;
 518   // get tclones arrays.
 519   fTrackParamAtHit = fTrack->GetTrackParamAtHit();
 520   fHitForRecAtHit = fTrack->GetHitForRecAtHit();
 521
 522   // get size of arrays
 523   Int_t nTrackParam = fTrackParamAtHit->GetEntries();
 524   Int_t nHitForRec = fHitForRecAtHit->GetEntries();
 525   AliInfo(Form("Number of track param entries : %i ", nTrackParam));
 526   AliInfo(Form("Number of hit for rec entries : %i ", nHitForRec));
 527
 528   for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
 529     fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
 530     fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
 531     if (!fRecHit || !fTrackParam) continue;
 532     // fill local variables for this position --> one measurement
 533     FillDetElemData();
 534     FillRecPointData();
 535     FillTrackParamData();
 536 //     if (fDetElemId<500) continue;
 537     fTrackPos0[0]      = fTrackPos[0];
 538     fTrackPos0[1]      = fTrackPos[1];
 539     fTrackPos0[2]      = fTrackPos[2];
 540     fTrackSlope0[0] = fTrackSlope[0];
 541     fTrackSlope0[1] = fTrackSlope[1];
 542     break;
 543   }
 544
 545   for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
 546     // and get new pointers
 547     fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
 548     fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
 549     if (!fRecHit || !fTrackParam) continue;
 550     // fill local variables for this position --> one measurement
 551     FillDetElemData();
 552     FillRecPointData();
 553     FillTrackParamData();
 554 //     if (fDetElemId<500) continue;
 555     AliDebug(1,Form("cluster: %i", iHit));
 556     AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
 557     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));
 558
 559     AliDebug(1,Form("Track Parameter: %i", iHit));
 560     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]));
 561     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]));
 562
 563     fCosPhi = TMath::Cos(fPhi);
 564     fSinPhi = TMath::Sin(fPhi);
 565     if (fBFieldOn){
 566       fMeas[0] = fTrackPos[0] - fClustPos[0];
 567       fMeas[1] = fTrackPos[1] - fClustPos[1];
 568     }
 569     else {
 570       fMeas[0] = - fClustPos[0];
 571       fMeas[1] = - fClustPos[1];
 572     }
 573     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]));
 574     // Set local equations
 575     LocalEquationX();
 576     LocalEquationY();
 577   }
 578 }
 579
 580 void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
 581   /// Call local fit for this tracks
 582   Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
 583   if (iRes && !lSingleFit) {
 584     fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
 585   }
 586 }
 587
 588 void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
 589   /// Call global fit; Global parameters are stored in parameters
 590   fMillepede->GlobalFit(parameters,errors,pulls);
 591   AliInfo("Done fitting global parameters!");
 592   for (int iGlob=0; iGlob<fgNDetElem; iGlob++){
 593     printf("%d\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+2]);
 594   }
 595 }
 596
 597 Double_t AliMUONAlignment::GetParError(Int_t iPar) {
 598   /// Get error of parameter iPar
 599   Double_t lErr = fMillepede->GetParError(iPar);
 600   return lErr;
 601 }
 602
 603 void AliMUONAlignment::PrintGlobalParameters() {
 604   /// Print global parameters
 605   fMillepede->PrintGlobalParameters();
 606 }
 607
 608 //_________________________________________________________________________
 609 TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *detElemMisAlignment) const
 610 {
 611   /// Realign given transformation by given misalignment and return the misaligned transformation
 612
 613   Double_t cartMisAlig[3] = {0,0,0};
 614   Double_t angMisAlig[3] = {0,0,0};
 615   const Double_t *trans = transform.GetTranslation();
 616   TGeoRotation *rot;
 617   // check if the rotation we obtain is not NULL
 618   if (transform.GetRotation()) {
 619     rot = transform.GetRotation();
 620   }
 621   else {
 622     rot = new TGeoRotation("rot");
 623   }                     // default constructor.
 624
 625   cartMisAlig[0] = -detElemMisAlignment[0];
 626   cartMisAlig[1] = -detElemMisAlignment[1];
 627   angMisAlig[2] = -detElemMisAlignment[2]*180./TMath::Pi();
 628
 629   TGeoTranslation newTrans(cartMisAlig[0] + trans[0], cartMisAlig[1] + trans[1], cartMisAlig[2] + trans[2]);
 630
 631   rot->RotateX(angMisAlig[0]);
 632   rot->RotateY(angMisAlig[1]);
 633   rot->RotateZ(angMisAlig[2]);
 634
 635   return TGeoCombiTrans(newTrans, *rot);
 636 }
 637
 638 //______________________________________________________________________
 639 AliMUONGeometryTransformer *
 640 AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
 641                             double *misAlignments, Bool_t verbose)
 642
 643 {
 644   /////////////////////////////////////////////////////////////////////
 645   //   Takes the internal geometry module transformers, copies them
 646   // and gets the Detection Elements from them.
 647   // Takes misalignment parameters and applies these
 648   // to the local transform of the Detection Element
 649   // Obtains the global transform by multiplying the module transformer
 650   // transformation with the local transformation
 651   // Applies the global transform to a new detection element
 652   // Adds the new detection element to a new module transformer
 653   // Adds the new module transformer to a new geometry transformer
 654   // Returns the new geometry transformer
 655
 656   Double_t lDetElemMisAlignment[3] = {0.,0.,0.};
 657
 658   AliMUONGeometryTransformer *newGeometryTransformer =
 659     new AliMUONGeometryTransformer(kTRUE);
 660   for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++)
 661     {                           // module transformers
 662
 663       const AliMUONGeometryModuleTransformer *kModuleTransformer =
 664         transformer->GetModuleTransformer(iMt, true);
 665
 666       AliMUONGeometryModuleTransformer *newModuleTransformer =
 667         new AliMUONGeometryModuleTransformer(iMt);
 668       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
 669
 670       TGeoCombiTrans moduleTransform =
 671         TGeoCombiTrans(*kModuleTransformer->GetTransformation());
 672       TGeoCombiTrans *newModuleTransform = new TGeoCombiTrans(moduleTransform);
 673               // same module transform as the previous one
 674               // no mis align object created
 675       newModuleTransformer->SetTransformation(moduleTransform);
 676
 677       AliMUONGeometryStore *detElements =
 678         kModuleTransformer->GetDetElementStore();
 679
 680       if (verbose)
 681         AliInfo(Form
 682                 ("%i DEs in old GeometryStore  %i",
 683                  detElements->GetNofEntries(), iMt));
 684       Int_t iBase = !iMt ? 0 : fgSNDetElemCh[iMt-1];
 685       for (Int_t iDe = 0; iDe < detElements->GetNofEntries(); iDe++)
 686         {                       // detection elements.
 687           AliMUONGeometryDetElement *detElement =
 688             (AliMUONGeometryDetElement *) detElements->GetEntry(iDe);
 689           if (!detElement)
 690             AliFatal("Detection element not found.");
 691
 692           /// make a new detection element
 693           AliMUONGeometryDetElement *newDetElement =
 694             new AliMUONGeometryDetElement(detElement->GetId(),
 695                                           detElement->GetVolumePath());
 696           for (int i=0; i<fgNParCh; i++) {
 697             AliInfo(Form("iMt %i, iBase %i, iDe %i, iPar %i",iMt, iBase, iDe, (iBase+iDe)*fgNParCh+i));
 698             lDetElemMisAlignment[i] =
 699               (iMt<fgNCh) ? misAlignments[(iBase+iDe)*fgNParCh+i] : 0.;
 700           }
 701           // local transformation of this detection element.
 702           TGeoCombiTrans localTransform
 703             = TGeoCombiTrans(*detElement->GetLocalTransformation());
 704           TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
 705           newDetElement->SetLocalTransformation(newLocalTransform);
 706
 707           // global transformation
 708           TGeoHMatrix newGlobalTransform =
 709             AliMUONGeometryBuilder::Multiply(*newModuleTransform,
 710                                               newLocalTransform);
 711           newDetElement->SetGlobalTransformation(newGlobalTransform);
 712
 713           // add this det element to module
 714           newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
 715                                                           newDetElement);
 716           // Get delta transformation:
 717           // Tdelta = Tnew * Told.inverse
 718           TGeoHMatrix  deltaTransform
 719             = AliMUONGeometryBuilder::Multiply(
 720                 newGlobalTransform,
 721                 detElement->GetGlobalTransformation()->Inverse());
 722
 723           // Create mis alignment matrix
 724           newGeometryTransformer
 725             ->AddMisAlignDetElement(detElement->GetId(), deltaTransform);
 726         }
 727       if (verbose)
 728         AliInfo(Form("Added module transformer %i to the transformer", iMt));
 729       newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
 730     }
 731   return newGeometryTransformer;
 732 }
 733