]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONAlignmentTask.cxx
o updates (Alla, Filip)
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignmentTask.cxx
index 661fd000c99dc40bc8977e31d4bc74120054c8ba..3fbbd6c64cd6a23ca6574526f6ad942b321d72a3 100644 (file)
 /// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
 //-----------------------------------------------------------------------------
 
-#include <fstream>
-
-#include <TString.h>
-#include <TError.h>
-#include <TGraphErrors.h>
-#include <TTree.h>
-#include <TChain.h>
-#include <TClonesArray.h>
-#include <TGeoGlobalMagField.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 "AliAODInputHandler.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODHeader.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 "AliMpCDB.h"
+#include "AliMUONCDB.h"
 #include "AliMUONAlignment.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackExtrap.h"
 #include "AliMUONGeometryTransformer.h"
 #include "AliMUONESDInterface.h"
 
-#include "AliMUONAlignmentTask.h"
+// system headers
+#include <fstream>
+#include <TString.h>
+#include <TError.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) 
-  : AliAnalysisTask(name, ""),
-    fESD(0x0),
-    fAlign(0x0),
-    fGeoFilename(geofilename),
-    fTransform(0x0),
-    fTrackTot(0),
-    fTrackOk(0),
-    fMSDEx(0x0), 
-    fMSDEy(0x0), 
-    fMSDEz(0x0),
-    fMSDEp(0x0),
-    fList(0x0)
+AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *newalignocdb, const char *oldalignocdb, const char *defaultocdb, const char *geofilename):
+  AliAnalysisTaskSE(name),
+  fReadRecords( kFALSE ),
+  fWriteRecords( kFALSE ),
+  fDoAlignment( kTRUE ),
+  fAlign(0x0),
+  fGeoFilename(geofilename),
+  fDefaultStorage(defaultocdb),
+  fOldAlignStorage(oldalignocdb),
+  fNewAlignStorage(newalignocdb),
+  fOldGeoTransformer(NULL),
+  fNewGeoTransformer(NULL),
+  fLoadOCDBOnce(kFALSE),
+  fOCDBLoaded(kFALSE),
+  fTrackTot(0),
+  fTrackOk(0),
+  fLastRunNumber(-1),
+  fMSDEx(0x0),
+  fMSDEy(0x0),
+  fMSDEz(0x0),
+  fMSDEp(0x0),
+  fList(0x0),
+  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());  
+  DefineOutput(1, TList::Class());
 
   // initialize parameters ...
-  for(Int_t k=0;k<4*156;k++) {
+  for(Int_t k=0;k<4*156;k++)
+  {
     fParameters[k]=0.;
     fErrors[k]=0.;
     fPulls[k]=0.;
   }
 
-  fAlign = new AliMUONAlignment();
-  fTransform = new AliMUONGeometryTransformer();  
-}
+  fOldGeoTransformer = new AliMUONGeometryTransformer();
 
+}
 
 //________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& obj) 
-  : AliAnalysisTask(obj),
-    fESD(0x0),
-    fAlign(0x0),
-    fGeoFilename(""),
-    fTransform(0x0),
-    fTrackTot(0),
-    fTrackOk(0),
-    fMSDEx(0x0), 
-    fMSDEy(0x0), 
-    fMSDEz(0x0),
-    fMSDEp(0x0),
-    fList(0x0)
+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 )
 {
-  /// 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;  
+
+  // 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) 
+AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other)
 {
   /// Assignment
-  AliAnalysisTask::operator=(other);
-  fESD = other.fESD;
+  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;
-  fTransform = other.fTransform;
-  fTrackTot = other.fTrackTot;  
-  fTrackOk = other.fTrackOk;  
-  fMSDEx = other.fMSDEx; 
-  fMSDEy = other.fMSDEy; 
+  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;  
-  
+  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];
+
+  }
+
   return *this;
 }
 
 //________________________________________________________________________
