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

// ---
// Macro for MUON alignment using physics tracks. The macro uses 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.
// Ny 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.
//
// Authors: B. Becker and J. Castillo
// ---

void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString fileName = "galice.root", TString fileList = "")
{
 
  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
  if (!gGeoManager) {
    TGeoManager::Import(geoFilename);
    if (!gGeoManager) {
      Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
      return;
    }
  }
  
  // set  mag field 
  // waiting for mag field in CDB 
  printf("Loading field map...\n");
  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
  AliTracker::SetFieldMap(field, kFALSE);
  // set the magnetic field for track extrapolations
  AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());

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

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

  Int_t iPar = 0;
  // Set initial values here, good guess may help convergence
  // St 1 
  //   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(true);
  transform->ReadGeometryData("volpath.dat", gGeoManager);
  alig->SetGeometryTransformer(transform);

  char cFileName[100];  
  AliMUONDataInterface amdi;

  Int_t lMaxFile = 1000;
  Int_t iFile = 0;
  bool bKeepLoop = kTRUE;
  while(bKeepLoop && iFile<lMaxFile){
    iFile++;
    if (bLoop) {
      sFileList.getline(cFileName,100);
      if (sFileList.eof()) bKeepLoop = kFALSE;
    }
    else {
      sprintf(cFileName,fileName.Data());
      bKeepLoop = kFALSE;
    }
    if (!strstr(cFileName,"galice.root")) continue;      
    cout << "Using file: " << cFileName << endl;
    amdi.SetFile(cFileName);
    Int_t nevents = amdi.NumberOfEvents();
    for(Int_t event = 0; event < nevents; event++) {
      amdi.GetEvent(event);
      Int_t ntracks = amdi.NumberOfRecTracks();
      cout << " there are " << ntracks << " tracks in event " << event << endl;
      Int_t iTrack=0;
      Int_t iTrackOk=0;
      AliMUONTrack* track = amdi.RecTrack(iTrack);
      while(track) {   
	AliMUONTrackParam trackParam(*((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First())));
	AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.);
	Double_t invBenMom = trackParam.GetInverseBendingMomentum();
	fInvBenMom->Fill(invBenMom);
	fBenMom->Fill(1./invBenMom);
	if (TMath::Abs(invBenMom)<=1.04) {
	  cout << "Track " << iTrack << endl;
	  alig->ProcessTrack(track);
	  cout << "Calling LocalFit" << endl;
	  alig->LocalFit(iTrackOk++,trackParams,0);
	}
	iTrack++;
	track = amdi.RecTrack(iTrack);
      }
    }
  }
  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 MSDExt[156] = {0};
  Double_t MSDEyt[156] = {0};
  Double_t DEidErr[156] = {0};
  Double_t MSDExErr[156] = {0};
  Double_t MSDEyErr[156] = {0};
  Double_t MSDExtErr[156] = {0};
  Double_t MSDEytErr[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[3*iDE+0];
    MSDEy[iDE]=parameters[3*iDE+1];
    MSDExt[iDE]=parameters[3*iDE+2];
    MSDEyt[iDE]=parameters[3*iDE+2];
    MSDExErr[iDE]=(Double_t)alig->GetParError(3*iDE+0);
    MSDEyErr[iDE]=(Double_t)alig->GetParError(3*iDE+1);
    MSDExtErr[iDE]=(Double_t)alig->GetParError(3*iDE+2);
    MSDEytErr[iDE]=(Double_t)alig->GetParError(3*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 *gMSDExt = new TGraphErrors(lNDetElem,DEid,MSDExt,DEidErr,MSDExtErr); 
  TGraphErrors *gMSDEyt = new TGraphErrors(lNDetElem,DEid,MSDEyt,DEidErr,MSDEytErr); 

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

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

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

  gMSDEx->Write("gMSDEx");
  gMSDEy->Write("gMSDEy");
  gMSDExt->Write("gMSDExt");
  gMSDEyt->Write("gMSDEyt");
  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
  TClonesArray* array = newTransform->GetMisAlignmentData();
   
  // 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, 0); 
  cdbManager->Put(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
	18	// ---
	19	// Macro for MUON alignment using physics tracks. The macro uses AliMUONAlignment
	20	// class to calculate the alignment parameters.
	21	// An array for the alignment parameters is created and can be filled with
	22	// initial values that will be used as starting values by the alignment
	23	// algorithm.
	24	// Ny default the macro run over galice.root in the working directory. If a file list
	25	// of galice.root is provided as third argument the macro will run over all files.
	26	// The macro loop over the files, events and tracks. For each track
	27	// AliMUONAlignment::ProcessTrack(AliMUONTrack * track) and then
	28	// AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t
	29	// lSingleFit) are called. After all tracks have been procesed and fitted
	30	// AliMUONAlignment::GlobalFit(Double_t parameters,Double_t errors,Double_t *pulls)
	31	// is called. The array parameters contains the obatained misalignement parameters.
	32	// A realigned geometry is generated in a local CDB.
	33	//
	34	// Authors: B. Becker and J. Castillo
	35	// ---
	36
59830a87	37	void MUONAlignment(Int_t nEvents = 100000, char* geoFilename = "geometry.root", TString fileName = "galice.root", TString fileList = "")
af9df62e	38	{
af9df62e	39
59830a87	40	// Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
	41	if (!gGeoManager) {
	42	TGeoManager::Import(geoFilename);
	43	if (!gGeoManager) {
	44	Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
	45	return;
	46	}
	47	}
	48
	49	// set mag field
	50	// waiting for mag field in CDB
	51	printf("Loading field map...\n");
	52	AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
	53	AliTracker::SetFieldMap(field, kFALSE);
	54	// set the magnetic field for track extrapolations
	55	AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
	56
af9df62e	57	Double_t parameters[3*156];
	58	Double_t errors[3*156];
	59	Double_t pulls[3*156];
	60	for(Int_t k=0;k<3*156;k++) {
	61	parameters[k]=0.;
	62	errors[k]=0.;
	63	pulls[k]=0.;
	64	}
	65
	66	Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
	67
	68	Int_t iPar = 0;
	69	// Set initial values here, good guess may help convergence
	70	// St 1
	71	// parameters[iPar++] = 0.010300 ; parameters[iPar++] = 0.010600 ; parameters[iPar++] = 0.000396 ;
	72
	73	bool bLoop = kFALSE;
	74	ifstream sFileList;
	75	if (fileList.Contains(".list")) {
	76	cout << "Reading file list: " << fileList.Data() << endl;
	77	bLoop = kTRUE;
	78
	79	TString fullListName(".");
	80	fullListName +="/";
	81	fullListName +=fileList;
	82	sFileList.open(fileList.Data());
	83	}
	84
	85	TH1F *fInvBenMom = new TH1F("fInvBenMom","fInvBenMom",200,-0.1,0.1);
	86	TH1F *fBenMom = new TH1F("fBenMom","fBenMom",200,-40,40);
	87
	88	AliMUONAlignment* alig = new AliMUONAlignment();
	89	alig->InitGlobalParameters(parameters);
	90
	91	AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer(true);
59830a87	92	transform->ReadGeometryData("volpath.dat", gGeoManager);
af9df62e	93	alig->SetGeometryTransformer(transform);
	94
	95	char cFileName[100];
	96	AliMUONDataInterface amdi;
	97
	98	Int_t lMaxFile = 1000;
	99	Int_t iFile = 0;
	100	bool bKeepLoop = kTRUE;
	101	while(bKeepLoop && iFile<lMaxFile){
	102	iFile++;
	103	if (bLoop) {
	104	sFileList.getline(cFileName,100);
	105	if (sFileList.eof()) bKeepLoop = kFALSE;
	106	}
	107	else {
	108	sprintf(cFileName,fileName.Data());
	109	bKeepLoop = kFALSE;
	110	}
	111	if (!strstr(cFileName,"galice.root")) continue;
	112	cout << "Using file: " << cFileName << endl;
	113	amdi.SetFile(cFileName);
	114	Int_t nevents = amdi.NumberOfEvents();
	115	for(Int_t event = 0; event < nevents; event++) {
	116	amdi.GetEvent(event);
	117	Int_t ntracks = amdi.NumberOfRecTracks();
	118	cout << " there are " << ntracks << " tracks in event " << event << endl;
	119	Int_t iTrack=0;
	120	Int_t iTrackOk=0;
	121	AliMUONTrack* track = amdi.RecTrack(iTrack);
	122	while(track) {
8cde4af5	123	AliMUONTrackParam trackParam(((AliMUONTrackParam)(track->GetTrackParamAtHit()->First())));
	124	AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.);
	125	Double_t invBenMom = trackParam.GetInverseBendingMomentum();
af9df62e	126	fInvBenMom->Fill(invBenMom);
	127	fBenMom->Fill(1./invBenMom);
	128	if (TMath::Abs(invBenMom)<=1.04) {
	129	cout << "Track " << iTrack << endl;
	130	alig->ProcessTrack(track);
	131	cout << "Calling LocalFit" << endl;
	132	alig->LocalFit(iTrackOk++,trackParams,0);
	133	}
	134	iTrack++;
	135	track = amdi.RecTrack(iTrack);
	136	}
	137	}
	138	}
	139	alig->GlobalFit(parameters,errors,pulls);
	140
	141	cout << "Done with GlobalFit " << endl;
	142
	143	// Store results
	144	Double_t DEid[156] = {0};
	145	Double_t MSDEx[156] = {0};
	146	Double_t MSDEy[156] = {0};
	147	Double_t MSDExt[156] = {0};
	148	Double_t MSDEyt[156] = {0};
	149	Double_t DEidErr[156] = {0};
	150	Double_t MSDExErr[156] = {0};
	151	Double_t MSDEyErr[156] = {0};
	152	Double_t MSDExtErr[156] = {0};
	153	Double_t MSDEytErr[156] = {0};
	154	Int_t lNDetElem = 42+42+182+262+26*2;
	155	Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
	156	Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
	157	Int_t idOffset = 0; // 400
	158	Int_t lSDetElemCh = 0;
	159	for(Int_t iDE=0; iDE<lNDetElem; iDE++){
	160	DEidErr[iDE] = 0.;
	161	DEid[iDE] = idOffset+100;
	162	DEid[iDE] += iDE;
	163	lSDetElemCh = 0;
	164	for(Int_t iCh=0; iCh<9; iCh++){
	165	lSDetElemCh += lNDetElemCh[iCh];
	166	if (iDE>=lSDetElemCh) {
	167	DEid[iDE] += 100;
	168	DEid[iDE] -= lNDetElemCh[iCh];
	169	}
	170	}
	171	MSDEx[iDE]=parameters[3*iDE+0];
	172	MSDEy[iDE]=parameters[3*iDE+1];
	173	MSDExt[iDE]=parameters[3*iDE+2];
	174	MSDEyt[iDE]=parameters[3*iDE+2];
	175	MSDExErr[iDE]=(Double_t)alig->GetParError(3*iDE+0);
	176	MSDEyErr[iDE]=(Double_t)alig->GetParError(3*iDE+1);
	177	MSDExtErr[iDE]=(Double_t)alig->GetParError(3*iDE+2);
	178	MSDEytErr[iDE]=(Double_t)alig->GetParError(3*iDE+2);
	179	}
	180
	181	cout << "Let's create graphs ... " << endl;
	182
	183	TGraphErrors *gMSDEx = new TGraphErrors(lNDetElem,DEid,MSDEx,DEidErr,MSDExErr);
	184	TGraphErrors *gMSDEy = new TGraphErrors(lNDetElem,DEid,MSDEy,DEidErr,MSDEyErr);
	185	TGraphErrors *gMSDExt = new TGraphErrors(lNDetElem,DEid,MSDExt,DEidErr,MSDExtErr);
	186	TGraphErrors *gMSDEyt = new TGraphErrors(lNDetElem,DEid,MSDEyt,DEidErr,MSDEytErr);
	187
	188	cout << "... graphs created, open file ... " << endl;
	189
190	TFile *hFile = new TFile("measShifts.root","RECREATE");
191
192	cout << "... file opened ... " << endl;
193
194	gMSDEx->Write("gMSDEx");
195	gMSDEy->Write("gMSDEy");
196	gMSDExt->Write("gMSDExt");
197	gMSDEyt->Write("gMSDEyt");
198	fInvBenMom->Write();
199	fBenMom->Write();
200	hFile->Close();
201
202	cout << "... and closed!" << endl;
203	// Re Align
204	AliMUONGeometryTransformer *newTransform = alig->ReAlign(transform,parameters,true);
205	newTransform->WriteTransformations("transform2ReAlign.dat");
206
207	// Generate realigned data in local cdb
208	TClonesArray* array = newTransform->GetMisAlignmentData();
209
210	// CDB manager
211	AliCDBManager* cdbManager = AliCDBManager::Instance();
212	cdbManager->SetDefaultStorage("local://ReAlignCDB");
213
214	AliCDBMetaData* cdbData = new AliCDBMetaData();
215	cdbData->SetResponsible("Dimuon Offline project");
216	cdbData->SetComment("MUON alignment objects with residual misalignment");
217	AliCDBId id("MUON/Align/Data", 0, 0);
218	cdbManager->Put(array, id, cdbData);
219
220	}