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