[u/mrichter/AliRoot.git] / MUON / AliMUONAlignmentTask.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

// $Id$

//-----------------------------------------------------------------------------
/// \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.
/// 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
//-----------------------------------------------------------------------------

#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>

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliMagF.h"
#include "AliCDBManager.h"
#include "AliGRPManager.h"
#include "AliCDBMetaData.h"
#include "AliCDBId.h"
#include "AliGeomManager.h"

#include "AliMpCDB.h"
#include "AliMUONAlignment.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackExtrap.h"
#include "AliMUONTrackParam.h"
#include "AliMUONGeometryTransformer.h"
#include "AliMUONESDInterface.h"

#include "AliMUONAlignmentTask.h"

///\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),
    fLastRunNumber(-1),
    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 AliMUONAlignmentTask& obj) 
  : AliAnalysisTask(obj),
    fESD(0x0),
    fAlign(0x0),
    fGeoFilename(""),
    fMisAlignOCDB(""),
    fDefaultOCDB(""),
    fTransform(0x0),
    fTrackTot(0),
    fTrackOk(0),
    fLastRunNumber(-1),
    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;  
  fLastRunNumber = obj.fLastRunNumber;
  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]=obj.fParameters[k];
    fErrors[k]=obj.fErrors[k];
    fPulls[k]=obj.fPulls[k];
  }

}

//________________________________________________________________________
AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other) 
{
  /// 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;  
  fLastRunNumber = other.fLastRunNumber;
  fMSDEx = other.fMSDEx; 
  fMSDEy = other.fMSDEy; 
  fMSDEz = other.fMSDEz;
  fMSDEp = other.fMSDEp;
  fList = other.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];
  }
  
  return *this;
}

//________________________________________________________________________
AliMUONAlignmentTask::~AliMUONAlignmentTask() 
{ 
  /// Destructor
  if (fAlign) delete fAlign;
  if (fTransform) delete fTransform;
}

//________________________________________________________________________
void AliMUONAlignmentTask::LocalInit() 
{
  /// Local initialization, called once per task on the client machine 
  /// where the analysis train is assembled
  fLastRunNumber = 0; 
  //  Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
  Prepare(fGeoFilename.Data(),"local://$ALICE_ROOT/OCDB",fMisAlignOCDB.Data());
  fLastRunNumber = -1;

  // 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 ;  


  fAlign->InitGlobalParameters(fParameters);


  fTransform->LoadGeometryData();
  fAlign->SetGeometryTransformer(fTransform);

  // 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);

}

//________________________________________________________________________
void AliMUONAlignmentTask::ConnectInputData(Option_t *) 
{
  /// 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();
    }
  }
}

//________________________________________________________________________
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 *) 
{
  /// Main loop, called for each event
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }

  Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
  if (fESD->GetRunNumber()!=fLastRunNumber){
    fLastRunNumber = fESD->GetRunNumber();
    Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
  }
 
  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);
    }
    fTrackTot++;
    cout << "Processed " << fTrackTot << " Tracks." << endl;
  }
  
  // Post final data. Write histo list to a file with option "RECREATE"
  PostData(0,fList);
}      

//________________________________________________________________________
void AliMUONAlignmentTask::FinishTaskOutput()
{
  /// Called once per task on the client machine at the end of the analysis.

  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; 
    lSDetElemCh = 0;
    for(Int_t iCh=0; iCh<9; iCh++){
      lSDetElemCh += 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));
  }

  // Post final data. Write histo list to a file with option "RECREATE"
  PostData(0,fList);

  // Re Align
  AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true); 
  newTransform->WriteTransformations("transform2ReAlign.dat");
  
  // Generate realigned data in local cdb
  const TClonesArray* array = newTransform->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(fDefaultOCDB.Data());
  cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.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);

}

//________________________________________________________________________
void AliMUONAlignmentTask::Terminate(const Option_t*)
{
  /// Called once per task on the client machine at the end of the analysis.

}

