rewrote AliMUONAlignment to - use AliMillePede2 instead of AliMillepede - use AliMill...
authorhupereir <pereira@hep.saclay.cea.fr>
Wed, 9 Apr 2014 07:19:14 +0000 (09:19 +0200)
committerlaphecet <laurent.aphecetche@subatech.in2p3.fr>
Mon, 7 Jul 2014 13:28:24 +0000 (15:28 +0200)
AddTaskMuonAlignment has been modified accordingly and the arguments simplified

MUON/AddTaskMuonAlignment.C
MUON/AliMUONAlignment.cxx
MUON/AliMUONAlignment.h
MUON/AliMUONAlignmentTask.cxx
MUON/AliMUONAlignmentTask.h

index 9e8dfb7..837b2d8 100644 (file)
@@ -3,65 +3,71 @@
 /// \brief Macro to add an AliMUONAlignmentTask to an analysis train
 ///
 /// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
+/// \author Hugo Pereira Da Costa, CEA/Saclay - Irfu/SPhN
 
-AliMUONAlignmentTask4 *AddTaskMuonAlignment(bool readrecs, bool doalign, bool writerecs, const char *name = "AliMUONAlignmentTask", const char *newalignocdb = "local://ReAlignOCDB", const char *oldalignocdb = "alien://Folder=/alice/data/2010/OCDB", const char *defaultocdb = "alien://Folder=/alice/data/2011/OCDB", const char *geofilename = "geometry.root")
+AliMUONAlignmentTask *AddTaskMuonAlignment(
+  TString oldAlignmentOCDB,
+  TString newAlignmentOCDB,
+  Bool_t doAlignment = kTRUE,
+  Bool_t writeRecords = kTRUE,
+  Bool_t readRecords = kFALSE
+ )
 {
-/// Creates a Muon Alignment task and adds it to the analysis manager.
 
-       // Get the pointer to the existing analysis manager via the static access method.
-       //==============================================================================
-       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-       if (!mgr) {
-               ::Error("AddTaskMuonAlignment", "No analysis manager to connect to.");
-               return NULL;
-       }   
-   
-       // This task requires an ESD input handler.
-       // Check this using the analysis manager.
-       //===============================================================================
-       TString type = mgr->GetInputEventHandler()->GetDataType();
-       if (readrecs) {         
-               if (!type.Contains("AOD")) {
-                       ::Error("AddTaskMuonAlignment", "Alignment task is asked to read records but input data type is not AOD! Game over!");
-                       return NULL;
-               } 
-       } else {
-               if (!type.Contains("ESD")) {
-                       ::Error("AddTaskMuonAlignment", "Alignment task is not asked to read records but input data type is not ESD! Game over!");
-                       return NULL;
-               } 
-       }
+  /// Creates a Muon Alignment task and adds it to the analysis manager.
 
-       // Create the task, add it to the manager and configure it.
-       //===========================================================================   
-       // Muons
-       AliMUONAlignmentTask *muonalign = new AliMUONAlignmentTask(name, newalignocdb, oldalignocdb,  defaultocdb, geofilename);
-       muonalign->SetLoadOCDBOnce(true);
-       muonalign->SetReadRecords(readrecs);
-  muonalign->SetDoAlignment(doalign);
-  muonalign->SetWriteRecords(writerecs);
-       mgr->AddTask(muonalign);
-       
-//    // Cuts on primary tracks
-//    AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard");
-//    esdTrackCutsL->SetMinNClustersTPC(50);
+  // Get the pointer to the existing analysis manager via the static access method.
+  AliAnalysisManager *analysisManager = AliAnalysisManager::GetAnalysisManager();
+  if( !analysisManager )
+  {
+    ::Error("AddTaskMuonAlignment", "No analysis manager to connect to.");
+    return NULL;
+  }
 
-//    AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");
-//    trackFilter->AddCuts(esdTrackCutsL);
+  // get input event handler and check type
+  TString type = analysisManager->GetInputEventHandler()->GetDataType();
+  if( readRecords )
+  {
 
-//    muonalign->SetTrackFilter(trackFilter);
+    // when reading records, AOD are required
+    if (!type.Contains( "AOD" ) )
+    {
+      Error("AddTaskMuonRefit", "AOD input handler needed!");
+      return NULL;
+    }
 
+  } else {
 
-       // Create ONLY the output containers for the data produced by the task.
-       // Get and connect other common input/output containers via the manager as below
-       //==============================================================================
-       mgr->ConnectInput(muonalign,  0, mgr->GetCommonInputContainer());
-       if (doalign) {
-               AliAnalysisDataContainer *listOut = mgr->CreateContainer("output", TList::Class(), AliAnalysisManager::kOutputContainer, "measShifts.root");
-               mgr->ConnectOutput(muonalign,  1, listOut);
-       }
-       if (writerecs) {
-               mgr->ConnectOutput(muonalign,  0, mgr->GetCommonOutputContainer());
-       }
-       return muonalign;
-}   
+    // ESDs are required otherwise
+    if( !type.Contains( "ESD" ) )
+    {
+      Error("AddTaskMuonRefit", "AOD input handler needed!");
+      return NULL;
+    }
+
+  }
+
+  // Create the task, add it to the manager and configure it.
+  AliMUONAlignmentTask *muonAlign = new AliMUONAlignmentTask( "AliMUONAlignmentTask" );
+  muonAlign->SetOldAlignStorage( oldAlignmentOCDB );
+  muonAlign->SetNewAlignStorage( newAlignmentOCDB );
+  muonAlign->SetLoadOCDBOnce( kTRUE );
+  muonAlign->SetReadRecords( readRecords );
+  muonAlign->SetDoAlignment( doAlignment );
+  muonAlign->SetWriteRecords( writeRecords );
+  muonAlign->SetMergeAlignmentCDBs( kTRUE );
+  muonAlign->SetUnbias( kFALSE );
+
+  analysisManager->AddTask(muonAlign);
+
+  // connect input
+  analysisManager->ConnectInput(muonAlign,  0, analysisManager->GetCommonInputContainer());
+
+  // when writting records, also connect output
+  if( writeRecords )
+  { analysisManager->ConnectOutput(muonAlign,  0, analysisManager->GetCommonOutputContainer()); }
+
+  // return created task
+  return muonAlign;
+
+}
index 595d470..d1c6546 100644 (file)
@@ -13,7 +13,7 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id$ */
+/* $Id: AliMUONAlignment.cxx 51000 2011-08-08 17:58:17Z ivana $ */
 
 //-----------------------------------------------------------------------------
 /// \class AliMUONAlignment
@@ -35,7 +35,7 @@
 #include "AliMUONGeometryModuleTransformer.h"
 #include "AliMUONGeometryDetElement.h"
 #include "AliMUONGeometryBuilder.h"
-#include "AliMillepede.h"
+#include "AliMillePede2.h"
 
 #include "AliMpExMap.h"
 #include "AliMpExMapIterator.h"
 #include "AliAlignObjMatrix.h"
 #include "AliLog.h"
 
-#include "TMath.h"
-#include "TMatrixDSym.h"
-#include "TClonesArray.h"
+#include <TMath.h>
+#include <TMatrixDSym.h>
+#include <TClonesArray.h>
+#include <TGraphErrors.h>
 
 /// \cond CLASSIMP
 ClassImp(AliMUONAlignment)
@@ -53,1279 +54,838 @@ ClassImp(AliMUONAlignment)
 
 //_____________________________________________________________________
 // static variables
-Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
-Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
-Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
-Int_t AliMUONAlignment::fgNParCh = 4;
-Int_t AliMUONAlignment::fgNTrkMod = 16;
-Int_t AliMUONAlignment::fgNCh = 10;
-Int_t AliMUONAlignment::fgNSt = 5;
+const Int_t AliMUONAlignment::fgNDetElemCh[AliMUONAlignment::fgNCh] = { 4, 4, 4, 4, 18, 18, 26, 26, 26, 26 };
+const Int_t AliMUONAlignment::fgSNDetElemCh[AliMUONAlignment::fgNCh+1] = { 0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156 };
+
+//_____________________________________________________________________
+/// self initialized array, used for adding constraints
+class Array
+{
+
+  public:
+
+  /// contructor
+  Array( void )
+  {
+    for( Int_t i=0; i < AliMUONAlignment::fNGlobal; ++i )
+    { values[i] = 0; }
+  }
+
+  /// array
+  Double_t values[AliMUONAlignment::fNGlobal];
+
+  private:
+
+  /// Not implemented
+  Array(const Array& );
+
+  /// Not implemented
+  Array&  operator = (const Array& );
+
+};
 
 //_____________________________________________________________________
 AliMUONAlignment::AliMUONAlignment()
   : TObject(),
-    fBFieldOn(kTRUE),
-    fStartFac(256.),
-    fResCutInitial(100.),
-    fResCut(100.),
-    fMillepede(0),
-    fTrack(0),
-    fCluster(0),
-    fTrackParam(0),
-    fNGlobal(fgNDetElem*fgNParCh),
-    fNLocal(4),
-    fNStdDev(3),
-    fDetElemId(0),
-    fDetElemNumber(0),
-    fPhi(0.),
-    fCosPhi(1.),
-    fSinPhi(0.),
+    fInitialized( kFALSE ),
+    fRunNumber( 0 ),
+    fBFieldOn( kTRUE ),
+    fStartFac( 256 ),
+    fResCutInitial( 100 ),
+    fResCut( 100 ),
+    fMillepede( 0L ),
+    fCluster( 0L ),
+    fNStdDev( 3 ),
+    fDetElemNumber( 0 ),
+    fUnbias( kFALSE ),
     fTrackRecord(),
-    fTransform(0)
+    fTransform( 0 ),
+    fGeoCombiTransInverse()
 {
-  /// Default constructor      
+  /// constructor
   fSigma[0] = 1.5e-1;
   fSigma[1] = 1.0e-2;
 
-  AliInfo(Form("fSigma[0]: %f\t fSigma[1]: %f",fSigma[0],fSigma[1]));
-
-  fDoF[0] = kTRUE;  fDoF[1] = kTRUE;  fDoF[2] = kTRUE;  fDoF[3] = kTRUE;
-  fAllowVar[0] = 0.05;  fAllowVar[1] = 0.05;  fAllowVar[2] = 0.001;  fAllowVar[3] = 0.5;
+  // default allowed variations
+  fAllowVar[0] = 0.5;  // x
+  fAllowVar[1] = 0.5;  // y
+  fAllowVar[2] = 0.01; // phi_z
+  fAllowVar[3] = 5;    // z
 
-  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));
+  // initialize millepede
+  fMillepede = new AliMillePede2();
 
-  fMillepede = new AliMillepede();
+  // initialize degrees of freedom
+  // by default all parameters are free
+  for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
+  { fGlobalParameterStatus[iPar] = kFreeParId; }
 
-  Init(fNGlobal, fNLocal, fNStdDev);
+  // initialize local equations
+  for(int i=0; i<fNLocal; ++i )
+  { fLocalDerivatives[i] = 0.0; }
 
-  ResetLocalEquation();
-  AliInfo("Parameters initialized to zero");
+  for(int i=0; i<fNGlobal; ++i )
+  { fGlobalDerivatives[i] = 0.0; }
 
 }
 
 //_____________________________________________________________________
-AliMUONAlignment::AliMUONAlignment(TRootIOCtor* /*dummy*/)
-: TObject(),
-fBFieldOn(kFALSE),
-fStartFac(0.),
-fResCutInitial(0.),
-fResCut(0.),
-fMillepede(0),
-fTrack(0),
-fCluster(0),
-fTrackParam(0),
-fNGlobal(0),
-fNLocal(0),
-fNStdDev(0),
-fDetElemId(0),
-fDetElemNumber(0),
-fPhi(0.),
-fCosPhi(0.),
-fSinPhi(0.),
-fTrackRecord(),
-fTransform(0)
+AliMUONAlignment::~AliMUONAlignment()
 {
-  /// Root IO constructor
-  ResetConstraints();
-  for (Int_t iCh=0; iCh<10; iCh++) { 
-    fChOnOff[iCh] = kFALSE; 
-  }
-  fSpecLROnOff[0] = kFALSE; fSpecLROnOff[1] = kFALSE;
-  for (Int_t iDoF=0; iDoF<4; iDoF++) {
-    fDoF[iDoF] = kFALSE;
-    fAllowVar[iDoF] = kFALSE;
-  }
-  for (Int_t i=0; i<3; i++) {
-    fClustPos[i] = 0;
-    fClustPosLoc[i] = 0;
-    fTrackPos0[i] = 0;  
-    fTrackPos[i] = 0;   
-    fTrackPosLoc[i] = 0;
-       fDetElemPos[i] = 0;
-  }
-  fTrackSlope0[0]=0;   fTrackSlope0[1]=0;
-  fTrackSlope[0]=0;    fTrackSlope[1]=0;
-  fMeas[0]=0;  fMeas[1]=0;
-  fSigma[0]=0; fSigma[1]=0;
-  for (Int_t iLPar=0; iLPar<4; iLPar++) {
-    fLocalDerivatives[iLPar]=0;        
-  }
-  for (Int_t iGPar=0; iGPar<624; iGPar++) {
-    fGlobalDerivatives[iGPar]=0;       
-  }
+  /// destructor
 }
 
 //_____________________________________________________________________
-AliMUONAlignment::~AliMUONAlignment()
-{}
-
-//_____________________________________________________________________
-void AliMUONAlignment::Init(
-  Int_t nGlobal,  /* number of global paramers */
-  Int_t nLocal,   /* number of local parameters */
-  Int_t nStdDev   /* std dev cut */ )
+void AliMUONAlignment::Init( void )
 {
-  /// Initialization of AliMillepede. Fix parameters, define constraints ...
-  fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
-
-  //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-  //  Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-  //  Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
 
-  //  AllowVariations(bChOnOff);
+  /// initialize
+  /**
+  initialize millipede
+  must be called after necessary detectors have been fixed,
+  but before constrains are added and before global parameters initial value are set
+  */
+  if( fInitialized )
+  { AliFatal( "Millepede already initialized" ); }
+
+  // assign proper groupID to free parameters
+  Int_t nGlobal = 0;
+  for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
+  {
 
-  // Fix parameters or add constraints here
-  //   for (Int_t iSt=0; iSt<5; iSt++)
-  //   { if (!bStOnOff[iSt]) FixStation(iSt+1); }
+    if( fGlobalParameterStatus[iPar] == kFixedParId )
+    {
+      // fixed parameters are left unchanged
+      continue;
 
-  //   for (Int_t iCh=0; iCh<10; iCh++)
-  //   { if (!bChOnOff[iCh]) FixChamber(iCh+1); }
+    } else if( fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId ) {
 
-//   FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
+      // free parameters or first element of group are assigned a new group id
+      fGlobalParameterStatus[iPar] = nGlobal++;
+      continue;
 
-  ResetConstraints();
+    } else if( fGlobalParameterStatus[iPar] < kGroupBaseId ) {
 
-  // Define global constrains to be applied
-  // X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
-  Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-  Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
-  //  AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
+      // get detector element id from status, get chamber parameter id
+      const Int_t iDeBase( kGroupBaseId - 1 - fGlobalParameterStatus[iPar] );
+      const Int_t iParBase = iPar%fgNParCh;
 
-  // Other possible way to add constrains
-  bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
-  bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
-  //  AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
+      // check
+      if( iDeBase < 0 || iDeBase >= iPar/fgNParCh )
+      { AliFatal( Form( "Group for parameter index %i has wrong base detector element: %i", iPar, iDeBase ) ); }
 
-  bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE;
-  //  AddConstraints(bStOnOff,bVarXYT);
+      // assign identical group id to current
+      fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase*fgNParCh + iParBase];
+      AliInfo( Form( "Parameter %i grouped to detector %i (%s)", iPar, iDeBase, GetParameterMaskString( 1<<iParBase ).Data() ) );
 
-  // Set iterations
-  if (fStartFac>1) fMillepede->SetIterations(fStartFac);
-}
+    } else AliFatal( Form( "Unrecognized parameter status for index %i: %i", iPar, fGlobalParameterStatus[iPar] ) );
 
-//_____________________________________________________________________
-void AliMUONAlignment::FixStation(Int_t iSt)
-{
-  /// Fix all detection elements of station iSt
-  Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
-  Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
-  for (Int_t i = iDetElemFirst; i < iDetElemLast; i++)
-  {
-    FixParameter(i*fgNParCh+0, 0.0);
-    FixParameter(i*fgNParCh+1, 0.0);
-    FixParameter(i*fgNParCh+2, 0.0);
-    FixParameter(i*fgNParCh+3, 0.0);
   }
 
-}
+  AliInfo( Form( "Free Parameters: %i out of %i", nGlobal, fNGlobal ) );
 
-//_____________________________________________________________________
-void AliMUONAlignment::FixChamber(Int_t iCh)
-{
-  /// Fix all detection elements of chamber iCh
-  Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0;
-  Int_t iDetElemLast = fgSNDetElemCh[iCh-1];
-  for (Int_t i = iDetElemFirst; i < iDetElemLast; i++)
-  {
-    FixParameter(i*fgNParCh+0, 0.0);
-    FixParameter(i*fgNParCh+1, 0.0);
-    FixParameter(i*fgNParCh+2, 0.0);
-    FixParameter(i*fgNParCh+3, 0.0);
-  }
-}
+  // initialize millepede
+  fMillepede->InitMille( fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus );
+  fInitialized = kTRUE;
 
-//_____________________________________________________________________
-void AliMUONAlignment::FixDetElem(Int_t iDetElemId, TString sVarXYT)
-{
+  // some debug output
+  for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
+  {  AliInfo( Form( "fAllowVar[%i]= %f", iPar, fAllowVar[iPar] ) ); }
 
-  /// Fix a given detection element
-  Int_t iDetElemNumber = iDetElemId%100;
-  for (int iCh=0; iCh<iDetElemId/100-1; iCh++)
+  // set allowed variations for all parameters
+  for( Int_t iDet = 0; iDet < fgNDetElem; ++iDet )
   {
-    iDetElemNumber += fgNDetElemCh[iCh];
+    for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
+    { fMillepede->SetParSigma( iDet*fgNParCh + iPar, fAllowVar[iPar] ); }
   }
 
-  if (sVarXYT.Contains("X"))
-  {
-    // X constraint
-    FixParameter(iDetElemNumber*fgNParCh+0, 0.0);
-  }
-
-  if (sVarXYT.Contains("Y"))
-  {
-    // Y constraint
-    FixParameter(iDetElemNumber*fgNParCh+1, 0.0);
-  }
-
-  if (sVarXYT.Contains("T"))
-  {
-    // T constraint
-    FixParameter(iDetElemNumber*fgNParCh+2, 0.0);
-  }
-
-  if (sVarXYT.Contains("Z"))
-  {
-    // T constraint
-    FixParameter(iDetElemNumber*fgNParCh+3, 0.0);
-  }
+  // Set iterations
+  if (fStartFac>1) fMillepede->SetIterations(fStartFac);
 
 }
 
