[u/mrichter/AliRoot.git] / MUON / MUONAlignment.C

/**************************************************************************
 * 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$ */

/// \ingroup macros
/// \file MUONAlignment.C
/// \brief Macro for MUON alignment using physics tracks. 
///
/// The macro uses the AliMUONAlignment class to calculate the alignment parameters.
/// An array for the alignment parameters is created and can be filled with
/// initial values that will be used as starting values by the alignment
/// algorithm.
///
/// By default the macro run over galice.root in the working directory. If a file list
/// of galice.root is provided as third argument the macro will run over all files.
/// The macro loop over the files, events and tracks. For each track
/// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then
/// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t
/// lSingleFit) are called. After all tracks have been procesed and fitted
/// AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
/// is called. The array parameters contains the obatained misalignement parameters.
/// A realigned geometry is generated in a local CDB.
///
/// \author: B. Becker and J. Castillo

#if !defined(__CINT__) || defined(__MAKECINT__)

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

#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliMagF.h"
#include "AliTracker.h"
#include "AliGRPManager.h"
#include "AliCDBManager.h"
#include "AliCDBMetaData.h"
#include "AliCDBId.h"
#include "AliGeomManager.h"

#include <TString.h>
#include <TError.h>
#include <TH1.h>
#include <TGraphErrors.h>
#include <TFile.h>
#include <TTree.h>
#include <TClonesArray.h>
#include <Riostream.h>

#include <fstream>

#endif

void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString esdFileName = "AliESDs.root", TString fileList = "")
{
 
  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
  if ( ! AliGeomManager::GetGeometry() ) {
    AliGeomManager::LoadGeometry(geoFilename);
    if (! AliGeomManager::GetGeometry() ) {
      Error("MUONAlignment", "getting geometry from file %s failed", geoFilename);
      return;
    }
  }
  
  // CDB manager
  AliCDBManager* cdbManager = AliCDBManager::Instance();
  cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  cdbManager->SetRun(0);

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

  Double_t parameters[4*156];
  Double_t errors[4*156];
  Double_t pulls[4*156];
  for(Int_t k=0;k<4*156;k++) {
    parameters[k]=0.;
    errors[k]=0.;
    pulls[k]=0.;
  }

  Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};

  // Set initial values here, good guess may help convergence
  // St 1 
  //  Int_t iPar = 0;
  //  parameters[iPar++] =  0.010300 ;  parameters[iPar++] =  0.010600 ;  parameters[iPar++] =  0.000396 ;  

  bool bLoop = kFALSE;
  ifstream sFileList;
  if (fileList.Contains(".list")) {
    cout << "Reading file list: " << fileList.Data() << endl;
    bLoop = kTRUE;

    TString fullListName(".");
    fullListName +="/";
    fullListName +=fileList;
    sFileList.open(fileList.Data());
  }
  
  TH1F *fInvBenMom = new TH1F("fInvBenMom","fInvBenMom",200,-0.1,0.1); 
  TH1F *fBenMom = new TH1F("fBenMom","fBenMom",200,-40,40); 

  AliMUONAlignment* alig = new AliMUONAlignment();
  alig->InitGlobalParameters(parameters);

  AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
  transform->LoadGeometryData();
  alig->SetGeometryTransformer(transform);

  // Do alignment with magnetic field off
  alig->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,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};

  // Set degrees of freedom to align (see AliMUONAlignment)
  alig->AllowVariations(bChOnOff);

  // Fix parameters or add constraints here