//-----------------------------------------------------------------------
void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
{
  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
      
  // Load mapping
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage(defaultOCDB);
  man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
  man->Print();
  man->SetRun(fLastRunNumber);
  if ( ! AliMpCDB::LoadDDLStore() ) {
    Error("MUONRefit","Could not access mapping from OCDB !");
    exit(-1);
  }

  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
  if (!gGeoManager) {
    AliGeomManager::LoadGeometry(geoFilename);
    if (!gGeoManager) {
      Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
      return;
    }
  }

  // set mag field
  if (!TGeoGlobalMagField::Instance()->GetField()) {
    printf("Loading field map...\n");
    AliGRPManager *grpMan = new AliGRPManager();
    grpMan->ReadGRPEntry();
    grpMan->SetMagField();
    delete grpMan;
  }
  // set the magnetic field for track extrapolations
  AliMUONTrackExtrap::SetField();
  
}
Commit	Line	Data
81f1d3ae	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// $Id$
	17
4d610fd5	18	//-----------------------------------------------------------------------------
	19	/// \class AliMUONAlignmentTask
	20	/// AliAnalysisTask to align the MUON spectrometer.
	21	/// The Task reads as input ESDs and feeds the MUONTracks to AliMUONAlignment.
	22	/// The alignment itself is performed by AliMillepede.
	23	/// A OCDB entry is written with the alignment parameters.
	24	/// A list of graph are written to monitor the alignment parameters.
	25	///
	26	/// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
	27	//-----------------------------------------------------------------------------
	28
	29	#include <fstream>
	30
	31	#include <TString.h>
2df5c2cf	32	#include <TError.h>
4d610fd5	33	#include <TGraphErrors.h>
	34	#include <TTree.h>
	35	#include <TChain.h>
	36	#include <TClonesArray.h>
	37	#include <TGeoGlobalMagField.h>
cd8521dd	38	#include <TGeoManager.h>
4d610fd5	39	#include <Riostream.h>
	40
	41	#include "AliAnalysisTask.h"
	42	#include "AliAnalysisManager.h"
	43	#include "AliESDInputHandler.h"
	44	#include "AliESDEvent.h"
	45	#include "AliESDMuonTrack.h"
	46	#include "AliMagF.h"
	47	#include "AliCDBManager.h"
cd8521dd	48	#include "AliGRPManager.h"
4d610fd5	49	#include "AliCDBMetaData.h"
	50	#include "AliCDBId.h"
	51	#include "AliGeomManager.h"
	52
	53	#include "AliMpCDB.h"
	54	#include "AliMUONAlignment.h"
	55	#include "AliMUONTrack.h"
	56	#include "AliMUONTrackExtrap.h"
	57	#include "AliMUONTrackParam.h"
	58	#include "AliMUONGeometryTransformer.h"
	59	#include "AliMUONESDInterface.h"
	60
	61	#include "AliMUONAlignmentTask.h"
	62
	63	///\cond CLASSIMP
	64	ClassImp(AliMUONAlignmentTask)
	65	///\endcond
	66
	67	// //________________________________________________________________________
	68	// AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name)
	69	// : AliAnalysisTask(name, ""),
	70	// fESD(0x0),
	71	// fAlign(0x0),
	72	// fGeoFilename(""),
	73	// fTransform(0x0),
	74	// fTrackTot(0),
	75	// fTrackOk(0),
	76	// fMSDEx(0x0),
	77	// fMSDEy(0x0),
	78	// fMSDEz(0x0),
	79	// fMSDEp(0x0),
	80	// fList(0x0)
	81	// {
	82	// /// Default Constructor
	83	// // Define input and output slots here
	84	// // Input slot #0 works with a TChain
	85	// DefineInput(0, TChain::Class());
	86	// // Output slot #0 writes NTuple/histos into a TList
	87	// DefineOutput(0, TList::Class());
	88
	89	// // initialize parameters ...
	90	// for(Int_t k=0;k<4*156;k++) {
	91	// fParameters[k]=0.;
	92	// fErrors[k]=0.;
	93	// fPulls[k]=0.;
	94	// }
	95
	96	// fAlign = new AliMUONAlignment();
	97	// fTransform = new AliMUONGeometryTransformer();
	98	// }
	99
	100	//________________________________________________________________________