-//_____________________________________________________________________
-void AliMUONAlignment::FixHalfSpectrometer(const Bool_t *lChOnOff, const Bool_t *lSpecLROnOff)
+//_____________________________________________________
+AliMillePedeRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight )
 {
+  /// process track for alignment minimization
+  /**
+  returns the alignment records for this track.
+  They can be stored in some output for later reprocessing.
+  */
 
-  /// Fix left or right detector
-  for (Int_t i = 0; i < fgNDetElem; i++)
-  {
+  // reset track records
+  fTrackRecord.Reset();
+  if( fMillepede->GetRecord() ) fMillepede->GetRecord()->Reset();
 
-    Int_t iCh=0;
-    for (iCh=1; iCh<=fgNCh; iCh++)
-    { if (i<fgSNDetElemCh[iCh-1]) break; }
+  // get number of track parameters
+  Int_t nTrackParam = track->GetTrackParamAtCluster()->GetEntries();
 
-    if (lChOnOff[iCh-1])
-    {
+  Bool_t first( kTRUE );
+  for( Int_t iTrackParam = 0; iTrackParam < nTrackParam; ++iTrackParam )
+  {
 
-      Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
-      if (iCh>=1 && iCh<=4)
-      {
+    // get new pointers
+    AliMUONTrackParam* trackParam( (AliMUONTrackParam *) track->GetTrackParamAtCluster()->At(iTrackParam) );
+    if( !trackParam ) continue;
 
-        if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0])
-        {
-          // From track crossings
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+    AliMUONVCluster* cluster = trackParam->GetClusterPtr();
+    if( !cluster ) continue;
 
-        if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1])
-        {
-          // From track crossings
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+    // fill local variables for this position --> one measurement
+    FillDetElemData( cluster );
+    FillRecPointData( cluster );
 
-      }
+    // unbias and store track parameters
+    if( fUnbias && !UnbiasTrackParamData( trackParam ) ) continue;
+    FillTrackParamData( trackParam );
 
-      if (iCh>=5 && iCh<=6)
-      {
+    if( first )
+    {
 
-        if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0])
-        {
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+      // for first valid cluster, save track position as "starting" values
+      first = kFALSE;
 
-        if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
-          (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1])
-        {
+      fTrackPos0[0] = fTrackPos[0];
+      fTrackPos0[1] = fTrackPos[1];
+      fTrackPos0[2] = fTrackPos[2];
+      fTrackSlope0[0] = fTrackSlope[0];
+      fTrackSlope0[1] = fTrackSlope[1];
 
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+    }
 
-      }
+    // 'inverse' (GlobalToLocal) rotation matrix
+    const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
 
-      if (iCh>=7 && iCh<=10)
-      {
+    // calculate measurements
+    if( fBFieldOn )
+    {
 
-        if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0])
-        {
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+      // use residuals (cluster - track) for measurement
+      fMeas[0] = r[0]*(fClustPos[0] - fTrackPos[0]) + r[1]*(fClustPos[1] - fTrackPos[1]);
+      fMeas[1] = r[3]*(fClustPos[0] - fTrackPos[0]) + r[4]*(fClustPos[1] - fTrackPos[1]);
 
-        if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
-          (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1])
-        {
-          FixParameter(i*fgNParCh+0, 0.0);
-          FixParameter(i*fgNParCh+1, 0.0);
-          FixParameter(i*fgNParCh+2, 0.0);
-          FixParameter(i*fgNParCh+3, 0.0);
-        }
+    } else {
 
-      }
+      // use cluster position for measurement
+      fMeas[0] = ( r[0]*fClustPos[0] + r[1]*fClustPos[1] );
+      fMeas[1] = ( r[3]*fClustPos[0] + r[4]*fClustPos[1] );
 
     }
 
-  }
+    // Set local equations
+    LocalEquationX();
+    LocalEquationY();
 
-}
+  }
 
-//______________________________________________________________________
-void AliMUONAlignment::SetNonLinear(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
-{
+  // copy track record
+  fMillepede->SetRecordRun(fRunNumber);
+  fMillepede->SetRecordWeight(weight);
+  fTrackRecord = *fMillepede->GetRecord();
 
-  /// Set non linear parameter flag selected chambers and degrees of freedom
-  for (Int_t i = 0; i < fgNDetElem; i++)
+  // save record data
+  if( doAlignment )
   {
+    fMillepede->SaveRecordData();
+  }
 
-    Int_t iCh=0;
-    for (iCh=1; iCh<=fgNCh; iCh++)
-    { if (i<fgSNDetElemCh[iCh-1]) break; }
-
-    if (lChOnOff[iCh-1])
-    {
+  // return record
+  return &fTrackRecord;
 
-      if (lVarXYT[0])
-      {
-        // X constraint
-        SetNonLinear(i*fgNParCh+0);
-      }
+}
 
-      if (lVarXYT[1])
-      {
-        // Y constraint
-        SetNonLinear(i*fgNParCh+1);
-      }
+//______________________________________________________________________________
+void AliMUONAlignment::ProcessTrack( AliMillePedeRecord* trackRecord )
+{
+  /// process track record
+  if( !trackRecord ) return;
 
-      if (lVarXYT[2])
-      {
-        // T constraint
-        SetNonLinear(i*fgNParCh+2);
-      }
+  // make sure record storage is initialized
+  if( !fMillepede->GetRecord() ) fMillepede->InitDataRecStorage();
 
-      if (lVarXYT[3])
-      {
-        // Z constraint
-        SetNonLinear(i*fgNParCh+3);
-      }
+  // copy content
+  *fMillepede->GetRecord() = *trackRecord;
 
-    }
+  // save record
+  fMillepede->SaveRecordData();
 
-  }
+  return;
 
 }
 
-//______________________________________________________________________
-void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
+//_____________________________________________________________________
+void AliMUONAlignment::FixAll( UInt_t mask )
 {
+  /// fix parameters matching mask, for all chambers
+  AliInfo( Form( "Fixing %s for all detector elements", GetParameterMaskString( mask ).Data() ) );
 
-  /// Add constraint equations for selected chambers and degrees of freedom
-  for (Int_t i = 0; i < fgNDetElem; i++)
+  // fix all stations
+  for( Int_t i = 0; i < fgNDetElem; ++i )
   {
+    if( mask & ParX )  FixParameter(i, 0);
+    if( mask & ParY )  FixParameter(i, 1);
+    if( mask & ParTZ ) FixParameter(i, 2);
+    if( mask & ParZ )  FixParameter(i, 3);
+  }
 
-    Int_t iCh=0;
-    for (iCh=1; iCh<=fgNCh; iCh++)
-    {
-      if (i<fgSNDetElemCh[iCh-1]) break;
-    }
+}
 
-    if (lChOnOff[iCh-1])
-    {
-      if (lVarXYT[0])
-      {
-        // X constraint
-        fConstraintX[i*fgNParCh+0]=1.0;
-      }
+//_____________________________________________________________________
+void AliMUONAlignment::FixChamber( Int_t iCh, UInt_t mask )
+{
+  /// fix parameters matching mask, for all detector elements in a given chamber, counting from 1
 
-      if (lVarXYT[1])
-      {
-        // Y constraint
-        fConstraintY[i*fgNParCh+1]=1.0;
-      }
+  // check boundaries
+  if( iCh < 1 || iCh > 10 )
+  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
 
-      if (lVarXYT[2])
-      {
-        // T constraint
-        fConstraintP[i*fgNParCh+2]=1.0;
-      }
-      //       if (lVarXYT[3]) { // Z constraint
-      //       fConstraintP[i*fgNParCh+3]=1.0;
-      //       }
+  // get first and last element
+  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
+  {
 
-    }
-  }
+    AliInfo( Form( "Fixing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
 
-  if (lVarXYT[0])
-  {
-    // X constraint
-    AddConstraint(fConstraintX,0.0);
-  }
+    if( mask & ParX )  FixParameter(i, 0);
+    if( mask & ParY )  FixParameter(i, 1);
+    if( mask & ParTZ ) FixParameter(i, 2);
+    if( mask & ParZ )  FixParameter(i, 3);
 
-  if (lVarXYT[1])
-  {
-    // Y constraint
-    AddConstraint(fConstraintY,0.0);
   }
+}
 
-  if (lVarXYT[2])
-  {
-    // T constraint
-    AddConstraint(fConstraintP,0.0);
-  }
+//_____________________________________________________________________
+void AliMUONAlignment::FixDetElem( Int_t iDetElemId, UInt_t mask )
+{
+  /// fix parameters matching mask, for a given detector element, counting from 0
+  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+  if ( mask & ParX )  FixParameter(iDet, 0);
+  if ( mask & ParY )  FixParameter(iDet, 1);
+  if ( mask & ParTZ ) FixParameter(iDet, 2);
+  if ( mask & ParZ )  FixParameter(iDet, 3);
 
-//   if (lVarXYT[3]) { // Z constraint
-//     AddConstraint(fConstraintP,0.0);
-//   }
 }
 
-//______________________________________________________________________
-void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, const Bool_t *lDetTLBR, const Bool_t *lSpecLROnOff)
+//_____________________________________________________________________
+void AliMUONAlignment::FixHalfSpectrometer( const Bool_t *lChOnOff, UInt_t sidesMask, UInt_t mask )
 {
-  /// Add constraint equations for selected chambers, degrees of freedom and detector half
-  Double_t lDetElemLocX = 0.;
-  Double_t lDetElemLocY = 0.;
-  Double_t lDetElemLocZ = 0.;
-  Double_t lDetElemGloX = 0.;
-  Double_t lDetElemGloY = 0.;
-  Double_t lDetElemGloZ = 0.;
-  Double_t lMeanY = 0.;
-  Double_t lSigmaY = 0.;
-  Double_t lMeanZ = 0.;
-  Double_t lSigmaZ = 0.;
-  Int_t lNDetElem = 0;
-  for (Int_t i = 0; i < fgNDetElem; i++)
+
+  /// Fix parameters matching mask for all detectors in selected chambers and selected sides of the spectrometer
+  for( Int_t i = 0; i < fgNDetElem; ++i )
   {
 
-    Int_t iCh=0;
-    for (iCh=1; iCh<=fgNCh; iCh++){
-      if (i<fgSNDetElemCh[iCh-1]) break;
-    }
-    if (lChOnOff[iCh-1]){
-      Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
-      Int_t lDetElemId = iCh*100+lDetElemNumber;
-      fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
-             lDetElemGloX,lDetElemGloY,lDetElemGloZ);
-      if (iCh>=1 && iCh<=4){
-  if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-  if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-      }
-      if (iCh>=5 && iCh<=6){
-  if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-  if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
-       (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-      }
-      if (iCh>=7 && iCh<=10){
-  if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-  if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
-       (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){
-    lMeanY += lDetElemGloY;
-    lSigmaY += lDetElemGloY*lDetElemGloY;
-    lMeanZ += lDetElemGloZ;
-    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
-    lNDetElem++;
-  }
-      }
-    }
-  }
-  if (lNDetElem) {
-    lMeanY /= lNDetElem;
-    lSigmaY /= lNDetElem;
-    lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
-    lMeanZ /= lNDetElem;
-    lSigmaZ /= lNDetElem;
-    lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
-     AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ));
-  } else {
-    AliError("No detection elements to constrain!!!");
-    return;    
-  }
-               
-  for (Int_t i = 0; i < fgNDetElem; i++){
-    Int_t iCh=0;
-    for (iCh=1; iCh<=fgNCh; iCh++){
-      if (i<fgSNDetElemCh[iCh-1]) break;
-    }
-    if (lChOnOff[iCh-1]){
-      Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
-      Int_t lDetElemId = iCh*100+lDetElemNumber;
-      fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
-             lDetElemGloX,lDetElemGloY,lDetElemGloZ);
-      if (lVarXYT[0]) { // X constraint
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
-      }
-      if (lVarXYT[1]) { // Y constraint
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
-      }
-      if (lVarXYT[2]) { // P constraint
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
-      }
-      if (lVarXYT[3]) { // X-Z shearing
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
-      }
-      if (lVarXYT[4]) { // Y-Z shearing
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
-      }
-      if (lVarXYT[5]) { // P-Z rotation
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
-      }
-      if (lVarXYT[6]) { // X-Y shearing
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
-      }
-      if (lVarXYT[7]) { // Y-Y scaling
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
-      }
-      if (lVarXYT[8]) { // P-Y rotation
-  if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
-  if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
-  if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
-  if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
-      }
-    }
-  }
-  if (lVarXYT[0]) { // X constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
-  }
-  if (lVarXYT[1]) { // Y constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
-  }
-  if (lVarXYT[2]) { // T constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
-  }
-  if (lVarXYT[3]) { // X-Z constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintXZT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintXZL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintXZB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintXZR,0.0); // Right half
-  }
-  if (lVarXYT[4]) { // Y-Z constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintYZT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintYZL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintYZB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintYZR,0.0); // Right half
-  }
-  if (lVarXYT[5]) { // P-Z constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintPZT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintPZL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintPZB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintPZR,0.0); // Right half
-  }
-  if (lVarXYT[6]) { // X-Y constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintXYT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintXYL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintXYB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintXYR,0.0); // Right half
-  }
-  if (lVarXYT[7]) { // Y-Y constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintYYT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintYYL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintYYB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintYYR,0.0); // Right half
-  }
-  if (lVarXYT[8]) { // P-Y constraint
-    if (lDetTLBR[0]) AddConstraint(fConstraintPYT,0.0); // Top half
-    if (lDetTLBR[1]) AddConstraint(fConstraintPYL,0.0); // Left half
-    if (lDetTLBR[2]) AddConstraint(fConstraintPYB,0.0); // Bottom half
-    if (lDetTLBR[3]) AddConstraint(fConstraintPYR,0.0); // Right half
-  }
-}
+    // get chamber matching detector
+    const Int_t iCh( GetChamberId(i) );
+    if( !lChOnOff[iCh-1] ) continue;
 
-void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t /*lWeight*/) const{
-  /// Set constrain equation for top half of spectrometer
-  Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
-  if (lCh>=1 && lCh<=4){
-    if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
-      lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
-    }
-  }
-  if (lCh>=5 && lCh<=6){
-    if (lDetElemNumber>=0&&lDetElemNumber<=9){
-      lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
-    }
-  }
-  if (lCh>=7 && lCh<=10){
-    if (lDetElemNumber>=0&&lDetElemNumber<=13){
-      lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
-    }
-  }
-}
+    // get detector element in chamber
+    Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
 
-//______________________________________________________________________
-void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight) const{
-  /// Set constrain equation for left half of spectrometer
-  Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
-  if (lCh>=1 && lCh<=4){
-    if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
-      lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
-    }
-  }
-  if (lCh>=5 && lCh<=6){
-    if (lDetElemNumber>=5&&lDetElemNumber<=13){
-      lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
-    }
-  }
-  if (lCh>=7 && lCh<=10){
-    if (lDetElemNumber>=7&&lDetElemNumber<=19){
-      lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
+    // skip detector if its side is off
+    // stations 1 and 2
+    if( iCh>=1 && iCh<=4 )
+    {
+      if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
     }
-  }
-}
 
-//______________________________________________________________________
-void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t /*lWeight*/) const{
-  /// Set constrain equation for bottom half of spectrometer
-  Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
-  if (lCh>=1 && lCh<=4){
-    if (lDetElemNumber==2 || lDetElemNumber==3){ // From track crossings
-      lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
-    }
-  }
-  if (lCh>=5 && lCh<=6){
-    if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
-  (lDetElemNumber==0)){
-      lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
-    }
-  }
-  if (lCh>=7 && lCh<=10){
-    if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
-  (lDetElemNumber==0)){
-      lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
+    // station 3
+    if (iCh>=5 && iCh<=6)
+    {
+      if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
     }
-  }
-}
 
-//______________________________________________________________________
-void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight) const{
-  /// Set constrain equation for right half of spectrometer
-  Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
-  if (lCh>=1 && lCh<=4){
-    if (lDetElemNumber==0 || lDetElemNumber==3){ // From track crossings
-      lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
-    }
-  }
-  if (lCh>=5 && lCh<=6){
-    if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
-  (lDetElemNumber>=14&&lDetElemNumber<=17)){
-      lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
-    }
-  }
-  if (lCh>=7 && lCh<=10){
-    if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
-  (lDetElemNumber>=20&&lDetElemNumber<=25)){
-      lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
+    // stations 4 and 5
+    if (iCh>=7 && iCh<=10)
+    {
+      if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
     }
-  }
-}
 