//   for (Int_t iSt=0; iSt<5; iSt++)
//     if (!bStOnOff[iSt]) alig->FixStation(iSt+1);
  for (Int_t iCh=0; iCh<10; iCh++)
    if (!bChOnOff[iCh]) alig->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};
  alig->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);

  alig->SetChOnOff(bChOnOff);
  alig->SetSpecLROnOff(bChOnOff);

  // Here we can fix some detection elements
  alig->FixDetElem(908);
  alig->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};
  //  alig->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);


  char cFileName[100];  
    
  Int_t lMaxFile = 1000;
  Int_t iFile = 0;
  Int_t iEvent = 0;
  bool bKeepLoop = kTRUE;
  Int_t iTrackTot=0;
  Int_t iTrackOk=0;

  while(bKeepLoop && iFile<lMaxFile){
    iFile++;
    if (bLoop) {
      sFileList.getline(cFileName,100);
      if (sFileList.eof()) bKeepLoop = kFALSE;
    }
    else {
      sprintf(cFileName,esdFileName.Data());
      bKeepLoop = kFALSE;
    }
    if (!strstr(cFileName,"AliESDs")) continue;      
    cout << "Using file: " << cFileName << endl;
    
    // load ESD event
    TFile* esdFile = TFile::Open(cFileName); // open the file
    if (!esdFile || !esdFile->IsOpen()) {
      cout << "opening ESD file " << cFileName << "failed" << endl;
      continue;
    }
    TTree* esdTree = (TTree*) esdFile->Get("esdTree"); // get the tree
    if (!esdTree) {
      cout << "no ESD tree found" << endl;
      esdFile->Close();
      continue;
    }
    AliESDEvent* esdEvent = new AliESDEvent(); // link ESD event to the tree
    esdEvent->ReadFromTree(esdTree);

    Int_t nevents = esdTree->GetEntries();
    cout << "... with " << nevents << endl;
    for(Int_t event = 0; event < nevents; event++) {
      if (iEvent >= nEvents){
	bKeepLoop = kFALSE;
	break;
      }
      iEvent++;

      if (esdTree->GetEvent(event) <= 0) {
	cout << "fails to read ESD object for event " << event << endl;
	continue;
      }

      Int_t nTracks = Int_t(esdEvent->GetNumberOfMuonTracks());
      if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
      for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
	AliESDMuonTrack* esdTrack = esdEvent->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);
	  alig->ProcessTrack(&track);
	  alig->LocalFit(iTrackOk++,trackParams,0);
	}
	iTrackTot++;
      }
    }
    delete esdEvent;
    esdFile->Close();
    cout << "Processed " << iTrackTot << " Tracks so far." << endl;
  }
  alig->GlobalFit(parameters,errors,pulls);

  cout << "Done with GlobalFit " << endl;

  // 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]=parameters[4*iDE+0];
    MSDEy[iDE]=parameters[4*iDE+1];
    MSDEz[iDE]=parameters[4*iDE+3];
    MSDEp[iDE]=parameters[4*iDE+2];
    MSDExErr[iDE]=(Double_t)alig->GetParError(4*iDE+0);
    MSDEyErr[iDE]=(Double_t)alig->GetParError(4*iDE+1);
    MSDEzErr[iDE]=(Double_t)alig->GetParError(4*iDE+3);
    MSDEpErr[iDE]=(Double_t)alig->GetParError(4*iDE+2);
  }

  cout << "Let's create graphs ...  " << endl;

  TGraphErrors *gMSDEx = new TGraphErrors(lNDetElem,DEid,MSDEx,DEidErr,MSDExErr); 
  TGraphErrors *gMSDEy = new TGraphErrors(lNDetElem,DEid,MSDEy,DEidErr,MSDEyErr); 
  TGraphErrors *gMSDEz = new TGraphErrors(lNDetElem,DEid,MSDEz,DEidErr,MSDEzErr); 
  TGraphErrors *gMSDEp = new TGraphErrors(lNDetElem,DEid,MSDEp,DEidErr,MSDEpErr); 

  cout << "... graphs created, open file ...  " << endl;

  TFile *hFile = new TFile("measShifts.root","RECREATE");

  cout << "... file opened ...  " << endl;

  gMSDEx->Write("gMSDEx");
  gMSDEy->Write("gMSDEy");
  gMSDEz->Write("gMSDEz");
  gMSDEp->Write("gMSDEp");
  fInvBenMom->Write();
  fBenMom->Write();
  hFile->Close();
  
  cout << "... and closed!" << endl;
  // Re Align
  AliMUONGeometryTransformer *newTransform = alig->ReAlign(transform,parameters,true); 
  newTransform->WriteTransformations("transform2ReAlign.dat");
  
  // Generate realigned data in local cdb
  const TClonesArray* array = newTransform->GetMisAlignmentData();

  // 100 mum residual resolution for chamber misalignments?
  alig->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);

  cdbManager->SetSpecificStorage("MUON/Align/Data","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
af9df62e	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
e54bf126	18	/// \ingroup macros
	19	/// \file MUONAlignment.C
	20	/// \brief Macro for MUON alignment using physics tracks.
	21	///
	22	/// The macro uses the AliMUONAlignment class to calculate the alignment parameters.
	23	/// An array for the alignment parameters is created and can be filled with
	24	/// initial values that will be used as starting values by the alignment
	25	/// algorithm.
	26	///
	27	/// By default the macro run over galice.root in the working directory. If a file list
	28	/// of galice.root is provided as third argument the macro will run over all files.
	29	/// The macro loop over the files, events and tracks. For each track
	30	/// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then
	31	/// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t
	32	/// lSingleFit) are called. After all tracks have been procesed and fitted
	33	/// AliMUONAlignment::GlobalFit(Double_t parameters,Double_t errors,Double_t *pulls)
	34	/// is called. The array parameters contains the obatained misalignement parameters.
	35	/// A realigned geometry is generated in a local CDB.
	36	///
	37	/// \author: B. Becker and J. Castillo
af9df62e	38
2688ea9b	39	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	40
	41	#include "AliMUONAlignment.h"
	42	#include "AliMUONTrack.h"
	43	#include "AliMUONTrackExtrap.h"
	44	#include "AliMUONTrackParam.h"
	45	#include "AliMUONGeometryTransformer.h"
103e6575	46	#include "AliMUONESDInterface.h"
103e6575	47
61fed964	48	#include "AliESDEvent.h"
61fed964	49	#include "AliESDMuonTrack.h"
f7a1cc68	50	#include "AliMagF.h"
2688ea9b	51	#include "AliTracker.h"
cd8521dd	52	#include "AliGRPManager.h"
2688ea9b	53	#include "AliCDBManager.h"
	54	#include "AliCDBMetaData.h"
	55	#include "AliCDBId.h"
a02fe82d	56	#include "AliGeomManager.h"
2688ea9b	57
2688ea9b	58	#include <TString.h>
2688ea9b	59	#include <TError.h>
	60	#include <TH1.h>
	61	#include <TGraphErrors.h>
	62	#include <TFile.h>
61fed964	63	#include <TTree.h>
2688ea9b	64	#include <TClonesArray.h>
	65	#include <Riostream.h>
	66
	67	#include <fstream>
	68
	69	#endif
	70
61fed964	71	void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString esdFileName = "AliESDs.root", TString fileList = "")
af9df62e	72	{
af9df62e	73
59830a87	74	// Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
a02fe82d	75	if ( ! AliGeomManager::GetGeometry() ) {
	76	AliGeomManager::LoadGeometry(geoFilename);
	77	if (! AliGeomManager::GetGeometry() ) {
	78	Error("MUONAlignment", "getting geometry from file %s failed", geoFilename);
59830a87	79	return;
	80	}
	81	}
	82
cd8521dd	83	// CDB manager
	84	AliCDBManager* cdbManager = AliCDBManager::Instance();
	85	cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
	86	cdbManager->SetRun(0);
	87
59830a87	88	// set mag field
f7a1cc68	89	if (!TGeoGlobalMagField::Instance()->GetField()) {
f7a1cc68	90	printf("Loading field map...\n");
cd8521dd	91	AliGRPManager *grpMan = new AliGRPManager();
	92	grpMan->ReadGRPEntry();
	93	grpMan->SetMagField();
	94	delete grpMan;
f7a1cc68	95	}
59830a87	96	// set the magnetic field for track extrapolations
f7a1cc68	97	AliMUONTrackExtrap::SetField();
59830a87	98
723b0b5b	99	Double_t parameters[4*156];
	100	Double_t errors[4*156];
	101	Double_t pulls[4*156];
	102	for(Int_t k=0;k<4*156;k++) {
af9df62e	103	parameters[k]=0.;
	104	errors[k]=0.;
	105	pulls[k]=0.;
	106	}
	107
	108	Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
	109
af9df62e	110	// Set initial values here, good guess may help convergence
af9df62e	111	// St 1
2688ea9b	112	// Int_t iPar = 0;
2688ea9b	113	// parameters[iPar++] = 0.010300 ; parameters[iPar++] = 0.010600 ; parameters[iPar++] = 0.000396 ;
af9df62e	114
	115	bool bLoop = kFALSE;
	116	ifstream sFileList;
	117	if (fileList.Contains(".list")) {
	118	cout << "Reading file list: " << fileList.Data() << endl;
	119	bLoop = kTRUE;
	120
	121	TString fullListName(".");
	122	fullListName +="/";
	123	fullListName +=fileList;
	124	sFileList.open(fileList.Data());
	125	}
	126
	127	TH1F *fInvBenMom = new TH1F("fInvBenMom","fInvBenMom",200,-0.1,0.1);
	128	TH1F *fBenMom = new TH1F("fBenMom","fBenMom",200,-40,40);
	129
	130	AliMUONAlignment* alig = new AliMUONAlignment();
	131	alig->InitGlobalParameters(parameters);
	132
a02fe82d	133	AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
a02fe82d	134	transform->LoadGeometryData();
af9df62e	135	alig->SetGeometryTransformer(transform);
7f8c1308	136
	137	// Do alignment with magnetic field off
	138	alig->SetBFieldOn(kFALSE);
587965bd	139
	140	// Set tracking station to use
	141	Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
d8ad38e5	142	Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
587965bd	143
f3e4edfb	144	// Set degrees of freedom to align (see AliMUONAlignment)
	145	alig->AllowVariations(bChOnOff);
	146
587965bd	147	// Fix parameters or add constraints here
d8ad38e5	148	// for (Int_t iSt=0; iSt<5; iSt++)
	149	// if (!bStOnOff[iSt]) alig->FixStation(iSt+1);
	150	for (Int_t iCh=0; iCh<10; iCh++)
	151	if (!bChOnOff[iCh]) alig->FixChamber(iCh+1);
587965bd	152
	153	// Left and right sides of the detector are independent, one can choose to align
	154	// only one side
	155	Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
d8ad38e5	156	alig->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
	157
	158	alig->SetChOnOff(bChOnOff);
	159	alig->SetSpecLROnOff(bChOnOff);
587965bd	160
723b0b5b	161	// Here we can fix some detection elements
	162	alig->FixDetElem(908);
	163	alig->FixDetElem(1020);
	164
587965bd	165	// Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
	166	Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
	167	Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
723b0b5b	168	// alig->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
587965bd	169
af9df62e	170
af9df62e	171	char cFileName[100];
7f8c1308	172
af9df62e	173	Int_t lMaxFile = 1000;
af9df62e	174	Int_t iFile = 0;
587965bd	175	Int_t iEvent = 0;
af9df62e	176	bool bKeepLoop = kTRUE;
587965bd	177	Int_t iTrackTot=0;
	178	Int_t iTrackOk=0;
	179
af9df62e	180	while(bKeepLoop && iFile<lMaxFile){
	181	iFile++;
	182	if (bLoop) {
	183	sFileList.getline(cFileName,100);
	184	if (sFileList.eof()) bKeepLoop = kFALSE;
	185	}
	186	else {
61fed964	187	sprintf(cFileName,esdFileName.Data());
af9df62e	188	bKeepLoop = kFALSE;
af9df62e	189	}
723b0b5b	190	if (!strstr(cFileName,"AliESDs")) continue;
af9df62e	191	cout << "Using file: " << cFileName << endl;
61fed964	192
	193	// load ESD event
	194	TFile* esdFile = TFile::Open(cFileName); // open the file
	195	if (!esdFile \|\| !esdFile->IsOpen()) {
	196	cout << "opening ESD file " << cFileName << "failed" << endl;
	197	continue;
	198	}
	199	TTree* esdTree = (TTree*) esdFile->Get("esdTree"); // get the tree
	200	if (!esdTree) {
	201	cout << "no ESD tree found" << endl;
	202	esdFile->Close();
	203	continue;
	204	}
	205	AliESDEvent* esdEvent = new AliESDEvent(); // link ESD event to the tree
	206	esdEvent->ReadFromTree(esdTree);
7f8c1308	207
61fed964	208	Int_t nevents = esdTree->GetEntries();
587965bd	209	cout << "... with " << nevents << endl;
af9df62e	210	for(Int_t event = 0; event < nevents; event++) {
587965bd	211	if (iEvent >= nEvents){
	212	bKeepLoop = kFALSE;
	213	break;
	214	}
	215	iEvent++;
	216
61fed964	217	if (esdTree->GetEvent(event) <= 0) {
	218	cout << "fails to read ESD object for event " << event << endl;
	219	continue;
	220	}
7f8c1308	221
61fed964	222	Int_t nTracks = Int_t(esdEvent->GetNumberOfMuonTracks());
	223	if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
	224	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
	225	AliESDMuonTrack* esdTrack = esdEvent->GetMuonTrack(iTrack);
	226	if (!esdTrack->ClustersStored()) continue;
	227	Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
af9df62e	228	fInvBenMom->Fill(invBenMom);
	229	fBenMom->Fill(1./invBenMom);
	230	if (TMath::Abs(invBenMom)<=1.04) {
103e6575	231	AliMUONTrack track;
ad3c6eda	232	AliMUONESDInterface::ESDToMUON(*esdTrack, track);
61fed964	233	alig->ProcessTrack(&track);
af9df62e	234	alig->LocalFit(iTrackOk++,trackParams,0);
af9df62e	235	}
587965bd	236	iTrackTot++;
af9df62e	237	}
af9df62e	238	}
61fed964	239	delete esdEvent;
61fed964	240	esdFile->Close();
587965bd	241	cout << "Processed " << iTrackTot << " Tracks so far." << endl;
af9df62e	242	}
	243	alig->GlobalFit(parameters,errors,pulls);
	244
	245	cout << "Done with GlobalFit " << endl;
	246
	247	// Store results
	248	Double_t DEid[156] = {0};
	249	Double_t MSDEx[156] = {0};
	250	Double_t MSDEy[156] = {0};
