]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONAlignmentTask.cxx
rewrote AliMUONAlignment to - use AliMillePede2 instead of AliMillepede - use AliMill...
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignmentTask.cxx
index a32e639fac466ea90732fa84e396c26109f01b94..370cdc89b62a6982673c93cea5a749b75f4bacd7 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
 //-----------------------------------------------------------------------------
 
-#include <fstream>
-
-#include <TString.h>
-#include <TGraphErrors.h>
-#include <TTree.h>
-#include <TChain.h>
-#include <TClonesArray.h>
-#include <TGeoGlobalMagField.h>
-#include <TGeoManager.h>
-#include <Riostream.h>
+// this class header must always come first
+#include "AliMUONAlignmentTask.h"
 
-#include "AliAnalysisTask.h"
+// 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 "AliCDBMetaData.h"
 #include "AliCDBId.h"
 #include "AliGeomManager.h"
-#include "AliLog.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 "AliMUONGeometryTransformer.h"
 #include "AliMUONESDInterface.h"
 
-#include "AliMUONAlignmentTask.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>
+#include <TClonesArray.h>
+#include <TGeoGlobalMagField.h>
+#include <TGeoManager.h>
+#include <Riostream.h>
 
-///\cond CLASSIMP  
+///\cond CLASSIMP
 ClassImp(AliMUONAlignmentTask)
 ///\endcond
 
-// //________________________________________________________________________
-// AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name) 
-//   : AliAnalysisTask(name, ""),
-//     fESD(0x0),
-//     fAlign(0x0),
-//     fGeoFilename(""),
-//     fTransform(0x0),
-//     fTrackTot(0),
-//     fTrackOk(0),
-//     fMSDEx(0x0), 
-//     fMSDEy(0x0), 
-//     fMSDEz(0x0),
-//     fMSDEp(0x0),
-//     fList(0x0)
-// {
-//   /// Default Constructor
-//   // Define input and output slots here
-//   // Input slot #0 works with a TChain
-//   DefineInput(0, TChain::Class());
-//   // Output slot #0 writes NTuple/histos into a TList
-//   DefineOutput(0, TList::Class());  
-
-//   // initialize parameters ...
-//   for(Int_t k=0;k<4*156;k++) {
-//     fParameters[k]=0.;
-//     fErrors[k]=0.;
-//     fPulls[k]=0.;
-//   }
-
-//   fAlign = new AliMUONAlignment();
-//   fTransform = new AliMUONGeometryTransformer();  
-// }
-
 //________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *geofilename, const char *defaultocdb, const char *misalignocdb) 
-  : AliAnalysisTask(name, ""),
-    fESD(0x0),
-    fAlign(0x0),
-    fGeoFilename(geofilename),
-    fMisAlignOCDB(misalignocdb),
-    fDefaultOCDB(defaultocdb),
-    fTransform(0x0),
-    fTrackTot(0),
-    fTrackOk(0),
-    fMSDEx(0x0), 
-    fMSDEy(0x0), 
-    fMSDEz(0x0),
-    fMSDEp(0x0),
-    fList(0x0)
+AliMUONAlignmentTask::AliMUONAlignmentTask( const char *name ):
+  AliAnalysisTaskSE(name),
+  fReadRecords( kFALSE ),
+  fWriteRecords( kFALSE ),
+  fDoAlignment( kTRUE ),
+  fMergeAlignmentCDBs( kTRUE ),
+  fForceBField( kFALSE ),
+  fBFieldOn( kFALSE ),
+  fUnbias( kFALSE ),
+  fAlign(0x0),
+  fNewAlignStorage( "local://ReAlignOCDB" ),
+  fOldGeoTransformer(0x0),
+  fNewGeoTransformer(0x0),
+  fLoadOCDBOnce(kFALSE),
+  fOCDBLoaded(kFALSE),
+  fEvent(0),
+  fTrackTot(0),
+  fTrackOk(0),
+  fRunNumberMin(0),
+  fRunNumberMax(0),
+  fRecords(0x0),
+  fRecordCount(0)
 {
   /// Default Constructor
   // Define input and output slots here
   // Input slot #0 works with a TChain
   DefineInput(0, TChain::Class());
-  // Output slot #0 writes NTuple/histos into a TList
-  DefineOutput(0, 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.;
   }
 
+  // create alignment object
   fAlign = new AliMUONAlignment();
-  fTransform = new AliMUONGeometryTransformer();  
-}
 