-//______________________________________________________________________
-void AliMUONAlignment::ResetConstraints(){
-  /// Reset all constraint equations
-  for (Int_t i = 0; i < fgNDetElem; i++){
-       fConstraintX3[i*fgNParCh+0]=0.0; 
-       fConstraintY3[i*fgNParCh+0]=0.0; 
-       fConstraintX4[i*fgNParCh+0]=0.0; 
-       fConstraintY4[i*fgNParCh+0]=0.0; 
-       fConstraintP4[i*fgNParCh+0]=0.0; 
-       fConstraintX5[i*fgNParCh+0]=0.0; 
-       fConstraintY5[i*fgNParCh+0]=0.0; 
-       fConstraintX[i*fgNParCh+0]=0.0;
-    fConstraintX[i*fgNParCh+1]=0.0;
-    fConstraintX[i*fgNParCh+2]=0.0;
-    fConstraintY[i*fgNParCh+0]=0.0;
-    fConstraintY[i*fgNParCh+1]=0.0;
-    fConstraintY[i*fgNParCh+2]=0.0;
-    fConstraintP[i*fgNParCh+0]=0.0;
-    fConstraintP[i*fgNParCh+1]=0.0;
-    fConstraintP[i*fgNParCh+2]=0.0;
-    fConstraintXT[i*fgNParCh+0]=0.0;
-    fConstraintXT[i*fgNParCh+1]=0.0;
-    fConstraintXT[i*fgNParCh+2]=0.0;
-    fConstraintYT[i*fgNParCh+0]=0.0;
-    fConstraintYT[i*fgNParCh+1]=0.0;
-    fConstraintYT[i*fgNParCh+2]=0.0;
-    fConstraintPT[i*fgNParCh+0]=0.0;
-    fConstraintPT[i*fgNParCh+1]=0.0;
-    fConstraintPT[i*fgNParCh+2]=0.0;
-    fConstraintXZT[i*fgNParCh+0]=0.0;
-    fConstraintXZT[i*fgNParCh+1]=0.0;
-    fConstraintXZT[i*fgNParCh+2]=0.0;
-    fConstraintYZT[i*fgNParCh+0]=0.0;
-    fConstraintYZT[i*fgNParCh+1]=0.0;
-    fConstraintYZT[i*fgNParCh+2]=0.0;
-    fConstraintPZT[i*fgNParCh+0]=0.0;
-    fConstraintPZT[i*fgNParCh+1]=0.0;
-    fConstraintPZT[i*fgNParCh+2]=0.0;
-    fConstraintXYT[i*fgNParCh+0]=0.0;
-    fConstraintXYT[i*fgNParCh+1]=0.0;
-    fConstraintXYT[i*fgNParCh+2]=0.0;
-    fConstraintYYT[i*fgNParCh+0]=0.0;
-    fConstraintYYT[i*fgNParCh+1]=0.0;
-    fConstraintYYT[i*fgNParCh+2]=0.0;
-    fConstraintPYT[i*fgNParCh+0]=0.0;
-    fConstraintPYT[i*fgNParCh+1]=0.0;
-    fConstraintPYT[i*fgNParCh+2]=0.0;
-    fConstraintXL[i*fgNParCh+0]=0.0;
-    fConstraintXL[i*fgNParCh+1]=0.0;
-    fConstraintXL[i*fgNParCh+2]=0.0;
-    fConstraintYL[i*fgNParCh+0]=0.0;
-    fConstraintYL[i*fgNParCh+1]=0.0;
-    fConstraintYL[i*fgNParCh+2]=0.0;
-    fConstraintPL[i*fgNParCh+0]=0.0;
-    fConstraintPL[i*fgNParCh+1]=0.0;
-    fConstraintPL[i*fgNParCh+2]=0.0;
-    fConstraintXZL[i*fgNParCh+0]=0.0;
-    fConstraintXZL[i*fgNParCh+1]=0.0;
-    fConstraintXZL[i*fgNParCh+2]=0.0;
-    fConstraintYZL[i*fgNParCh+0]=0.0;
-    fConstraintYZL[i*fgNParCh+1]=0.0;
-    fConstraintYZL[i*fgNParCh+2]=0.0;
-    fConstraintPZL[i*fgNParCh+0]=0.0;
-    fConstraintPZL[i*fgNParCh+1]=0.0;
-    fConstraintPZL[i*fgNParCh+2]=0.0;
-    fConstraintXYL[i*fgNParCh+0]=0.0;
-    fConstraintXYL[i*fgNParCh+1]=0.0;
-    fConstraintXYL[i*fgNParCh+2]=0.0;
-    fConstraintYYL[i*fgNParCh+0]=0.0;
-    fConstraintYYL[i*fgNParCh+1]=0.0;
-    fConstraintYYL[i*fgNParCh+2]=0.0;
-    fConstraintPYL[i*fgNParCh+0]=0.0;
-    fConstraintPYL[i*fgNParCh+1]=0.0;
-    fConstraintPYL[i*fgNParCh+2]=0.0;
-    fConstraintXB[i*fgNParCh+0]=0.0;
-    fConstraintXB[i*fgNParCh+1]=0.0;
-    fConstraintXB[i*fgNParCh+2]=0.0;
-    fConstraintYB[i*fgNParCh+0]=0.0;
-    fConstraintYB[i*fgNParCh+1]=0.0;
-    fConstraintYB[i*fgNParCh+2]=0.0;
-    fConstraintPB[i*fgNParCh+0]=0.0;
-    fConstraintPB[i*fgNParCh+1]=0.0;
-    fConstraintPB[i*fgNParCh+2]=0.0;
-    fConstraintXZB[i*fgNParCh+0]=0.0;
-    fConstraintXZB[i*fgNParCh+1]=0.0;
-    fConstraintXZB[i*fgNParCh+2]=0.0;
-    fConstraintYZB[i*fgNParCh+0]=0.0;
-    fConstraintYZB[i*fgNParCh+1]=0.0;
-    fConstraintYZB[i*fgNParCh+2]=0.0;
-    fConstraintPZB[i*fgNParCh+0]=0.0;
-    fConstraintPZB[i*fgNParCh+1]=0.0;
-    fConstraintPZB[i*fgNParCh+2]=0.0;
-    fConstraintXYB[i*fgNParCh+0]=0.0;
-    fConstraintXYB[i*fgNParCh+1]=0.0;
-    fConstraintXYB[i*fgNParCh+2]=0.0;
-    fConstraintYYB[i*fgNParCh+0]=0.0;
-    fConstraintYYB[i*fgNParCh+1]=0.0;
-    fConstraintYYB[i*fgNParCh+2]=0.0;
-    fConstraintPYB[i*fgNParCh+0]=0.0;
-    fConstraintPYB[i*fgNParCh+1]=0.0;
-    fConstraintPYB[i*fgNParCh+2]=0.0;
-    fConstraintXR[i*fgNParCh+0]=0.0;
-    fConstraintXR[i*fgNParCh+1]=0.0;
-    fConstraintXR[i*fgNParCh+2]=0.0;
-    fConstraintYR[i*fgNParCh+0]=0.0;
-    fConstraintYR[i*fgNParCh+1]=0.0;
-    fConstraintYR[i*fgNParCh+2]=0.0;
-    fConstraintPR[i*fgNParCh+0]=0.0;
-    fConstraintPR[i*fgNParCh+1]=0.0;
-    fConstraintPR[i*fgNParCh+2]=0.0;
-    fConstraintXZR[i*fgNParCh+0]=0.0;
-    fConstraintXZR[i*fgNParCh+1]=0.0;
-    fConstraintXZR[i*fgNParCh+2]=0.0;
-    fConstraintYZR[i*fgNParCh+0]=0.0;
-    fConstraintYZR[i*fgNParCh+1]=0.0;
-    fConstraintYZR[i*fgNParCh+2]=0.0;
-    fConstraintPZR[i*fgNParCh+0]=0.0;
-    fConstraintPZR[i*fgNParCh+1]=0.0;
-    fConstraintPZR[i*fgNParCh+2]=0.0;
-    fConstraintPZR[i*fgNParCh+0]=0.0;
-    fConstraintPZR[i*fgNParCh+1]=0.0;
-    fConstraintPZR[i*fgNParCh+2]=0.0;
-    fConstraintXYR[i*fgNParCh+0]=0.0;
-    fConstraintXYR[i*fgNParCh+1]=0.0;
-    fConstraintXYR[i*fgNParCh+2]=0.0;
-    fConstraintYYR[i*fgNParCh+0]=0.0;
-    fConstraintYYR[i*fgNParCh+1]=0.0;
-    fConstraintYYR[i*fgNParCh+2]=0.0;
-    fConstraintPYR[i*fgNParCh+0]=0.0;
-    fConstraintPYR[i*fgNParCh+1]=0.0;
-    fConstraintPYR[i*fgNParCh+2]=0.0;
+    // detector is accepted, fix it
+    FixDetElem( i, mask );
+
   }
-}
 
-//______________________________________________________________________
-void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
-  /// Constrain equation defined by par to value
-  fMillepede->SetGlobalConstraint(par, value);
-  AliInfo("Adding constraint");
 }
 
 //______________________________________________________________________
-void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
-  /// Initialize global parameters with par array
-  fMillepede->SetGlobalParameters(par);
-  AliInfo("Init Global Parameters");
-}
+void AliMUONAlignment::FixParameter( Int_t iPar )
+{
 
-//______________________________________________________________________
-void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
-  /// Parameter iPar is encourage to vary in [-value;value].
-  /// If value == 0, parameter is fixed
-  fMillepede->SetParSigma(iPar, value);
-  if (TMath::Abs(value)<1e-4) AliInfo(Form("Parameter %i Fixed", iPar));
-}
+  /// fix a given parameter, counting from 0
+  if( fInitialized )
+  { AliFatal( "Millepede already initialized" ); }
+
+  fGlobalParameterStatus[iPar] = kFixedParId;
 
-//______________________________________________________________________
-void AliMUONAlignment::ResetLocalEquation()
-{
-  /// Reset the derivative vectors
-  for(int i=0; i<fNLocal; i++) {
-    fLocalDerivatives[i] = 0.0;
-  }
-  for(int i=0; i<fNGlobal; i++) {
-    fGlobalDerivatives[i] = 0.0;
-  }
 }
 
-//______________________________________________________________________
-void AliMUONAlignment::AllowVariations(const Bool_t *bChOnOff)
+
+//_____________________________________________________________________
+void AliMUONAlignment::ReleaseChamber( Int_t iCh, UInt_t mask )
 {
+  /// release parameters matching mask, for all detector elements in a given chamber, counting from 1
 
-  /// Set allowed variation for selected chambers based on fDoF and fAllowVar
-  for (Int_t iCh=1; iCh<=10; iCh++)
-  {
-    if (bChOnOff[iCh-1])
-    {
+  // check boundaries
+  if( iCh < 1 || iCh > 10 )
+  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
 
-      Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0;
-      Int_t iDetElemLast = fgSNDetElemCh[iCh-1];
-      for (int i=0; i<fgNParCh; i++)
-      {
-        AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
-        if (fDoF[i])
-        {
-          for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
-            FixParameter(j*fgNParCh+i, fAllowVar[i]);
-          }
-        }
+  // get first and last element
+  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
+  {
 
-      }
+    AliInfo( Form( "Releasing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
 
-    }
+    if( mask & ParX )  ReleaseParameter(i, 0);
+    if( mask & ParY )  ReleaseParameter(i, 1);
+    if( mask & ParTZ ) ReleaseParameter(i, 2);
+    if( mask & ParZ )  ReleaseParameter(i, 3);
 
   }
-
 }
 
-//______________________________________________________________________
-void AliMUONAlignment::SetNonLinear(Int_t iPar  /* set non linear flag */ )
+//_____________________________________________________________________
+void AliMUONAlignment::ReleaseDetElem( Int_t iDetElemId, UInt_t mask )
 {
-  /// Set nonlinear flag for parameter iPar
-  fMillepede->SetNonLinear(iPar);
-  AliInfo(Form("Parameter %i set to non linear", iPar));
+  /// release parameters matching mask, for a given detector element, counting from 0
+  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+  if ( mask & ParX )  ReleaseParameter(iDet, 0);
+  if ( mask & ParY )  ReleaseParameter(iDet, 1);
+  if ( mask & ParTZ ) ReleaseParameter(iDet, 2);
+  if ( mask & ParZ )  ReleaseParameter(iDet, 3);
+
 }
 
 //______________________________________________________________________
-void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
+void AliMUONAlignment::ReleaseParameter( Int_t iPar )
 {
 
-  /// Set expected measurement resolution
-  fSigma[0] = sigmaX;   fSigma[1] = sigmaY;
-  AliInfo(Form("Using fSigma[0]=%f and fSigma[1]=%f",fSigma[0],fSigma[1]));
-}
+  /// release a given parameter, counting from 0
+  if( fInitialized )
+  { AliFatal( "Millepede already initialized" ); }
 
+  fGlobalParameterStatus[iPar] = kFreeParId;
 
-//______________________________________________________________________
-void AliMUONAlignment::LocalEquationX( Bool_t doAlignment )
+}
+
+//_____________________________________________________________________
+void AliMUONAlignment::GroupChamber( Int_t iCh, UInt_t mask )
 {
+  /// group parameters matching mask for all detector elements in a given chamber, counting from 1
+  if( iCh < 1 || iCh > 10 )
+  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
 
+  const Int_t detElemMin = 100*iCh;
+  const Int_t detElemMax = 100*iCh + fgNDetElemCh[iCh]-1;
+  GroupDetElems( detElemMin, detElemMax, mask );
 
-  // local cluster record
-  AliMUONAlignmentClusterRecord clusterRecord;
+}
 
-  // store detector and measurement
-  clusterRecord.SetDetElemId( fDetElemId );
-  clusterRecord.SetDetElemNumber( fDetElemNumber );
-  clusterRecord.SetMeas( fMeas[0] );
-  clusterRecord.SetSigma( fSigma[0] );
+//_____________________________________________________________________
+void AliMUONAlignment::GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t mask )
+{
+  /// group parameters matching mask for all detector elements between min and max
+  // check number of detector elements
+  const Int_t nDetElem = detElemMax - detElemMin + 1;
+  if( nDetElem<2 )
+  { AliFatal( Form( "Requested group of DEs %d-%d contains less than 2 DE's", detElemMin, detElemMax ) ); }
 
-  // store local derivatives
-  clusterRecord.SetLocalDerivative( 0, fCosPhi );
-  clusterRecord.SetLocalDerivative( 1, fCosPhi*(fTrackPos[2] - fTrackPos0[2]) );
-  clusterRecord.SetLocalDerivative( 2, fSinPhi );
-  clusterRecord.SetLocalDerivative( 3, fSinPhi*(fTrackPos[2] - fTrackPos0[2]) );
+  // create list
+  Int_t* detElemList = new int[nDetElem];
+  for( Int_t i = 0; i < nDetElem; ++i )
+  { detElemList[i] = detElemMin+i; }
+
+  // group
+  GroupDetElems( detElemList, nDetElem, mask );
+  delete[] detElemList;
 
-  // store global derivatives
-  clusterRecord.SetGlobalDerivative( 0, -fCosPhi );
-  clusterRecord.SetGlobalDerivative( 1, -fSinPhi );
+}
 
-  if (fBFieldOn)
+//_____________________________________________________________________
+void AliMUONAlignment::GroupDetElems( Int_t* detElemList, Int_t nDetElem, UInt_t mask )
+{
+  /// group parameters matching mask for all detector elements in list
+  if( fInitialized )
+  { AliFatal( "Millepede already initialized" ); }
+
+  const Int_t iDeBase( GetDetElemNumber( detElemList[0] ) );
+  for( Int_t i = 0; i < nDetElem; ++i )
   {
+    const Int_t iDeCurrent( GetDetElemNumber( detElemList[i] ) );
+    if( mask & ParX ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 0] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
+    if( mask & ParY ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 1] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
+    if( mask & ParTZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 2] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
+    if( mask & ParZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 3] = (i==0) ?  kGroupBaseId : (kGroupBaseId-iDeBase-1);
+
+    if( i== 0 ) AliInfo( Form( "Creating new group for detector %i and variable %s", detElemList[i], GetParameterMaskString( mask ).Data() ) );
+    else AliInfo( Form( "Adding detector element %i to current group", detElemList[i] ) );
+  }
 
-    clusterRecord.SetGlobalDerivative(
-      2,
-      -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
-      +fCosPhi*(fTrackPos[1]-fDetElemPos[1]) );
+}
 
-  } else {
+//______________________________________________________________________
+void AliMUONAlignment::SetChamberNonLinear( Int_t iCh, UInt_t mask )
+{
+  /// Set parameters matching mask as non linear, for all detector elements in a given chamber, counting from 1
+  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
+  {
 
-    clusterRecord.SetGlobalDerivative(
-      2,
-      -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
-      +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]) );
+      if( mask & ParX ) SetParameterNonLinear(i, 0);
+      if( mask & ParY ) SetParameterNonLinear(i, 1);
+      if( mask & ParTZ ) SetParameterNonLinear(i, 2);
+      if( mask & ParZ ) SetParameterNonLinear(i, 3);
 
   }
 
-  clusterRecord.SetGlobalDerivative( 3, fCosPhi*fTrackSlope0[0]+fSinPhi*fTrackSlope0[1] );
-
-  // append to trackRecord
-  fTrackRecord.AddClusterRecord( clusterRecord );
+}
 
-  // store local equation
-  if( doAlignment ) LocalEquation( clusterRecord );
+//_____________________________________________________________________
+void AliMUONAlignment::SetDetElemNonLinear( Int_t iDetElemId, UInt_t mask )
+{
+  /// Set parameters matching mask as non linear, for a given detector element, counting from 0
+  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+  if ( mask & ParX )  SetParameterNonLinear(iDet, 0);
+  if ( mask & ParY )  SetParameterNonLinear(iDet, 1);
+  if ( mask & ParTZ ) SetParameterNonLinear(iDet, 2);
+  if ( mask & ParZ )  SetParameterNonLinear(iDet, 3);
 
 }
 
 //______________________________________________________________________
-void AliMUONAlignment::LocalEquationY(Bool_t doAlignment )
+void AliMUONAlignment::SetParameterNonLinear( Int_t iPar )
 {
+  /// Set nonlinear flag for parameter iPar
+  if( !fInitialized )
+  { AliFatal( "Millepede not initialized" ); }
 
-  // local cluster record
-  AliMUONAlignmentClusterRecord clusterRecord;
-
-  // store detector and measurement
-  clusterRecord.SetDetElemId( fDetElemId );
-  clusterRecord.SetDetElemNumber( fDetElemNumber );
-  clusterRecord.SetMeas( fMeas[1] );
-  clusterRecord.SetSigma( fSigma[1] );
+  fMillepede->SetNonLinear( iPar );
+  AliInfo( Form( "Parameter %i set to non linear", iPar ) );
+}
 
