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 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 }