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