pi0 Re/Mi histograms added, pi0 parameterization set to PHOS13bcdef
[u/mrichter/AliRoot.git] / HMPID / HMPIDFindAlign.C
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 = (p1*xAl.Y()+p2*yAl.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.5*pParam->SizeAllX(),0.5*pParam->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.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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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 }