-AliMUONAlignmentTask::~AliMUONAlignmentTask() 
-{ 
-  /// Destructor
-  if (fAlign) delete fAlign;
-  if (fTransform) delete fTransform;
+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.
+  */
+  delete fAlign;
+  delete fOldGeoTransformer;
+  delete fNewGeoTransformer;
 }
 
 //________________________________________________________________________
-void AliMUONAlignmentTask::LocalInit() 
+void AliMUONAlignmentTask::LocalInit()
 {
-  /// Local initialization, called once per task on the client machine 
+  /// Local initialization, called once per task on the client machine
   /// where the analysis train is assembled
-  AliMpCDB::LoadMpSegmentation();
 
-  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
-  if ( ! AliGeomManager::GetGeometry() ) {
-    AliGeomManager::LoadGeometry(fGeoFilename.Data());
-    if (! AliGeomManager::GetGeometry() ) {
-      Error("MUONAlignment", "getting geometry from file %s failed", fGeoFilename.Data());
-      return;
-    }
-  }
-  
-  // set  mag field 
-  // waiting for mag field in CDB 
-  if (!TGeoGlobalMagField::Instance()->GetField()) {
-    printf("Loading field map...\n");
-    AliMagF* field = new AliMagF("Maps","Maps",2,0.,0., 10.,AliMagF::k5kG);
-    TGeoGlobalMagField::Instance()->SetField(field);
+  // print configuration
+  AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
+  AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
+  AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "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" ); }
+
+  /* must run alignment when reading records */
+  if( fReadRecords && !fDoAlignment )
+  {
+
+    AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
+    SetDoAlignment( kTRUE );
+
   }
-  // set the magnetic field for track extrapolations
-  AliMUONTrackExtrap::SetField();
 
   // Set initial values here, good guess may help convergence
-  // St 1 
+  // St 1
   //  Int_t iPar = 0;
-  //  fParameters[iPar++] =  0.010300 ;  fParameters[iPar++] =  0.010600 ;  fParameters[iPar++] =  0.000396 ;  
-
+  //  fParameters[iPar++] =  0.010300 ;  fParameters[iPar++] =  0.010600 ;  fParameters[iPar++] =  0.000396 ;
 
+  fAlign = new AliMUONAlignment();
   fAlign->InitGlobalParameters(fParameters);
 
-
-  fTransform->LoadGeometryData();
-  fAlign->SetGeometryTransformer(fTransform);
+//  AliCDBManager::Instance()->Print();
+//
+//  fAlign->SetGeometryTransformer(fOldGeoTransformer);
 
   // 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};
+  Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
 
   // Set degrees of freedom to align (see AliMUONAlignment)
   fAlign->AllowVariations(bChOnOff);
 
+  // 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 
+  // 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);
@@ -249,8 +289,9 @@ void AliMUONAlignmentTask::LocalInit()
   fAlign->SetSpecLROnOff(bChOnOff);
 
     // Here we can fix some detection elements
-  fAlign->FixDetElem(908);
-  fAlign->FixDetElem(1020);
+  //  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};
@@ -260,102 +301,219 @@ void AliMUONAlignmentTask::LocalInit()
 }
 
 //________________________________________________________________________
-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();
-    }
+
+  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 )
+  {
+    AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
+    if( handler )
+    {
+
+      // get AOD output handler and add Branch
+      fRecords = new TClonesArray( "AliMUONAlignmentTrackRecord", 0 );
+      fRecords->SetName( "records" );
+      handler->AddBranch("TClonesArray", &fRecords);
+
+      fRecordCount = 0;
+
+    } else AliInfo( "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;
+  // print configuration
+//  AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
+//  AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
+//  AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
+  // clear array
+  if( fWriteRecords && fRecords )
+  {
+    fRecords->Clear();
+    fRecordCount = 0;
   }
 
+  // local track parameters
   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,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( AliMUONAlignmentTrackRecord* record = dynamic_cast<AliMUONAlignmentTrackRecord*>( records->UncheckedAt(index) ) )
+      {
+
+        fAlign->ProcessTrack( record, fDoAlignment );
+        if( fDoAlignment )
+        { fAlign->LocalFit( fTrackOk++, trackParams, 0 ); }
+
+      } else AliInfo( Form( "Invalid record at %i", index ) );
+
+    }
+
+  } else {
+
+    /// Main loop, called for each event
+    AliESDEvent* lESD( dynamic_cast<AliESDEvent*>(InputEvent()) );
 
-  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);
+    // check ESD
+    if( !lESD )
+    {
+      AliInfo("Error: ESD event not available");
+      return;
     }