723b0b5b	251	Double_t MSDEz[156] = {0};
723b0b5b	252	Double_t MSDEp[156] = {0};
af9df62e	253	Double_t DEidErr[156] = {0};
	254	Double_t MSDExErr[156] = {0};
	255	Double_t MSDEyErr[156] = {0};
723b0b5b	256	Double_t MSDEzErr[156] = {0};
723b0b5b	257	Double_t MSDEpErr[156] = {0};
af9df62e	258	Int_t lNDetElem = 42+42+182+262+26*2;
af9df62e	259	Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
2688ea9b	260	// Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
af9df62e	261	Int_t idOffset = 0; // 400
	262	Int_t lSDetElemCh = 0;
	263	for(Int_t iDE=0; iDE<lNDetElem; iDE++){
	264	DEidErr[iDE] = 0.;
	265	DEid[iDE] = idOffset+100;
	266	DEid[iDE] += iDE;
	267	lSDetElemCh = 0;
	268	for(Int_t iCh=0; iCh<9; iCh++){
	269	lSDetElemCh += lNDetElemCh[iCh];
	270	if (iDE>=lSDetElemCh) {
	271	DEid[iDE] += 100;
	272	DEid[iDE] -= lNDetElemCh[iCh];
	273	}
	274	}
723b0b5b	275	MSDEx[iDE]=parameters[4*iDE+0];
	276	MSDEy[iDE]=parameters[4*iDE+1];
	277	MSDEz[iDE]=parameters[4*iDE+3];
	278	MSDEp[iDE]=parameters[4*iDE+2];
	279	MSDExErr[iDE]=(Double_t)alig->GetParError(4*iDE+0);
	280	MSDEyErr[iDE]=(Double_t)alig->GetParError(4*iDE+1);
	281	MSDEzErr[iDE]=(Double_t)alig->GetParError(4*iDE+3);
	282	MSDEpErr[iDE]=(Double_t)alig->GetParError(4*iDE+2);
af9df62e	283	}
	284
	285	cout << "Let's create graphs ... " << endl;
	286
	287	TGraphErrors *gMSDEx = new TGraphErrors(lNDetElem,DEid,MSDEx,DEidErr,MSDExErr);
	288	TGraphErrors *gMSDEy = new TGraphErrors(lNDetElem,DEid,MSDEy,DEidErr,MSDEyErr);