-  // store local derivatives
-  clusterRecord.SetLocalDerivative( 0, -fSinPhi );
-  clusterRecord.SetLocalDerivative( 1, -fSinPhi*(fTrackPos[2] - fTrackPos0[2] ) );
-  clusterRecord.SetLocalDerivative( 2, fCosPhi );
-  clusterRecord.SetLocalDerivative( 3, fCosPhi*(fTrackPos[2] - fTrackPos0[2] ) );
+//______________________________________________________________________
+void AliMUONAlignment::AddConstraints( const Bool_t *lChOnOff, UInt_t mask )
+{
+  /// Add constraint equations for selected chambers and degrees of freedom
 
-  // set global derivatives
-  clusterRecord.SetGlobalDerivative( 0,  fSinPhi);
-  clusterRecord.SetGlobalDerivative( 1, -fCosPhi);
+  Array fConstraintX;
+  Array fConstraintY;
+  Array fConstraintTZ;
+  Array fConstraintZ;
 
-  if (fBFieldOn)
+  for( Int_t i = 0; i < fgNDetElem; ++i )
   {
 
-    clusterRecord.SetGlobalDerivative(
-       2,
-      -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
-      -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
+    // get chamber matching detector
+    const Int_t iCh( GetChamberId(i) );
+    if (lChOnOff[iCh-1])
+    {
 
-  } else {
+      if( mask & ParX ) fConstraintX.values[i*fgNParCh+0]=1.0;
+      if( mask & ParY ) fConstraintY.values[i*fgNParCh+1]=1.0;
+      if( mask & ParTZ ) fConstraintTZ.values[i*fgNParCh+2]=1.0;
+      if( mask & ParZ ) fConstraintTZ.values[i*fgNParCh+3]=1.0;
 
-    clusterRecord.SetGlobalDerivative(
-       2,
-      -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
-      -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
+    }
   }
 
-  clusterRecord.SetGlobalDerivative( 3, -fSinPhi*fTrackSlope0[0]+fCosPhi*fTrackSlope0[1]);
-
-  // append to trackRecord
-  fTrackRecord.AddClusterRecord( clusterRecord );
-
-  // store local equation
-  if( doAlignment ) LocalEquation( clusterRecord );
+  if( mask & ParX ) AddConstraint(fConstraintX.values,0.0);
+  if( mask & ParY ) AddConstraint(fConstraintY.values,0.0);
+  if( mask & ParTZ ) AddConstraint(fConstraintTZ.values,0.0);
+  if( mask & ParZ ) AddConstraint(fConstraintZ.values,0.0);
 
 }
 
-//_____________________________________________________
-void AliMUONAlignment::LocalEquation( const AliMUONAlignmentClusterRecord& clusterRecord )
+//______________________________________________________________________
+void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask )
 {
+  /*
+  questions:
+  - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
+  - why is weight ignored for ConstrainT and ConstrainB
+  - why is there no constrain on z
+  */
 
-  // copy to local variables
-  for( Int_t index = 0; index < 4; ++index )
-  {
-    SetLocalDerivative( index, clusterRecord.GetLocalDerivative( index ) );
-    SetGlobalDerivative( clusterRecord.GetDetElemNumber()*fgNParCh + index, clusterRecord.GetGlobalDerivative( index ) );
-  }
+  /// Add constraint equations for selected chambers, degrees of freedom and detector half
+  Double_t lMeanY = 0.;
+  Double_t lSigmaY = 0.;
+  Double_t lMeanZ = 0.;
+  Double_t lSigmaZ = 0.;
+  Int_t lNDetElem = 0;
 
-  // pass equation parameters to millepede
-  fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, clusterRecord.GetMeas(), clusterRecord.GetSigma() );
+  for( Int_t i = 0; i < fgNDetElem; ++i )
+  {
 
-}
+    // get chamber matching detector
+    const Int_t iCh( GetChamberId(i) );
 
-//______________________________________________________________________
-void AliMUONAlignment::FillRecPointData()
-{
+    // skip detector if chamber is off
+    if( lChOnOff[iCh-1] ) continue;
 
-  /// Get information of current cluster
-  fClustPos[0] = fCluster->GetX();
-  fClustPos[1] = fCluster->GetY();
-  fClustPos[2] = fCluster->GetZ();
-  fTransform->Global2Local(
-    fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
-    fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
+    // get detector element id from detector element number
+    const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
+    const Int_t lDetElemId = iCh*100+lDetElemNumber;
 
-}
+    // skip detector if its side is off
+    // stations 1 and 2
+    if( iCh>=1 && iCh<=4 )
+    {
+      if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
+    }
 
-//______________________________________________________________________
-void AliMUONAlignment::FillTrackParamData()
-{
+    // station 3
+    if (iCh>=5 && iCh<=6)
+    {
+      if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
+    }
 
-  /// Get information of current track at current cluster
-  fTrackPos[0] = fTrackParam->GetNonBendingCoor();
-  fTrackPos[1] = fTrackParam->GetBendingCoor();
-  fTrackPos[2] = fTrackParam->GetZ();
-  fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
-  fTrackSlope[1] = fTrackParam->GetBendingSlope();
-  fTransform->Global2Local(
-    fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
-    fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
+    // stations 4 and 5
+    if (iCh>=7 && iCh<=10)
+    {
+      if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
+      if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
+      if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
+      if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
+    }
 
-}
+    // get global x, y and z position
+    Double_t lDetElemGloX = 0.;
+    Double_t lDetElemGloY = 0.;
+    Double_t lDetElemGloZ = 0.;
+    fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
 
-//_____________________________________________________
-void AliMUONAlignment::FillDetElemData()
-{
+    // increment mean Y, mean Z, sigmas and number of accepted detectors
+    lMeanY += lDetElemGloY;
+    lSigmaY += lDetElemGloY*lDetElemGloY;
+    lMeanZ += lDetElemGloZ;
+    lSigmaZ += lDetElemGloZ*lDetElemGloZ;
+    lNDetElem++;
 
-  /// Get information of current detection element
-  Double_t lDetElemLocX = 0.;
-  Double_t lDetElemLocY = 0.;
-  Double_t lDetElemLocZ = 0.;
-  fDetElemId = fCluster->GetDetElemId();
-  fDetElemNumber = fDetElemId%100;
-  for (int iCh=0; iCh<fDetElemId/100-1; iCh++)
-  { fDetElemNumber += fgNDetElemCh[iCh]; }
-
-  fTransform->Local2Global(
-    fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
-    fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
+  }
 
-}
+  // calculate mean values
+  lMeanY /= lNDetElem;
+  lSigmaY /= lNDetElem;
+  lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
+  lMeanZ /= lNDetElem;
+  lSigmaZ /= lNDetElem;
+  lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
+  AliInfo( Form( "Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ ) );
+
+  // create all possible arrays
+  Array fConstraintX[4];  //Array for constraint equation X
+  Array fConstraintY[4];  //Array for constraint equation Y
+  Array fConstraintP[4];  //Array for constraint equation P
+  Array fConstraintXZ[4];  //Array for constraint equation X vs Z
+  Array fConstraintYZ[4];  //Array for constraint equation Y vs Z
+  Array fConstraintPZ[4];  //Array for constraint equation P vs Z
+
+  // do we really need these ?
+  Array fConstraintXY[4];  //Array for constraint equation X vs Y
+  Array fConstraintYY[4];  //Array for constraint equation Y vs Y
+  Array fConstraintPY[4];  //Array for constraint equation P vs Y
+
+  // fill Bool_t sides array based on masks, for convenience
+  Bool_t lDetTLBR[4];
+  lDetTLBR[0] = sidesMask & SideTop;
+  lDetTLBR[1] = sidesMask & SideLeft;
+  lDetTLBR[2] = sidesMask & SideBottom;
+  lDetTLBR[3] = sidesMask & SideRight;
+
+  for( Int_t i = 0; i < fgNDetElem; ++i )
+  {
 
-//_____________________________________________________
-AliMUONAlignmentTrackRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment )
-{
-  /// Process track and set Local Equations
-  // store current track in running member.
-  fTrack = track;
+    // get chamber matching detector
+    const Int_t iCh( GetChamberId(i) );
 
-  // clear track record
-  fTrackRecord.Clear();
+    // skip detector if chamber is off
+    if( !lChOnOff[iCh-1] ) continue;
 
-  // get number of tracks
-  Int_t nTrackParam = fTrack->GetTrackParamAtCluster()->GetEntries();
-  AliDebug(1,Form("Number of track param entries : %i ", nTrackParam));
+    // get detector element id from detector element number
+    const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
+    const Int_t lDetElemId = iCh*100+lDetElemNumber;
 
-  Bool_t first( kTRUE );
-  for(Int_t iCluster=0; iCluster<nTrackParam; iCluster++)
-  {
+    // get global x, y and z position
+    Double_t lDetElemGloX = 0.;
+    Double_t lDetElemGloY = 0.;
+    Double_t lDetElemGloZ = 0.;
+    fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
 
-    // and get new pointers
-    fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtCluster()->At(iCluster);
-    if ( ! fTrackParam ) continue;
-    fCluster = fTrackParam->GetClusterPtr();
-    if ( ! fCluster ) continue;
-    // if (fDetElemId<500) continue;
+    // loop over sides
+    for( Int_t iSide = 0; iSide < 4; iSide++ )
+    {
 
-    // fill local variables for this position --> one measurement
-    FillDetElemData();
-    FillRecPointData();
-    FillTrackParamData();
+      // skip if side is not selected
+      if( !lDetTLBR[iSide] ) continue;
 
-    if( first )
-    {
+      // skip detector if it is not in the selected side
+      // stations 1 and 2
+      if( iCh>=1 && iCh<=4 )
+      {
+        if( lDetElemNumber == 0 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+        if( lDetElemNumber == 1 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+        if( lDetElemNumber == 2 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+        if( lDetElemNumber == 3 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
+      }
 
-      // for first valid cluster, save track position as "starting" values
-      first = kFALSE;
+      // station 3
+      if (iCh>=5 && iCh<=6)
+      {
+        if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+        if( lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+        if( lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+        if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
+      }
 
-      fTrackPos0[0] = fTrackPos[0];
-      fTrackPos0[1] = fTrackPos[1];
-      fTrackPos0[2] = fTrackPos[2];
-      fTrackSlope0[0] = fTrackSlope[0];
-      fTrackSlope0[1] = fTrackSlope[1];
+      // stations 4 and 5
+      if (iCh>=7 && iCh<=10)
+      {
+        if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+        if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+        if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+        if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
+      }
 
-    }
+      // constrain x
+      if( lVarXYT[0] ) fConstraintX[iSide].values[i*fgNParCh+0] = 1;
 
-    // calculate measurements
-    fCosPhi = TMath::Cos(fPhi);
-    fSinPhi = TMath::Sin(fPhi);
-    if( fBFieldOn )
-    {
+      // constrain y
+      if( lVarXYT[1] ) fConstraintY[iSide].values[i*fgNParCh+1] = 1;
 
-      fMeas[0] = fTrackPos[0] - fClustPos[0];
-      fMeas[1] = fTrackPos[1] - fClustPos[1];
+      // constrain phi (rotation around z)
+      if( lVarXYT[2] ) fConstraintP[iSide].values[i*fgNParCh+2] = 1;
 
-    } else {
+      // x-z shearing
+      if( lVarXYT[3] ) fConstraintXZ[iSide].values[i*fgNParCh+0] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
 
-      fMeas[0] = - fClustPos[0];
-      fMeas[1] = - fClustPos[1];
+      // y-z shearing
+      if( lVarXYT[4] ) fConstraintYZ[iSide].values[i*fgNParCh+1] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
 
-    }
+      // phi-z shearing
+      if( lVarXYT[5] ) fConstraintPZ[iSide].values[i*fgNParCh+2] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
 
-    // soùe debugging
-    AliDebug(1,Form("cluster: %i", iCluster));
-    AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
-    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));
+      // x-y shearing
+      if( lVarXYT[6] ) fConstraintXY[iSide].values[i*fgNParCh+0] = (lDetElemGloY-lMeanY)/lSigmaY;
 
-    AliDebug(1,Form("Track Parameter: %i", iCluster));
-    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]));
-    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]));
+      // y-y shearing
+      if( lVarXYT[7] ) fConstraintYY[iSide].values[i*fgNParCh+1] = (lDetElemGloY-lMeanY)/lSigmaY;
 
-    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]));
+      // phi-y shearing
+      if( lVarXYT[8] ) fConstraintPY[iSide].values[i*fgNParCh+2] = (lDetElemGloY-lMeanY)/lSigmaY;
 
-    // Set local equations
-    LocalEquationX( doAlignment );
-    LocalEquationY( doAlignment );
+    }
 
   }
 
-  return &fTrackRecord;
+  // pass constraints to millepede
+  for( Int_t iSide = 0; iSide < 4; iSide++ )
+  {
+    // skip if side is not selected
+    if( !lDetTLBR[iSide] ) continue;
+
+    if( lVarXYT[0] ) AddConstraint(fConstraintX[iSide].values,0.0);
+    if( lVarXYT[1] ) AddConstraint(fConstraintY[iSide].values,0.0);
+    if( lVarXYT[2] ) AddConstraint(fConstraintP[iSide].values,0.0);
+    if( lVarXYT[3] ) AddConstraint(fConstraintXZ[iSide].values,0.0);
+    if( lVarXYT[4] ) AddConstraint(fConstraintYZ[iSide].values,0.0);
+    if( lVarXYT[5] ) AddConstraint(fConstraintPZ[iSide].values,0.0);
+    if( lVarXYT[6] ) AddConstraint(fConstraintXY[iSide].values,0.0);
+    if( lVarXYT[7] ) AddConstraint(fConstraintYY[iSide].values,0.0);
+    if( lVarXYT[8] ) AddConstraint(fConstraintPY[iSide].values,0.0);
+  }
 
 }
 
-//______________________________________________________________________________
-void AliMUONAlignment::ProcessTrack( AliMUONAlignmentTrackRecord* track, Bool_t doAlignment )
+//______________________________________________________________________
+void AliMUONAlignment::InitGlobalParameters(Double_t *par)
 {
-  /// Process track (from record) and set Local Equations
-  if( !( track && doAlignment ) ) return;
+  /// Initialize global parameters with par array
+  if( !fInitialized )
+  { AliFatal( "Millepede is not initialized" ); }
 
-  // loop over clusters
-  for( Int_t index = 0; index < track->GetNRecords(); ++index )
-  { if( AliMUONAlignmentClusterRecord* record = track->GetRecord( index ) ) LocalEquation( *record ); }
+  fMillepede->SetGlobalParameters(par);
+}
 
-  return;
+//______________________________________________________________________
+void AliMUONAlignment::SetAllowedVariation( Int_t iPar, Double_t value )
+{
+  /// "Encouraged" variation for degrees of freedom
+  // check initialization
+  if( fInitialized )
+  { AliFatal( "Millepede already initialized" ); }
+
+  // check initialization
+  if( !(iPar >= 0 && iPar < fgNParCh ) )
+  { AliFatal( Form( "Invalid index: %i", iPar ) ); }
 
+  fAllowVar[iPar] = value;
 }
 
-//_____________________________________________________
-void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit)
+//______________________________________________________________________
+void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
 {
 
-  /// Call local fit for this tracks
-  Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
-  if (iRes && !lSingleFit)
-  { fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1); }
+  /// Set expected measurement resolution
+  fSigma[0] = sigmaX;
+  fSigma[1] = sigmaY;
+
+  // print
+  for( Int_t i=0; i<2; ++i )
+  { AliInfo( Form( "fSigma[%i]=%f", i, fSigma[i] ) ); }
 
 }
 
 //_____________________________________________________
-void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
+void AliMUONAlignment::GlobalFit( Double_t *parameters, Double_t *errors, Double_t *pulls )
 {
 
   /// Call global fit; Global parameters are stored in parameters
-  fMillepede->GlobalFit(parameters,errors,pulls);
-
-  AliInfo("Done fitting global parameters!");
-  for (int iGlob=0; iGlob<fgNDetElem; iGlob++)
-  { printf("%d\t %f\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+3],parameters[iGlob*fgNParCh+2]); }
+  fMillepede->GlobalFit( parameters, errors, pulls );
 
+  AliInfo( "Done fitting global parameters" );
+  for( int iDet=0; iDet<fgNDetElem; ++iDet )
+  {
+    AliInfo( Form( "%d\t %f\t %f\t %f\t %f",
+      iDet,
+      parameters[iDet*fgNParCh+0],parameters[iDet*fgNParCh+1],
+      parameters[iDet*fgNParCh+3],parameters[iDet*fgNParCh+2]
+      ) );
+  }
 
 }
 
 //_____________________________________________________
-Double_t AliMUONAlignment::GetParError(Int_t iPar)
-{
-  /// Get error of parameter iPar
-  Double_t lErr = fMillepede->GetParError(iPar);
-  return lErr;
-}
+void AliMUONAlignment::PrintGlobalParameters() const
+{ fMillepede->PrintGlobalParameters(); }
 
 //_____________________________________________________
-void AliMUONAlignment::PrintGlobalParameters()
-{
-  /// Print global parameters
-  fMillepede->PrintGlobalParameters();
-}
-
-//_________________________________________________________________________
-TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, const double *lMisAlignment) const
-{
-  /// Realign given transformation by given misalignment and return the misaligned transformation
-
-  Double_t cartMisAlig[3] = {0,0,0};
-  Double_t angMisAlig[3] = {0,0,0};
-//   const Double_t *trans = transform.GetTranslation();
-//   TGeoRotation *rot;
-//   // check if the rotation we obtain is not NULL
-//   if (transform.GetRotation()) {
-//     rot = transform.GetRotation();
-//   }
-//   else {
-//     rot = new TGeoRotation("rot");
-//   }                 // default constructor.
-
-  cartMisAlig[0] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0])*lMisAlignment[0];
-  cartMisAlig[1] = -TMath::Sign(1.0,transform.GetRotationMatrix()[4])*lMisAlignment[1];
-  cartMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[8])*lMisAlignment[3];
-  angMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0]*transform.GetRotationMatrix()[4])*lMisAlignment[2]*180./TMath::Pi();
-
-  TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
-  TGeoRotation deltaRot;
-  deltaRot.RotateX(angMisAlig[0]);
-  deltaRot.RotateY(angMisAlig[1]);
-  deltaRot.RotateZ(angMisAlig[2]);
-
-  TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
-  TGeoHMatrix newTransfMat = transform * deltaTransf;
-
-  return TGeoCombiTrans(newTransfMat);
-}
+Double_t AliMUONAlignment::GetParError(Int_t iPar) const
+{ return fMillepede->GetParError(iPar); }
 
 //______________________________________________________________________
-AliMUONGeometryTransformer *
-AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
-          const double *misAlignments, Bool_t verbose)
-
+AliMUONGeometryTransformer* AliMUONAlignment::ReAlign(
+  const AliMUONGeometryTransformer * transformer,
+  const double *misAlignments, Bool_t )
 {
+
   /// Returns a new AliMUONGeometryTransformer with the found misalignments
   /// applied.
 
@@ -1340,42 +900,64 @@ AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
   // Adds the new module transformer to a new geometry transformer
   // Returns the new geometry transformer
 
-  Double_t lModuleMisAlignment[4] = {0.,0.,0.,0.};
-  Double_t lDetElemMisAlignment[4] = {0.,0.,0.,0.};
-  Int_t iDetElemId = 0;
-  Int_t iDetElemNumber = 0;
+  Double_t lModuleMisAlignment[fgNParCh] = {0};
+  Double_t lDetElemMisAlignment[fgNParCh] = {0};
+  const TClonesArray* oldMisAlignArray( transformer->GetMisAlignmentData() );
+
+  AliMUONGeometryTransformer *newGeometryTransformer = new AliMUONGeometryTransformer();
+  for( Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt )
+  {
 
-  AliMUONGeometryTransformer *newGeometryTransformer =
-    new AliMUONGeometryTransformer();
-  for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
     // module transformers
-    const AliMUONGeometryModuleTransformer *kModuleTransformer =
-      transformer->GetModuleTransformer(iMt, true);
+    const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
 
-    AliMUONGeometryModuleTransformer *newModuleTransformer =
-      new AliMUONGeometryModuleTransformer(iMt);
+    AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
 
-    TGeoCombiTrans moduleTransform =
-      TGeoCombiTrans(*kModuleTransformer->GetTransformation());
-    // New module transformation
-    TGeoCombiTrans newModuleTransform = ReAlign(moduleTransform,lModuleMisAlignment);
+    // get transformation
+    TGeoHMatrix deltaModuleTransform( DeltaTransform( lModuleMisAlignment ) );
+
+    // update module
+    TGeoHMatrix moduleTransform( *kModuleTransformer->GetTransformation() );
+    TGeoHMatrix newModuleTransform( AliMUONGeometryBuilder::Multiply( deltaModuleTransform, moduleTransform ) );
     newModuleTransformer->SetTransformation(newModuleTransform);
 
-    // Get delta transformation:
-    // Tdelta = Tnew * Told.inverse
-    TGeoHMatrix deltaModuleTransform =
-      AliMUONGeometryBuilder::Multiply(newModuleTransform,
-               kModuleTransformer->GetTransformation()->Inverse());
+    // Get matching old alignment and update current matrix accordingly
+    if( oldMisAlignArray )
+    {
+
+      const AliAlignObjMatrix* oldAlignObj(0);
+      const Int_t moduleId( kModuleTransformer->GetModuleId() );
+      const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId );
+      for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
+      {
+
+        const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
+        if( localAlignObj && localAlignObj->GetVolUID() == volId )
+        {
+          oldAlignObj = localAlignObj;
+          break;
+        }
+
+      }
+
+      // multiply
+      if( oldAlignObj )
+      {
+
+        TGeoHMatrix oldMatrix;
+        oldAlignObj->GetMatrix( oldMatrix );
+        deltaModuleTransform.Multiply( &oldMatrix );
+
+      }
+
+    }
+
     // Create module mis alignment matrix
-    newGeometryTransformer
-      ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
+    newGeometryTransformer ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
 
     AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
 
-    if (verbose)
-      AliInfo(Form("%i DEs in old GeometryStore  %i",detElements->GetSize(), iMt));
-
     TIter next(detElements->CreateIterator());
     AliMUONGeometryDetElement* detElement;
     Int_t iDe(-1);
@@ -1383,114 +965,432 @@ AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
     {
       ++iDe;
       // make a new detection element
-      AliMUONGeometryDetElement *newDetElement =
-  new AliMUONGeometryDetElement(detElement->GetId(),
-              detElement->GetVolumePath());
+      AliMUONGeometryDetElement *newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
       TString lDetElemName(detElement->GetDEName());
       lDetElemName.ReplaceAll("DE","");
-      iDetElemId = lDetElemName.Atoi();
-      iDetElemNumber = iDetElemId%100;
-      for (int iCh=0; iCh<iDetElemId/100-1; iCh++){
-  iDetElemNumber += fgNDetElemCh[iCh];
+
+      // store detector element id and number
+      const Int_t iDetElemId = lDetElemName.Atoi();
+      if( !DetElemIsValid( iDetElemId ) )
+      {
+        AliInfo( Form( "Skipping invalid detector element %i", iDetElemId ) );
+        continue;
       }
-      for (int i=0; i<fgNParCh; i++) {
-  lDetElemMisAlignment[i] = 0.0;
-  if (iMt<fgNTrkMod) {
-    AliInfo(Form("iMt %i, iCh %i, iDe %i, iDeId %i, iDeNb %i, iPar %i",iMt, iDetElemId/100, iDe, iDetElemId, iDetElemNumber, iDetElemNumber*fgNParCh+i));
-    lDetElemMisAlignment[i] =  misAlignments[iDetElemNumber*fgNParCh+i];
-  }
+
+      const Int_t iDetElemNumber( GetDetElemNumber( iDetElemId ) );
+
+      for( int i=0; i<fgNParCh; ++i )
+      {
+        lDetElemMisAlignment[i] = 0.0;
+        if( iMt<fgNTrkMod ) { lDetElemMisAlignment[i] =  misAlignments[iDetElemNumber*fgNParCh+i]; }
       }
-      // local transformation of this detection element.
-      TGeoCombiTrans localTransform
-  = TGeoCombiTrans(*detElement->GetLocalTransformation());
-      TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
-      newDetElement->SetLocalTransformation(newLocalTransform);
-
-      // global transformation
-      TGeoHMatrix newGlobalTransform =
-  AliMUONGeometryBuilder::Multiply(newModuleTransform,
-           newLocalTransform);
-      newDetElement->SetGlobalTransformation(newGlobalTransform);
-
-      // add this det element to module
-      newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
-                  newDetElement);
-
-      // In the Alice Alignment Framework misalignment objects store
-      // global delta transformation
-      // Get detection "intermediate" global transformation
-      TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform;
-      // Get detection element global delta transformation:
-      // Tdelta = Tnew * Told.inverse
-      TGeoHMatrix  deltaGlobalTransform
-  = AliMUONGeometryBuilder::Multiply(newGlobalTransform,
-             newOldGlobalTransform.Inverse());
-
-      // Create mis alignment matrix
-      newGeometryTransformer
-  ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
+
+      // get transformation
+      TGeoHMatrix deltaGlobalTransform( DeltaTransform( lDetElemMisAlignment ) );
+
+      // update module
+      TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
+      TGeoHMatrix newGlobalTransform( AliMUONGeometryBuilder::Multiply( deltaGlobalTransform, globalTransform ) );
+      newDetElement->SetGlobalTransformation( newGlobalTransform );
+      newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
+
+      // Get matching old alignment and update current matrix accordingly
+      if( oldMisAlignArray )
+      {
+
+        const AliAlignObjMatrix* oldAlignObj(0);
+        const int detElemId( detElement->GetId() );
+        const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
+        for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
+        {
+
+          const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
+          if( localAlignObj && localAlignObj->GetVolUID() == volId )
+          {
+            oldAlignObj = localAlignObj;
+            break;
+          }
+
+        }
+
+        // multiply
+        if( oldAlignObj )
+        {
+
+          TGeoHMatrix oldMatrix;
+          oldAlignObj->GetMatrix( oldMatrix );
+          deltaGlobalTransform.Multiply( &oldMatrix );
+
+        }
+
+      }
+
+      // Create misalignment matrix
+      newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
+
     }
 
-    if (verbose)
-      AliInfo(Form("Added module transformer %i to the transformer", iMt));
     newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
   }
+
   return newGeometryTransformer;
+
 }
 
 //______________________________________________________________________