-    fTrackTot++;
+
+// 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++ )
+    {
+
+      AliESDMuonTrack* esdTrack = lESD->GetMuonTrack(iTrack);
+      if (!esdTrack->ClustersStored()) continue;
+      if (!esdTrack->ContainTriggerData()) 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);
+
+        // 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++;
+
+        // store in array
+        if( fWriteRecords && fRecords )
+        {                                      
+          new((*fRecords)[fRecordCount]) AliMUONAlignmentTrackRecord( *lRecords );
+          ++fRecordCount;
+        }
+
+      }                        
+                       
+      ++fTrackTot;
+
+      if(!(fTrackTot%1000))
+      { 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);
+
+    }
+
+    // save AOD
+    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" );
+    }
+
   }
-  
-  // Post final data. Write histo list to a file with option "RECREATE"
-  PostData(0,fList);
-}      
+
+}
 
 //________________________________________________________________________
 void AliMUONAlignmentTask::Terminate(const Option_t*)
+{ return; }
+
+//________________________________________________________________________
+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);
+//   // 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};
@@ -373,59 +531,140 @@ void AliMUONAlignmentTask::Terminate(const Option_t*)
   // 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<lNDetElem; iDE++)
+  {
+
     DEidErr[iDE] = 0.;
     DEid[iDE] = idOffset+100;
-    DEid[iDE] += iDE; 
+    DEid[iDE] += iDE;
     lSDetElemCh = 0;
-    for(Int_t iCh=0; iCh<9; iCh++){
+    for(Int_t iCh=0; iCh<9; iCh++)
+    {
       lSDetElemCh += lNDetElemCh[iCh];
-      if (iDE>=lSDetElemCh) {
-       DEid[iDE] += 100;
-       DEid[iDE] -= lNDetElemCh[iCh];
+
+      if (iDE>=lSDetElemCh)
+      {
+        DEid[iDE] += 100;
+        DEid[iDE] -= lNDetElemCh[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));
+
+    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(0,fList);
+  PostData(1,fList);
 
+  // HUGO: stop here to test reproducibility
+       //  return;
 
   // Re Align
-  AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true); 
-  newTransform->WriteTransformations("transform2ReAlign.dat");
-  
+  fNewGeoTransformer = fAlign->ReAlign(fOldGeoTransformer,fParameters,true);
+  //  newTransform->WriteTransformations("transform2ReAlign.dat");
+
   // 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);
-   
+
   // CDB manager
-  AliCDBManager* cdbManager = AliCDBManager::Instance();
-  cdbManager->SetDefaultStorage("local://ReAlignCDB");
-  
+  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());
+
+       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());
+
+
   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()); 
-  cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
+  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");
 
 }
 
+//________________________________________________________________________
+void AliMUONAlignmentTask::NotifyRun()
+{
+  
+  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");
+
+  // 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);
+
+  if (!AliMUONCDB::LoadField()) { AliError("Problem field"); return;}
+
+  // set the magnetic field for track extrapolations
+  AliMUONTrackExtrap::SetField();
+
+  if (!AliMUONCDB::LoadMapping(kTRUE)) { AliError("Problem 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());
+
+       // 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")
+  {
+
+               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());
+               
+    AliGeomManager::ApplyAlignObjsFromCDB("MUON");
+
+  }
+
+  // fOldGeoTransformer = new AliMUONGeometryTransformer();
+       fOldGeoTransformer->LoadGeometryData();
+       fAlign->SetGeometryTransformer(fOldGeoTransformer);
+
+       fOCDBLoaded = kTRUE;
+
+}