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

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

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliMagF.h"
#include "AliCDBManager.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) 
  : AliAnalysisTask(name, ""),
    fESD(0x0),
    fAlign(0x0),
    fGeoFilename(geofilename),
    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 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)
{
  /// 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) 
{
  /// Assignment
  AliAnalysisTask::operator=(other);
  fESD = other.fESD;
  fAlign = other.fAlign;
  fGeoFilename = other.fGeoFilename;
  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;
}

//________________________________________________________________________
void AliMUONAlignmentTask::LocalInit() 
{
  /// 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);
  }
  // set the magnetic field for track extrapolations
  AliMUONTrackExtrap::SetField();

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


  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++;
  }
  
  // Post final data. Write histo list to a file with option "RECREATE"
  PostData(0,fList);
}      

//________________________________________________________________________
void AliMUONAlignmentTask::Terminate(const Option_t*)
{
  /// 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("local://ReAlignCDB");
  
  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);

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