cd8521dd	101	AliMUONAlignmentTask::AliMUONAlignmentTask(const char name, const char geofilename, const char defaultocdb, const char misalignocdb)
4d610fd5	102	: AliAnalysisTask(name, ""),
	103	fESD(0x0),
	104	fAlign(0x0),
	105	fGeoFilename(geofilename),
cd8521dd	106	fMisAlignOCDB(misalignocdb),
cd8521dd	107	fDefaultOCDB(defaultocdb),
4d610fd5	108	fTransform(0x0),
	109	fTrackTot(0),
	110	fTrackOk(0),
2df5c2cf	111	fLastRunNumber(-1),
4d610fd5	112	fMSDEx(0x0),
	113	fMSDEy(0x0),
	114	fMSDEz(0x0),
	115	fMSDEp(0x0),
	116	fList(0x0)
	117	{
	118	/// Default Constructor
	119	// Define input and output slots here
	120	// Input slot #0 works with a TChain
	121	DefineInput(0, TChain::Class());
	122	// Output slot #0 writes NTuple/histos into a TList
	123	DefineOutput(0, TList::Class());
	124
	125	// initialize parameters ...
	126	for(Int_t k=0;k<4*156;k++) {
	127	fParameters[k]=0.;
	128	fErrors[k]=0.;
	129	fPulls[k]=0.;
	130	}
	131
	132	fAlign = new AliMUONAlignment();
	133	fTransform = new AliMUONGeometryTransformer();
	134	}
	135
	136
	137	//________________________________________________________________________
	138	AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& obj)
	139	: AliAnalysisTask(obj),
	140	fESD(0x0),
	141	fAlign(0x0),
	142	fGeoFilename(""),
cd8521dd	143	fMisAlignOCDB(""),
cd8521dd	144	fDefaultOCDB(""),
4d610fd5	145	fTransform(0x0),
	146	fTrackTot(0),
	147	fTrackOk(0),