-void AliMUONAlignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY)
+void AliMUONAlignment::SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY )
 {
 
   /// Set alignment resolution to misalign objects to be stored in CDB
-  Int_t chIdMin = (rChId<0)? 0 : rChId;
-  Int_t chIdMax = (rChId<0)? 9 : rChId;
-  Double_t chResX = rChResX;
-  Double_t chResY = rChResY;
-  Double_t deResX = rDeResX;
-  Double_t deResY = rDeResY;
-
+  /// if rChId is > 0 set parameters for this chamber only, counting from 1
   TMatrixDSym mChCorrMatrix(6);
   mChCorrMatrix[0][0]=chResX*chResX;
   mChCorrMatrix[1][1]=chResY*chResY;
-  //  mChCorrMatrix.Print();
 
   TMatrixDSym mDECorrMatrix(6);
   mDECorrMatrix[0][0]=deResX*deResX;
   mDECorrMatrix[1][1]=deResY*deResY;
-  //  mDECorrMatrix.Print();
 
   AliAlignObjMatrix *alignMat = 0x0;
 
-  for(Int_t chId=chIdMin; chId<=chIdMax; chId++) {
+  for( Int_t chId = 0; chId <= 9; ++chId )
+  {
+
+    // skip chamber if selection is valid, and does not match
+    if( rChId > 0 && chId+1 != rChId ) continue;
+
     TString chName1;
     TString chName2;
-    if (chId<4){
+    if (chId<4)
+    {
+
       chName1 = Form("GM%d",chId);
       chName2 = Form("GM%d",chId);
+
     } else {
+
       chName1 = Form("GM%d",4+(chId-4)*2);
       chName2 = Form("GM%d",4+(chId-4)*2+1);
+
     }
 
-    for (int i=0; i<misAlignArray->GetEntries(); i++) {
+    for( int i=0; i<misAlignArray->GetEntries(); ++i )
+    {
+
       alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
       TString volName(alignMat->GetSymName());
       if((volName.Contains(chName1)&&
-    ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
-     (volName.Length()==volName.Index(chName1)+chName1.Length())))||
-   (volName.Contains(chName2)&&
-    ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
-     (volName.Length()==volName.Index(chName2)+chName2.Length())))){
-  volName.Remove(0,volName.Last('/')+1);
-  if (volName.Contains("GM")) {
-    // alignMat->Print("NULL");
-    alignMat->SetCorrMatrix(mChCorrMatrix);
-  } else if (volName.Contains("DE")) {
-    // alignMat->Print("NULL");
-    alignMat->SetCorrMatrix(mDECorrMatrix);
-  }
+        ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
+        (volName.Length()==volName.Index(chName1)+chName1.Length())))||
+        (volName.Contains(chName2)&&
+        ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
+        (volName.Length()==volName.Index(chName2)+chName2.Length()))))
+      {
+
+        volName.Remove(0,volName.Last('/')+1);
+        if (volName.Contains("GM")) alignMat->SetCorrMatrix(mChCorrMatrix);
+        else if (volName.Contains("DE")) alignMat->SetCorrMatrix(mDECorrMatrix);
+
       }
+
     }
+
+  }
+
+}
+
+
+//_____________________________________________________
+void AliMUONAlignment::FillDetElemData( AliMUONVCluster* cluster )
+{
+
+  /// Get information of current detection element
+  // get detector element number from Alice ID
+  const Int_t detElemId = cluster->GetDetElemId();
+  fDetElemNumber = GetDetElemNumber( detElemId );
+
+  // get detector element
+  const AliMUONGeometryDetElement* detElement = fTransform->GetDetElement( detElemId );
+
+  /*
+  get the global transformation matrix and store its inverse, in order to manually perform
+  the global to Local transformations needed to calculate the derivatives
+  */
+  fGeoCombiTransInverse = detElement->GetGlobalTransformation()->Inverse();
+
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::FillRecPointData( AliMUONVCluster* cluster )
+{
+
+  /// Get information of current cluster
+  fClustPos[0] = cluster->GetX();
+  fClustPos[1] = cluster->GetY();
+  fClustPos[2] = cluster->GetZ();
+
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::FillTrackParamData( AliMUONTrackParam* trackParam )
+{
+
+  /// Get information of current track at current cluster
+  fTrackPos[0] = trackParam->GetNonBendingCoor();
+  fTrackPos[1] = trackParam->GetBendingCoor();
+  fTrackPos[2] = trackParam->GetZ();
+  fTrackSlope[0] = trackParam->GetNonBendingSlope();
+  fTrackSlope[1] = trackParam->GetBendingSlope();
+
+}
+
+//______________________________________________________________________
+Bool_t AliMUONAlignment::UnbiasTrackParamData( AliMUONTrackParam* trackParam ) const
+{
+
+  /**
+  calculate unbiased track parameters at a given detector, that is,
+  after taking out the contribution of the detector's cluster from the track
+  */
+
+  // check track parameters
+  if( !trackParam ) return kFALSE;
+
+  // Remove cluster contibution from smoothed track param
+  TMatrixD smoothParameters( trackParam->GetSmoothParameters() );
+  TMatrixD smoothCovariances( trackParam->GetSmoothCovariances() );
+
+  AliMUONVCluster* cluster = trackParam->GetClusterPtr();
+  // p' = p + K(m - H*p)
+  // K  = C H (-V + H C H^t)^-1
+  // C' = C - K H C
+  // where p,C are smoothed param,cov at cluster position
+  // H is the matrix converting the state vector to measurement
+  // V is the measurement cov.matrix
+  // m is the measurement vector
+  static TMatrixD H(2,5);
+  H.Zero();
+  H(0,0)=H(1,2) = 1.;
+
+  // (-Vk+H_k C^n_k H_k^T)^-1
+  static TMatrixD df(2,2);
+  df.Zero();
+  df(0,0) = smoothCovariances(0,0) - cluster->GetErrX2();
+  df(1,1) = smoothCovariances(2,2) - cluster->GetErrY2();
+  df(0,1) = smoothCovariances(0,2);
+  df(1,0) = smoothCovariances(2,0);
+
+  if (df.Determinant() != 0) df.Invert();
+  else {
+    AliInfo( "Determinant = 0\n" );
+    return kFALSE;
+  }
+
+  // gain matrix
+  TMatrixD kTmp( smoothCovariances, TMatrixD::kMultTranspose, H );
+  TMatrixD K(kTmp, TMatrixD::kMult, df);
+
+  TMatrixD dfc(2,1);
+  dfc.Zero();
+  dfc(0,0) = cluster->GetX()-smoothParameters(0,0);
+  dfc(1,0) = cluster->GetY()-smoothParameters(2,0);
+  TMatrixD tmp0(K,TMatrixD::kMult, dfc);
+  smoothParameters += tmp0;
+
+  TMatrixD tmp1(K,   TMatrixD::kMult, H);
+  TMatrixD tmp2(tmp1,TMatrixD::kMult, smoothCovariances);
+  smoothCovariances -= tmp2;
+
+  // update track parameters
+  trackParam->SetParameters( smoothParameters );
+  trackParam->SetCovariances( smoothCovariances );
+
+  return kTRUE;
+
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::LocalEquationX( void )
+{
+  /// local equation along X
+
+  // 'inverse' (GlobalToLocal) rotation matrix
+  const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
+
+  // local derivatives
+  SetLocalDerivative( 0, r[0] );
+  SetLocalDerivative( 1, r[0]*(fTrackPos[2] - fTrackPos0[2]) );
+  SetLocalDerivative( 2, r[1] );
+  SetLocalDerivative( 3, r[1]*(fTrackPos[2] - fTrackPos0[2]) );
+
+  // global derivatives
+  /*
+  alignment parameters are
+  0: delta_x
+  1: delta_y
+  2: delta_phiz
+  3: delta_z
+  */
+
+  SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[0] );
+  SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[1] );
+
+  if( fBFieldOn )
+  {
+
+    // use local position for derivatives vs 'delta_phi_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*fTrackPos[0] + r[0]*fTrackPos[1] );
+
+    // use local slopes for derivatives vs 'delta_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope[0] + r[1]*fTrackSlope[1] );
+
+  } else {
+
+    // local copy of extrapolated track positions
+    const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
+    const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
+
+    // use properly extrapolated position for derivatives vs 'delta_phi_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*trackPosX + r[0]*trackPosY );
+
+    // use slopes at origin for derivatives vs 'delta_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope0[0] + r[1]*fTrackSlope0[1] );
+
+  }
+
+  // store local equation
+  fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0] );
+
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::LocalEquationY( void )
+{
+  /// local equation along Y
+
+  // 'inverse' (GlobalToLocal) rotation matrix
+  const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
+
+  // store local derivatives
+  SetLocalDerivative( 0, r[3] );
+  SetLocalDerivative( 1, r[3]*(fTrackPos[2] - fTrackPos0[2] ) );
+  SetLocalDerivative( 2, r[4] );
+  SetLocalDerivative( 3, r[4]*(fTrackPos[2] - fTrackPos0[2] ) );
+
+  // set global derivatives
+  SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[3]);
+  SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[4]);
+
+  if( fBFieldOn )
+  {
+
+    // use local position for derivatives vs 'delta_phi'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*fTrackPos[0] + r[3]*fTrackPos[1]);
+
+    // use local slopes for derivatives vs 'delta_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope[0]+r[4]*fTrackSlope[1] );
+
+  } else {
+
+    // local copy of extrapolated track positions
+    const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
+    const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
+
+    // use properly extrapolated position for derivatives vs 'delta_phi'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*trackPosX + r[3]*trackPosY );
+
+    // use slopes at origin for derivatives vs 'delta_z'
+    SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope0[0]+r[4]*fTrackSlope0[1] );
+
   }
+
+  // store local equation
+  fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1] );
+
+}
+
+//_________________________________________________________________________
+TGeoCombiTrans AliMUONAlignment::DeltaTransform( const double *lMisAlignment) const
+{
+  /// Get Delta Transformation, based on alignment parameters
+
+  // translation
+  const TGeoTranslation deltaTrans( lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
+
+  // rotation
+  TGeoRotation deltaRot;
+  deltaRot.RotateZ(lMisAlignment[2]*180./TMath::Pi());
+
+  // combined rotation and translation.
+  return TGeoCombiTrans(deltaTrans,deltaRot);
+
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value)
+{
+  /// Constrain equation defined by par to value
+  if( !fInitialized )
+  { AliFatal( "Millepede is not initialized" ); }
+
+  fMillepede->SetGlobalConstraint(par, value);
+}
+
+//______________________________________________________________________
+Bool_t AliMUONAlignment::DetElemIsValid( Int_t iDetElemId ) const
+{
+  /// return true if given detector element is valid (and belongs to muon tracker)
+  const Int_t iCh = iDetElemId/100;
+  const Int_t iDet = iDetElemId%100;
+  return ( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] );
+}
+
+//______________________________________________________________________
+Int_t AliMUONAlignment::GetDetElemNumber( Int_t iDetElemId ) const
+{
+  /// get det element number from ID
+  // get chamber and element number in chamber
+  const Int_t iCh = iDetElemId/100;
+  const Int_t iDet = iDetElemId%100;
+
+  // make sure detector index is valid
+  if( !( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] ) )
+  { AliFatal( Form( "Invalid detector element id: %i", iDetElemId ) ); }
+
+  // add number of detectors up to this chamber
+  return iDet + fgSNDetElemCh[iCh-1];
+
+}
+
+//______________________________________________________________________
+Int_t AliMUONAlignment::GetChamberId( Int_t iDetElemNumber ) const
+{
+  /// get chamber (counting from 1) matching a given detector element id
+  Int_t iCh( 0 );
+  for( iCh=0; iCh<fgNCh; iCh++ )
+  { if( iDetElemNumber < fgSNDetElemCh[iCh] ) break; }
+
+  return iCh;
+}
+
+//______________________________________________________________________
+TString AliMUONAlignment::GetParameterMaskString( UInt_t mask ) const
+{
+  TString out;
+  if( mask & ParX ) out += "X";
+  if( mask & ParY ) out += "Y";
+  if( mask & ParZ ) out += "Z";
+  if( mask & ParTZ ) out += "T";
+  return out;
+}
+
+//______________________________________________________________________
+TString AliMUONAlignment::GetSidesMaskString( UInt_t mask ) const
+{
+  TString out;
+  if( mask & SideTop ) out += "T";
+  if( mask & SideLeft ) out += "L";
+  if( mask & SideBottom ) out += "B";
+  if( mask & SideRight ) out += "R";
+  return out;
 }
index f398cca..4756956 100644 (file)
@@ -1,24 +1,24 @@
 #ifndef ALIMUONALIGNMENT_H
 #define ALIMUONALIGNMENT_H
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
+* See cxx source for full Copyright notice                               */
 
-/* $Id$ */
+/* $Id: AliMUONAlignment.h 50858 2011-07-29 10:58:20Z ivana $ */
 
 /// \ingroup rec
 /// \class AliMUONAlignment
 /// \brief Class for alignment of muon spectrometer
 //
-// Authors: Bruce Becker, Javier Castillo
+// Authors: Bruce Becker, Javier Castillo, Hugo Pereira Da Costa
+
+#include "AliMillePedeRecord.h"
 
 #include <TObject.h>
 #include <TString.h>
+#include <TGeoMatrix.h>
 
-#include "AliMUONAlignmentRecord.h"
-
-class TGeoCombiTrans;
 class TClonesArray;
-class AliMillepede;
+class AliMillePede2;
 class AliMUONGeometryTransformer;
 class AliMUONTrack;
 class AliMUONTrackParam;
@@ -29,98 +29,182 @@ class AliMUONAlignment:public TObject
 
   public:
 
-
-  /// constructor
   AliMUONAlignment();
-       /// rootio constructor
-       AliMUONAlignment(TRootIOCtor* dummy);
 
-  /// destructor
   virtual ~AliMUONAlignment();
 
-  /// process track for alignment minimization
-  /**
-  returns the alignment records for this track.
-  They can be stored in some output for later reprocessing.
-  */
-  AliMUONAlignmentTrackRecord* ProcessTrack(AliMUONTrack *track, Bool_t doAlignment );
+  void Init( void );
 
-  /// process track record
-  void ProcessTrack( AliMUONAlignmentTrackRecord* track, Bool_t doAlignment = kTRUE );
+  // array dimendions
+  enum
+  {
+    /// Number tracking stations
+    fgNSt = 5,
 
-  /// returns current track record
-  AliMUONAlignmentTrackRecord* GetTrackRecord( void )
-  { return &fTrackRecord; }
+    /// Number tracking chambers
+    fgNCh = 10,
 
-  /// Set geometry transformer
-  void SetGeometryTransformer(AliMUONGeometryTransformer * transformer)
-  { fTransform = transformer; }
+    /// Number of tracking modules (4 ch + 6*2 half-ch)
+    fgNTrkMod = 16,
 
-  /// Set flag for Magnetic field On/Off
-  void SetBFieldOn(Bool_t bBFieldOn)
-  { fBFieldOn =  bBFieldOn; }
+    /// Total number of detection elements
+    /// (4*2 + 4*2 + 18*2 + 26*2 + 26*2)
+    fgNDetElem = 156,
+
+    /// Number of local parameters
+    fNLocal = 4,
+
+    /// Number of degrees of freedom per chamber
+    fgNParCh = 4,
+
+    /// Number of global parameters
+    fNGlobal = fgNParCh*fgNDetElem
+  };
 
-  /// Define chambers to align
-  void SetChOnOff(Bool_t *bChOnOff)
+  /// Number of detection elements per chamber
+  static const Int_t fgNDetElemCh[fgNCh];
+
+  /// Sum of detection elements up to this chamber
+  static const Int_t fgSNDetElemCh[fgNCh+1];
+
+  /// global parameter bit set, used for masks
+  enum ParameterMask
   {
-    for(int iCh=0; iCh<10; iCh++)
-    { fChOnOff[iCh] =  bChOnOff[iCh]; }
-  }
+    ParX = 1<<0,
+    ParY = 1<<1,
+    ParZ = 1<<2,
+    ParTZ = 1<<3,
+
+    ParAllTranslations = ParX|ParY|ParZ,
+    ParAllRotations = ParTZ,
+    ParAll = ParAllTranslations|ParAllRotations
 