+  // create old geometry transformer
+  fOldGeoTransformer = new AliMUONGeometryTransformer();
 
-//________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& obj) 
-  : AliAnalysisTask(obj),
-    fESD(0x0),
-    fAlign(0x0),
-    fGeoFilename(""),
-    fMisAlignOCDB(""),
-    fDefaultOCDB(""),
-    fTransform(0x0),
-    fTrackTot(0),
-    fTrackOk(0),
-    fMSDEx(0x0), 
-    fMSDEy(0x0), 
-    fMSDEz(0x0),
-    fMSDEp(0x0),
-    fList(0x0)
-{
-  /// Copy constructor
-  fESD = obj.fESD;
-  fAlign = obj.fAlign;
-  fGeoFilename = obj.fGeoFilename;
-  fTransform = obj.fTransform;
-  fTrackTot = obj.fTrackTot;  
-  fTrackOk = obj.fTrackOk;  
-  fMSDEx = obj.fMSDEx; 
-  fMSDEy = obj.fMSDEy; 
-  fMSDEz = obj.fMSDEz;
-  fMSDEp = obj.fMSDEp;
-  fList = obj.fList;  
 }
 
 //________________________________________________________________________
-AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other) 
+AliMUONAlignmentTask::~AliMUONAlignmentTask()
 {
-  /// Assignment
-  AliAnalysisTask::operator=(other);
-  fESD = other.fESD;
-  fAlign = other.fAlign;
-  fGeoFilename = other.fGeoFilename;
-  fMisAlignOCDB = other.fMisAlignOCDB;
-  fDefaultOCDB = other.fDefaultOCDB;
-  fTransform = other.fTransform;
-  fTrackTot = other.fTrackTot;  
-  fTrackOk = other.fTrackOk;  
-  fMSDEx = other.fMSDEx; 
-  fMSDEy = other.fMSDEy; 
-  fMSDEz = other.fMSDEz;
-  fMSDEp = other.fMSDEp;
-  fList = other.fList;  
-  
-  return *this;
-}
-
-//________________________________________________________________________
-AliMUONAlignmentTask::~AliMUONAlignmentTask() 
-{ 
-  /// Destructor
-  if (fAlign) delete fAlign;
-  if (fTransform) delete fTransform;
+  /// destructor
+  delete fAlign;
+  delete fOldGeoTransformer;
+  delete fNewGeoTransformer;
 }
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::LocalInit() 
+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
+  **/
 
-  Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
+  /* must run alignment when reading records */
+  if( fReadRecords && !fDoAlignment )
+  {
 
-  // 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( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
+    SetDoAlignment( kTRUE );
 
+  }
 
-  fAlign->InitGlobalParameters(fParameters);
+  // 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 ) )
+  { AliFatal( "Need at least one of the three flags fReadRecords, fWriteRecords and fDoAlignment set to kTRUE" ); }
+
+  /* cannot read and write records at the same time */
+  if( fReadRecords && fWriteRecords )
+  { AliFatal( "Cannot set both fReadRecords and fWriteRecords to kTRUE at the same time" ); }
+
+  /*
+  fix detectors
+  warning, counting chambers from 1.
+  this must be done before calling the Init method
+  */
+  if( fDoAlignment )
+  {
+
+    fAlign->FixChamber(1);
+    fAlign->FixChamber(10);
+
+  } else {
+
+    AliInfo( "Not fixing detector elements, since alignment is not required" );
+
+  }
 
+  // initialize
+  fAlign->Init();
 
