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="Vertex.Performance.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 xtrue,ytrue,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 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
 296   nt->SetBranchAddress("xtrue",&xtrue);
 297   nt->SetBranchAddress("ytrue",&ytrue);
 298   nt->SetBranchAddress("ztrue",&ztrue);
 299
 300   nt->SetBranchAddress("triggered",&triggered);
 301
 302   nt->SetBranchAddress("SPD3D",&vtx3D);
 303
 304
 305   // loop on events
 306   for(Int_t iev=0;iev<nnnev;iev++){
 307     nt->GetEvent(iev);
 308
 309     xdiffVtx=10000.*(xVtx-xtrue);
 310     ydiffVtx=10000.*(yVtx-ytrue);
 311     zdiffVtx=10000.*(zVtx-ztrue);
 312
 313
 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   }
 387   totev+=nnnev;
 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 //----------------------------------------------------------------------------
1023 void ComputeVtxMean(TString vtxtype="TRK",
1024                     TString fname="Vertex.Performance.root",
1025                     TString ntname="fNtupleVertexESD",
1026                     Int_t nEventsToUse=10000,
1027                     Int_t mincontr=1) {
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
1057   TH1F* hx = new TH1F("hx","",200,-0.1,0.1);
1058   hx->SetXTitle("vertex x [#mu m]");
1059   hx->SetYTitle("events");
1060   TH1F* hy = new TH1F("hy","",200,-0.1,0.1);
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());
1071   Int_t nnnev=nt->GetEntries();
1072   printf("Events = %d\n",nnnev);
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
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
1113   Int_t divider=(Int_t)(nnnev/nEventsToUse);
1114   // first loop on events
1115   for(Int_t iev=0;iev<nnnev;iev++) {
1116     if(iev%divider!=0) continue;
1117     total++;
1118     nt->GetEvent(iev);
1119     if(!vtxtype.Contains("SPD")) vtx3D=1.;
1120     //if(vtx3D<0.5) continue;
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
1154   for(Int_t iev=0;iev<nnnev;iev++){
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 }