-  /// Possibility to align only one side of the detector
-  void SetSpecLROnOff(Bool_t *bSpecLROnOff)
+  };
+
+  /// detector sides bit set, used for selecting sides in constrains
+  enum SidesMask
   {
-    fSpecLROnOff[0] =  bSpecLROnOff[0];
-    fSpecLROnOff[1] =  bSpecLROnOff[1];
-  }
+    SideTop = 1<<0,
+    SideLeft = 1<<1,
+    SideBottom = 1<<2,
+    SideRight = 1<<3,
+    SideTopLeft = SideTop|SideLeft,
+    SideTopRight = SideTop|SideRight,
+    SideBottomLeft = SideBottom|SideLeft,
+    SideBottomRight = SideBottom|SideRight,
+    AllSides = SideTop|SideBottom|SideLeft|SideRight
+  };
 
-  void FixStation(Int_t iSt);
+  AliMillePedeRecord* ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight = 1 );
 
-  void FixChamber(Int_t iCh);
+  void ProcessTrack( AliMillePedeRecord* );
 
-  void FixDetElem(Int_t iDetElemId, TString sVarXYT = "XYTZ");
+  //@name modifiers
+  //@{
 
-  void FixHalfSpectrometer(const Bool_t* bChOnOff, const Bool_t* bSpecLROnOff);
+  /// run number
+  void SetRunNumber( Int_t id )
+  {fRunNumber = id;}
 
-  void AllowVariations(const Bool_t* bChOnOff);
+  /// Set flag for Magnetic field On/Off
+  void SetBFieldOn( Bool_t value )
+  { fBFieldOn = value; }
 
-  void SetNonLinear(const Bool_t *bChOnOff, const Bool_t *bVarXYT);
+  /// use unbiased residuals
+  void SetUnbias( Bool_t value )
+  { fUnbias = value; }
 
-  void AddConstraints(const Bool_t *bChOnOff, const Bool_t *bVarXYT);
+  void SetAllowedVariation( Int_t iPar, Double_t value );
 
-  void AddConstraints(const Bool_t *bChOnOff, const Bool_t *bVarXYT, const Bool_t *bDetTLBR, const Bool_t *bSpecLROnOff);
+  void SetSigmaXY(Double_t sigmaX, Double_t sigmaY);
 
-  void ResetConstraints();
+  /// Set geometry transformer
+  void SetGeometryTransformer( AliMUONGeometryTransformer * transformer )
+  { fTransform = transformer; }
 
-  void FixParameter(Int_t param, Double_t value);
+  //@}
 
-  void SetNonLinear(Int_t param);
+  //@name fixing detectors
+  //@{
 
-  void AddConstraint(Double_t *factor, Double_t value );
+  void FixAll( UInt_t parameterMask = ParAll );
 
-  void InitGlobalParameters(Double_t *par);
+  void FixChamber( Int_t iCh, UInt_t parameterMask = ParAll );
 
-  void SetSigmaXY(Double_t sigmaX, Double_t sigmaY);
+  void FixDetElem( Int_t iDetElemId, UInt_t parameterMask = ParAll );
 
-  /// Set array of local derivatives
-  void SetLocalDerivative(Int_t index, Double_t value)
-  { fLocalDerivatives[index] = value; }
+  void FixHalfSpectrometer( const Bool_t* bChOnOff, UInt_t sidesMask = AllSides, UInt_t parameterMask = ParAll );
 
-  /// Set array of global derivatives
-  void SetGlobalDerivative(Int_t index, Double_t value)
-  { fGlobalDerivatives[index] = value; }
+  void FixParameter( Int_t iPar );
+
+  void FixParameter( Int_t iDetElem, Int_t iPar )
+  { FixParameter( iDetElem*fgNParCh + iPar ); }
+
+  //@}
+
+  //@name releasing detectors
+  //@{
+
+  void ReleaseChamber( Int_t iCh, UInt_t parameterMask = ParAll );
+
+  void ReleaseDetElem( Int_t iDetElemId, UInt_t parameterMask = ParAll );
+
+  void ReleaseParameter( Int_t iPar );
+
+  void ReleaseParameter( Int_t iDetElem, Int_t iPar )
+  { ReleaseParameter( iDetElem*fgNParCh + iPar ); }
+
+  //@}
+
+  //@name grouping detectors
+  //@{
+
+  void GroupChamber( Int_t iCh, UInt_t parameterMask = ParAll );
+
+  void GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t parameterMask = ParAll );
+
+  void GroupDetElems( Int_t *detElemList, Int_t nDetElem, UInt_t parameterMask = ParAll );
+
+  //@}
+
+  //@name define non linearity
+  //@{
 
-  void LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit);
+  void SetChamberNonLinear( Int_t iCh, UInt_t parameterMask);
 
-  void GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls);
+  void SetDetElemNonLinear( Int_t iSt, UInt_t parameterMask);
 
-  void PrintGlobalParameters();
+  void SetParameterNonLinear( Int_t iPar );
 
-  Double_t GetParError(Int_t iPar);
+  void SetParameterNonLinear( Int_t iDetElem, Int_t iPar )
+  { SetParameterNonLinear( iDetElem*fgNParCh + iPar ); }
 
-  AliMUONGeometryTransformer* ReAlign(const AliMUONGeometryTransformer * transformer, const double *misAlignments, Bool_t verbose);
+  //@}
 
-  void SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t chId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY);
+  //@name constraints
+  //@{
+
+  void AddConstraints( const Bool_t *bChOnOff, UInt_t parameterMask );
+
+  void AddConstraints( const Bool_t *bChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask = AllSides );
+
+  //@}
+
+  /// initialize global parameters to a give set of values
+  void InitGlobalParameters( Double_t *par );
+
+  /// perform global fit
+  void GlobalFit( Double_t *parameters,Double_t *errors,Double_t *pulls );
+
+  /// print global parameters
+  void PrintGlobalParameters( void ) const;
+
+  /// get error on a given parameter
+  Double_t GetParError( Int_t iPar ) const;
+
+  AliMUONGeometryTransformer* ReAlign( const AliMUONGeometryTransformer * transformer, const double *misAlignments, Bool_t verbose );
+
+  void SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t chId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY );
 
   private:
 
@@ -130,141 +214,127 @@ class AliMUONAlignment:public TObject
   /// Not implemented
   AliMUONAlignment&  operator = (const AliMUONAlignment& right);
 
-  void Init(Int_t nGlobal, Int_t nLocal, Int_t nStdDev);
+  /// Set array of local derivatives
+  void SetLocalDerivative(Int_t index, Double_t value)
+  { fLocalDerivatives[index] = value; }
+
+  /// Set array of global derivatives
+  void SetGlobalDerivative(Int_t index, Double_t value)
+  { fGlobalDerivatives[index] = value; }
+
+  void FillDetElemData( AliMUONVCluster* );
+
+  void FillRecPointData( AliMUONVCluster* );
 
-  void ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t lWeight=1.0) const;
+  void FillTrackParamData( AliMUONTrackParam* );
 
-  void ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight=1.0) const;
+  Bool_t UnbiasTrackParamData( AliMUONTrackParam* ) const;
 
-  void ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t lWeight=1.0) const;
+  void LocalEquationX( void );
 
-  void ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight=1.0) const;
+  void LocalEquationY( void );
 
-  void FillDetElemData();
+  TGeoCombiTrans DeltaTransform(const double *detElemMisAlignment) const;
 
-  void FillRecPointData();
+  ///@name utilities
+  //@{
 
-  void FillTrackParamData();
+  void AddConstraint(Double_t* parameters, Double_t value );
 
-  void ResetLocalEquation();
+  Int_t GetChamberId( Int_t iDetElemNumber ) const;
 
-  /// local equation along X
-  void LocalEquationX( Bool_t doAlignment = kTRUE );
+  Bool_t DetElemIsValid( Int_t iDetElemId ) const;
 
-  /// local equation along Y
-  void LocalEquationY( Bool_t doAlignment = kTRUE );
+  Int_t GetDetElemNumber( Int_t iDetElemId ) const;
 
-  /// local equation using cluster alignment record
-  void LocalEquation( const AliMUONAlignmentClusterRecord& clusterRecord);
+  TString GetParameterMaskString( UInt_t parameterMask ) const;
 
-  TGeoCombiTrans ReAlign(const TGeoCombiTrans& transform, const double *detElemMisAlignment) const;
+  TString GetSidesMaskString( UInt_t sidesMask ) const;
 
-  Bool_t fBFieldOn; ///< Flag for Magnetic filed On/Off
+  //@}
 
-  Bool_t fChOnOff[10];  ///< Flags for chamber On/Off
+  /// true when initialized
+  Bool_t fInitialized;
 
-  Bool_t fSpecLROnOff[2];  ///< Flags for left right On/Off
+  /// current run id
+  Int_t fRunNumber;
 
-  Bool_t fDoF[4];  ///< Flags degrees of freedom to align (x,y,phi)
+  /// Flag for Magnetic filed On/Off
+  Bool_t fBFieldOn;
 
-  Double_t fAllowVar[4];  ///< "Encouraged" variation for degrees of freedom
+  /// "Encouraged" variation for degrees of freedom
+  Double_t fAllowVar[fgNParCh];
 
+  /// Initial value for chi2 cut
   /** if > 1 Iterations in AliMillepede are turned on */
-  Double_t fStartFac;   ///< Initial value for chi2 cut
-
-  Double_t fResCutInitial;  ///< Cut on residual for first iteration
-
-  Double_t fResCut;  ///< Cut on residual for other iterations
-
-  AliMillepede *fMillepede;   ///< Detector independent alignment class
-
-  AliMUONTrack *fTrack;   //!< AliMUONTrack
-
-  AliMUONVCluster *fCluster;  //!< AliMUONVCluster
-
-  AliMUONTrackParam *fTrackParam;  //!< Track parameters
-
-  Int_t fNGlobal;  ///< Number of global parameters
-  Int_t fNLocal;   ///< Number of local parameters
-  Int_t fNStdDev;  ///< Number of standard deviations for chi2 cut
-  Double_t fClustPos[3];    ///< Cluster position
-  Double_t fClustPosLoc[3]; ///< Cluster position in local coordinates
-  Double_t fTrackSlope0[2]; ///< Track slope at reference point
-  Double_t fTrackSlope[2];  ///< Track slope at current point
-  Double_t fTrackPos0[3];   ///< Track intersection at reference point
-  Double_t fTrackPos[3];    ///< Track intersection at current point
-  Double_t fTrackPosLoc[3]; ///< Track intersection at current point in local coordinates
-  Double_t fMeas[2];        ///< Current measurement (depend on B field On/Off)
-  Double_t fSigma[2];       ///< Estimated resolution on measurement
-
-  Double_t fGlobalDerivatives[624]; ///< Array of global derivatives
-  Double_t fLocalDerivatives[4];    ///< Array of local derivatives
-
-  Double_t fConstraintX[624];   ///< Array for constraint equation all X
-  Double_t fConstraintY[624];   ///< Array for constraint equation all Y
-  Double_t fConstraintP[624];   ///< Array for constraint equation all P
-  Double_t fConstraintXT[624];  ///< Array for constraint equation X Top half
-  Double_t fConstraintYT[624];  ///< Array for constraint equation Y Top half
-  Double_t fConstraintPT[624];  ///< Array for constraint equation P Top half
-  Double_t fConstraintXZT[624];  ///< Array for constraint equation X vs Z Top half
-  Double_t fConstraintYZT[624];  ///< Array for constraint equation Y vs Z Top half
-  Double_t fConstraintPZT[624];  ///< Array for constraint equation P vs Z Top half
-  Double_t fConstraintXYT[624];  ///< Array for constraint equation X vs Y Top half
-  Double_t fConstraintYYT[624];  ///< Array for constraint equation Y vs Y Top half
-  Double_t fConstraintPYT[624];  ///< Array for constraint equation P vs Y Top half
-  Double_t fConstraintXB[624];  ///< Array for constraint equation X Bottom half
-  Double_t fConstraintYB[624];  ///< Array for constraint equation Y Bottom half
-  Double_t fConstraintPB[624];  ///< Array for constraint equation P Bottom half
-  Double_t fConstraintXZB[624];  ///< Array for constraint equation X vs Z Bottom half
-  Double_t fConstraintYZB[624];  ///< Array for constraint equation Y vs Z Bottom half
-  Double_t fConstraintPZB[624];  ///< Array for constraint equation P vs Z Bottom half
-  Double_t fConstraintXYB[624];  ///< Array for constraint equation X vs Y Bottom half
-  Double_t fConstraintYYB[624];  ///< Array for constraint equation Y vs Y Bottom half
-  Double_t fConstraintPYB[624];  ///< Array for constraint equation P vs Y Bottom half
-  Double_t fConstraintXR[624];  ///< Array for constraint equation X Right half
-  Double_t fConstraintYR[624];  ///< Array for constraint equation Y Right half
-  Double_t fConstraintPR[624];  ///< Array for constraint equation P Right half
-  Double_t fConstraintXZR[624];  ///< Array for constraint equation X vs Z Right half
-  Double_t fConstraintYZR[624];  ///< Array for constraint equation Y vs Z Right half
-  Double_t fConstraintPZR[624];  ///< Array for constraint equation P vs Z Right half
-  Double_t fConstraintXYR[624];  ///< Array for constraint equation X vs Y Right half
-  Double_t fConstraintYYR[624];  ///< Array for constraint equation Y vs Y Right half
-  Double_t fConstraintPYR[624];  ///< Array for constraint equation P vs Y Right half
-  Double_t fConstraintXL[624];  ///< Array for constraint equation X Left half
-  Double_t fConstraintYL[624];  ///< Array for constraint equation Y Left half
-  Double_t fConstraintPL[624];  ///< Array for constraint equation P Left half
-  Double_t fConstraintXZL[624];  ///< Array for constraint equation X vs Z Left half
-  Double_t fConstraintYZL[624];  ///< Array for constraint equation Y vs Z Left half
-  Double_t fConstraintPZL[624];  ///< Array for constraint equation P vs Z Left half
-  Double_t fConstraintXYL[624];  ///< Array for constraint equation X vs Y Left half
-  Double_t fConstraintYYL[624];  ///< Array for constraint equation Y vs Y Left half
-  Double_t fConstraintPYL[624];  ///< Array for constraint equation P vs Y Left half
-  Double_t fConstraintX3[624];  ///< Array for constraint equation St3 X
-  Double_t fConstraintY3[624];  ///< Array for constraint equation St3 Y
-  Double_t fConstraintX4[624];  ///< Array for constraint equation St4 X
-  Double_t fConstraintY4[624];  ///< Array for constraint equation St4 Y
-  Double_t fConstraintP4[624];  ///< Array for constraint equation St4 P
-  Double_t fConstraintX5[624];  ///< Array for constraint equation St5 X
-  Double_t fConstraintY5[624];  ///< Array for constraint equation St5 Y
-
-  Int_t fDetElemId;        ///< Detection element id
-  Int_t fDetElemNumber;    ///< Detection element number
-  Double_t fPhi;           ///< Azimuthal tilt of detection element
-  Double_t fCosPhi;        ///< Cosine of fPhi
-  Double_t fSinPhi;        ///< Sine of fPhi
-  Double_t fDetElemPos[3]; ///< Position of detection element
-
-  AliMUONAlignmentTrackRecord fTrackRecord;  //!< running Track record
-
-  AliMUONGeometryTransformer *fTransform;  ///< Geometry transformation
-
-  static Int_t fgNSt;            ///< Number tracking stations
-  static Int_t fgNCh;            ///< Number tracking chambers
-  static Int_t fgNTrkMod;        ///< Number of tracking modules (4 ch + 6*2 half-ch)
-  static Int_t fgNParCh;         ///< Number of degrees of freedom per chamber
-  static Int_t fgNDetElem;       ///< Total number of detection elements
-  static Int_t fgNDetElemCh[10]; ///< Number of detection elements per chamber
-  static Int_t fgSNDetElemCh[10];///< Sum of detection elements up to this chamber (inc)
+  Double_t fStartFac;
+
+  /// Cut on residual for first iteration
+  Double_t fResCutInitial;
+
+  /// Cut on residual for other iterations
+  Double_t fResCut;
+
+  /// Detector independent alignment class
+  AliMillePede2* fMillepede;
+
+  /// running AliMUONVCluster
+  AliMUONVCluster* fCluster;
+
+  /// Number of standard deviations for chi2 cut
+  Int_t fNStdDev;
+
+  /// Cluster (global) position
+  Double_t fClustPos[3];
+
+  /// Track slope at reference point
+  Double_t fTrackSlope0[2];
+
+  /// Track slope at current point
+  Double_t fTrackSlope[2];
+
+  /// Track intersection at reference point
+  Double_t fTrackPos0[3];
+
+  /// Track intersection at current point
+  Double_t fTrackPos[3];
+
+  /// Current measurement (depend on B field On/Off)
+  Double_t fMeas[2];
+
+  /// Estimated resolution on measurement
+  Double_t fSigma[2];
+
+  /// degrees of freedom
+  enum
+  {
+    kFixedParId = -1,
+    kFreeParId = kFixedParId-1,
+    kGroupBaseId = -10
+  };
+
+  /// Array of effective degrees of freedom
+  /// it is used to fix detectors, parameters, etc.
+  Int_t fGlobalParameterStatus[fNGlobal];
+
+  /// Array of global derivatives
+  Double_t fGlobalDerivatives[fNGlobal];
+
+  /// Array of local derivatives
+  Double_t fLocalDerivatives[fNLocal];
+
+  /// current detection element number
+  Int_t fDetElemNumber;
+
+  /// unbias
+  Bool_t fUnbias;
+
+  /// running Track record
+  AliMillePedeRecord fTrackRecord;
+
+  /// Geometry transformation
+  AliMUONGeometryTransformer *fTransform;
+  TGeoCombiTrans fGeoCombiTransInverse;
 
   ClassDef(AliMUONAlignment, 2)
 
index 80de748..370cdc8 100644 (file)
 /// \class AliMUONAlignmentTask
 /// AliAnalysisTask to align the MUON spectrometer.
 /// The Task reads as input ESDs and feeds the MUONTracks to AliMUONAlignment.
-/// The alignment itself is performed by AliMillepede.
+/// The alignment itself is performed by AliMillePede2.
 /// A OCDB entry is written with the alignment parameters.
-/// A list of graph are written to monitor the alignment parameters.
 ///
 /// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
+/// \author Hugo Pereira Da Costa, CEA/Saclay - Irfu/SPhN
 //-----------------------------------------------------------------------------
 
 // this class header must always come first
 
 // local header must come before system headers
 #include "AliAnalysisManager.h"
+#include "AliAlignObjMatrix.h"
 #include "AliAODInputHandler.h"
 #include "AliAODHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
+#include "AliCDBEntry.h"
 #include "AliESDInputHandler.h"
 #include "AliESDEvent.h"
 #include "AliESDMuonTrack.h"
-#include "AliMagF.h"
 #include "AliCDBStorage.h"
 #include "AliCDBManager.h"
 #include "AliGRPManager.h"
 #include "AliCDBId.h"
 #include "AliGeomManager.h"
 
+#include "AliMagF.h"
+#include "AliMillePedeRecord.h"
 #include "AliMpCDB.h"
 #include "AliMUONCDB.h"
 #include "AliMUONAlignment.h"
+#include "AliMUONConstants.h"
+#include "AliMUONRecoParam.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackExtrap.h"
 #include "AliMUONTrackParam.h"
 #include "AliMUONESDInterface.h"
 
 // system headers
+#include <cassert>
 #include <fstream>
+
+// root headers
 #include <TString.h>
 #include <TError.h>