2df5c2cf	148	fLastRunNumber(-1),
4d610fd5	149	fMSDEx(0x0),
	150	fMSDEy(0x0),
	151	fMSDEz(0x0),
	152	fMSDEp(0x0),
	153	fList(0x0)
	154	{
	155	/// Copy constructor
	156	fESD = obj.fESD;
	157	fAlign = obj.fAlign;
	158	fGeoFilename = obj.fGeoFilename;
	159	fTransform = obj.fTransform;
	160	fTrackTot = obj.fTrackTot;
	161	fTrackOk = obj.fTrackOk;
2df5c2cf	162	fLastRunNumber = obj.fLastRunNumber;
4d610fd5	163	fMSDEx = obj.fMSDEx;
	164	fMSDEy = obj.fMSDEy;
	165	fMSDEz = obj.fMSDEz;
	166	fMSDEp = obj.fMSDEp;
	167	fList = obj.fList;
f4c2989f	168
	169	// initialize parameters ...
	170	for(Int_t k=0;k<4*156;k++) {
	171	fParameters[k]=obj.fParameters[k];
	172	fErrors[k]=obj.fErrors[k];
	173	fPulls[k]=obj.fPulls[k];
	174	}
	175
4d610fd5	176	}
	177
	178	//________________________________________________________________________
	179	AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other)
	180	{
	181	/// Assignment
	182	AliAnalysisTask::operator=(other);
	183	fESD = other.fESD;
	184	fAlign = other.fAlign;
	185	fGeoFilename = other.fGeoFilename;
cd8521dd	186	fMisAlignOCDB = other.fMisAlignOCDB;
cd8521dd	187	fDefaultOCDB = other.fDefaultOCDB;
4d610fd5	188	fTransform = other.fTransform;
	189	fTrackTot = other.fTrackTot;
	190	fTrackOk = other.fTrackOk;
2df5c2cf	191	fLastRunNumber = other.fLastRunNumber;
4d610fd5	192	fMSDEx = other.fMSDEx;
	193	fMSDEy = other.fMSDEy;
	194	fMSDEz = other.fMSDEz;
	195	fMSDEp = other.fMSDEp;
	196	fList = other.fList;
f4c2989f	197
	198	// initialize parameters ...
	199	for(Int_t k=0;k<4*156;k++) {
	200	fParameters[k]=other.fParameters[k];
	201	fErrors[k]=other.fErrors[k];
	202	fPulls[k]=other.fPulls[k];
	203	}
4d610fd5	204
	205	return *this;
	206	}
	207
	208	//________________________________________________________________________
	209	AliMUONAlignmentTask::~AliMUONAlignmentTask()
	210	{
	211	/// Destructor
	212	if (fAlign) delete fAlign;
	213	if (fTransform) delete fTransform;
	214	}
	215
	216	//________________________________________________________________________
	217	void AliMUONAlignmentTask::LocalInit()
	218	{
	219	/// Local initialization, called once per task on the client machine
	220	/// where the analysis train is assembled
2df5c2cf	221	fLastRunNumber = 0;
	222	// Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
	223	Prepare(fGeoFilename.Data(),"local://$ALICE_ROOT/OCDB",fMisAlignOCDB.Data());
	224	fLastRunNumber = -1;
4d610fd5	225
	226	// Set initial values here, good guess may help convergence
	227	// St 1
	228	// Int_t iPar = 0;
	229	// fParameters[iPar++] = 0.010300 ; fParameters[iPar++] = 0.010600 ; fParameters[iPar++] = 0.000396 ;
	230
	231
	232	fAlign->InitGlobalParameters(fParameters);
	233
	234
	235	fTransform->LoadGeometryData();
	236	fAlign->SetGeometryTransformer(fTransform);
	237
	238	// Do alignment with magnetic field off
	239	fAlign->SetBFieldOn(kFALSE);
	240
	241	// Set tracking station to use
	242	// Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
	243	Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kFALSE,kTRUE,kTRUE,kTRUE,kTRUE};
	244
	245	// Set degrees of freedom to align (see AliMUONAlignment)
	246	fAlign->AllowVariations(bChOnOff);
	247
	248	// Fix parameters or add constraints here
	249	// for (Int_t iSt=0; iSt<5; iSt++)
	250	// if (!bStOnOff[iSt]) fAlign->FixStation(iSt+1);
	251	for (Int_t iCh=0; iCh<10; iCh++)
	252	if (!bChOnOff[iCh]) fAlign->FixChamber(iCh+1);
	253
	254	// Left and right sides of the detector are independent, one can choose to align
	255	// only one side
	256	Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
	257	fAlign->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
	258
	259	fAlign->SetChOnOff(bChOnOff);
	260	fAlign->SetSpecLROnOff(bChOnOff);
	261
	262	// Here we can fix some detection elements
	263	fAlign->FixDetElem(908);
	264	fAlign->FixDetElem(1020);
	265
	266	// Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
	267	// Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
	268	// Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
	269	// fAlign->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
	270
	271	}
	272
	273	//________________________________________________________________________
	274	void AliMUONAlignmentTask::ConnectInputData(Option_t *)
	275	{
	276	/// Connect ESD here. Called on each input data change.
	277
	278	// Connect ESD here
	279	TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
	280	if (!esdTree) {
	281	Printf("ERROR: Could not read chain from input slot 0");
	282	}
	283	else {
	284	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	285	if (!esdH) {
	286	Printf("ERROR: Could not get ESDInputHandler");
	287	}
	288	else {
289	fESD = esdH->GetEvent();
290	}
291	}
292	}
293
294	//________________________________________________________________________
295	void AliMUONAlignmentTask::CreateOutputObjects()
296	{
297	/// Executed once on each worker (machine actually running the analysis code)
298	//
299	// This method has to be called INSIDE the user redefined CreateOutputObjects
300	// method, before creating each object corresponding to the output containers
301	// that are to be written to a file. This need to be done in general for the big output
302	// objects that may not fit memory during processing.
303	// OpenFile(0);
304
305	// Creating graphs
306	fMSDEx = new TGraphErrors(156);
307	fMSDEy = new TGraphErrors(156);
308	fMSDEz = new TGraphErrors(156);
309	fMSDEp = new TGraphErrors(156);
310
311	// Add Ntuples to the list
312	fList = new TList();
313	fList->Add(fMSDEx);
314	fList->Add(fMSDEy);
315	fList->Add(fMSDEz);
316	fList->Add(fMSDEp);
317	}
318
319	//________________________________________________________________________
320	void AliMUONAlignmentTask::Exec(Option_t *)
321	{
322	/// Main loop, called for each event
323	if (!fESD) {
324	Printf("ERROR: fESD not available");
325	return;
326	}
327
328	Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
2df5c2cf	329	if (fESD->GetRunNumber()!=fLastRunNumber){
	330	fLastRunNumber = fESD->GetRunNumber();
	331	Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
	332	}
4d610fd5	333
	334	Int_t nTracks = Int_t(fESD->GetNumberOfMuonTracks());
	335	// if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
	336	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
	337	AliESDMuonTrack* esdTrack = fESD->GetMuonTrack(iTrack);
	338	if (!esdTrack->ClustersStored()) continue;
	339	Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
	340	// fInvBenMom->Fill(invBenMom);
	341	// fBenMom->Fill(1./invBenMom);
	342	if (TMath::Abs(invBenMom)<=1.04) {
	343	AliMUONTrack track;
	344	AliMUONESDInterface::ESDToMUON(*esdTrack, track);
	345	fAlign->ProcessTrack(&track);
	346	fAlign->LocalFit(fTrackOk++,trackParams,0);
	347	}
	348	fTrackTot++;
2df5c2cf	349	cout << "Processed " << fTrackTot << " Tracks." << endl;
4d610fd5	350	}
	351
	352	// Post final data. Write histo list to a file with option "RECREATE"
	353	PostData(0,fList);
	354	}
	355
	356	//________________________________________________________________________