-  fTransform->LoadGeometryData();
-  fAlign->SetGeometryTransformer(fTransform);
+  // use unbiased residuals
+  fAlign->SetUnbias( fUnbias );
 
   // Do alignment with magnetic field off
-  fAlign->SetBFieldOn(kFALSE);
-  
-  // Set tracking station to use
-  //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-  Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kFALSE,kTRUE,kTRUE,kTRUE,kTRUE};
-
-  // Set degrees of freedom to align (see AliMUONAlignment)
-  fAlign->AllowVariations(bChOnOff);
-
-  // 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(1020);
-
-  // 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);
+  fAlign->SetBFieldOn( fBFieldOn );
+
+  // Set expected resolution (see AliMUONAlignment)
+  fAlign->SetSigmaXY(0.15,0.10);
+
+  // initialize global parameters to provided values
+  fAlign->InitGlobalParameters(fParameters);
 
 }
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::ConnectInputData(Option_t *) 
+void AliMUONAlignmentTask::UserCreateOutputObjects()
 {
-  /// Connect ESD here. Called on each input data change.
 
-  // Connect ESD here
-  TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
-  if (!esdTree) {
-    Printf("ERROR: Could not read chain from input slot 0");
-  } 
-  else {
-    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());   
-    if (!esdH) {
-      Printf("ERROR: Could not get ESDInputHandler");
-    } 
-    else {
-      fESD = esdH->GetEvent();
-    }
+  // connect AOD output
+  if( fWriteRecords )
+  {
+    AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
+    if( handler )
+    {
+
+      // get AOD output handler and add Branch
+      fRecords = new TClonesArray( "AliMillePedeRecord", 0 );
+      fRecords->SetName( "records" );
+      handler->AddBranch("TClonesArray", &fRecords);
+
+      fRecordCount = 0;
+
+    } else AliFatal( "Error: invalid output event handler" );
+
   }
-}
 