+#include <TFile.h>
 #include <TGraphErrors.h>
 #include <TTree.h>
 #include <TChain.h>
@@ -72,28 +81,26 @@ ClassImp(AliMUONAlignmentTask)
 ///\endcond
 
 //________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *newalignocdb, const char *oldalignocdb, const char *defaultocdb, const char *geofilename):
+AliMUONAlignmentTask::AliMUONAlignmentTask( const char *name ):
   AliAnalysisTaskSE(name),
   fReadRecords( kFALSE ),
   fWriteRecords( kFALSE ),
   fDoAlignment( kTRUE ),
+  fMergeAlignmentCDBs( kTRUE ),
+  fForceBField( kFALSE ),
+  fBFieldOn( kFALSE ),
+  fUnbias( kFALSE ),
   fAlign(0x0),
-  fGeoFilename(geofilename),
-  fDefaultStorage(defaultocdb),
-  fOldAlignStorage(oldalignocdb),
-  fNewAlignStorage(newalignocdb),
-  fOldGeoTransformer(NULL),
-  fNewGeoTransformer(NULL),
+  fNewAlignStorage( "local://ReAlignOCDB" ),
+  fOldGeoTransformer(0x0),
+  fNewGeoTransformer(0x0),
   fLoadOCDBOnce(kFALSE),
   fOCDBLoaded(kFALSE),
+  fEvent(0),
   fTrackTot(0),
   fTrackOk(0),
-  fLastRunNumber(-1),
-  fMSDEx(0x0),
-  fMSDEy(0x0),
-  fMSDEz(0x0),
-  fMSDEp(0x0),
-  fList(0x0),
+  fRunNumberMin(0),
+  fRunNumberMax(0),
   fRecords(0x0),
   fRecordCount(0)
 {
@@ -102,119 +109,26 @@ AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *newalig
   // Input slot #0 works with a TChain
   DefineInput(0, TChain::Class());
 
-  // Output slot #0 writes NTuple/histos into a TList
-  DefineOutput(1, TList::Class());
-
   // initialize parameters ...
-  for(Int_t k=0;k<4*156;k++)
+  for(Int_t k=0;k<AliMUONAlignment::fNGlobal;k++)
   {
     fParameters[k]=0.;
     fErrors[k]=0.;
     fPulls[k]=0.;
   }
 
-  fOldGeoTransformer = new AliMUONGeometryTransformer();
-
-}
-
-//________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& other):
-  AliAnalysisTaskSE(other),
-  fReadRecords( other.fReadRecords ),
-  fWriteRecords( other.fWriteRecords ),
-  fDoAlignment( other.fDoAlignment ),
-  fAlign( other.fAlign ),
-  fGeoFilename( other.fGeoFilename ),
-  fDefaultStorage( other.fDefaultStorage ),
-  fOldAlignStorage( other.fOldAlignStorage ),
-  fNewAlignStorage( other.fNewAlignStorage ),
-  fOldGeoTransformer( other.fOldGeoTransformer ),
-  fNewGeoTransformer( other.fNewGeoTransformer ),
-  fLoadOCDBOnce( other.fLoadOCDBOnce ),
-  fOCDBLoaded( other.fOCDBLoaded ),
-  fTrackTot( other.fTrackTot ),
-  fTrackOk( other.fTrackOk ),
-  fLastRunNumber( other.fLastRunNumber ),
-  fMSDEx( other.fMSDEx ),
-  fMSDEy( other.fMSDEy ),
-  fMSDEz( other.fMSDEz ),
-  fMSDEp( other.fMSDEp ),
-  fList( other.fList ),
-  fRecords( other.fRecords ),
-  fRecordCount( other.fRecordCount )
-{
-
-  // initialize parameters
-  for(Int_t k=0;k<4*156;k++)
-  {
-
-    fParameters[k]=other.fParameters[k];
-    fErrors[k]=other.fErrors[k];
-    fPulls[k]=other.fPulls[k];
-
-  }
-
-}
-
-//________________________________________________________________________
-AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other)
-{
-  /// Assignment
-  if(&other == this) return *this;
-  AliAnalysisTaskSE::operator=(other);
-
-  fReadRecords = other.fReadRecords;
-  fWriteRecords = other.fWriteRecords;
-  fDoAlignment = other.fDoAlignment;
-
-  // this breaks in destructor
-  fAlign = other.fAlign;
-
-  fGeoFilename = other.fGeoFilename;
-  fDefaultStorage = other.fDefaultStorage;
-  fOldAlignStorage = other.fOldAlignStorage;
-  fNewAlignStorage = other.fNewAlignStorage;
-
-  // this breaks in destructor
-  fOldGeoTransformer = other.fOldGeoTransformer;
-  fNewGeoTransformer = other.fNewGeoTransformer;
-
-  fLoadOCDBOnce = other.fLoadOCDBOnce;
-  fOCDBLoaded = other.fOCDBLoaded;
-  fTrackTot = other.fTrackTot;
-  fTrackOk = other.fTrackOk;
-  fLastRunNumber = other.fLastRunNumber;
-  fMSDEx = other.fMSDEx;
-  fMSDEy = other.fMSDEy;
-  fMSDEz = other.fMSDEz;
-  fMSDEp = other.fMSDEp;
-  fList = other.fList;
-  fRecords = other.fRecords;
-  fRecordCount = other.fRecordCount;
-
-  // initialize parameters ...
-  for( Int_t k=0; k<4*156; ++k)
-  {
-
-    fParameters[k]=other.fParameters[k];
-    fErrors[k]=other.fErrors[k];
-    fPulls[k]=other.fPulls[k];
+  // create alignment object
+  fAlign = new AliMUONAlignment();
 
-  }
+  // create old geometry transformer
+  fOldGeoTransformer = new AliMUONGeometryTransformer();
 
-  return *this;
 }
 
 //________________________________________________________________________
 AliMUONAlignmentTask::~AliMUONAlignmentTask()
 {
-  /*
-  it is safe to delete NULL pointers, so
-  there is no need to test their validity here.
-  However, it crashes here, because of incorrect assignment operator
-  and copy constructor, resulting in double deletion of objects.
-  Would require deep copy instead.
-  */
+  /// destructor
   delete fAlign;
   delete fOldGeoTransformer;
   delete fNewGeoTransformer;
@@ -223,14 +137,53 @@ AliMUONAlignmentTask::~AliMUONAlignmentTask()
 //________________________________________________________________________
 void AliMUONAlignmentTask::LocalInit()
 {
-  /// Local initialization, called once per task on the client machine
-  /// where the analysis train is assembled
+  /**
+  Local initialization, called once per task on the client machine
+  where the analysis train is assembled
+  **/
+
+  /* must run alignment when reading records */
+  if( fReadRecords && !fDoAlignment )
+  {
+
+    AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
+    SetDoAlignment( kTRUE );
+
+  }
 
   // print configuration
+  if( fRunNumberMin > 0 || fRunNumberMax > 0 )
+  { AliInfo( Form( "run range: %i - %i", fRunNumberMin, fRunNumberMax ) ); }
+
   AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
   AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
   AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
 
+  if( fDoAlignment )
+  {
+    // merge alignment DB flag is irrelevant if not actually performing the alignemnt
+    AliInfo( Form( "fMergeAlignmentCDBs: %s", (fMergeAlignmentCDBs ? "kTRUE":"kFALSE" ) ) );
+  }
+
+  // storage elements
+  if( !fDefaultStorage.IsNull() ) AliInfo( Form( "fDefaultStorage: %s", fDefaultStorage.Data() ) );
+  if( !fOldAlignStorage.IsNull() ) AliInfo( Form( "fOldAlignStorage: %s", fOldAlignStorage.Data() ) );
+
+  if( fDoAlignment )
+  {
+    // new alignment storage is irrelevant if not actually performing the alignemnt
+    if( !fNewAlignStorage.IsNull() ) AliInfo( Form( "fNewAlignStorage: %s", fNewAlignStorage.Data() ) );
+    else AliFatal( "Invalid new alignment storage path" );
+  }
+
+  if( !fReadRecords )
+  {
+    // following flags are only relevant if not reading records
+    if( fForceBField ) AliInfo( Form( "fBFieldOn: %s", (fBFieldOn ? "kTRUE":"kFALSE" ) ) );
+    else AliInfo( "fBFieldOn: from GRP" );
+    AliInfo( Form( "fUnbias: %s", (fUnbias ? "kTRUE":"kFALSE" ) ) );
+  }
+
   // consistency checks between flags
   /* need at least one of the flags set to true */
   if( !( fReadRecords || fWriteRecords || fDoAlignment ) )
@@ -240,63 +193,37 @@ void AliMUONAlignmentTask::LocalInit()
   if( fReadRecords && fWriteRecords )
   { AliFatal( "Cannot set both fReadRecords and fWriteRecords to kTRUE at the same time" ); }
 
-  /* must run alignment when reading records */
-  if( fReadRecords && !fDoAlignment )
+  /*
+  fix detectors
+  warning, counting chambers from 1.
+  this must be done before calling the Init method
+  */
+  if( fDoAlignment )
   {
 
-    AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
-    SetDoAlignment( kTRUE );
+    fAlign->FixChamber(1);
+    fAlign->FixChamber(10);
 
-  }
+  } else {
 
-  // Set initial values here, good guess may help convergence
-  // St 1
-  //  Int_t iPar = 0;
-  //  fParameters[iPar++] =  0.010300 ;  fParameters[iPar++] =  0.010600 ;  fParameters[iPar++] =  0.000396 ;
+    AliInfo( "Not fixing detector elements, since alignment is not required" );
 
-  fAlign = new AliMUONAlignment();
-  fAlign->InitGlobalParameters(fParameters);
-
-//  AliCDBManager::Instance()->Print();
-//
-//  fAlign->SetGeometryTransformer(fOldGeoTransformer);
+  }
 
-  // Do alignment with magnetic field off
-  fAlign->SetBFieldOn(kFALSE);
+  // initialize
+  fAlign->Init();
 
-  // Set tracking station to use
-  //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-  Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
+  // use unbiased residuals
+  fAlign->SetUnbias( fUnbias );
 
-  // Set degrees of freedom to align (see AliMUONAlignment)
-  fAlign->AllowVariations(bChOnOff);
+  // Do alignment with magnetic field off
+  fAlign->SetBFieldOn( fBFieldOn );
 
   // Set expected resolution (see AliMUONAlignment)
   fAlign->SetSigmaXY(0.15,0.10);
 
-  // Fix parameters or add constraints here
-  //   for (Int_t iSt=0; iSt<5; iSt++)
-  //     if (!bStOnOff[iSt]) fAlign->FixStation(iSt+1);
-  for (Int_t iCh=0; iCh<10; iCh++)
-    if (!bChOnOff[iCh]) fAlign->FixChamber(iCh+1);
-
-  // Left and right sides of the detector are independent, one can choose to align
-  // only one side
-  Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
-  fAlign->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
-
-  fAlign->SetChOnOff(bChOnOff);
-  fAlign->SetSpecLROnOff(bChOnOff);
-
-    // Here we can fix some detection elements
-  //  fAlign->FixDetElem(908);
-       //  fAlign->FixDetElem(1012);
-       fAlign->FixDetElem(608);
-
-  // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
-//   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-//   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
-  //  fAlign->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
+  // initialize global parameters to provided values
+  fAlign->InitGlobalParameters(fParameters);
 
 }
 
@@ -304,26 +231,6 @@ void AliMUONAlignmentTask::LocalInit()
 void AliMUONAlignmentTask::UserCreateOutputObjects()
 {
 
-  if( fDoAlignment )
-  {
-
-    // Creating graphs
-    fMSDEx = new TGraphErrors(156);
-    fMSDEy = new TGraphErrors(156);
-    fMSDEz = new TGraphErrors(156);
-    fMSDEp = new TGraphErrors(156);
-
-    // Add Ntuples to the list
-    fList = new TList();
-    fList->Add(fMSDEx);
-    fList->Add(fMSDEy);
-    fList->Add(fMSDEz);
-    fList->Add(fMSDEp);
-
-    fList->SetOwner(kTRUE);
-    PostData(1, fList);
-  }
-
   // connect AOD output
   if( fWriteRecords )
   {
@@ -332,13 +239,13 @@ void AliMUONAlignmentTask::UserCreateOutputObjects()
     {
 
       // get AOD output handler and add Branch
-      fRecords = new TClonesArray( "AliMUONAlignmentTrackRecord", 0 );
+      fRecords = new TClonesArray( "AliMillePedeRecord", 0 );
       fRecords->SetName( "records" );
       handler->AddBranch("TClonesArray", &fRecords);
 
       fRecordCount = 0;
 
-    } else AliInfo( "Error: invalid output event handler" );
+    } else AliFatal( "Error: invalid output event handler" );
 
   }
 
@@ -347,10 +254,14 @@ void AliMUONAlignmentTask::UserCreateOutputObjects()
 //________________________________________________________________________
 void AliMUONAlignmentTask::UserExec(Option_t *)
 {
-  // print configuration
-//  AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
-//  AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
-//  AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
+
+  // do nothing if run number not in range
+  if( fRunNumberMin > 0 && fCurrentRunNumber < fRunNumberMin ) return;
+  if( fRunNumberMax > 0 && fCurrentRunNumber > fRunNumberMax ) return;
+
+  // increase event count
+  ++fEvent;
+
   // clear array
   if( fWriteRecords && fRecords )
   {
@@ -358,10 +269,8 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
     fRecordCount = 0;
   }
 
-  // local track parameters
-  Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
-
-  if( fReadRecords ) {
+  if( fReadRecords )
+  {
 
     AliAODEvent* lAOD( dynamic_cast<AliAODEvent*>(InputEvent() ) );
 
@@ -388,12 +297,14 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
     for( Int_t index = 0; index < lRecordsCount; ++index )
     {
 
-      if( AliMUONAlignmentTrackRecord* record = dynamic_cast<AliMUONAlignmentTrackRecord*>( records->UncheckedAt(index) ) )
+      if( AliMillePedeRecord* record = dynamic_cast<AliMillePedeRecord*>( records->UncheckedAt(index) ) )
       {
 
-        fAlign->ProcessTrack( record, fDoAlignment );
-        if( fDoAlignment )
-        { fAlign->LocalFit( fTrackOk++, trackParams, 0 ); }
+        fAlign->ProcessTrack( record );
+        ++fTrackOk;
+
+        if(!(fTrackTot%100))
+        { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
 
       } else AliInfo( Form( "Invalid record at %i", index ) );
 
@@ -411,20 +322,7 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
       return;
     }
 
-// HUGO: Comment out check on run number, to be able to run on MC
-//     if (!lESD->GetRunNumber())
-//     {
-//       AliInfo( Form( "Current Run Number: %i", lESD->GetRunNumber() ) );
-//       return;
-//     }
-
-    //  if (lESD->GetRunNumber()!=fLastRunNumber){
-    //    fLastRunNumber = lESD->GetRunNumber();
-    //    Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
-    //  }
-
     Int_t nTracks = Int_t(lESD->GetNumberOfMuonTracks());
-//    cout << " there are " << nTracks << " tracks" << endl;
     for( Int_t iTrack = 0; iTrack < nTracks; iTrack++ )
     {
 
@@ -433,8 +331,6 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
       if (!esdTrack->ContainTriggerData()) continue;
 
       Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
-      //     fInvBenMom->Fill(invBenMom);
-      //     fBenMom->Fill(1./invBenMom);
       if (TMath::Abs(invBenMom)<=1.04)
       {
 
@@ -442,46 +338,41 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
         AliMUONESDInterface::ESDToMUON(*esdTrack, track);
 
         // process track and retrieve corresponding records, for storage
-        const AliMUONAlignmentTrackRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
-
-        // do the fit, if required
-        if( fDoAlignment ) fAlign->LocalFit(fTrackOk++,trackParams,0);
-        else fTrackOk++;
+        const AliMillePedeRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
+        ++fTrackOk;
 
         // store in array
         if( fWriteRecords && fRecords )
-        {                                      
-          new((*fRecords)[fRecordCount]) AliMUONAlignmentTrackRecord( *lRecords );
+        {
+          new((*fRecords)[fRecordCount]) AliMillePedeRecord( *lRecords );
           ++fRecordCount;
         }
 
-      }                        
-                       
+      }
+
       ++fTrackTot;
 
-      if(!(fTrackTot%1000))
+      if(!(fTrackTot%100))
       { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
 
-//      cout << "Processed " << fTrackTot << " Tracks." << endl;
       // Post final data. Write histo list to a file with option "RECREATE"
-      PostData(1,fList);
+      // PostData(1,fList);
 
     }
 
     // save AOD
-    if( fWriteRecords && fRecordCount > 0 ) { 
+    if( fWriteRecords && fRecordCount > 0 )
+    {
       AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
       if( handler )
-       {
-//                     printf("handler: %p\n",handler);
-                       AliAODEvent* aod = handler->GetAOD();
-//                     printf("aod: %p\n",aod);
-                       AliAODHeader* header = aod->GetHeader();
-//                     printf("header: %p\n",header);
-                       header->SetRunNumber(lESD->GetRunNumber());
-//                     printf("RunNumber: %d\n",lESD->GetRunNumber());
-                       AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
-       } else AliInfo( "Error: invalid output event handler" );
+      {
+        AliAODEvent* aod = handler->GetAOD();
+        AliAODHeader* header = aod->GetHeader();
+        header->SetRunNumber(lESD->GetRunNumber());
+        AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
+
+      } else AliInfo( "Error: invalid output event handler" );
+
     }
 
   }
@@ -489,10 +380,6 @@ void AliMUONAlignmentTask::UserExec(Option_t *)
 }
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::Terminate(const Option_t*)
-{ return; }
-
-//________________________________________________________________________
 void AliMUONAlignmentTask::FinishTaskOutput()
 {
 
@@ -507,164 +394,283 @@ void AliMUONAlignmentTask::FinishTaskOutput()
 
   // Perform global fit
   fAlign->GlobalFit(fParameters,fErrors,fPulls);
-
-//   // Update pointers reading them from the output slot
-//   fList = (TList*)GetOutputData(0);
-//   fMSDEx = (TGraphErrors*)fList->At(0);
-//   fMSDEy = (TGraphErrors*)fList->At(1);
-//   fMSDEz = (TGraphErrors*)fList->At(2);
-//   fMSDEp = (TGraphErrors*)fList->At(3);
-
-  // Store results
-  Double_t DEid[156] = {0};
-  Double_t MSDEx[156] = {0};
-  Double_t MSDEy[156] = {0};
-  Double_t MSDEz[156] = {0};
-  Double_t MSDEp[156] = {0};
-  Double_t DEidErr[156] = {0};
-  Double_t MSDExErr[156] = {0};
-  Double_t MSDEyErr[156] = {0};
-  Double_t MSDEzErr[156] = {0};
-  Double_t MSDEpErr[156] = {0};
-  Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
-  Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
-  // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
   Int_t idOffset = 0; // 400
   Int_t lSDetElemCh = 0;
-  for(Int_t iDE=0; iDE<lNDetElem; iDE++)
+
+  for( Int_t iDE=0; iDE<AliMUONAlignment::fgNDetElem; iDE++ )
   {
 
-    DEidErr[iDE] = 0.;
-    DEid[iDE] = idOffset+100;
-    DEid[iDE] += iDE;
+    Int_t detElemId = idOffset+100;
+    detElemId += iDE;
     lSDetElemCh = 0;
     for(Int_t iCh=0; iCh<9; iCh++)
     {
-      lSDetElemCh += lNDetElemCh[iCh];
 
+      lSDetElemCh += AliMUONAlignment::fgNDetElemCh[iCh];
       if (iDE>=lSDetElemCh)
       {
-        DEid[iDE] += 100;
-        DEid[iDE] -= lNDetElemCh[iCh];
+        detElemId += 100;
+        detElemId -= AliMUONAlignment::fgNDetElemCh[iCh];
       }
 
     }
 
-    MSDEx[iDE]=fParameters[4*iDE+0];
-    MSDEy[iDE]=fParameters[4*iDE+1];
-    MSDEz[iDE]=fParameters[4*iDE+3];
-    MSDEp[iDE]=fParameters[4*iDE+2];
-    MSDExErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+0);
-    MSDEyErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+1);
-    MSDEzErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+3);
-    MSDEpErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+2);
-    fMSDEx->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+0]);
-    fMSDEx->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+0));
-    fMSDEy->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+1]);
-    fMSDEy->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+1));
-    fMSDEz->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+3]);
-    fMSDEz->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+3));
-    fMSDEp->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+2]);
-    fMSDEp->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+2));
-
   }
 
   // Post final data. Write histo list to a file with option "RECREATE"
