new functionality and new class added
[u/mrichter/AliRoot.git] / HMPID / HMPIDFindAlign.C
CommitLineData
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
18TVector3 AlignAxis(Int_t ch,Int_t axis);
19Double_t FindChCenterD(Int_t ch);
20void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist);
21TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc);
22void help();
23
24//--------------------------------------------------
25void 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 = (p1*xAl.Y()+p2*yAl.Y())/(p1+p2);
62
63 FindAligners(ich,xAl.X(),yAl.X(),distMean);
64 }
65}
66//--------------------------------------------------
67TVector3 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//------------------------------------------------------------------
180Double_t FindChCenterD(Int_t ch)
181{
182
183 AliHMPIDParam *pParam = AliHMPIDParam::Instance();
184 TVector3 xyz = pParam->Lors2Mars(ch,0.5*pParam->SizeAllX(),0.5*pParam->SizeAllY());
185 return xyz.Mag();
186}
187//------------------------------------------------------------------
188void 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//------------------------------------------------------------------
219TVector3 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.5*pParam->SizeAllX()+xShLoc,0.5*pParam->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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
232void 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}