]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/macros/PlotVertexESD.C
correct also the eta position in A side
[u/mrichter/AliRoot.git] / PWG1 / macros / PlotVertexESD.C
CommitLineData
388ca814 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <Riostream.h>
3#include <TCanvas.h>
4#include <TLegend.h>
5#include <TLegendEntry.h>
6#include <TFile.h>
7#include <TLine.h>
8#include <TNtuple.h>
9#include <TStyle.h>
10#include <TString.h>
11#include <TSystem.h>
12#include <TH1F.h>
13#include <TGraphErrors.h>
14#include <TF1.h>
15#endif
16
17void PlotVertexESD(TString vtxtype="SPD",
b583315d 18 TString fname="Vertex.Performance.root",
388ca814 19 TString ntname="fNtupleVertexESD",
20 Bool_t useztrue=kTRUE,
21 Int_t optgif=0) {
22 //-------------------------------------------------------------------------
23 //
24 // Plot output of AliAnalysisTaskVertexESD.
25 // Origin: francesco.prino@to.infn.it
26 //
27 //-------------------------------------------------------------------------
28
29 // ranges for ideal geometry
30 /*
31 Float_t rangeHPull=40.;
32 Float_t rangeHRes=2500.;
33 Int_t binsHRes=250;
34 Float_t rangeGrAve = 100; // micron
35 Float_t rangeGrRms = 900; // micron
36 Float_t rangeGrPull = 5;
37 */
38
39 // ranges for misaligned geometry
40 Float_t rangeHPull=40.;
41 Float_t rangeHRes=10000.*0.5;
42 Int_t binsHRes=250*0.5;
43 Float_t rangeGrAve = 10000; // micron
44 Float_t rangeGrRms = 5000; // micron
45 Float_t rangeGrPull = 15;
46
47
48 Float_t limcont[10]={0.,2.5,3.5,4.5,6.5,9.5,14.5,20.5,99999.};
49 const Int_t nbinsmult=7;
50 Float_t limmult[nbinsmult+1]={0.,10.,14.,24.,30.,44.,60.,99999.};
51 //const Int_t nbinseff=11;
52 //Float_t limeff[nbinseff+1]={0.,0.5,1.5,2.5,3.5,4.5,5.5,7.5,10.5,13.,15.,999999.};
53 const Int_t nbinseff=9;
54 Float_t limeff[nbinseff+1]={0.,0.5,1.5,2.5,3.5,4.5,5.5,7.5,10.5,999999.};
55 const Int_t nbinz = 14;
56 Float_t limitzt[nbinz+1]={-20,-18,-15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18,20};
57 //Float_t limitzt[nbinz+1]={-3,0.,3.,6.,9.,11.,13.,15,17,19,21};
58
59 TH1F **hxm=new TH1F*[nbinsmult];
60 TH1F **hym=new TH1F*[nbinsmult];
61 TH1F **hzm=new TH1F*[nbinsmult];
62 TH1F **hxpullm=new TH1F*[nbinsmult];
63 TH1F **hypullm=new TH1F*[nbinsmult];
64 TH1F **hzpullm=new TH1F*[nbinsmult];
65 TH1F **hmult=new TH1F*[nbinsmult];
66
67 TH1F **hxc=new TH1F*[8];
68 TH1F **hyc=new TH1F*[8];
69 TH1F **hzc=new TH1F*[8];
70 TH1F **hxpullc=new TH1F*[8];
71 TH1F **hypullc=new TH1F*[8];
72 TH1F **hzpullc=new TH1F*[8];
73 TH1F **hcont=new TH1F*[8];
74
75 TH1F **hxz=new TH1F*[nbinz];
76 TH1F **hyz=new TH1F*[nbinz];
77 TH1F **hzz=new TH1F*[nbinz];
78 TH1F **hz=new TH1F*[nbinz];
79
80 TH1F **hmulteff=new TH1F*[nbinseff];
81 TH1F **haux=new TH1F*[nbinseff];
82 TH1F *htot=new TH1F("htot","",100,-1000,1000);
83
84 TH1F **htrkseff=new TH1F*[nbinseff];
85 TH1F **htrksaux=new TH1F*[nbinseff];
86 TH1F **htrks3Daux=new TH1F*[nbinseff];
87
88 Char_t hisnam[10];
89 TCanvas *c1=new TCanvas("c1","Residuals",1000,700);
90 c1->Divide(nbinsmult,3,0.001,0.001);
91 TCanvas *c1z=new TCanvas("c1z","Residuals z",1000,700);
92 c1z->Divide(nbinz,3,0.001,0.001);
93 TCanvas *cp=new TCanvas("cp","Pulls",1000,700);
94 cp->Divide(nbinsmult,3,0.001,0.001);
95
96 TGraphErrors *gavexm=new TGraphErrors(0);
97 TGraphErrors *gaveym=new TGraphErrors(0);
98 TGraphErrors *gavezm=new TGraphErrors(0);
99 TGraphErrors *grmsxm=new TGraphErrors(0);
100 TGraphErrors *grmsym=new TGraphErrors(0);
101 TGraphErrors *grmszm=new TGraphErrors(0);
102 TGraphErrors *gavexmg=new TGraphErrors(0);
103 TGraphErrors *gaveymg=new TGraphErrors(0);
104 TGraphErrors *gavezmg=new TGraphErrors(0);
105 TGraphErrors *grmsxmg=new TGraphErrors(0);
106 TGraphErrors *grmsymg=new TGraphErrors(0);
107 TGraphErrors *grmszmg=new TGraphErrors(0);
108 TGraphErrors *gpullxm=new TGraphErrors(0);
109 TGraphErrors *gpullym=new TGraphErrors(0);
110 TGraphErrors *gpullzm=new TGraphErrors(0);
111 TGraphErrors *gpullxmg=new TGraphErrors(0);
112 TGraphErrors *gpullymg=new TGraphErrors(0);
113 TGraphErrors *gpullzmg=new TGraphErrors(0);
114 TGraphErrors *geffm=new TGraphErrors(0);
115 TGraphErrors *gefftrks=new TGraphErrors(0);
116 TGraphErrors *geff3Dtrks=new TGraphErrors(0);
117
118 TGraphErrors *gavexc=new TGraphErrors(0);
119 TGraphErrors *gaveyc=new TGraphErrors(0);
120 TGraphErrors *gavezc=new TGraphErrors(0);
121 TGraphErrors *grmsxc=new TGraphErrors(0);
122 TGraphErrors *grmsyc=new TGraphErrors(0);
123 TGraphErrors *grmszc=new TGraphErrors(0);
124 TGraphErrors *gpullxc=new TGraphErrors(0);
125 TGraphErrors *gpullyc=new TGraphErrors(0);
126 TGraphErrors *gpullzc=new TGraphErrors(0);
127
128 TGraph * gbeamxz=new TGraph(0);
129 TGraph * gbeamyz=new TGraph(0);
130 TGraph * gbeamxy=new TGraph(0);
131 TH2F *hbeamxz = new TH2F("hbeamxz","",100,-14.,14.,100,-1.5,1.5);
132 TH2F *hbeamyz = new TH2F("hbeamyz","",100,-14.,14.,100,-1.5,1.5);
133 TH1F *hbeamx = new TH1F("hbeamx","",1000,-1.5,1.5);
134 TH1F *hbeamy = new TH1F("hbeamy","",1000,-1.5,1.5);
135 TH2F *hbeamxy = new TH2F("hbeamxy","",100,-1.5,1.5,100,-1.5,1.5);
136
137
138
139 gavexm->SetName("gavexm");
140 grmsxm->SetName("grmsxm");
141 gavexmg->SetName("gavexmg");
142 grmsxmg->SetName("grmsxmg");
143 gpullxm->SetName("gpullxm");
144 gpullxmg->SetName("gpullxmg");
145 gaveym->SetName("gaveym");
146 grmsym->SetName("grmsym");
147 gaveymg->SetName("gaveymg");
148 grmsymg->SetName("grmsymg");
149 gpullym->SetName("gpullym");
150 gpullymg->SetName("gpullymg");
151 gavezm->SetName("gavezm");
152 grmszm->SetName("grmszm");
153 gavezmg->SetName("gavezmg");
154 grmszmg->SetName("grmszmg");
155 gpullzm->SetName("gpullzm");
156 gpullzmg->SetName("gpullzmg");
157 geffm->SetName("geffm");
158 gefftrks->SetName("gefftrks");
159 geff3Dtrks->SetName("geff3Dtrks");
160 gavexc->SetName("gavexc");
161 grmsxc->SetName("grmsxc");
162 gpullxc->SetName("gpullxc");
163 gaveyc->SetName("gaveyc");
164 grmsyc->SetName("grmsyc");
165 gpullyc->SetName("gpullyc");
166 gavezc->SetName("gavezc");
167 grmszc->SetName("grmszc");
168 gpullzc->SetName("gpullzc");
169
170 TGraphErrors *gavexz=new TGraphErrors(0);
171 TGraphErrors *gaveyz=new TGraphErrors(0);
172 TGraphErrors *gavezz=new TGraphErrors(0);
173 TGraphErrors *grmsxz=new TGraphErrors(0);
174 TGraphErrors *grmsyz=new TGraphErrors(0);
175 TGraphErrors *grmszz=new TGraphErrors(0);
176 TGraphErrors *gavexzg=new TGraphErrors(0);
177 TGraphErrors *gaveyzg=new TGraphErrors(0);
178 TGraphErrors *gavezzg=new TGraphErrors(0);
179 TGraphErrors *grmsxzg=new TGraphErrors(0);
180 TGraphErrors *grmsyzg=new TGraphErrors(0);
181 TGraphErrors *grmszzg=new TGraphErrors(0);
182 TGraphErrors *geffz=new TGraphErrors(0);
183
184 gavexz->SetName("gavexz");
185 grmsxz->SetName("grmsxz");
186 gavexzg->SetName("gavexzg");
187 grmsxzg->SetName("grmsxzg");
188 gaveyz->SetName("gaveyz");
189 grmsyz->SetName("grmsyz");
190 gaveyzg->SetName("gaveyzg");
191 grmsyzg->SetName("grmsyzg");
192 gavezz->SetName("gavezz");
193 grmszz->SetName("grmszz");
194 gavezzg->SetName("gavezzg");
195 grmszzg->SetName("grmszzg");
196 geffz->SetName("geffz");
197
198 TF1 *fitf;
199
200 // histogram creation
201 for(Int_t khis=0;khis<nbinseff;khis++){
202 sprintf(hisnam,"hmeff%d",khis);
203 hmulteff[khis]=new TH1F(hisnam,"",100,0.,200.);
204 sprintf(hisnam,"htrkseff%d",khis);
205 htrkseff[khis]=new TH1F(hisnam,"",100,0.,200.);
206 sprintf(hisnam,"haux%d",khis);
207 haux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
208 sprintf(hisnam,"htrksaux%d",khis);
209 htrksaux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
210 sprintf(hisnam,"htrks3Daux%d",khis);
211 htrks3Daux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
212 }
213 for(Int_t khis=0;khis<nbinsmult;khis++){
214 sprintf(hisnam,"hxm%d",khis);
215 hxm[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
216 sprintf(hisnam,"hym%d",khis);
217 hym[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
218 sprintf(hisnam,"hzm%d",khis);
219 hzm[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
220 sprintf(hisnam,"hxpm%d",khis);
221 hxpullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
222 sprintf(hisnam,"hypm%d",khis);
223 hypullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
224 sprintf(hisnam,"hzpm%d",khis);
225 hzpullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
226 sprintf(hisnam,"hmult%d",khis);
227 hmult[khis]=new TH1F(hisnam,"",100,0.,200.);
228 }
229
230 for(Int_t khis=0;khis<8;khis++){
231 sprintf(hisnam,"hxc%d",khis);
232 hxc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
233 sprintf(hisnam,"hyc%d",khis);
234 hyc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
235 sprintf(hisnam,"hzc%d",khis);
236 hzc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
237 sprintf(hisnam,"hxpc%d",khis);
238 hxpullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
239 sprintf(hisnam,"hypc%d",khis);
240 hypullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
241 sprintf(hisnam,"hzpc%d",khis);
242 hzpullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
243 sprintf(hisnam,"hcont%d",khis);
244 hcont[khis]=new TH1F(hisnam,"",100,0.,200.);
245 }
246
247 for(Int_t khis=0;khis<nbinz;khis++){
248 sprintf(hisnam,"hxz%d",khis);
249 hxz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
250 sprintf(hisnam,"hyz%d",khis);
251 hyz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
252 sprintf(hisnam,"hzz%d",khis);
253 hzz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
254 sprintf(hisnam,"hz%d",khis);
255 hz[khis]=new TH1F(hisnam,"",100,-20.,20.);
256 }
257
258 Float_t totev=0,totevtriggered=0,nvtx3D=0,nvtxZ=0;
259 TFile *f=new TFile(fname.Data());
260 TList *cOutput = (TList*)f->Get("cOutput");
261 TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
97805f7d 262 Int_t nnnev=nt->GetEntries();
263 printf("Events = %d\n",nnnev);
388ca814 264 Float_t xVtx,xdiffVtx,xerrVtx;
265 Float_t yVtx,ydiffVtx,yerrVtx;
266 Float_t zVtx,zdiffVtx,zerrVtx;
267 Float_t ntrklets,ncontrVtx,dndy,triggered,vtx3D;
97805f7d 268 Float_t xtrue,ytrue,ztrue,zref;
388ca814 269
270 TString sxx="x"; sxx.Append(vtxtype.Data());
271 nt->SetBranchAddress(sxx.Data(),&xVtx);
272 TString syy="y"; syy.Append(vtxtype.Data());
273 nt->SetBranchAddress(syy.Data(),&yVtx);
274 TString szz="z"; szz.Append(vtxtype.Data());
275 nt->SetBranchAddress(szz.Data(),&zVtx);
276
388ca814 277 TString xerr="xerr"; xerr.Append(vtxtype.Data());
278 nt->SetBranchAddress(xerr.Data(),&xerrVtx);
279 TString yerr="yerr"; yerr.Append(vtxtype.Data());
280 nt->SetBranchAddress(yerr.Data(),&yerrVtx);
281 TString zerr="zerr"; zerr.Append(vtxtype.Data());
282 nt->SetBranchAddress(zerr.Data(),&zerrVtx);
283
284 TString trkstitle;
285 if(vtxtype.Contains("TPC")) {
286 nt->SetBranchAddress("nTPCin",&ntrklets);
287 trkstitle="TPC tracks pointing to beam pipe";
288 } else {
289 nt->SetBranchAddress("ntrklets",&ntrklets);
290 trkstitle="SPD tracklets";
291 }
292 TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
293 nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
294 nt->SetBranchAddress("dndygen",&dndy);
295
97805f7d 296 nt->SetBranchAddress("xtrue",&xtrue);
297 nt->SetBranchAddress("ytrue",&ytrue);
388ca814 298 nt->SetBranchAddress("ztrue",&ztrue);
299
300 nt->SetBranchAddress("triggered",&triggered);
301
302 nt->SetBranchAddress("SPD3D",&vtx3D);
303
304
305 // loop on events
97805f7d 306 for(Int_t iev=0;iev<nnnev;iev++){
388ca814 307 nt->GetEvent(iev);
97805f7d 308
309 xdiffVtx=10000.*(xVtx-xtrue);
310 ydiffVtx=10000.*(yVtx-ytrue);
311 zdiffVtx=10000.*(zVtx-ztrue);
312
313
388ca814 314 zref = (useztrue ? ztrue : zVtx);
315 if(!vtxtype.Contains("SPD")) vtx3D=1.;
316 if(triggered<0.5) continue; // not triggered
317 totevtriggered += 1.;
318 if(ncontrVtx>0) {
319 htot->Fill(zdiffVtx);
320 if(vtx3D>0.5) {
321 nvtx3D += 1.;
322 } else {
323 nvtxZ += 1.;
324 }
325 }
326
327 if(ncontrVtx>0 && vtx3D>0.5) {
328 gbeamxz->SetPoint(gbeamxz->GetN(),zVtx,xVtx);
329 gbeamyz->SetPoint(gbeamxz->GetN(),zVtx,yVtx);
330 gbeamxy->SetPoint(gbeamxz->GetN(),xVtx,yVtx);
331
332 hbeamx->Fill(xVtx);
333 hbeamy->Fill(yVtx);
334 hbeamxz->Fill(zVtx,xVtx);
335 hbeamyz->Fill(zVtx,yVtx);
336 hbeamxy->Fill(xVtx,yVtx);
337 }
338
339 for(Int_t khis=0;khis<nbinseff;khis++){
340 if(dndy>=limeff[khis] && dndy<limeff[khis+1]){
341 hmulteff[khis]->Fill(dndy);
342 if(ncontrVtx>0) haux[khis]->Fill(zdiffVtx);
343 }
344 if(ntrklets>=limeff[khis] && ntrklets<limeff[khis+1]){
345 htrkseff[khis]->Fill(ntrklets);
346 if(ncontrVtx>0) htrksaux[khis]->Fill(zdiffVtx);
347 if(ncontrVtx>0 && vtx3D>0.5) htrks3Daux[khis]->Fill(zdiffVtx);
348 }
349 }
350 for(Int_t khis=0;khis<nbinsmult;khis++){
351 if(ntrklets>=limmult[khis] && ntrklets<limmult[khis+1]){
352 hmult[khis]->Fill(ntrklets);
353 if(ncontrVtx>0){
354 if(vtx3D>0.5) hxm[khis]->Fill(xdiffVtx);
355 if(vtx3D>0.5) hym[khis]->Fill(ydiffVtx);
356 hzm[khis]->Fill(zdiffVtx);
357 if(vtx3D>0.5) hxpullm[khis]->Fill(xdiffVtx/10000./xerrVtx);
358 if(vtx3D>0.5) hypullm[khis]->Fill(ydiffVtx/10000./yerrVtx);
359 hzpullm[khis]->Fill(zdiffVtx/10000./zerrVtx);
360 }
361 }
362 }
363 for(Int_t khis=0;khis<8;khis++){
364 if(ncontrVtx>=limcont[khis]&&ncontrVtx<limcont[khis+1]){
365 hcont[khis]->Fill(ncontrVtx);
366 if(ncontrVtx>0){
367 if(vtx3D>0.5)hxc[khis]->Fill(xdiffVtx);
368 if(vtx3D>0.5)hyc[khis]->Fill(ydiffVtx);
369 hzc[khis]->Fill(zdiffVtx);
370 if(vtx3D>0.5)hxpullc[khis]->Fill(xdiffVtx/10000./xerrVtx);
371 if(vtx3D>0.5)hypullc[khis]->Fill(ydiffVtx/10000./yerrVtx);
372 hzpullc[khis]->Fill(zdiffVtx/10000./zerrVtx);
373 }
374 }
375 }
376 for(Int_t khis=0;khis<nbinz;khis++){
377 if(zref>=limitzt[khis]&&zref<limitzt[khis+1]){
378 hz[khis]->Fill(zref);
379 if(ncontrVtx>0){
380 if(vtx3D>0.5)hxz[khis]->Fill(xdiffVtx);
381 if(vtx3D>0.5)hyz[khis]->Fill(ydiffVtx);
382 hzz[khis]->Fill(zdiffVtx);
383 }
384 }
385 }
386 }
97805f7d 387 totev+=nnnev;
388ca814 388
389 if(totev==0){
390 printf("Total number of events = 0\n");
391 return;
392 }
393 for(Int_t khis=0;khis<nbinseff;khis++){
394 Double_t x=hmulteff[khis]->GetMean();
395 Double_t ex=hmulteff[khis]->GetRMS();;
396 Float_t nEv=(Float_t)(hmulteff[khis]->GetEntries());
397 Float_t trkeff=-1.;
398 cout<<"Eff. dNch/dy bin "<<khis<<" ("<<limeff[khis]<<"-"<<limeff[khis+1]<<") # Events ="<<nEv<<" with vertex ="<<haux[khis]->GetEntries()<<endl;
399 if(nEv>0) trkeff=(Float_t)(haux[khis]->GetEntries())/nEv;
400 geffm->SetPoint(khis,x,trkeff);
401 geffm->SetPointError(khis,ex,0.);
402 }
403 for(Int_t khis=0;khis<nbinseff;khis++){
404 Double_t x=htrkseff[khis]->GetMean();
405 Double_t ex=htrkseff[khis]->GetRMS();;
406 Float_t nEv=(Float_t)(htrkseff[khis]->GetEntries());
407 Float_t trkeff=-1.;
408 cout<<"Eff. trks bin "<<khis<<" ("<<limeff[khis]<<"-"<<limeff[khis+1]<<") # Events ="<<nEv<<" with vertex ="<<htrksaux[khis]->GetEntries()<<endl;
409 if(nEv>0) trkeff=(Float_t)(htrksaux[khis]->GetEntries())/nEv;
410 gefftrks->SetPoint(khis,x,trkeff);
411 gefftrks->SetPointError(khis,ex,0.);
412 if(nEv>0) trkeff=(Float_t)(htrks3Daux[khis]->GetEntries())/nEv;
413 geff3Dtrks->SetPoint(khis,x,trkeff);
414 geff3Dtrks->SetPointError(khis,ex,0.);
415 }
416
417 for(Int_t khis=0;khis<nbinsmult;khis++){
418
419 c1->cd(khis+1);
420 if(hxm[khis]->GetEntries()<10) continue;
421 hxm[khis]->Draw();
422 hxm[khis]->Fit("gaus","Q0");
423 fitf= hxm[khis]->GetFunction("gaus");
424 Double_t avexg=fitf->GetParameter(1);
425 Double_t eavexg=fitf->GetParError(1);
426 Double_t rmsxg=fitf->GetParameter(2);
427 Double_t ermsxg=fitf->GetParError(2);
428 c1->cd(nbinsmult+khis+1);
429 hym[khis]->Draw();
430 hym[khis]->Fit("gaus","Q0");
431 fitf= hym[khis]->GetFunction("gaus");
432 Double_t aveyg=fitf->GetParameter(1);
433 Double_t eaveyg=fitf->GetParError(1);
434 Double_t rmsyg=fitf->GetParameter(2);
435 Double_t ermsyg=fitf->GetParError(2);
436 c1->cd(2*nbinsmult+khis+1);
437 hzm[khis]->Draw();
438 hzm[khis]->Fit("gaus","Q0");
439 fitf= hzm[khis]->GetFunction("gaus");
440 Double_t avezg=fitf->GetParameter(1);
441 Double_t eavezg=fitf->GetParError(1);
442 Double_t rmszg=fitf->GetParameter(2);
443 Double_t ermszg=fitf->GetParError(2);
444
445 cp->cd(khis+1);
446 hxpullm[khis]->Draw();
447 hxpullm[khis]->Fit("gaus","Q0");
448 fitf= hxpullm[khis]->GetFunction("gaus");
449 Double_t pullxg=fitf->GetParameter(2);
450 Double_t epullxg=fitf->GetParError(2);
451 cp->cd(nbinsmult+khis+1);
452 hypullm[khis]->Draw();
453 hypullm[khis]->Fit("gaus","Q0");
454 fitf= hypullm[khis]->GetFunction("gaus");
455 Double_t pullyg=fitf->GetParameter(2);
456 Double_t epullyg=fitf->GetParError(2);
457 cp->cd(2*nbinsmult+khis+1);
458 hzpullm[khis]->Draw();
459 hzpullm[khis]->Fit("gaus","Q0");
460 fitf= hzpullm[khis]->GetFunction("gaus");
461 Double_t pullzg=fitf->GetParameter(2);
462 Double_t epullzg=fitf->GetParError(2);
463
464
465 Double_t rmsxt=hxm[khis]->GetRMS();
466 Double_t rmsyt=hym[khis]->GetRMS();
467 Double_t rmszt=hzm[khis]->GetRMS();
468 Double_t ermsxt=hxm[khis]->GetRMSError();
469 Double_t ermsyt=hym[khis]->GetRMSError();
470 Double_t ermszt=hzm[khis]->GetRMSError();
471 Double_t avext=hxm[khis]->GetMean();
472 Double_t aveyt=hym[khis]->GetMean();
473 Double_t avezt=hzm[khis]->GetMean();
474 Double_t eavext=hxm[khis]->GetMeanError();
475 Double_t eaveyt=hym[khis]->GetMeanError();
476 Double_t eavezt=hzm[khis]->GetMeanError();
477 Double_t pullxt=hxpullm[khis]->GetRMS();
478 Double_t pullyt=hypullm[khis]->GetRMS();
479 Double_t pullzt=hzpullm[khis]->GetRMS();
480 Double_t epullxt=hxpullm[khis]->GetRMSError();
481 Double_t epullyt=hypullm[khis]->GetRMSError();
482 Double_t epullzt=hzpullm[khis]->GetRMSError();
483
484 Double_t x=hmult[khis]->GetMean();
485 if(hmult[khis]->GetEntries()==0) x=-1;
486 Double_t ex=hmult[khis]->GetRMS();;
487
488 gavexm->SetPoint(khis,x,avext);
489 gavexm->SetPointError(khis,ex,eavext);
490 gaveym->SetPoint(khis,x,aveyt);
491 gaveym->SetPointError(khis,ex,eaveyt);
492 gavezm->SetPoint(khis,x,avezt);
493 gavezm->SetPointError(khis,ex,eavezt);
494 grmsxm->SetPoint(khis,x,rmsxt);
495 grmsxm->SetPointError(khis,ex,ermsxt);
496 grmsym->SetPoint(khis,x,rmsyt);
497 grmsym->SetPointError(khis,ex,ermsyt);
498 grmszm->SetPoint(khis,x,rmszt);
499 grmszm->SetPointError(khis,ex,ermszt);
500
501 gavexmg->SetPoint(khis,x,avexg);
502 gavexmg->SetPointError(khis,ex,eavexg);
503 gaveymg->SetPoint(khis,x,aveyg);
504 gaveymg->SetPointError(khis,ex,eaveyg);
505 gavezmg->SetPoint(khis,x,avezg);
506 gavezmg->SetPointError(khis,ex,eavezg);
507 grmsxmg->SetPoint(khis,x,rmsxg);
508 grmsxmg->SetPointError(khis,ex,ermsxg);
509 grmsymg->SetPoint(khis,x,rmsyg);
510 grmsymg->SetPointError(khis,ex,ermsyg);
511 grmszmg->SetPoint(khis,x,rmszg);
512 grmszmg->SetPointError(khis,ex,ermszg);
513
514 gpullxm->SetPoint(khis,x,pullxt);
515 gpullxm->SetPointError(khis,ex,epullxt);
516 gpullym->SetPoint(khis,x,pullyt);
517 gpullym->SetPointError(khis,ex,epullyt);
518 gpullzm->SetPoint(khis,x,pullzt);
519 gpullzm->SetPointError(khis,ex,epullzt);
520
521 gpullxmg->SetPoint(khis,x,pullxg);
522 gpullxmg->SetPointError(khis,ex,epullxg);
523 gpullymg->SetPoint(khis,x,pullyg);
524 gpullymg->SetPointError(khis,ex,epullyg);
525 gpullzmg->SetPoint(khis,x,pullzg);
526 gpullzmg->SetPointError(khis,ex,epullzg);
527
528 Float_t nEv=hmult[khis]->GetEntries();
529 cout<<"Mult. bin "<<khis<<" # Events ="<<nEv<<endl;
530 }
531
532 for(Int_t khis=0;khis<8;khis++){
533
534 Double_t rmsxt=hxc[khis]->GetRMS();
535 Double_t rmsyt=hyc[khis]->GetRMS();
536 Double_t rmszt=hzc[khis]->GetRMS();
537 Double_t ermsxt=hxc[khis]->GetRMSError();
538 Double_t ermsyt=hyc[khis]->GetRMSError();
539 Double_t ermszt=hzc[khis]->GetRMSError();
540 Double_t avext=hxc[khis]->GetMean();
541 Double_t aveyt=hyc[khis]->GetMean();
542 Double_t avezt=hzc[khis]->GetMean();
543 Double_t eavext=hxc[khis]->GetMeanError();
544 Double_t eaveyt=hyc[khis]->GetMeanError();
545 Double_t eavezt=hzc[khis]->GetMeanError();
546 Double_t pullxt=hxpullc[khis]->GetRMS();
547 Double_t pullyt=hypullc[khis]->GetRMS();
548 Double_t pullzt=hzpullc[khis]->GetRMS();
549 Double_t epullxt=hxpullc[khis]->GetRMSError();
550 Double_t epullyt=hypullc[khis]->GetRMSError();
551 Double_t epullzt=hzpullc[khis]->GetRMSError();
552
553 Double_t x=hcont[khis]->GetMean();
554 Double_t ex=hcont[khis]->GetRMS();;
555
556 gavexc->SetPoint(khis,x,avext);
557 gavexc->SetPointError(khis,ex,eavext);
558 gaveyc->SetPoint(khis,x,aveyt);
559 gaveyc->SetPointError(khis,ex,eaveyt);
560 gavezc->SetPoint(khis,x,avezt);
561 gavezc->SetPointError(khis,ex,eavezt);
562 grmsxc->SetPoint(khis,x,rmsxt);
563 grmsxc->SetPointError(khis,ex,ermsxt);
564 grmsyc->SetPoint(khis,x,rmsyt);
565 grmsyc->SetPointError(khis,ex,ermsyt);
566 grmszc->SetPoint(khis,x,rmszt);
567 grmszc->SetPointError(khis,ex,ermszt);
568
569 gpullxc->SetPoint(khis,x,pullxt);
570 gpullxc->SetPointError(khis,ex,epullxt);
571 gpullyc->SetPoint(khis,x,pullyt);
572 gpullyc->SetPointError(khis,ex,epullyt);
573 gpullzc->SetPoint(khis,x,pullzt);
574 gpullzc->SetPointError(khis,ex,epullzt);
575
576 Float_t nEv=hcont[khis]->GetEntries();
577 cout<<"Contrib. bin "<<khis<<" # Events ="<<nEv<<endl;
578 }
579
580 for(Int_t khis=0; khis<nbinz; khis++){
581
582 Double_t rmsxt=hxz[khis]->GetRMS();
583 Double_t rmsyt=hyz[khis]->GetRMS();
584 Double_t rmszt=hzz[khis]->GetRMS();
585 Double_t ermsxt=hxz[khis]->GetRMSError();
586 Double_t ermsyt=hyz[khis]->GetRMSError();
587 Double_t ermszt=hzz[khis]->GetRMSError();
588 Double_t avext=hxz[khis]->GetMean();
589 Double_t aveyt=hyz[khis]->GetMean();
590 Double_t avezt=hzz[khis]->GetMean();
591 Double_t eavext=hxz[khis]->GetMeanError();
592 Double_t eaveyt=hyz[khis]->GetMeanError();
593 Double_t eavezt=hzz[khis]->GetMeanError();
594
595 Float_t nEv=hz[khis]->GetEntries();
596 Double_t x=-999.;
597 if(nEv>0) x=hz[khis]->GetMean();
598 Double_t ex=hz[khis]->GetRMS();;
599
600 gavexz->SetPoint(khis,x,avext);
601 gavexz->SetPointError(khis,ex,eavext);
602 gaveyz->SetPoint(khis,x,aveyt);
603 gaveyz->SetPointError(khis,ex,eaveyt);
604 gavezz->SetPoint(khis,x,avezt);
605 gavezz->SetPointError(khis,ex,eavezt);
606 grmsxz->SetPoint(khis,x,rmsxt);
607 grmsxz->SetPointError(khis,ex,ermsxt);
608 grmsyz->SetPoint(khis,x,rmsyt);
609 grmsyz->SetPointError(khis,ex,ermsyt);
610 grmszz->SetPoint(khis,x,rmszt);
611 grmszz->SetPointError(khis,ex,ermszt);
612
613 c1z->cd(khis+1);
614 hxz[khis]->Draw();
615 hxz[khis]->Fit("gaus","Q0");
616 fitf= hxz[khis]->GetFunction("gaus");
617 Double_t avexg=fitf->GetParameter(1);
618 Double_t eavexg=fitf->GetParError(1);
619 Double_t rmsxg=fitf->GetParameter(2);
620 Double_t ermsxg=fitf->GetParError(2);
621 c1z->cd(1*nbinz+khis+1);
622 hyz[khis]->Draw();
623 hyz[khis]->Fit("gaus","Q0");
624 fitf= hyz[khis]->GetFunction("gaus");
625 Double_t aveyg=fitf->GetParameter(1);
626 Double_t eaveyg=fitf->GetParError(1);
627 Double_t rmsyg=fitf->GetParameter(2);
628 Double_t ermsyg=fitf->GetParError(2);
629 c1z->cd(2*nbinz+khis+1);
630 hzz[khis]->Draw();
631 hzz[khis]->Fit("gaus","Q0");
632 fitf= hzz[khis]->GetFunction("gaus");
633 Double_t avezg=fitf->GetParameter(1);
634 Double_t eavezg=fitf->GetParError(1);
635 Double_t rmszg=fitf->GetParameter(2);
636 Double_t ermszg=fitf->GetParError(2);
637
638 gavexzg->SetPoint(khis,x,avexg);
639 gavexzg->SetPointError(khis,ex,eavexg);
640 gaveyzg->SetPoint(khis,x,aveyg);
641 gaveyzg->SetPointError(khis,ex,eaveyg);
642 gavezzg->SetPoint(khis,x,avezg);
643 gavezzg->SetPointError(khis,ex,eavezg);
644 grmsxzg->SetPoint(khis,x,rmsxg);
645 grmsxzg->SetPointError(khis,ex,ermsxg);
646 grmsyzg->SetPoint(khis,x,rmsyg);
647 grmsyzg->SetPointError(khis,ex,ermsyg);
648 grmszzg->SetPoint(khis,x,rmszg);
649 grmszzg->SetPointError(khis,ex,ermszg);
650
651 Float_t zeff=-999.;
652 if(nEv>0) zeff=hzz[khis]->GetEntries()/nEv;
653 geffz->SetPoint(khis,x,zeff);
654 geffz->SetPointError(khis,ex,0.);
655
656 cout<<"Z bin "<<khis<<" # Events ="<<nEv<<endl;
657 }
658
659 Double_t efftrk=htot->GetEntries()/totevtriggered;
660
661 printf("EVENTS STATISTICS:\n Total: %d\n Triggered (MB1): %d\n Triggered and with vertex %d\n",totev,totevtriggered,htot->GetEntries());
662 if(vtxtype.Contains("SPD")) printf(" %d with Vertexer3D, %d with VertexerZ\n",nvtx3D,nvtxZ);
663 printf("Overall efficiency (for triggered evts) Vertexer%s = %f / %f = %f\n",vtxtype.Data(),htot->GetEntries(),totevtriggered,efftrk);
664
665 TFile* in = new TFile("vert-graphs.root","recreate");
666 gbeamxz->Write("gbeamxz");
667 gbeamyz->Write("gbeamyz");
668 grmsxm->Write();
669 grmsxmg->Write();
670 gpullxm->Write();
671 gpullxmg->Write();
672 gavexm->Write();
673 gavexmg->Write();
674 grmsym->Write();
675 grmsymg->Write();
676 gpullym->Write();
677 gpullymg->Write();
678 gaveym->Write();
679 gaveymg->Write();
680 grmszm->Write();
681 grmszmg->Write();
682 gpullzm->Write();
683 gpullzmg->Write();
684 gavezm->Write();
685 gavezmg->Write();
686 geffm->Write();
687 gefftrks->Write();
688 geff3Dtrks->Write();
689 grmsxc->Write();
690 gpullxc->Write();
691 gavexc->Write();
692 grmsyc->Write();
693 gpullyc->Write();
694 gaveyc->Write();
695 grmszc->Write();
696 gpullzc->Write();
697 gavezc->Write();
698 gavexz->Write();
699 gaveyz->Write();
700 gavezz->Write();
701 grmsxz->Write();
702 grmsyz->Write();
703 grmszz->Write();
704 gavexzg->Write();
705 gaveyzg->Write();
706 gavezzg->Write();
707 grmsxzg->Write();
708 grmsyzg->Write();
709 grmszzg->Write();
710
711 in->Close();
712
713 gStyle->SetOptTitle(0);
714
715 Char_t outgif[100];
716 TLine *lin0=new TLine(0,0,60,0);
717 lin0->SetLineStyle(2);
718
719 TCanvas *cg1=new TCanvas("cg1","Histo mean");
720 cg1->SetBottomMargin(0.14);
721 cg1->SetTopMargin(0.08);
722 cg1->SetLeftMargin(0.14);
723 cg1->SetRightMargin(0.08);
724 gavexm->SetMarkerStyle(20);
725 gavexm->SetMarkerColor(1);
726 gaveym->SetMarkerStyle(21);
727 gaveym->SetMarkerColor(2);
728 gavezm->SetMarkerStyle(22);
729 gavezm->SetMarkerColor(4);
730 TLegend *leg=new TLegend(0.18,0.70,0.25,0.90);
731 leg->SetBorderSize(0);
732 leg->SetFillStyle(0);
733 TLegendEntry *ent=leg->AddEntry(gavexm,"x","P");
734 ent->SetTextColor(1);
735 ent=leg->AddEntry(gaveym,"y","P");
736 ent->SetTextColor(2);
737 ent=leg->AddEntry(gavezm,"z","P");
738 ent->SetTextColor(4);
739 gavezm->GetXaxis()->SetLimits(0.,60.);
740 gavezm->SetMinimum(-rangeGrAve);
741 gavezm->SetMaximum(rangeGrAve);
742 gavezm->Draw("AP");
743 gavezm->GetXaxis()->SetTitle(trkstitle.Data());
744 gavezm->GetXaxis()->SetTitleSize(0.05);
745 gavezm->GetYaxis()->SetTitle("<Pos_{found}-Pos_{true}> [#mum]");
746 gavezm->GetYaxis()->SetTitleSize(0.05);
747 gavexm->Draw("PSAME");
748 gaveym->Draw("PSAME");
749 leg->Draw();
750 lin0->Draw();
751 sprintf(outgif,"vert%s-ave-mult.gif",vtxtype.Data());
752 if(optgif) cg1->SaveAs(outgif);
753
754
755 TCanvas *cg2=new TCanvas("cg2","Histo RMS");
756 cg2->SetBottomMargin(0.14);
757 cg2->SetTopMargin(0.08);
758 cg2->SetLeftMargin(0.14);
759 cg2->SetRightMargin(0.08);
760 grmsxm->SetMarkerStyle(20);
761 grmsxm->SetMarkerColor(1);
762 grmsym->SetMarkerStyle(21);
763 grmsym->SetMarkerColor(2);
764 grmszm->SetMarkerStyle(22);
765 grmszm->SetMarkerColor(4);
766 grmszm->SetMinimum(0);
767 grmszm->SetMaximum(rangeGrRms);
768 grmszm->GetXaxis()->SetLimits(0,60);
769 grmszm->Draw("AP");
770 grmszm->GetXaxis()->SetTitle(trkstitle.Data());
771 grmszm->GetXaxis()->SetTitleSize(0.05);
772 grmszm->GetYaxis()->SetTitle("Resolution [#mum]");
773 grmszm->GetYaxis()->SetTitleSize(0.05);
774 grmsym->Draw("PSAME");
775 grmsxm->Draw("PSAME");
776 grmsxmg->SetMarkerStyle(24);
777 grmsxmg->SetMarkerColor(1);
778 grmsymg->SetMarkerStyle(25);
779 grmsymg->SetMarkerColor(2);
780 grmszmg->SetMarkerStyle(26);
781 grmszmg->SetMarkerColor(4);
782 grmszmg->SetMinimum(0);
783 grmszmg->SetMaximum(rangeGrRms);
784 grmszmg->GetXaxis()->SetLimits(0,60);
785 grmszmg->Draw("PSAME");
786 grmszmg->GetXaxis()->SetTitle(trkstitle.Data());
787 grmszmg->GetXaxis()->SetTitleSize(0.05);
788 grmszmg->GetYaxis()->SetTitle("Resolution [#mum]");
789 grmszmg->GetYaxis()->SetTitleSize(0.05);
790 grmsymg->Draw("PSAME");
791 grmsxmg->Draw("PSAME");
792 leg->Draw();
793 sprintf(outgif,"vert%s-rms-mult.gif",vtxtype.Data());
794 if(optgif) cg2->SaveAs(outgif);
795
796
797
798
799 TCanvas *cg3=new TCanvas("cg3","Efficiency vs dNch/dy");
800 cg3->SetBottomMargin(0.14);
801 cg3->SetTopMargin(0.08);
802 cg3->SetLeftMargin(0.14);
803 cg3->SetRightMargin(0.08);
804 geffm->SetMarkerStyle(22);
805 geffm->SetMarkerColor(1);
806 geffm->GetXaxis()->SetLimits(0.,40.);
807 geffm->SetMinimum(0.);
808 geffm->SetMaximum(1.2);
809 geffm->Draw("AP");
810 geffm->GetXaxis()->SetTitle("MC dN_{ch}/dy in |y|<1");
811 geffm->GetXaxis()->SetTitleSize(0.05);
812 geffm->GetYaxis()->SetTitle("efficiency");
813 geffm->GetYaxis()->SetTitleSize(0.05);
814 sprintf(outgif,"vert%s-eff-mult.gif",vtxtype.Data());
815 if(optgif) cg3->SaveAs(outgif);
816
817
818 TCanvas *cg3b=new TCanvas("cg3b","Efficiency vs tracks");
819 cg3b->SetBottomMargin(0.14);
820 cg3b->SetTopMargin(0.08);
821 cg3b->SetLeftMargin(0.14);
822 cg3b->SetRightMargin(0.08);
823 gefftrks->SetMarkerStyle(22);
824 gefftrks->SetMarkerColor(1);
825 gefftrks->GetXaxis()->SetLimits(0.,40.);
826 gefftrks->SetMinimum(0.);
827 gefftrks->SetMaximum(1.2);
828 gefftrks->Draw("AP");
829 gefftrks->GetXaxis()->SetTitle(trkstitle.Data());
830 gefftrks->GetXaxis()->SetTitleSize(0.05);
831 gefftrks->GetYaxis()->SetTitle("efficiency");
832 gefftrks->GetYaxis()->SetTitleSize(0.05);
833 if(vtxtype.Contains("SPD")) {
834 geff3Dtrks->SetMarkerStyle(26);
835 geff3Dtrks->SetMarkerColor(1);
836 geff3Dtrks->Draw("P");
837 }
838 sprintf(outgif,"vert%s-eff-mult.gif",vtxtype.Data());
839 if(optgif) cg3b->SaveAs(outgif);
840
841
842 TCanvas *cg4=new TCanvas("cg4","Pulls");
843 cg4->SetBottomMargin(0.14);
844 cg4->SetTopMargin(0.08);
845 cg4->SetLeftMargin(0.14);
846 cg4->SetRightMargin(0.08);
847 gpullxm->SetMarkerStyle(20);
848 gpullxm->SetMarkerColor(1);
849 gpullym->SetMarkerStyle(21);
850 gpullym->SetMarkerColor(2);
851 gpullzm->SetMarkerStyle(22);
852 gpullzm->SetMarkerColor(4);
853 gpullzm->GetXaxis()->SetLimits(0,60);
854 gpullzm->SetMinimum(0);
855 gpullzm->SetMaximum(rangeGrPull);
856 gpullzm->Draw("AP");
857 gpullzm->GetXaxis()->SetTitle(trkstitle.Data());
858 gpullzm->GetXaxis()->SetTitleSize(0.05);
859 gpullzm->GetYaxis()->SetTitle("PULL");
860 gpullzm->GetYaxis()->SetTitleSize(0.05);
861 gpullxm->Draw("PSAME");
862 gpullym->Draw("PSAME");
863 gpullxmg->SetMarkerStyle(24);
864 gpullxmg->SetMarkerColor(1);
865 gpullymg->SetMarkerStyle(25);
866 gpullymg->SetMarkerColor(2);
867 gpullzmg->SetMarkerStyle(26);
868 gpullzmg->SetMarkerColor(4);
869 gpullzmg->GetXaxis()->SetLimits(0,60);
870 gpullzmg->SetMinimum(0);
871 gpullzmg->SetMaximum(rangeGrPull);
872 gpullzmg->Draw("PSAME");
873 gpullzmg->GetXaxis()->SetTitle(trkstitle.Data());
874 gpullzmg->GetXaxis()->SetTitleSize(0.05);
875 gpullzmg->GetYaxis()->SetTitle("PULL");
876 gpullzmg->GetYaxis()->SetTitleSize(0.05);
877 gpullxmg->Draw("PSAME");
878 gpullymg->Draw("PSAME");
879 TLine *lin=new TLine(0,1,60,1);
880 lin->SetLineStyle(2);
881 lin->Draw();
882 leg->Draw();
883 sprintf(outgif,"vert%s-pull-mult.gif",vtxtype.Data());
884 if(optgif) cg4->SaveAs(outgif);
885
886
887
888
889 TCanvas *cz1=new TCanvas("cz1","Efficiency vs. Z");
890 cz1->SetBottomMargin(0.14);
891 cz1->SetTopMargin(0.08);
892 cz1->SetLeftMargin(0.14);
893 cz1->SetRightMargin(0.08);
894 geffz->SetMarkerStyle(22);
895 geffz->SetMarkerColor(1);
896 geffz->GetXaxis()->SetLimits(-20,20.);
897 geffz->SetMinimum(0.);
898 geffz->SetMaximum(1.);
899 geffz->Draw("AP");
900 geffz->GetXaxis()->SetTitle("Z [cm]");
901 geffz->GetXaxis()->SetTitleSize(0.05);
902 geffz->GetYaxis()->SetTitle("efficiency");
903 geffz->GetYaxis()->SetTitleSize(0.05);
904 sprintf(outgif,"vert%s-eff-z.gif",vtxtype.Data());
905 if(optgif) cz1->SaveAs(outgif);
906
907
908 TLine *lin0z=new TLine(-20,0,20,0);
909 lin0z->SetLineStyle(2);
910 TCanvas *cz2=new TCanvas("cz2","Mean vs. Z");
911 cz2->SetBottomMargin(0.14);
912 cz2->SetTopMargin(0.08);
913 cz2->SetLeftMargin(0.14);
914 cz2->SetRightMargin(0.08);
915 gavexz->SetMarkerStyle(20);
916 gavexz->SetMarkerColor(1);
917 gaveyz->SetMarkerStyle(21);
918 gaveyz->SetMarkerColor(2);
919 gavezz->SetMarkerStyle(22);
920 gavezz->SetMarkerColor(4);
921 gavezz->GetXaxis()->SetLimits(-20,20.);
922 gavezz->SetMinimum(-rangeGrAve);
923 gavezz->SetMaximum(rangeGrAve);
924 gavezz->Draw("AP");
925 gavezz->GetXaxis()->SetTitle("Z [cm]");
926 gavezz->GetXaxis()->SetTitleSize(0.05);
927 gavezz->GetYaxis()->SetTitle("<Pos_{found}-Pos_{true}> [#mum]");
928 gavezz->GetYaxis()->SetTitleSize(0.05);
929 gavexz->Draw("PSAME");
930 gaveyz->Draw("PSAME");
931 lin0z->Draw();
932 gavexzg->SetMarkerStyle(24);
933 gavexzg->SetMarkerColor(1);
934 gavexzg->Draw("P");
935 gaveyzg->SetMarkerStyle(25);
936 gaveyzg->SetMarkerColor(2);
937 gaveyzg->Draw("P");
938 gavezzg->SetMarkerStyle(26);
939 gavezzg->SetMarkerColor(4);
940 gavezzg->Draw("P");
941 leg->Draw();
942 sprintf(outgif,"vert%s-ave-z.gif",vtxtype.Data());
943 if(optgif) cz2->SaveAs(outgif);
944
945
946 TCanvas *cz3=new TCanvas("cz3","Resolution vs. Z");
947 cz3->SetBottomMargin(0.14);
948 cz3->SetTopMargin(0.08);
949 cz3->SetLeftMargin(0.14);
950 cz3->SetRightMargin(0.08);
951 grmsxz->SetMarkerStyle(20);
952 grmsxz->SetMarkerColor(1);
953 grmsyz->SetMarkerStyle(21);
954 grmsyz->SetMarkerColor(2);
955 grmszz->SetMarkerStyle(22);
956 grmszz->SetMarkerColor(4);
957 grmszz->SetMinimum(0);
958 grmszz->SetMaximum(rangeGrRms);
959 grmszz->GetXaxis()->SetLimits(-20,20);
960 grmszz->Draw("AP");
961 grmszz->GetXaxis()->SetTitle("Z [cm]");
962 grmszz->GetXaxis()->SetTitleSize(0.05);
963 grmszz->GetYaxis()->SetTitle("Resolution [#mum]");
964 grmszz->GetYaxis()->SetTitleSize(0.05);
965 grmsxz->Draw("PSAME");
966 grmsyz->Draw("PSAME");
967 leg->Draw();
968 sprintf(outgif,"vert%s-rms-z.gif",vtxtype.Data());
969 if(optgif) cz3->SaveAs(outgif);
970
971 gStyle->SetPalette(1);
972 TCanvas *cbeam = new TCanvas("cbeam","Beam Long",800,800);
973 cbeam->Divide(2,2);
974 cbeam->cd(1);
975 hbeamx->GetYaxis()->SetTitle("X [cm]");
976 hbeamx->Draw();
977 cbeam->cd(3);
978 hbeamx->GetYaxis()->SetTitle("Y [cm]");
979 hbeamy->Draw();
980 cbeam->cd(2);
981 cbeam_2->SetLogz();
982 //gbeamxz->SetMarkerStyle(7);
983 hbeamxz->Draw("colz");
984 hbeamxz->GetXaxis()->SetTitle("Z [cm]");
985 hbeamxz->GetYaxis()->SetTitle("X [cm]");
986 cbeam_1->Update();
987 TPaveStats *st1=(TPaveStats*)hbeamxz->GetListOfFunctions()->FindObject("stats");
988 st1->SetX1NDC(0.13);
989 st1->SetX2NDC(0.33);
990 cbeam_2->Modified();
991 cbeam_2->Update();
992 cbeam->cd(4);
993 cbeam_4->SetLogz();
994 //gbeamyz->SetMarkerStyle(7);
995 hbeamyz->Draw("colz");
996 hbeamyz->GetXaxis()->SetTitle("Z [cm]");
997 hbeamyz->GetYaxis()->SetTitle("Y [cm]");
998 cbeam_4->Update();
999 TPaveStats *st2=(TPaveStats*)hbeamyz->GetListOfFunctions()->FindObject("stats");
1000 st2->SetX1NDC(0.13);
1001 st2->SetX2NDC(0.33);
1002 cbeam_4->Modified();
1003 cbeam_4->Update();
1004 cbeam->Update();
1005
1006 TCanvas *cbeam2 = new TCanvas("cbeam2","Beam Transv",500,500);
1007 cbeam2->SetLogz();
1008 cbeam2->SetRightMargin(0.14);
1009 //gbeamxy->SetMarkerStyle(7);
1010 hbeamxy->Draw("colz");
1011 hbeamxy->GetXaxis()->SetTitle("X [cm]");
1012 hbeamxy->GetYaxis()->SetTitle("Y [cm]");
1013 cbeam2->Update();
1014 TPaveStats *st3=(TPaveStats*)hbeamxy->GetListOfFunctions()->FindObject("stats");
1015 st3->SetX1NDC(0.13);
1016 st3->SetX2NDC(0.33);
1017 cbeam2->Modified();
1018 cbeam2->Update();
1019
1020 return;
1021}
1022//----------------------------------------------------------------------------
97805f7d 1023void ComputeVtxMean(TString vtxtype="TRK",
1024 TString fname="Vertex.Performance.root",
388ca814 1025 TString ntname="fNtupleVertexESD",
1026 Int_t nEventsToUse=10000,
97805f7d 1027 Int_t mincontr=1) {
388ca814 1028 //-----------------------------------------------------------------------
1029 // Compute weighted mean and cov. matrix from the ntuple
1030 //-----------------------------------------------------------------------
1031 gStyle->SetOptStat(0);
1032 //gStyle->SetOptFit(0);
1033
1034 Double_t diamondx=0.0200.,diamondy=0.0200.,diamondz=7.5.;
1035
1036 Double_t avx=0.;
1037 Double_t avy=0.;
1038 Double_t avz=0.;
1039 Double_t wgtavx=0.;
1040 Double_t wgtavy=0.;
1041 Double_t wgtavz=0.;
1042 Double_t sum=0.;
1043 Double_t sumwgtx=0.;
1044 Double_t sumwgty=0.;
1045 Double_t sumwgtz=0.;
1046 Double_t rmsx=0;
1047 Double_t rmsy=0;
1048 Double_t rmsz=0;
1049 Double_t varx=0.;
1050 Double_t vary=0.;
1051 Double_t varz=0.;
1052 Double_t covxy=0.;
1053 Double_t covxz=0.;
1054 Double_t covyz=0.;
1055 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1056
97805f7d 1057 TH1F* hx = new TH1F("hx","",200,-0.1,0.1);
388ca814 1058 hx->SetXTitle("vertex x [#mu m]");
1059 hx->SetYTitle("events");
97805f7d 1060 TH1F* hy = new TH1F("hy","",200,-0.1,0.1);
388ca814 1061 hy->SetXTitle("vertex y [#mu m]");
1062 hy->SetYTitle("events");
1063 TH1F* hz = new TH1F("hz","",200,-20,20);
1064 hz->SetXTitle("vertex z [cm]");
1065 hz->SetYTitle("events");
1066
1067
1068 TFile *f=new TFile(fname.Data());
1069 TList *cOutput = (TList*)f->Get("cOutput");
1070 TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
97805f7d 1071 Int_t nnnev=nt->GetEntries();
1072 printf("Events = %d\n",nnnev);
388ca814 1073 Float_t xVtx,xdiffVtx,xerrVtx;
1074 Float_t yVtx,ydiffVtx,yerrVtx;
1075 Float_t zVtx,zdiffVtx,zerrVtx;
1076 Float_t ntrklets,ncontrVtx,dndy,triggered,vtx3D;
1077 Float_t ztrue,zref;
1078
1079 TString sxx="x"; sxx.Append(vtxtype.Data());
1080 nt->SetBranchAddress(sxx.Data(),&xVtx);
1081 TString syy="y"; syy.Append(vtxtype.Data());
1082 nt->SetBranchAddress(syy.Data(),&yVtx);
1083 TString szz="z"; szz.Append(vtxtype.Data());
1084 nt->SetBranchAddress(szz.Data(),&zVtx);
1085
388ca814 1086 TString xerr="xerr"; xerr.Append(vtxtype.Data());
1087 nt->SetBranchAddress(xerr.Data(),&xerrVtx);
1088 TString yerr="yerr"; yerr.Append(vtxtype.Data());
1089 nt->SetBranchAddress(yerr.Data(),&yerrVtx);
1090 TString zerr="zerr"; zerr.Append(vtxtype.Data());
1091 nt->SetBranchAddress(zerr.Data(),&zerrVtx);
1092
1093 TString trkstitle;
1094 if(vtxtype.Contains("TPC")) {
1095 nt->SetBranchAddress("nTPCin",&ntrklets);
1096 trkstitle="TPC tracks pointing to beam pipe";
1097 } else {
1098 nt->SetBranchAddress("ntrklets",&ntrklets);
1099 trkstitle="SPD tracklets";
1100 }
1101 TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
1102 nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
1103 nt->SetBranchAddress("dndygen",&dndy);
1104
1105 nt->SetBranchAddress("ztrue",&ztrue);
1106
1107 nt->SetBranchAddress("triggered",&triggered);
1108
1109 nt->SetBranchAddress("SPD3D",&vtx3D);
1110
1111 Int_t total=0;
1112
97805f7d 1113 Int_t divider=(Int_t)(nnnev/nEventsToUse);
388ca814 1114 // first loop on events
97805f7d 1115 for(Int_t iev=0;iev<nnnev;iev++) {
388ca814 1116 if(iev%divider!=0) continue;
1117 total++;
1118 nt->GetEvent(iev);
1119 if(!vtxtype.Contains("SPD")) vtx3D=1.;
97805f7d 1120 //if(vtx3D<0.5) continue;
388ca814 1121 if(triggered<0.5) continue; // not triggered
1122 if(ncontrVtx<=0) continue; // no vertex
1123
1124 if(ncontrVtx<mincontr) continue;
1125
1126 avx += xVtx;
1127 avy += yVtx;
1128 avz += zVtx;
1129 sum += 1.;
1130 wgtavx += xVtx/xerrVtx/xerrVtx;
1131 wgtavy += yVtx/yerrVtx/yerrVtx;
1132 wgtavz += zVtx/zerrVtx/zerrVtx;
1133 sumwgtx += 1./xerrVtx/xerrVtx;
1134 sumwgty += 1./yerrVtx/yerrVtx;
1135 sumwgtz += 1./zerrVtx/zerrVtx;
1136
1137 hx->Fill(xVtx);
1138 hy->Fill(yVtx);
1139 hz->Fill(zVtx);
1140 }
1141
1142 avx /= sum;
1143 avy /= sum;
1144 avz /= sum;
1145 wgtavx /= sumwgtx;
1146 wgtavy /= sumwgty;
1147 wgtavz /= sumwgtz;
1148 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1149 ewgtavy = 1./TMath::Sqrt(sumwgty);
1150 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1151
1152
1153 // second loop on events
97805f7d 1154 for(Int_t iev=0;iev<nnnev;iev++){
388ca814 1155 if(iev%divider!=0) continue;
1156 nt->GetEvent(iev);
1157 if(!vtxtype.Contains("SPD")) vtx3D=1.;
1158 if(vtx3D<0.5) continue;
1159 if(triggered<0.5) continue; // not triggered
1160 if(ncontrVtx<=0) continue; // no vertex
1161
1162 if(ncontrVtx<mincontr) continue;
1163
1164 varx += (xVtx-avx)*(xVtx-avx);
1165 vary += (yVtx-avy)*(yVtx-avy);
1166 varz += (zVtx-avz)*(zVtx-avz);
1167 covxy += (xVtx-avx)*(yVtx-avy);
1168 covxz += (xVtx-avx)*(zVtx-avz);
1169 covyz += (yVtx-avy)*(zVtx-avz);
1170 }
1171
1172 varx /= sum;
1173 vary /= sum;
1174 varz /= sum;
1175 covxy /= sum;
1176 covxz /= sum;
1177 covyz /= sum;
1178 rmsx = TMath::Sqrt(varx);
1179 rmsy = TMath::Sqrt(vary);
1180 rmsz = TMath::Sqrt(varz);
1181 eavx = rmsx/TMath::Sqrt(sum);
1182 eavy = rmsy/TMath::Sqrt(sum);
1183 eavz = rmsz/TMath::Sqrt(sum);
1184
1185
1186 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1187 printf("Minimum number of contributors: %d\n",mincontr);
1188 printf("Average:\n x = (%f +- %f) cm\n y = (%f +- %f) cm\n z = (%f +- %f) cm\n",avx,eavx,avy,eavy,avz,eavz);
1189 printf("Weighted Average:\n x = (%f +- %f) cm\n y = (%f +- %f) cm\n z = (%f +- %f) cm\n",wgtavx,ewgtavx,wgtavy,ewgtavy,wgtavz,ewgtavz);
1190 printf("RMS:\n x = %f cm\n y = %f cm\n z = %f cm\n",rmsx,rmsy,rmsz);
1191
1192 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1193 c->Divide(3,1);
1194 c->cd(1);
1195 hx->Draw();
1196 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1197 gx->SetLineColor(2);
1198 hx->Fit(gx,"Q");
1199 TF1 *gxx = (TF1*)gx->Clone("gxx");
1200 gxx->FixParameter(2,diamondx);
1201 gxx->SetLineStyle(2);
1202 gxx->Draw("same");
1203 c->cd(2);
1204 hy->Draw();
1205 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1206 gy->SetLineColor(2);
1207 hy->Fit(gy,"Q");
1208 TF1 *gyy = (TF1*)gy->Clone("gyy");
1209 gyy->FixParameter(2,diamondy);
1210 gyy->SetLineStyle(2);
1211 gyy->Draw("same");
1212 c->cd(3);
1213 hz->Draw();
1214 TF1 *gz = new TF1("gz","gaus",-10,10);
1215 gz->SetLineColor(2);
1216 hz->Fit(gz,"Q");
1217 TF1 *gzz = (TF1*)gz->Clone("gzz");
1218 gzz->FixParameter(2,diamondz);
1219 gzz->SetLineStyle(2);
1220 gzz->Draw("same");
1221
1222
1223 return;
1224}