2df5c2cf	357	void AliMUONAlignmentTask::FinishTaskOutput()
4d610fd5	358	{
	359	/// Called once per task on the client machine at the end of the analysis.
	360
	361	cout << "Processed " << fTrackTot << " Tracks." << endl;
	362	// Perform global fit
	363	fAlign->GlobalFit(fParameters,fErrors,fPulls);
	364
	365	cout << "Done with GlobalFit " << endl;
	366
2df5c2cf	367	// // Update pointers reading them from the output slot
	368	// fList = (TList*)GetOutputData(0);
	369	// fMSDEx = (TGraphErrors*)fList->At(0);
	370	// fMSDEy = (TGraphErrors*)fList->At(1);
	371	// fMSDEz = (TGraphErrors*)fList->At(2);
	372	// fMSDEp = (TGraphErrors*)fList->At(3);
4d610fd5	373
4d610fd5	374	// Store results
9ee1d6ff	375	Double_t deId[156] = {0};
	376	Double_t msdEx[156] = {0};
	377	Double_t msdEy[156] = {0};
	378	Double_t msdEz[156] = {0};
	379	Double_t msdEp[156] = {0};
	380	Double_t deIdErr[156] = {0};
	381	Double_t msdExErr[156] = {0};
	382	Double_t msdEyErr[156] = {0};
	383	Double_t msdEzErr[156] = {0};
	384	Double_t msdEpErr[156] = {0};
4d610fd5	385	Int_t lNDetElem = 42+42+182+262+26*2;
	386	Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
	387	// Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
	388	Int_t idOffset = 0; // 400
	389	Int_t lSDetElemCh = 0;
	390	for(Int_t iDE=0; iDE<lNDetElem; iDE++){
9ee1d6ff	391	deIdErr[iDE] = 0.;
	392	deId[iDE] = idOffset+100;
	393	deId[iDE] += iDE;
4d610fd5	394	lSDetElemCh = 0;
	395	for(Int_t iCh=0; iCh<9; iCh++){
	396	lSDetElemCh += lNDetElemCh[iCh];
	397	if (iDE>=lSDetElemCh) {
9ee1d6ff	398	deId[iDE] += 100;
9ee1d6ff	399	deId[iDE] -= lNDetElemCh[iCh];
4d610fd5	400	}
4d610fd5	401	}
9ee1d6ff	402	msdEx[iDE]=fParameters[3*iDE+0];
	403	msdEy[iDE]=fParameters[3*iDE+1];
	404	msdEz[iDE]=fParameters[3*iDE+3];
	405	msdEp[iDE]=fParameters[3*iDE+2];
	406	msdExErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+0);
	407	msdEyErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+1);
	408	msdEzErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+3);
	409	msdEpErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+2);
	410	fMSDEx->SetPoint(iDE,deId[iDE],fParameters[3*iDE+0]);
	411	fMSDEx->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+0));
	412	fMSDEy->SetPoint(iDE,deId[iDE],fParameters[3*iDE+1]);
	413	fMSDEy->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+1));
	414	fMSDEp->SetPoint(iDE,deId[iDE],fParameters[3*iDE+2]);
	415	fMSDEp->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+2));
	416	fMSDEz->SetPoint(iDE,deId[iDE],fParameters[3*iDE+3]);
	417	fMSDEz->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+3));
