[u/mrichter/AliRoot.git] / HMPID / HMPIDFindAlign.C

#include <TSystem.h>
#include <TVector3.h>
#include <TVector2.h>
#include <TNtuple.h>
#include <TFile.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TPaveText.h>
#include <TText.h>
#include <TF1.h>
#include <AliHMPIDParam.h>
#include <AliGeomManager.h>
#include <AliCDBEntry.h>
#include <TClonesArray.h>

TVector3 AlignAxis(Int_t ch,Int_t axis);
Double_t FindChCenterD(Int_t ch);
void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist);
TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc);
void help();

//--------------------------------------------------
void HMPIDFindAlign(Int_t cham=-1)
{
  Printf("");
  Printf("******************************************");
  Printf("* You can type help() for instructions...*");
  Printf("******************************************");
  Printf("");
  
  gSystem->Sleep(3000);
  
  AliGeomManager::LoadGeometry("geometry.root");

  TFile *file = TFile::Open("real_geometry.root");
  if(file) {
    AliCDBEntry *entry = (AliCDBEntry*)file->Get("AliCDBEntry");
    TClonesArray *array = (TClonesArray*)entry->GetObject();
    AliGeomManager::ApplyAlignObjsToGeom(*array);
}      
  AliHMPIDParam::Instance(); //just to print at the beginning...

  Int_t ichMin = 0;
  Int_t ichMax = 6;         
  if(cham>=0&&cham<=6) { ichMin = cham; ichMax = cham;}
  
  for(Int_t ich=ichMin;ich<=ichMax;ich++) {
    
    Printf("");
    Printf("*********SUMMARY FOR CHAMBER %i **********",ich);
    Printf ("Distance from IP of chamber %i = %f (cm)",ich,FindChCenterD(ich));
    
    TVector3 xAl = AlignAxis(ich,0);
    TVector3 yAl = AlignAxis(ich,1);
    
    Double_t p1 = 1/(xAl.Z()*xAl.Z());
    Double_t p2 = 1/(yAl.Z()*yAl.Z());
    
    Double_t distMean = (p1*xAl.Y()+p2*yAl.Y())/(p1+p2);
  
    FindAligners(ich,xAl.X(),yAl.X(),distMean); 
  }
}      
//--------------------------------------------------
TVector3 AlignAxis(Int_t ch,Int_t axis)
{
  
// Axis = 0  X Local Axis;
// Axis = 1  Y Local Axis;
  
  
  TH1F* delta[6];
  
  TFile* file = TFile::Open("histos.root","old"); 
  TNtuple *pNtuple;
  if(!(TNtuple*)file->Get("nTupla")) {
// new version because Giacomo uses now a Task to strio the data...
  TList *hmpoutput = (TList*)(file->FindObjectAny("hmpoutput"));
  pNtuple = (TNtuple*)hmpoutput->At(0);
} else {
  pNtuple = (TNtuple*)file->Get("nTupla");
}
  
  for(Int_t j=0;j<6;j++) {
    if(axis == 0) {
      delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
    } else {
      delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
    }
  }
  Double_t mean[6];
  Double_t sigma[6];
  Double_t siz[6];
  Double_t bin[6];
  
  TCanvas *c = new TCanvas();

  Int_t index = 0;  
  for(Int_t j=0;j<6;j++) {
    Double_t xMin = j*20;
    Double_t xMax = xMin+20;
    if(axis == 0) {
      pNtuple->Draw(Form("Xtrk-Xmip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Xmip>%f&&Xmip<%f",ch,xMin,xMax),"0");
    } else {
      pNtuple->Draw(Form("Ytrk-Ymip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Ymip>%f&&Ymip<%f",ch,xMin,xMax),"0");
    }
    if(delta[j]->GetEntries()<50)  continue;

    Double_t dIntFit = 5.;
    siz[index] = 0.5*(xMax-xMin);
    bin[index]= 0.5*(xMax+xMin);
    Double_t xMaxPos = delta[j]->GetBinCenter(delta[j]->FindFirstBinAbove(delta[j]->GetMaximum()-1));
    delta[j]->Fit("gaus","Q0","",xMaxPos-dIntFit,xMaxPos+dIntFit);
    mean[index]  = (delta[j]->GetFunction("gaus"))->GetParameter(1);
    sigma[index] = (delta[j]->GetFunction("gaus"))->GetParError(1);
    index++;
  }
  c->Clear();
  c->Divide(3,3);

  for(Int_t j=0;j<6;j++) {
    c->cd(j+1);delta[j]->Draw();
  }
      
  c->cd(8);
  
  TGraphErrors *pGr = new TGraphErrors(index,bin,mean,siz,sigma);
  
  pGr->SetTitle(Form("LRS for chamber %i",ch));
  if(axis == 0) {
    pGr->GetXaxis()->SetTitle("Xmip (cm)");
    pGr->GetYaxis()->SetTitle("deltaX (cm)");
  } else {  
    pGr->GetXaxis()->SetTitle("Ymip (cm)");
    pGr->GetYaxis()->SetTitle("DeltaY (cm)");
  }
  pGr->SetMinimum(-5.);
  pGr->SetMaximum( 5.);
  //pGr->SetMinimum(TMath::MinElement(index,mean)-1.);
  //pGr->SetMaximum(TMath::MaxElement(index,mean)+1.);
  pGr->Draw("ALP");
  pGr->Fit("pol1","Q");  
  Double_t  shift      = (pGr->GetFunction("pol1"))->GetParameter(0);
  Double_t  errshift   = (pGr->GetFunction("pol1"))->GetParError(0);
  Double_t  coefang    = (pGr->GetFunction("pol1"))->GetParameter(1);
  Double_t  errcoefang = (pGr->GetFunction("pol1"))->GetParError(1);
  Double_t  prob       = (pGr->GetFunction("pol1"))->GetProb();
  Double_t dist = FindChCenterD(ch);
  
  Double_t rfShift    = -coefang/(1.+coefang)*dist;
  Double_t errrfShift = TMath::Abs(dist/((1.+coefang)*(1+coefang))*errcoefang);
  if (axis == 0) {
    Printf(" x shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
  } else {
    Printf(" y shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
  }

  c->cd(7);
  
  TPaveText *pPa = new TPaveText();
  pPa->DrawPave(0.2,0.2,0.8,0.8);
  TText *t = new TText();
  t->SetTextSize(0.08);
  t->DrawText(0.21,0.70,Form(" shift %5.3f +/- %5.3f",shift,errshift));
  t->DrawText(0.21,0.50,Form(" dist shift %5.3f +/- %5.3f",rfShift,errrfShift));
  t->DrawText(0.21,0.30,Form(" dist from IP %5.3f ",dist));
  
  if (axis == 0) {
    c->SaveAs(Form("ShiftXChamber%i.gif",ch));
  } else {
    c->SaveAs(Form("ShiftYChamber%i.gif",ch));
  }
    
  TVector3 align(shift,rfShift,errrfShift);
  return align;
}
//------------------------------------------------------------------
Double_t FindChCenterD(Int_t ch)
{
  
  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
  TVector3 xyz = pParam->Lors2Mars(ch,0.5*pParam->SizeAllX(),0.5*pParam->SizeAllY());
  return xyz.Mag();
}
//------------------------------------------------------------------
void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist) 
{
  TVector3 old = FindNewCenter(ch,0,0,0);
    
  Printf("");
  Printf("Old X Y Z of the chamber n. %i",ch);
  Printf("XC = %f",old.X());
  Printf("YC = %f",old.Y());
  Printf("ZC = %f",old.Z());
  Printf("Dist. from IP = %f",old.Mag());
  Printf("");
  

  Printf("Local X shift = %f",xShift);
  Printf("Local Y shift = %f",yShift);
  Printf("RPhi    shift = %f",dist);
    
  TVector3 shift = FindNewCenter(ch,xShift,yShift,dist);

  Printf("");
  Printf("New X Y Z of the chamber n. %i",ch);
  Printf("");
  Printf("XC = %f YC = %f ZC = %f",shift.X(),shift.Y(),shift.Z());
  Printf("Dist. from IP = %f",shift.Mag());
  Printf("");
  Printf("(next line are the values to be inserted in the Align Object)");
  Printf("");
  Printf("XC shift = %5.2f YC shift = %5.2f ZC shift = %5.2f  for chamber %i",shift.X()-old.X(),shift.Y()-old.Y(),shift.Z()-old.Z(),ch);
  Printf("");
}  
//------------------------------------------------------------------
TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc) 
{
//First shift in (LocX,LocY) and go in MRS
//Then shift radially...
  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
  TVector3 xyz = pParam->Lors2Mars(ch,0.5*pParam->SizeAllX()+xShLoc,0.5*pParam->SizeAllY()+yShLoc,-1);
  TVector3 xyzSh(1,1,1);
  xyzSh.SetTheta(xyz.Theta());
  xyzSh.SetPhi(xyz.Phi());
  xyzSh.SetMag(xyz.Mag()+zShLoc);
  return xyzSh;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void help()
{
  Printf("Instructions to run this macros:");
  Printf("    1 - In the same dir where the macro runs it must exist the file geometry.root");
  Printf("    2 - assign symbolic links as it follows:");
  Printf("       a -  file from OCDB/HMPID/Align/Data linked to real_geometry.root");
  Printf("       b -  file with the NTupla of ESD Data (produced from AliHMPIDAnalysisTask.C) linked as histos.root");
  Printf("    ");
  Printf("   The macro have to be run in the following way:");
  Printf("    ");
  Printf("    > aliroot");
  Printf("    ");
  Printf("    AliRoot> .L HMPIDFindAlign.C");
  Printf("    AliRoot> HMPIDFindAlign(); summary.txt");
  Printf("    AliRoot> .q");
  Printf("    ");
  Printf("The results of the alignment procedures is in the text file summary.txt");
  Printf("In addition 14 Canvases (7 chamb. x X and Y) are produced to show the residual distrib. in Local Coord.");
  Printf("if some plots are funny, please try to run for a signle chamber modifying some parameters in the macro:");
  Printf("until you find a reasonable result.");
  Printf("    ");
  Printf("AliRoot> HMPIDFindAlign(chamber)");
  Printf("    ");
  Printf("Ciao!."); 
}
Commit	Line	Data
c0a66394	1	#include <TSystem.h>
	2	#include <TVector3.h>
	3	#include <TVector2.h>
	4	#include <TNtuple.h>
	5	#include <TFile.h>
	6	#include <TH2F.h>
	7	#include <TH1F.h>
	8	#include <TCanvas.h>
	9	#include <TGraphErrors.h>
	10	#include <TPaveText.h>
	11	#include <TText.h>
	12	#include <TF1.h>
	13	#include <AliHMPIDParam.h>
	14	#include <AliGeomManager.h>
	15	#include <AliCDBEntry.h>
	16	#include <TClonesArray.h>
	17
	18	TVector3 AlignAxis(Int_t ch,Int_t axis);
	19	Double_t FindChCenterD(Int_t ch);
	20	void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist);
	21	TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc);
	22	void help();
	23
	24	//--------------------------------------------------
	25	void HMPIDFindAlign(Int_t cham=-1)
	26	{
	27	Printf("");
	28	Printf("******************************************");
	29	Printf("* You can type help() for instructions...*");
	30	Printf("******************************************");
	31	Printf("");
	32
	33	gSystem->Sleep(3000);
	34
	35	AliGeomManager::LoadGeometry("geometry.root");
	36
	37	TFile *file = TFile::Open("real_geometry.root");
	38	if(file) {
	39	AliCDBEntry entry = (AliCDBEntry)file->Get("AliCDBEntry");
	40	TClonesArray array = (TClonesArray)entry->GetObject();
	41	AliGeomManager::ApplyAlignObjsToGeom(*array);
	42	}
	43	AliHMPIDParam::Instance(); //just to print at the beginning...
	44
	45	Int_t ichMin = 0;
	46	Int_t ichMax = 6;
	47	if(cham>=0&&cham<=6) { ichMin = cham; ichMax = cham;}
	48
	49	for(Int_t ich=ichMin;ich<=ichMax;ich++) {
	50
	51	Printf("");
	52	Printf("*******SUMMARY FOR CHAMBER %i ********",ich);
	53	Printf ("Distance from IP of chamber %i = %f (cm)",ich,FindChCenterD(ich));
	54
	55	TVector3 xAl = AlignAxis(ich,0);
	56	TVector3 yAl = AlignAxis(ich,1);
	57
	58	Double_t p1 = 1/(xAl.Z()*xAl.Z());
	59	Double_t p2 = 1/(yAl.Z()*yAl.Z());
	60
	61	Double_t distMean = (p1xAl.Y()+p2yAl.Y())/(p1+p2);
	62
	63	FindAligners(ich,xAl.X(),yAl.X(),distMean);
	64	}
65	}
66	//--------------------------------------------------
67	TVector3 AlignAxis(Int_t ch,Int_t axis)
68	{
69
70	// Axis = 0 X Local Axis;
71	// Axis = 1 Y Local Axis;
72
73
74	TH1F* delta[6];
75
76	TFile* file = TFile::Open("histos.root","old");
77	TNtuple *pNtuple;
78	if(!(TNtuple*)file->Get("nTupla")) {
79	// new version because Giacomo uses now a Task to strio the data...
80	TList hmpoutput = (TList)(file->FindObjectAny("hmpoutput"));
81	pNtuple = (TNtuple*)hmpoutput->At(0);
82	} else {
83	pNtuple = (TNtuple*)file->Get("nTupla");
84	}
85
86	for(Int_t j=0;j<6;j++) {
87	if(axis == 0) {
88	delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
89	} else {
90	delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
91	}
92	}
93	Double_t mean[6];
94	Double_t sigma[6];
95	Double_t siz[6];
96	Double_t bin[6];
97
98	TCanvas *c = new TCanvas();
99
100	Int_t index = 0;
101	for(Int_t j=0;j<6;j++) {
102	Double_t xMin = j*20;
103	Double_t xMax = xMin+20;
104	if(axis == 0) {
105	pNtuple->Draw(Form("Xtrk-Xmip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Xmip>%f&&Xmip<%f",ch,xMin,xMax),"0");
106	} else {
107	pNtuple->Draw(Form("Ytrk-Ymip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Ymip>%f&&Ymip<%f",ch,xMin,xMax),"0");
108	}
109	if(delta[j]->GetEntries()<50) continue;
110
111	Double_t dIntFit = 5.;
112	siz[index] = 0.5*(xMax-xMin);
113	bin[index]= 0.5*(xMax+xMin);
114	Double_t xMaxPos = delta[j]->GetBinCenter(delta[j]->FindFirstBinAbove(delta[j]->GetMaximum()-1));
115	delta[j]->Fit("gaus","Q0","",xMaxPos-dIntFit,xMaxPos+dIntFit);
116	mean[index] = (delta[j]->GetFunction("gaus"))->GetParameter(1);
117	sigma[index] = (delta[j]->GetFunction("gaus"))->GetParError(1);
118	index++;
119	}
120	c->Clear();
121	c->Divide(3,3);
122
123	for(Int_t j=0;j<6;j++) {
124	c->cd(j+1);delta[j]->Draw();
125	}
126
127	c->cd(8);
128
129	TGraphErrors *pGr = new TGraphErrors(index,bin,mean,siz,sigma);
130
131	pGr->SetTitle(Form("LRS for chamber %i",ch));
132	if(axis == 0) {
133	pGr->GetXaxis()->SetTitle("Xmip (cm)");
134	pGr->GetYaxis()->SetTitle("deltaX (cm)");
135	} else {
136	pGr->GetXaxis()->SetTitle("Ymip (cm)");
137	pGr->GetYaxis()->SetTitle("DeltaY (cm)");
138	}
139	pGr->SetMinimum(-5.);
140	pGr->SetMaximum( 5.);
141	//pGr->SetMinimum(TMath::MinElement(index,mean)-1.);
142	//pGr->SetMaximum(TMath::MaxElement(index,mean)+1.);
143	pGr->Draw("ALP");
144	pGr->Fit("pol1","Q");
145	Double_t shift = (pGr->GetFunction("pol1"))->GetParameter(0);
146	Double_t errshift = (pGr->GetFunction("pol1"))->GetParError(0);
147	Double_t coefang = (pGr->GetFunction("pol1"))->GetParameter(1);
148	Double_t errcoefang = (pGr->GetFunction("pol1"))->GetParError(1);
149	Double_t prob = (pGr->GetFunction("pol1"))->GetProb();
150	Double_t dist = FindChCenterD(ch);
151
152	Double_t rfShift = -coefang/(1.+coefang)*dist;
153	Double_t errrfShift = TMath::Abs(dist/((1.+coefang)(1+coefang))errcoefang);
154	if (axis == 0) {
155	Printf(" x shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
156	} else {
157	Printf(" y shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
158	}
159
160	c->cd(7);
161
162	TPaveText *pPa = new TPaveText();
163	pPa->DrawPave(0.2,0.2,0.8,0.8);
164	TText *t = new TText();
165	t->SetTextSize(0.08);
166	t->DrawText(0.21,0.70,Form(" shift %5.3f +/- %5.3f",shift,errshift));
167	t->DrawText(0.21,0.50,Form(" dist shift %5.3f +/- %5.3f",rfShift,errrfShift));
168	t->DrawText(0.21,0.30,Form(" dist from IP %5.3f ",dist));
169
170	if (axis == 0) {
171	c->SaveAs(Form("ShiftXChamber%i.gif",ch));
172	} else {
173	c->SaveAs(Form("ShiftYChamber%i.gif",ch));
174	}
175
176	TVector3 align(shift,rfShift,errrfShift);
177	return align;
178	}
179	//------------------------------------------------------------------
180	Double_t FindChCenterD(Int_t ch)
181	{
182
183	AliHMPIDParam *pParam = AliHMPIDParam::Instance();
184	TVector3 xyz = pParam->Lors2Mars(ch,0.5pParam->SizeAllX(),0.5pParam->SizeAllY());
185	return xyz.Mag();
186	}
187	//------------------------------------------------------------------
188	void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist)
189	{
190	TVector3 old = FindNewCenter(ch,0,0,0);
191
192	Printf("");
193	Printf("Old X Y Z of the chamber n. %i",ch);
194	Printf("XC = %f",old.X());
195	Printf("YC = %f",old.Y());
196	Printf("ZC = %f",old.Z());
197	Printf("Dist. from IP = %f",old.Mag());
198	Printf("");
199
200
201	Printf("Local X shift = %f",xShift);
202	Printf("Local Y shift = %f",yShift);
203	Printf("RPhi shift = %f",dist);
204
205	TVector3 shift = FindNewCenter(ch,xShift,yShift,dist);
206
207	Printf("");
208	Printf("New X Y Z of the chamber n. %i",ch);
209	Printf("");
210	Printf("XC = %f YC = %f ZC = %f",shift.X(),shift.Y(),shift.Z());
211	Printf("Dist. from IP = %f",shift.Mag());
212	Printf("");
213	Printf("(next line are the values to be inserted in the Align Object)");
214	Printf("");
215	Printf("XC shift = %5.2f YC shift = %5.2f ZC shift = %5.2f for chamber %i",shift.X()-old.X(),shift.Y()-old.Y(),shift.Z()-old.Z(),ch);
216	Printf("");
217	}
218	//------------------------------------------------------------------
219	TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc)
220	{
221	//First shift in (LocX,LocY) and go in MRS
222	//Then shift radially...
223	AliHMPIDParam *pParam = AliHMPIDParam::Instance();
224	TVector3 xyz = pParam->Lors2Mars(ch,0.5pParam->SizeAllX()+xShLoc,0.5pParam->SizeAllY()+yShLoc,-1);
225	TVector3 xyzSh(1,1,1);
226	xyzSh.SetTheta(xyz.Theta());
227	xyzSh.SetPhi(xyz.Phi());
228	xyzSh.SetMag(xyz.Mag()+zShLoc);
229	return xyzSh;
230	}
231	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
232	void help()
233	{
234	Printf("Instructions to run this macros:");
235	Printf(" 1 - In the same dir where the macro runs it must exist the file geometry.root");
236	Printf(" 2 - assign symbolic links as it follows:");
237	Printf(" a - file from OCDB/HMPID/Align/Data linked to real_geometry.root");
238	Printf(" b - file with the NTupla of ESD Data (produced from AliHMPIDAnalysisTask.C) linked as histos.root");
239	Printf(" ");
240	Printf(" The macro have to be run in the following way:");
241	Printf(" ");
242	Printf(" > aliroot");
243	Printf(" ");
244	Printf(" AliRoot> .L HMPIDFindAlign.C");
245	Printf(" AliRoot> HMPIDFindAlign(); summary.txt");
246	Printf(" AliRoot> .q");
247	Printf(" ");
248	Printf("The results of the alignment procedures is in the text file summary.txt");
249	Printf("In addition 14 Canvases (7 chamb. x X and Y) are produced to show the residual distrib. in Local Coord.");
250	Printf("if some plots are funny, please try to run for a signle chamber modifying some parameters in the macro:");
251	Printf("until you find a reasonable result.");
252	Printf(" ");
253	Printf("AliRoot> HMPIDFindAlign(chamber)");
254	Printf(" ");
255	Printf("Ciao!.");
256	}