723b0b5b	289	TGraphErrors *gMSDEz = new TGraphErrors(lNDetElem,DEid,MSDEz,DEidErr,MSDEzErr);
723b0b5b	290	TGraphErrors *gMSDEp = new TGraphErrors(lNDetElem,DEid,MSDEp,DEidErr,MSDEpErr);
af9df62e	291
	292	cout << "... graphs created, open file ... " << endl;
	293
	294	TFile *hFile = new TFile("measShifts.root","RECREATE");
	295
	296	cout << "... file opened ... " << endl;
	297
	298	gMSDEx->Write("gMSDEx");
	299	gMSDEy->Write("gMSDEy");
723b0b5b	300	gMSDEz->Write("gMSDEz");
723b0b5b	301	gMSDEp->Write("gMSDEp");
af9df62e	302	fInvBenMom->Write();
	303	fBenMom->Write();
	304	hFile->Close();
	305
	306	cout << "... and closed!" << endl;
	307	// Re Align
	308	AliMUONGeometryTransformer *newTransform = alig->ReAlign(transform,parameters,true);
	309	newTransform->WriteTransformations("transform2ReAlign.dat");
	310
	311	// Generate realigned data in local cdb
2688ea9b	312	const TClonesArray* array = newTransform->GetMisAlignmentData();
4818a9b7	313
	314	// 100 mum residual resolution for chamber misalignments?
	315	alig->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
cd8521dd	316
cd8521dd	317	cdbManager->SetSpecificStorage("MUON/Align/Data","local://ReAlignCDB");
af9df62e	318
	319	AliCDBMetaData* cdbData = new AliCDBMetaData();
	320	cdbData->SetResponsible("Dimuon Offline project");
	321	cdbData->SetComment("MUON alignment objects with residual misalignment");
d8ad38e5	322	AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
2688ea9b	323	cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
af9df62e	324
61fed964	325	}
61fed964	326