4d610fd5	418	}
	419
	420	// Post final data. Write histo list to a file with option "RECREATE"
	421	PostData(0,fList);
	422
4d610fd5	423	// Re Align
	424	AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true);
	425	newTransform->WriteTransformations("transform2ReAlign.dat");
	426
	427	// Generate realigned data in local cdb
	428	const TClonesArray* array = newTransform->GetMisAlignmentData();
	429
	430	// 100 mum residual resolution for chamber misalignments?
	431	fAlign->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
	432
	433	// CDB manager
	434	AliCDBManager* cdbManager = AliCDBManager::Instance();
cd8521dd	435	cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
cd8521dd	436	cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.Data());
4d610fd5	437
	438	AliCDBMetaData* cdbData = new AliCDBMetaData();
	439	cdbData->SetResponsible("Dimuon Offline project");
	440	cdbData->SetComment("MUON alignment objects with residual misalignment");
	441	AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
	442	cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
	443
	444	}
	445
2df5c2cf	446	//________________________________________________________________________
	447	void AliMUONAlignmentTask::Terminate(const Option_t*)
	448	{
	449	/// Called once per task on the client machine at the end of the analysis.
	450
	451	}
	452
cd8521dd	453	//-----------------------------------------------------------------------
	454	void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
	455	{
	456	/// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
2df5c2cf	457
cd8521dd	458	// Load mapping
	459	AliCDBManager* man = AliCDBManager::Instance();
	460	man->SetDefaultStorage(defaultOCDB);
	461	man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
	462	man->Print();
2df5c2cf	463	man->SetRun(fLastRunNumber);
cd8521dd	464	if ( ! AliMpCDB::LoadDDLStore() ) {
2df5c2cf	465	Error("MUONRefit","Could not access mapping from OCDB !");
cd8521dd	466	exit(-1);
	467	}
	468
2df5c2cf	469	// Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
	470	if (!gGeoManager) {
	471	AliGeomManager::LoadGeometry(geoFilename);
	472	if (!gGeoManager) {
	473	Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
	474	return;
	475	}
	476	}
	477
cd8521dd	478	// set mag field
	479	if (!TGeoGlobalMagField::Instance()->GetField()) {
	480	printf("Loading field map...\n");
	481	AliGRPManager *grpMan = new AliGRPManager();
	482	grpMan->ReadGRPEntry();
	483	grpMan->SetMagField();
	484	delete grpMan;
	485	}
	486	// set the magnetic field for track extrapolations
	487	AliMUONTrackExtrap::SetField();
	488
	489	}