-  PostData(1,fList);
+  // PostData(1,fList);
 
-  // HUGO: stop here to test reproducibility
-       //  return;
+  // store misalignments from OCDB into old transformers
+  if( fMergeAlignmentCDBs )
+  { SaveMisAlignmentData( fOldGeoTransformer ); }
 
   // Re Align
-  fNewGeoTransformer = fAlign->ReAlign(fOldGeoTransformer,fParameters,true);
-  //  newTransform->WriteTransformations("transform2ReAlign.dat");
+  fNewGeoTransformer = fAlign->ReAlign( fOldGeoTransformer, fParameters, true );
 
   // Generate realigned data in local cdb
   const TClonesArray* array = fNewGeoTransformer->GetMisAlignmentData();
 
   // 100 mum residual resolution for chamber misalignments?
-  fAlign->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
+  fAlign->SetAlignmentResolution( array, -1, 0.01, 0.01, 0.004, 0.003 );
 
   // CDB manager
   AliLog::SetGlobalDebugLevel(2);
-  AliCDBManager* cdbm = AliCDBManager::Instance();
-  // cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
 
-       // recover default storage full name (raw:// cannot be used to set specific storage)
-       TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
-       if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
-       else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+  // recover default storage full name (raw:// cannot be used to set specific storage)
+  AliCDBManager* cdbManager = AliCDBManager::Instance();
+
+  // unload old alignment path
+  if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
+  { cdbManager->UnloadFromCache("MUON/Align/Data"); }
 
-       if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
-       if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
-       //      else cdbm->SetSpecificStorage("MUON/Align/Data",fDefaultStorage.Data());
-       else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+  // load new alignment path
+  if( !fNewAlignStorage.IsNull() ) cdbManager->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
+  else {
 
+    TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
+    if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+    else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+    cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+
+  }
 
+  // create new DB entry
   AliCDBMetaData* cdbData = new AliCDBMetaData();
   cdbData->SetResponsible("Dimuon Offline project");
   cdbData->SetComment("MUON alignment objects with residual misalignment");
   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
-  cdbm->Put(const_cast<TClonesArray*>(array), id, cdbData);
-
-       gSystem->Exec("cp -a ReAlignOCDB/MUON/Align/Data/Run0_999999999_v0_s0.root Run0_999999999_v0_s0.root");
-       gSystem->Exec("ls -l");
+  cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 
 }
 
 //________________________________________________________________________
 void AliMUONAlignmentTask::NotifyRun()
 {
-  
-  if (fOCDBLoaded && fLoadOCDBOnce) { AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
+
+  /// run number (re)initialization
+
+  // propagate run number to fAlign
+  if( fAlign ) fAlign->SetRunNumber( fCurrentRunNumber );
+
+  // update ocdb
+  if (fOCDBLoaded && fLoadOCDBOnce)
+  {
+    AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
     return;
   }
 
-  AliCDBManager* cdbm = AliCDBManager::Instance();
-  cdbm->SetDefaultStorage(fDefaultStorage.Data());
-//     cdbm->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  AliCDBManager* cdbManager = AliCDBManager::Instance();
+
+  // Initialize default storage
+  if( !( cdbManager->IsDefaultStorageSet() || fDefaultStorage.IsNull() ) )
+  {
+
+    AliInfo( Form( "Initializing default storage: %s", fDefaultStorage.Data() ) );
+    cdbManager->SetDefaultStorage(fDefaultStorage.Data());
+
+  } else if( !fDefaultStorage.IsNull() ) {
+
+    AliInfo( "Default storage already set. Ignoring fDefaultStorage" );
+
+  } else {
+
+    AliInfo( "Default storage already set" );
+
+  }
+
+  // Initialize run number
+  if( cdbManager->GetRun() < 0 )
+  {
+    AliInfo( Form( "Setting run number: %i", fCurrentRunNumber ) );
+    cdbManager->SetRun(fCurrentRunNumber);
+  }
+
+  // following initialization is not needed when reading records
+  if( !fReadRecords )
+  {
+
+    // load magnetic field if needed
+    if( !TGeoGlobalMagField::Instance()->IsLocked() )
+    {
 
-  // HUGO: add to comment out the specific settings below in order to be able to run.
-  //cdbm->SetSpecificStorage("GRP/GRP/Data",fDefaultStorage.Data());
-       //cdbm->SetSpecificStorage("GRP/Geometry/data",fDefaultStorage.Data());
-       //cdbm->SetSpecificStorage("MUON/Calib/MappingData",fDefaultStorage.Data());
-       cdbm->SetRun(fCurrentRunNumber);
+      AliInfo( "Loading magnetic field" );
+      if( !AliMUONCDB::LoadField() )
+      {
+        AliError( "Failed to load magnetic field" );
+        return;
+      }
 
-  if (!AliMUONCDB::LoadField()) { AliError("Problem field"); return;}
+    } else { AliInfo( "Magnetic field already locked" ); }
 
-  // set the magnetic field for track extrapolations
-  AliMUONTrackExtrap::SetField();
+    // checking magnitic field
+    if( !fForceBField )
+    {
+      AliInfo( "Reading magnetic field setup from GRP" );
+
+      // decide bFieldOn value base on dipole current
+      // propagete to Alignment class
+      // and printout
+      AliMagF* magF = dynamic_cast<AliMagF*>( TGeoGlobalMagField::Instance()->GetField() );
+      fBFieldOn = TMath::Abs( magF->GetFactorDip() ) > 1e-5;
+      fAlign->SetBFieldOn( fBFieldOn );
+      AliInfo( Form( "Dipole magnetic field factor: %.2f", magF->GetFactorDip() ) );
+      AliInfo( Form( "fBFieldOn = %s", (fBFieldOn ? "kTRUE":"kFALSE") ) );
+    }
 
-  if (!AliMUONCDB::LoadMapping(kTRUE)) { AliError("Problem mapping"); return;}
+    AliInfo( "Loading muon mapping" );
+    if( !AliMUONCDB::LoadMapping(kFALSE) )
+    {
+      AliError( "Failed to load muon mapping" );
+      return;
+    }
 
-       // recover default storage full name (raw:// cannot be used to set specific storage)
-       TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
-       if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
-       else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+    AliInfo( "Assigning field to Track extrapolator" );
+    AliMUONTrackExtrap::SetField();
 
-       // reset existing geometry/alignment if any
-       if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
-       if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
-       if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
+  }
 
-       // get original geometry transformer
-       AliGeomManager::LoadGeometry();
-       if (!AliGeomManager::GetGeometry()) { AliError("Problem geometry"); return;}
-       if (fOldAlignStorage != "none")
+  // load geometry if needed
+  if( !AliGeomManager::GetGeometry() )
   {
 
-               if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
-               else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
-               //              else cdbm->SetSpecificStorage("MUON/Align/Data",fDefaultStorage.Data());
-               
+    // reset existing geometry/alignment if any
+    if( cdbManager->GetEntryCache()->Contains("GRP/Geometry/Data") )
+    {
+      AliInfo( "Unloading GRP/Geometry/Data" );
+      cdbManager->UnloadFromCache("GRP/Geometry/Data");
+    }
+
+    if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
+    {
+      AliInfo( Form( "Unloading MUON/Align/Data from %s", cdbManager->GetSpecificStorage( "MUON/Align/Data" )->GetBaseFolder().Data() ) );
+      cdbManager->UnloadFromCache("MUON/Align/Data");
+    }
+
+    // get original geometry transformer
+    AliInfo( "Loading geometry" );
+    AliGeomManager::LoadGeometry();
+    if (!AliGeomManager::GetGeometry())
+    {
+      AliError("Failed to load geometry");
+      return;
+    }
+
+    if (!fOldAlignStorage.IsNull())
+    {
+
+      AliInfo( Form( "Initializing MUON/Align/Data using: %s", fOldAlignStorage.Data() ) );
+      cdbManager->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
+
+    } else {
+
+      // recover default storage full name (raw:// cannot be used to set specific storage)
+      TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
+      if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+      else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+
+      AliInfo( Form( "Re-initializing MUON/Align/Data using: %s", defaultStorage.Data() ) );
+      cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+
+    }
+
+    AliInfo("Loading muon Alignment objects");
     AliGeomManager::ApplyAlignObjsFromCDB("MUON");
 
+  } else { AliInfo( "Geometry already initialized - fOldAlignStorage ignored" ); }
+
+  // update geometry transformer and pass to Alignment object
+  AliInfo("Loading muon geometry in transformer");
+  fOldGeoTransformer->LoadGeometryData();
+  fAlign->SetGeometryTransformer(fOldGeoTransformer);
+
+  fOCDBLoaded = kTRUE;
+
+}
+
+//_____________________________________________________________________________________
+void AliMUONAlignmentTask::SaveMisAlignmentData( AliMUONGeometryTransformer* transformer ) const
+{
+
+  // clear transformer
+  transformer->ClearMisAlignmentData();
+
+  // load MUON/Align/Data from OCDB
+  AliCDBManager* cdbManager = AliCDBManager::Instance();
+  AliCDBEntry* cdbEntry = cdbManager->Get( "MUON/Align/Data" );
+  if (!cdbEntry)
+  {
+
+    AliError( "unable to load entry for path MUON/Align/Data" );
+    return;
+
+  }
+
+  // get TClonesArray and check
+  TClonesArray* misArray = (TClonesArray*)cdbEntry->GetObject();
+  if( !misArray )
+  {
+
+    AliError( "unable to load old misalignment array" );
+    return;
+
   }
 
-  // fOldGeoTransformer = new AliMUONGeometryTransformer();
-       fOldGeoTransformer->LoadGeometryData();
-       fAlign->SetGeometryTransformer(fOldGeoTransformer);
+  // loop over stored entries
+  for (Int_t index=0; index<misArray->GetEntriesFast(); ++index )
+  {
+
+    // load matrix and check
+    AliAlignObjMatrix* matrix = (AliAlignObjMatrix*) misArray->At( index );
+    if( !matrix )
+    {
+      AliError( Form( "unable to load matrix for index %i", index ) );
+      continue;
+    }
+
+    // get volume ID
+    AliGeomManager::ELayerID layerId;
+    Int_t moduleId;
+    matrix->GetVolUID( layerId, moduleId);
+
+    // make sure ELayerID is correct
+    assert( layerId == AliGeomManager::kMUON );
 
-       fOCDBLoaded = kTRUE;
+    // printout
+    // AliInfo( Form( "Found matrix for %s %i", matrix->GetSymName(), moduleId ) );
 
+    // get matrix
+    TGeoHMatrix misMatrix;
+    matrix->GetMatrix(misMatrix);
+
+    // add to geometry transformer
+    // need the detector element
+    // "detElement->GetId()"
+    // see fOldGeoTransformer->GetMisAlignmentData( ... )
+
+    if( TString( matrix->GetSymName() ).Contains( "DE" ) ) transformer->AddMisAlignDetElement( moduleId,  misMatrix);
+    else transformer->AddMisAlignModule( moduleId,  misMatrix);
+  }
+
+  return;
 }
+
index c921a76..abf94f1 100644 (file)
 ///
 //  Author Javier Castillo, CEA/Saclay - Irfu/SPhN
 
-class TList;
+#include "AliAnalysisTaskSE.h"
+#include "AliMUONAlignment.h"
+
 class TString;
-class TGraphErrors;
 class TClonesArray;
-class AliMUONAlignment;
 class AliMUONGeometryTransformer;
 
-#include "AliAnalysisTaskSE.h"
-
 class AliMUONAlignmentTask : public AliAnalysisTaskSE
 {
 
   public:
 
-
   /// constructor
-  AliMUONAlignmentTask(const char *name = "AliMUONAlignmentTask", const char *newalignocdb = "local://ReAlignOCDB", const char *oldalignocdb = "none", const char *defaultocdb = "raw://", const char *geofilename = "geometry.root");
-
-  /// copy constructor
-  AliMUONAlignmentTask(const AliMUONAlignmentTask& obj);
-
-  /// asignment operator
-  AliMUONAlignmentTask& operator=(const AliMUONAlignmentTask& other);
+  AliMUONAlignmentTask( const char *name = "AliMUONAlignmentTask" );
 
   /// destructor
   virtual ~AliMUONAlignmentTask();
 
-  virtual void   LocalInit();
+  /// get pointer to alignment class
+  AliMUONAlignment* alignment( void ) const
+  { return fAlign; }
+
+  virtual void LocalInit();
 
   ///@name flags
   //@{
@@ -55,43 +50,80 @@ class AliMUONAlignmentTask : public AliAnalysisTaskSE
   void SetDoAlignment( Bool_t value )
   { fDoAlignment = value; }
 
+  /// merge old and new Align CDBs into the new one.
+  void SetMergeAlignmentCDBs( Bool_t value )
+  { fMergeAlignmentCDBs = value; }
+
+  /// field on alignment
+  void SetBFieldOn( Bool_t value )
+  {
+    fForceBField = kTRUE;
+    fBFieldOn = value;
+  }
+
+  /// run range
+  void SetRunRange( Int_t runNumberMin, Int_t runNumberMax )
+  {
+    fRunNumberMin = runNumberMin;
+    fRunNumberMax = runNumberMax;
+  }
+
+  /// use unbiased residuals
+  void SetUnbias(Bool_t value )
+  {
+    fUnbias = value;
+    if( fAlign ) fAlign->SetUnbias( value );
+  }
+
+  /// use unbiased residuals
+  Bool_t GetUnbias() const
+  { return fUnbias; }
+
   //@}
 
   /// output data
   virtual void UserCreateOutputObjects();
 
   /// per-event method
-  virtual void UserExec(Option_t *option);
+  virtual void UserExec( Option_t* );
   virtual void NotifyRun();
 
   /// termination cleanup
-  virtual void Terminate(const Option_t*);
+  virtual void Terminate( const Option_t* )
+  {}
 
   /// end of task execution
   virtual void FinishTaskOutput();
 
-  /// Set geoemetry file name
-  void SetGeoFilename(const char* geoFilename)
-  {fGeoFilename = geoFilename;}
-
   /// Set default ocdb
-  void SetDefaultStorage(const char* defaultOCDB)
-  {fDefaultStorage = defaultOCDB;}
+  void SetDefaultStorage( TString defaultOCDB )
+  { fDefaultStorage = defaultOCDB; }
 
-  /// Set mis align ocdb
-  void SetOldAlignStorage(const char* oldalignOCDB)
-  {fOldAlignStorage = oldalignOCDB;}
+  /// Set old (misaligned) alignment path for ocdb
+  void SetOldAlignStorage( TString oldalignOCDB )
+  { fOldAlignStorage = oldalignOCDB; }
 
-  /// Set mis align ocdb
-  void SetNewAlignStorage(const char* newalignOCDB)
-  {fNewAlignStorage = newalignOCDB;}
+  /// Set new (realigned) alignment path for ocdb
+  void SetNewAlignStorage( TString newalignOCDB )
+  { fNewAlignStorage = newalignOCDB; }
 
-  /// Flag to set OCDB once at first notify
-  void SetLoadOCDBOnce(Bool_t loadOCDBOnce = kTRUE)
-  {fLoadOCDBOnce = loadOCDBOnce;}
+  /// Flag to set OCDB once at first run notify
+  void SetLoadOCDBOnce( Bool_t loadOCDBOnce = kTRUE )
+  { fLoadOCDBOnce = loadOCDBOnce; }
+
+  protected:
+
+  /// store misalignment matrices from OCDB into geometry transformer
+  void SaveMisAlignmentData( AliMUONGeometryTransformer* ) const;
 
   private:
 
+  /// copy constructor, not implemented
+  AliMUONAlignmentTask(const AliMUONAlignmentTask& obj);
+
+  /// asignment operator, not implemented
+  AliMUONAlignmentTask& operator=(const AliMUONAlignmentTask& other);
+
   ///@name flags
   //@{
 
@@ -104,14 +136,23 @@ class AliMUONAlignmentTask : public AliAnalysisTaskSE
   /// perform alignment (from either tracks or records depending on fReadRecords)
   Bool_t fDoAlignment;
 
+  /// merge old and new Align CDBs into the new one.
+  Bool_t fMergeAlignmentCDBs;
+
+  /// true if magnetic field was forced to value, instead of reading from GRP
+  Bool_t fForceBField;
+
+  /// Flag for Magnetic field On/Off
+  Bool_t fBFieldOn;
+
+  //! use unbiased residuals
+  Bool_t fUnbias;
+
   //@}
 
   /// The MUON alignment object
   AliMUONAlignment *fAlign;
 
-  /// Geometry file name
-  TString fGeoFilename;
-
   /// location of the default OCDB storage
   TString fDefaultStorage;
 
@@ -127,41 +168,33 @@ class AliMUONAlignmentTask : public AliAnalysisTaskSE
   /// new geometry transformer containing the new alignment to be applied
   AliMUONGeometryTransformer* fNewGeoTransformer;
 
+  /// set to true if not willing to re-initialize OCDB at every new run
   Bool_t fLoadOCDBOnce;
+
+  /// set to true when OCDB was loaded at least once
   Bool_t fOCDBLoaded;
 
+  //! event number (internal counter)
+  Int_t fEvent;
+
   /// Total number of track read
   Int_t fTrackTot;
 
   /// Number of tracks used for alignment
   Int_t fTrackOk;
 
-  /// Last run number
-  Int_t fLastRunNumber;
+  /// run range
+  Int_t fRunNumberMin;
+  Int_t fRunNumberMax;
 
   /// Array of alignment parameters
-  Double_t fParameters[4*156];
+  Double_t fParameters[AliMUONAlignment::fNGlobal];
 
   /// Array of alignment parameters errors
-  Double_t fErrors[4*156];
+  Double_t fErrors[AliMUONAlignment::fNGlobal];
 
   /// Array of alignment parameters pulls
-  Double_t fPulls[4*156];
-
-  /// Graph of translations along x
-  TGraphErrors *fMSDEx ;
-
-  /// Graph of translations along y
-  TGraphErrors *fMSDEy ;
-
-  /// Graph of translations along z
-  TGraphErrors *fMSDEz ;
-
-  /// Graph of rotation about z
-  TGraphErrors *fMSDEp;
-
-  /// list of graphs
-  TList *fList;
+  Double_t fPulls[AliMUONAlignment::fNGlobal];
 
   /// list of track records
   TClonesArray *fRecords;