-//________________________________________________________________________
-void AliMUONAlignmentTask::CreateOutputObjects()
-{
-  /// Executed once on each worker (machine actually running the analysis code)
-  //
-  // This method has to be called INSIDE the user redefined CreateOutputObjects
-  // method, before creating each object corresponding to the output containers
-  // that are to be written to a file. This need to be done in general for the big output
-  // objects that may not fit memory during processing. 
-  //  OpenFile(0); 
-
-  // 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);
 }
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::Exec(Option_t *) 
+void AliMUONAlignmentTask::UserExec(Option_t *)
 {
-  /// Main loop, called for each event
-  if (!fESD) {
-    Printf("ERROR: fESD not available");
-    return;
+
+  // 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 )
+  {
+    fRecords->Clear();
+    fRecordCount = 0;
   }
 
-  Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
-
-
-  Int_t nTracks = Int_t(fESD->GetNumberOfMuonTracks());
-  //  if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
-  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
-    AliESDMuonTrack* esdTrack = fESD->GetMuonTrack(iTrack);
-    if (!esdTrack->ClustersStored()) continue;
-    Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
-//     fInvBenMom->Fill(invBenMom);
-//     fBenMom->Fill(1./invBenMom);
-    if (TMath::Abs(invBenMom)<=1.04) {
-      AliMUONTrack track;
-      AliMUONESDInterface::ESDToMUON(*esdTrack, track);
-      fAlign->ProcessTrack(&track);
-      fAlign->LocalFit(fTrackOk++,trackParams,0);
+  if( fReadRecords )
+  {
+
+    AliAODEvent* lAOD( dynamic_cast<AliAODEvent*>(InputEvent() ) );
+
+    // check AOD
+    if( !lAOD )
+    {
+      AliInfo("Error: AOD event not available");
+      return;
+    }
+
+    // read records
+    TClonesArray* records = static_cast<TClonesArray*>( lAOD->FindListObject( "records" ) );
+    if( !records )
+    {
+      AliInfo( "Unable to read object name \"records\"" );
+      return;
+    }
+
+    // get number of records and print
+    const Int_t lRecordsCount( records->GetEntriesFast() );
+    fTrackTot += lRecordsCount;
+
+    // loop over records
+    for( Int_t index = 0; index < lRecordsCount; ++index )
+    {
+
+      if( AliMillePedeRecord* record = dynamic_cast<AliMillePedeRecord*>( records->UncheckedAt(index) ) )
+      {
+
+        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 ) );
+
+    }
+
+  } else {
+
+    /// Main loop, called for each event
+    AliESDEvent* lESD( dynamic_cast<AliESDEvent*>(InputEvent()) );
+
+    // check ESD
+    if( !lESD )
+    {
+      AliInfo("Error: ESD event not available");
+      return;
+    }
+
+    Int_t nTracks = Int_t(lESD->GetNumberOfMuonTracks());
+    for( Int_t iTrack = 0; iTrack < nTracks; iTrack++ )
+    {
+
+      AliESDMuonTrack* esdTrack = lESD->GetMuonTrack(iTrack);
+      if (!esdTrack->ContainTrackerData()) continue;
+      if (!esdTrack->ContainTriggerData()) continue;
+
+      Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
+      if (TMath::Abs(invBenMom)<=1.04)
+      {
+
+        AliMUONTrack track;
+        AliMUONESDInterface::ESDToMUON(*esdTrack, track);
+
+        // process track and retrieve corresponding records, for storage
+        const AliMillePedeRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
+        ++fTrackOk;
+
+        // store in array
+        if( fWriteRecords && fRecords )
+        {
+          new((*fRecords)[fRecordCount]) AliMillePedeRecord( *lRecords );
+          ++fRecordCount;
+        }
+
+      }
+
+      ++fTrackTot;
+
+      if(!(fTrackTot%100))
+      { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
+
+      // Post final data. Write histo list to a file with option "RECREATE"
+      // PostData(1,fList);
+
     }
-    fTrackTot++;
+
+    // save AOD
+    if( fWriteRecords && fRecordCount > 0 )
+    {
+      AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
+      if( 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" );
+
+    }
+
   }
-  
-  // Post final data. Write histo list to a file with option "RECREATE"
-  PostData(0,fList);
-}      
+
+}
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::Terminate(const Option_t*)
+void AliMUONAlignmentTask::FinishTaskOutput()
 {
+
   /// Called once per task on the client machine at the end of the analysis.
+  AliInfo( Form( "Processed %i tracks.", fTrackTot ) );
+  AliInfo( Form( "Accepted %i tracks.", fTrackOk ) );
+
+  // stop here if no alignment is to be performed
+  if( !fDoAlignment ) return;
+
+  AliLog::SetGlobalLogLevel(AliLog::kInfo);
 
-  cout << "Processed " << fTrackTot << " Tracks." << endl;
   // Perform global fit
   fAlign->GlobalFit(fParameters,fErrors,fPulls);
-
-  cout << "Done with GlobalFit " << endl;
-
-  // 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++){
-    DEidErr[iDE] = 0.;
-    DEid[iDE] = idOffset+100;
-    DEid[iDE] += iDE; 
+
+  for( Int_t iDE=0; iDE<AliMUONAlignment::fgNDetElem; iDE++ )
+  {
+
+    Int_t detElemId = idOffset+100;
+    detElemId += iDE;
     lSDetElemCh = 0;
-    for(Int_t iCh=0; iCh<9; iCh++){
-      lSDetElemCh += lNDetElemCh[iCh];
-      if (iDE>=lSDetElemCh) {
-       DEid[iDE] += 100;
-       DEid[iDE] -= lNDetElemCh[iCh];
+    for(Int_t iCh=0; iCh<9; iCh++)
+    {
+
+      lSDetElemCh += AliMUONAlignment::fgNDetElemCh[iCh];
+      if (iDE>=lSDetElemCh)
+      {
+        detElemId += 100;
+        detElemId -= AliMUONAlignment::fgNDetElemCh[iCh];
       }
+
     }
-    MSDEx[iDE]=fParameters[3*iDE+0];
-    MSDEy[iDE]=fParameters[3*iDE+1];
-    MSDEz[iDE]=fParameters[3*iDE+3];
-    MSDEp[iDE]=fParameters[3*iDE+2];
-    MSDExErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+0);
-    MSDEyErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+1);
-    MSDEzErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+3);
-    MSDEpErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+2);
-    fMSDEx->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+0]);
-    fMSDEx->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+0));
-    fMSDEy->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+1]);
-    fMSDEy->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+1));
-    fMSDEp->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+2]);
-    fMSDEp->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+2));
-    fMSDEz->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+3]);
-    fMSDEz->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+3));
+
   }
 
   // Post final data. Write histo list to a file with option "RECREATE"
-  PostData(0,fList);
+  // PostData(1,fList);
 
+  // store misalignments from OCDB into old transformers
+  if( fMergeAlignmentCDBs )
+  { SaveMisAlignmentData( fOldGeoTransformer ); }
 
   // Re Align
-  AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true); 
-  newTransform->WriteTransformations("transform2ReAlign.dat");
-  
+  fNewGeoTransformer = fAlign->ReAlign( fOldGeoTransformer, fParameters, true );
+
   // Generate realigned data in local cdb
-  const TClonesArray* array = newTransform->GetMisAlignmentData();
+  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);
+
+  // recover default storage full name (raw:// cannot be used to set specific storage)
   AliCDBManager* cdbManager = AliCDBManager::Instance();
-  cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
-  cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.Data());
-  
+
+  // unload old alignment path
+  if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
+  { cdbManager->UnloadFromCache("MUON/Align/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()); 
+  AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 
 }
 
-//-----------------------------------------------------------------------
-void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
+//________________________________________________________________________
+void AliMUONAlignmentTask::NotifyRun()
 {
-  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
-  
-  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
-  if (!gGeoManager) {
-    AliGeomManager::LoadGeometry(geoFilename);
-    if (!gGeoManager) {
-      AliError(Form("Getting geometry from file %s failed", "generated/galice.root"));
+
+  /// 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* 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() )
+    {
+
+      AliInfo( "Loading magnetic field" );
+      if( !AliMUONCDB::LoadField() )
+      {
+        AliError( "Failed to load magnetic field" );
+        return;
+      }
+
+    } else { AliInfo( "Magnetic field already locked" ); }
+
+    // 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") ) );
+    }
+
+    AliInfo( "Loading muon mapping" );
+    if( !AliMUONCDB::LoadMapping(kFALSE) )
+    {
+      AliError( "Failed to load muon mapping" );
+      return;
+    }
+
+    AliInfo( "Assigning field to Track extrapolator" );
+    AliMUONTrackExtrap::SetField();
+
+  }
+
+  // load geometry if needed
+  if( !AliGeomManager::GetGeometry() )
+  {
+
+    // 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;
+
   }
-    
-  // Load mapping
-  AliCDBManager* man = AliCDBManager::Instance();
-  man->SetDefaultStorage(defaultOCDB);
-  man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
-  man->Print();
-  man->SetRun(0);
-  if ( ! AliMpCDB::LoadDDLStore() ) {
-    AliError("Could not access mapping from OCDB !");
-    exit(-1);
+
+  // get TClonesArray and check
+  TClonesArray* misArray = (TClonesArray*)cdbEntry->GetObject();
+  if( !misArray )
+  {
+
+    AliError( "unable to load old misalignment array" );
+    return;
+
   }
 
-  // set mag field
-  if (!TGeoGlobalMagField::Instance()->GetField()) {
-    printf("Loading field map...\n");
-    AliGRPManager *grpMan = new AliGRPManager();
-    grpMan->ReadGRPEntry();
-    grpMan->SetMagField();
-    delete grpMan;
+  // 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 );
+
+    // 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);
   }
-  // set the magnetic field for track extrapolations
-  AliMUONTrackExtrap::SetField();
-  
+
+  return;
 }
+