New task for primary vertex analysis
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 9 Nov 2009 18:20:33 +0000 (18:20 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 9 Nov 2009 18:20:33 +0000 (18:20 +0000)
PWG1/CMake_libPWG1.txt
PWG1/PWG1LinkDef.h
PWG1/global/AliAnalysisTaskVertexESD.cxx [new file with mode: 0644]
PWG1/global/AliAnalysisTaskVertexESD.h [new file with mode: 0644]
PWG1/macros/AddTaskVertexESD.C [new file with mode: 0644]
PWG1/macros/PlotVertexESD.C [new file with mode: 0644]

index 72ee7eb..317274b 100644 (file)
@@ -30,6 +30,7 @@ set(SRCS
        AliPerformanceDEdx.cxx
        AliPerformanceDCA.cxx
        AliPerformanceTPC.cxx
+       AliAnalysisTaskVertexESD.cxx
        AliAlignmentDataFilterITS.cxx
        AliTrackMatchingTPCITSCosmics.cxx
         AliRelAlignerKalmanArray.cxx
index 94649ac..36571fb 100644 (file)
@@ -52,6 +52,7 @@
 #pragma link C++ class AliPerformancePtCalibMC+;
 #pragma link C++ class AliPerfAnalyzeInvPt+;
 
+#pragma link C++ class AliAnalysisTaskVertexESD+;
 #pragma link C++ class AliAlignmentDataFilterITS+;
 #pragma link C++ class AliTrackMatchingTPCITSCosmics+;
 #pragma link C++ class AliAnalysisTaskV0QA+;
diff --git a/PWG1/global/AliAnalysisTaskVertexESD.cxx b/PWG1/global/AliAnalysisTaskVertexESD.cxx
new file mode 100644 (file)
index 0000000..5cdbb34
--- /dev/null
@@ -0,0 +1,384 @@
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//*************************************************************************
+// Class AliAnalysisTaskVertexESD
+// AliAnalysisTask to extract from ESD the information for the analysis
+// of the primary vertex reconstruction efficiency and resolution
+// (for MC events) and distributions (for real data). Three vertices:
+// - SPD tracklets
+// - ITS+TPC tracks
+// - TPC-only tracks
+//
+// Author: A.Dainese, andrea.dainese@pd.infn.it
+//*************************************************************************
+
+#include <TChain.h>
+#include <TTree.h>
+#include <TNtuple.h>
+#include <TBranch.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
+#include <TH1F.h>
+#include <TH2F.h>  
+#include <TCanvas.h>
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliMultiplicity.h"
+#include "AliVertexerTracks.h"
+#include "AliESDtrack.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDfriend.h"
+#include "AliESDInputHandler.h"
+#include "AliTrackPointArray.h"
+#include "AliTrackReference.h"
+
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliLog.h"
+
+#include "AliGenEventHeader.h" 
+#include "AliAnalysisTaskVertexESD.h"
+
+
+ClassImp(AliAnalysisTaskVertexESD)
+
+//________________________________________________________________________
+AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : 
+AliAnalysisTask(name, "VertexESDTask"), 
+fReadMC(kFALSE),
+fESD(0), 
+fESDfriend(0), 
+fOutput(0), 
+fNtupleVertexESD(0),
+fhTrackPoints(0),
+fhTrackRefs(0)
+{
+  // Constructor
+
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TList container
+  DefineOutput(0, TList::Class());  //My private output
+}
+//________________________________________________________________________
+AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
+{
+  // Destructor
+
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *) 
+{
+  // Connect ESD or AOD here
+  // Called once
+
+  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+  if (!tree) {
+    Printf("ERROR: Could not read chain from input slot 0");
+  } else {
+    // Disable all branches and enable only the needed ones
+    // The next two lines are different when data produced as AliESDEvent is read
+    //    tree->SetBranchStatus("*", kFALSE);
+    tree->SetBranchStatus("fTriggerMask", 1);
+    tree->SetBranchStatus("fSPDVertex*", 1);
+    tree->SetBranchStatus("fSPDMult*", 1);
+    
+    tree->SetBranchStatus("ESDfriend*", 1);
+    tree->SetBranchAddress("ESDfriend.",&fESDfriend);
+
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if (!esdH) {
+      Printf("ERROR: Could not get ESDInputHandler");
+    } else fESD = esdH->GetEvent();
+  }
+  
+  return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskVertexESD::CreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+
+  // Several histograms are more conveniently managed in a TList
+  fOutput = new TList;
+  fOutput->SetOwner();
+
+  fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","xtrue:ytrue:ztrue:xSPD:xdiffSPD:xerrSPD:ySPD:ydiffSPD:yerrSPD:zSPD:zdiffSPD:zerrSPD:ntrksSPD:xTPC:xdiffTPC:xerrTPC:yTPC:ydiffTPC:yerrTPC:zTPC:zdiffTPC:zerrTPC:ntrksTPC:xTRK:xdiffTRK:xerrTRK:yTRK:ydiffTRK:yerrTRK:zTRK:zdiffTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK");
+
+  fOutput->Add(fNtupleVertexESD);
+
+  fhTrackPoints = new TH2F("fhTrackPoints","Track points; x; y",1000,-4,4,1000,-4,4);
+  fOutput->Add(fhTrackPoints);
+  fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
+  fOutput->Add(fhTrackRefs);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskVertexESD::Exec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  
+  if (!fESD) {
+    Printf("ERROR: fESD not available");
+    return;
+  }
+
+  
+  TArrayF mcVertex(3);
+  mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
+  Float_t dNchdy=-999.;
+
+  // ***********  MC info ***************
+  if (fReadMC) {
+    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+    if (!eventHandler) {
+      Printf("ERROR: Could not retrieve MC event handler");
+      return;
+    }
+    
+    AliMCEvent* mcEvent = eventHandler->MCEvent();
+    if (!mcEvent) {
+      Printf("ERROR: Could not retrieve MC event");
+      return;
+    }
+    
+    AliStack* stack = mcEvent->Stack();
+    if (!stack) {
+      AliDebug(AliLog::kError, "Stack not available");
+      return;
+    }
+    
+    AliHeader* header = mcEvent->Header();
+    if (!header) {
+      AliDebug(AliLog::kError, "Header not available");
+      return;
+    }
+    AliGenEventHeader* genHeader = header->GenEventHeader();
+    genHeader->PrimaryVertex(mcVertex);
+
+    Int_t ngenpart = (Int_t)stack->GetNtrack();
+    //printf("# generated particles = %d\n",ngenpart);
+    dNchdy=0;
+    for(Int_t ip=0; ip<ngenpart; ip++) {
+      TParticle* part = (TParticle*)stack->Particle(ip);
+      // keep only electorns, muons, pions, kaons and protons
+      Int_t apdg = TMath::Abs(part->GetPdgCode());
+      if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      
+      // reject secondaries
+      if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
+      // reject incoming protons
+      Double_t energy  = part->Energy();
+      if(energy>900.) continue;
+      Double_t pz = part->Pz();
+      Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
+      if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
+      // tracks refs
+      TClonesArray *trefs=0;
+      Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
+      if(ntrefs<0) continue;
+      for(Int_t iref=0; iref<ntrefs; iref++) {
+       AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
+       if(tref->R()>10.) continue;
+       fhTrackRefs->Fill(tref->X(),tref->Y());
+      }
+    }
+    //printf("# primary particles = %7.1f\n",dNchdy);
+  } 
+  // ***********  MC info ***************
+
+  // ***********  ESD friends ***********
+  fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD
+  // ***********  ESD friends ***********
+
+  //
+    
+  // Trigger
+  ULong64_t triggerMask;
+  ULong64_t spdFO = (1 << 14);
+  ULong64_t v0left = (1 << 11);
+  ULong64_t v0right = (1 << 12);
+  
+  triggerMask=fESD->GetTriggerMask();
+  // MB1: SPDFO || V0L || V0R
+  Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
+  //MB2: GFO && V0R
+  //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
+
+  Int_t ntracks = fESD->GetNumberOfTracks();
+  Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
+  //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
+  for(Int_t itr=0; itr<ntracks; itr++) {
+    AliESDtrack *t = fESD->GetTrack(itr);
+    if(t->GetNcls(0)>=5) nITS5or6++;
+    Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
+    if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
+      nTPCin++;
+      if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
+    }
+    // tracks with two SPD points
+    if(!TESTBIT(t->GetITSClusterMap(),0) ||
+       !TESTBIT(t->GetITSClusterMap(),1)) continue; 
+    // AliTrackPoints
+    const AliTrackPointArray *array = t->GetTrackPointArray();
+    AliTrackPoint point;
+    for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) {
+      array->GetPoint(point,ipt);
+      fhTrackPoints->Fill(point.GetX(),point.GetY());
+    }
+  }
+
+    
+  const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
+  const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
+  const AliESDVertex *trkv=fESD->GetPrimaryVertex();
+
+  //Float_t tpccontrorig=tpcv->GetNContributors();
+  //tpcv = 0;
+  //tpcv = ReconstructPrimaryVertexTPC();
+
+  const AliMultiplicity *alimult = fESD->GetMultiplicity();
+  Int_t ntrklets=0,spd0cls=0;
+  if(alimult) {
+    ntrklets = alimult->GetNumberOfTracklets();
+    for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
+      if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
+    }
+    spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
+  }
+
+
+  Float_t xnt[43];
+  
+  xnt[0]=mcVertex[0];
+  xnt[1]=mcVertex[1];
+  xnt[2]=mcVertex[2];
+  
+  xnt[3]=spdv->GetXv();
+  xnt[4]=(spdv->GetXv()-mcVertex[0])*10000.;
+  xnt[5]=spdv->GetXRes();
+  xnt[6]=spdv->GetYv();
+  xnt[7]=(spdv->GetYv()-mcVertex[1])*10000.;
+  xnt[8]=spdv->GetYRes();
+  xnt[9]=spdv->GetZv();
+  xnt[10]=(spdv->GetZv()-mcVertex[2])*10000.;
+  xnt[11]=spdv->GetZRes();
+  xnt[12]=spdv->GetNContributors();
+  
+  xnt[13]=tpcv->GetXv();
+  xnt[14]=(tpcv->GetXv()-mcVertex[0])*10000.;
+  xnt[15]=tpcv->GetXRes();
+  xnt[16]=tpcv->GetYv();
+  xnt[17]=(tpcv->GetYv()-mcVertex[1])*10000.;
+  xnt[18]=tpcv->GetYRes();
+  xnt[19]=tpcv->GetZv();
+  xnt[20]=(tpcv->GetZv()-mcVertex[2])*10000.;
+  xnt[21]=tpcv->GetZRes();
+  xnt[22]=tpcv->GetNContributors();
+  
+  xnt[23]=trkv->GetXv();
+  xnt[24]=(trkv->GetXv()-mcVertex[0])*10000.;
+  xnt[25]=trkv->GetXRes();
+  xnt[26]=trkv->GetYv();
+  xnt[27]=(trkv->GetYv()-mcVertex[1])*10000.;
+  xnt[28]=trkv->GetYRes();
+  xnt[29]=trkv->GetZv();
+  xnt[30]=(trkv->GetZv()-mcVertex[2])*10000.;
+  xnt[31]=trkv->GetZRes();
+  xnt[32]=trkv->GetNContributors();// tpccontrorig;
+  
+
+  xnt[33]=float(ntrklets);
+  xnt[34]=float(ntracks);
+  xnt[35]=float(nITS5or6);
+  xnt[36]=float(nTPCin);
+  xnt[37]=float(nTPCinEta09);
+
+  xnt[38]=float(dNchdy);
+
+  xnt[39]=0.; 
+  if(eventTriggered) xnt[39]=1.;
+
+  xnt[40]=0.; 
+  TString spdtitle = spdv->GetTitle();
+  if(spdtitle.Contains("vertexer: 3D")) xnt[40]=1.;
+
+  xnt[42]=0.; 
+  TString trktitle = trkv->GetTitle();
+  if(trktitle.Contains("WithConstraint")) xnt[42]=1.;
+
+  xnt[41]=spd0cls;
+
+  fNtupleVertexESD->Fill(xnt);
+    
+  // Post the data already here
+  PostData(0, fOutput);
+
+  return;
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+  fOutput = dynamic_cast<TList*> (GetOutputData(0));
+  if (!fOutput) {     
+    Printf("ERROR: fOutput not available");
+    return;
+  }
+
+  fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
+
+  return;
+}
+
+//_________________________________________________________________________
+AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() {
+  // On the fly reco of TPC vertex from ESD
+
+  AliVertexerTracks vertexer(fESD->GetMagneticField());
+  vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
+  //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
+  Double_t pos[3]={+0.0220,-0.0340,+0.270}; 
+  Double_t err[3]={0.0200,0.0200,7.5};
+  AliESDVertex *initVertex = new AliESDVertex(pos,err);
+  vertexer.SetVtxStart(initVertex);
+  delete initVertex;
+  //vertexer.SetConstraintOff();
+
+  return vertexer.FindPrimaryVertex(fESD);
+}
diff --git a/PWG1/global/AliAnalysisTaskVertexESD.h b/PWG1/global/AliAnalysisTaskVertexESD.h
new file mode 100644 (file)
index 0000000..2d35575
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef AliAnalysisTaskVertexESD_cxx
+#define AliAnalysisTaskVertexESD_cxx
+
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//*************************************************************************
+// Class AliAnalysisTaskVertexESD
+// AliAnalysisTask to extract from ESD the information for the analysis
+// of the primary vertex reconstruction efficiency and resolution
+// (for MC events) and distributions (for real data). Three vertices:
+// - SPD tracklets
+// - ITS+TPC tracks
+// - TPC-only tracks
+//
+// Author: A.Dainese, andrea.dainese@pd.infn.it
+//*************************************************************************
+
+class TNtuple;
+class TH2F;
+class AliESDEvent;
+class AliESDfriend;
+class AliESDVertex;
+
+#include "AliAnalysisTask.h"
+
+class AliAnalysisTaskVertexESD : public AliAnalysisTask 
+{
+ public:
+
+  AliAnalysisTaskVertexESD(const char *name = "AliAnalysisTaskVertexESD");
+  virtual ~AliAnalysisTaskVertexESD(); 
+  
+  virtual void   ConnectInputData(Option_t *);
+  virtual void   CreateOutputObjects();
+  virtual void   Exec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  Bool_t         GetReadMC() const { return fReadMC; }
+  void           SetReadMC(Bool_t flag=kTRUE) { fReadMC=flag; }
+  
+ protected:
+  Bool_t       fReadMC;         // read Monte Carlo
+  AliESDEvent *fESD;            // ESD object
+  AliESDfriend *fESDfriend;     // ESD friend object
+  TList       *fOutput;         //! list send on output slot 0
+  TNtuple     *fNtupleVertexESD;//! output ntuple
+  TH2F        *fhTrackPoints;   //! output histo
+  TH2F        *fhTrackRefs;     //! output histo
+
+ private:    
+
+  AliAnalysisTaskVertexESD(const AliAnalysisTaskVertexESD&); // not implemented
+  AliAnalysisTaskVertexESD& operator=(const AliAnalysisTaskVertexESD&); // not implemented
+  AliESDVertex* ReconstructPrimaryVertexTPC();
+  
+  ClassDef(AliAnalysisTaskVertexESD,1); // primary vertex analysis
+};
+
+#endif
diff --git a/PWG1/macros/AddTaskVertexESD.C b/PWG1/macros/AddTaskVertexESD.C
new file mode 100644 (file)
index 0000000..56f78bb
--- /dev/null
@@ -0,0 +1,37 @@
+AliAnalysisTaskVertexESD *AddTaskVertexESD() 
+{
+  //
+  // Task for validation of the primary vertices (SPD,TPC,ITS+TPC)
+  //
+  // andrea.dainese@pd.infn.it
+  //
+
+
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTask", "No analysis manager to connect to.");
+    return NULL;
+  }   
+
+  // Create the task
+  AliAnalysisTaskVertexESD *taskVtxESD = new AliAnalysisTaskVertexESD("VertexESD");
+  AliLog::SetClassDebugLevel("AliAnalysisTaskVertexESD",10);
+  // Add to the manager
+  mgr->AddTask(taskVtxESD);
+
+  //
+  // Create containers for input/output
+  AliAnalysisDataContainer *cInputVtxESD = mgr->CreateContainer("cInputVtxESD",TChain::Class(),AliAnalysisManager::kInputContainer);
+
+  AliAnalysisDataContainer *cOutputVtxESD = mgr->CreateContainer("cOutputVtxESD",TList::Class(),AliAnalysisManager::kOutputContainer,"VertexESD.root");
+
+
+  // Attach input
+  mgr->ConnectInput(taskVtxESD,0,mgr->GetCommonInputContainer());
+  // Attach output
+  mgr->ConnectOutput(taskVtxESD,0,cOutputVtxESD);
+  
+  return taskVtxESD;
+}
diff --git a/PWG1/macros/PlotVertexESD.C b/PWG1/macros/PlotVertexESD.C
new file mode 100644 (file)
index 0000000..e658e62
--- /dev/null
@@ -0,0 +1,1230 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TFile.h>
+#include <TLine.h>
+#include <TNtuple.h>
+#include <TStyle.h>
+#include <TString.h>
+#include <TSystem.h>
+#include <TH1F.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+#endif
+
+void PlotVertexESD(TString vtxtype="SPD",
+                  TString fname="VertexESD.root",
+                  TString ntname="fNtupleVertexESD",
+                  Bool_t useztrue=kTRUE,
+                  Int_t optgif=0) {
+  //-------------------------------------------------------------------------
+  // 
+  // Plot output of AliAnalysisTaskVertexESD.
+  // Origin: francesco.prino@to.infn.it
+  //
+  //-------------------------------------------------------------------------
+
+  // ranges for ideal geometry
+  /*
+  Float_t rangeHPull=40.;
+  Float_t rangeHRes=2500.;
+  Int_t binsHRes=250;
+  Float_t rangeGrAve = 100; // micron
+  Float_t rangeGrRms = 900; // micron
+  Float_t rangeGrPull = 5; 
+  */
+
+  // ranges for misaligned geometry
+  Float_t rangeHPull=40.;
+  Float_t rangeHRes=10000.*0.5;
+  Int_t binsHRes=250*0.5;
+  Float_t rangeGrAve = 10000; // micron
+  Float_t rangeGrRms = 5000; // micron
+  Float_t rangeGrPull = 15; 
+
+
+  Float_t limcont[10]={0.,2.5,3.5,4.5,6.5,9.5,14.5,20.5,99999.};
+  const Int_t nbinsmult=7;
+  Float_t limmult[nbinsmult+1]={0.,10.,14.,24.,30.,44.,60.,99999.};
+  //const Int_t nbinseff=11;
+  //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.};
+  const Int_t nbinseff=9;
+  Float_t limeff[nbinseff+1]={0.,0.5,1.5,2.5,3.5,4.5,5.5,7.5,10.5,999999.};
+  const Int_t nbinz = 14;
+  Float_t limitzt[nbinz+1]={-20,-18,-15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18,20};
+  //Float_t limitzt[nbinz+1]={-3,0.,3.,6.,9.,11.,13.,15,17,19,21};
+
+  TH1F **hxm=new TH1F*[nbinsmult];
+  TH1F **hym=new TH1F*[nbinsmult];
+  TH1F **hzm=new TH1F*[nbinsmult];
+  TH1F **hxpullm=new TH1F*[nbinsmult];
+  TH1F **hypullm=new TH1F*[nbinsmult];
+  TH1F **hzpullm=new TH1F*[nbinsmult];
+  TH1F **hmult=new TH1F*[nbinsmult];
+
+  TH1F **hxc=new TH1F*[8];
+  TH1F **hyc=new TH1F*[8];
+  TH1F **hzc=new TH1F*[8];
+  TH1F **hxpullc=new TH1F*[8];
+  TH1F **hypullc=new TH1F*[8];
+  TH1F **hzpullc=new TH1F*[8];
+  TH1F **hcont=new TH1F*[8];
+
+  TH1F **hxz=new TH1F*[nbinz];
+  TH1F **hyz=new TH1F*[nbinz];
+  TH1F **hzz=new TH1F*[nbinz];
+  TH1F **hz=new TH1F*[nbinz];
+
+  TH1F **hmulteff=new TH1F*[nbinseff];
+  TH1F **haux=new TH1F*[nbinseff];
+  TH1F *htot=new TH1F("htot","",100,-1000,1000);
+
+  TH1F **htrkseff=new TH1F*[nbinseff];
+  TH1F **htrksaux=new TH1F*[nbinseff];
+  TH1F **htrks3Daux=new TH1F*[nbinseff];
+
+  Char_t hisnam[10];
+  TCanvas *c1=new TCanvas("c1","Residuals",1000,700);
+  c1->Divide(nbinsmult,3,0.001,0.001);
+  TCanvas *c1z=new TCanvas("c1z","Residuals z",1000,700);
+  c1z->Divide(nbinz,3,0.001,0.001);
+  TCanvas *cp=new TCanvas("cp","Pulls",1000,700);
+  cp->Divide(nbinsmult,3,0.001,0.001);
+
+  TGraphErrors *gavexm=new TGraphErrors(0);
+  TGraphErrors *gaveym=new TGraphErrors(0);
+  TGraphErrors *gavezm=new TGraphErrors(0);
+  TGraphErrors *grmsxm=new TGraphErrors(0);
+  TGraphErrors *grmsym=new TGraphErrors(0);
+  TGraphErrors *grmszm=new TGraphErrors(0);
+  TGraphErrors *gavexmg=new TGraphErrors(0);
+  TGraphErrors *gaveymg=new TGraphErrors(0);
+  TGraphErrors *gavezmg=new TGraphErrors(0);
+  TGraphErrors *grmsxmg=new TGraphErrors(0);
+  TGraphErrors *grmsymg=new TGraphErrors(0);
+  TGraphErrors *grmszmg=new TGraphErrors(0);
+  TGraphErrors *gpullxm=new TGraphErrors(0);
+  TGraphErrors *gpullym=new TGraphErrors(0);
+  TGraphErrors *gpullzm=new TGraphErrors(0);
+  TGraphErrors *gpullxmg=new TGraphErrors(0);
+  TGraphErrors *gpullymg=new TGraphErrors(0);
+  TGraphErrors *gpullzmg=new TGraphErrors(0);
+  TGraphErrors *geffm=new TGraphErrors(0);
+  TGraphErrors *gefftrks=new TGraphErrors(0);
+  TGraphErrors *geff3Dtrks=new TGraphErrors(0);
+
+  TGraphErrors *gavexc=new TGraphErrors(0);
+  TGraphErrors *gaveyc=new TGraphErrors(0);
+  TGraphErrors *gavezc=new TGraphErrors(0);
+  TGraphErrors *grmsxc=new TGraphErrors(0);
+  TGraphErrors *grmsyc=new TGraphErrors(0);
+  TGraphErrors *grmszc=new TGraphErrors(0);
+  TGraphErrors *gpullxc=new TGraphErrors(0);
+  TGraphErrors *gpullyc=new TGraphErrors(0);
+  TGraphErrors *gpullzc=new TGraphErrors(0);
+
+  TGraph * gbeamxz=new TGraph(0);
+  TGraph * gbeamyz=new TGraph(0);
+  TGraph * gbeamxy=new TGraph(0);
+  TH2F *hbeamxz = new TH2F("hbeamxz","",100,-14.,14.,100,-1.5,1.5);
+  TH2F *hbeamyz = new TH2F("hbeamyz","",100,-14.,14.,100,-1.5,1.5);
+  TH1F *hbeamx = new TH1F("hbeamx","",1000,-1.5,1.5);
+  TH1F *hbeamy = new TH1F("hbeamy","",1000,-1.5,1.5);
+  TH2F *hbeamxy = new TH2F("hbeamxy","",100,-1.5,1.5,100,-1.5,1.5);
+
+
+
+  gavexm->SetName("gavexm");
+  grmsxm->SetName("grmsxm");
+  gavexmg->SetName("gavexmg");
+  grmsxmg->SetName("grmsxmg");
+  gpullxm->SetName("gpullxm");
+  gpullxmg->SetName("gpullxmg");
+  gaveym->SetName("gaveym");
+  grmsym->SetName("grmsym");
+  gaveymg->SetName("gaveymg");
+  grmsymg->SetName("grmsymg");
+  gpullym->SetName("gpullym");
+  gpullymg->SetName("gpullymg");
+  gavezm->SetName("gavezm");
+  grmszm->SetName("grmszm");
+  gavezmg->SetName("gavezmg");
+  grmszmg->SetName("grmszmg");
+  gpullzm->SetName("gpullzm");
+  gpullzmg->SetName("gpullzmg");
+  geffm->SetName("geffm");
+  gefftrks->SetName("gefftrks");
+  geff3Dtrks->SetName("geff3Dtrks");
+  gavexc->SetName("gavexc");
+  grmsxc->SetName("grmsxc");
+  gpullxc->SetName("gpullxc");
+  gaveyc->SetName("gaveyc");
+  grmsyc->SetName("grmsyc");
+  gpullyc->SetName("gpullyc");
+  gavezc->SetName("gavezc");
+  grmszc->SetName("grmszc");
+  gpullzc->SetName("gpullzc");
+
+  TGraphErrors *gavexz=new TGraphErrors(0);
+  TGraphErrors *gaveyz=new TGraphErrors(0);
+  TGraphErrors *gavezz=new TGraphErrors(0);
+  TGraphErrors *grmsxz=new TGraphErrors(0);
+  TGraphErrors *grmsyz=new TGraphErrors(0);
+  TGraphErrors *grmszz=new TGraphErrors(0);
+  TGraphErrors *gavexzg=new TGraphErrors(0);
+  TGraphErrors *gaveyzg=new TGraphErrors(0);
+  TGraphErrors *gavezzg=new TGraphErrors(0);
+  TGraphErrors *grmsxzg=new TGraphErrors(0);
+  TGraphErrors *grmsyzg=new TGraphErrors(0);
+  TGraphErrors *grmszzg=new TGraphErrors(0);
+  TGraphErrors *geffz=new TGraphErrors(0);
+
+  gavexz->SetName("gavexz");
+  grmsxz->SetName("grmsxz");
+  gavexzg->SetName("gavexzg");
+  grmsxzg->SetName("grmsxzg");
+  gaveyz->SetName("gaveyz");
+  grmsyz->SetName("grmsyz");
+  gaveyzg->SetName("gaveyzg");
+  grmsyzg->SetName("grmsyzg");
+  gavezz->SetName("gavezz");
+  grmszz->SetName("grmszz");
+  gavezzg->SetName("gavezzg");
+  grmszzg->SetName("grmszzg");
+  geffz->SetName("geffz");
+
+  TF1 *fitf;
+
+  // histogram creation
+  for(Int_t khis=0;khis<nbinseff;khis++){
+    sprintf(hisnam,"hmeff%d",khis);
+    hmulteff[khis]=new TH1F(hisnam,"",100,0.,200.);
+    sprintf(hisnam,"htrkseff%d",khis);
+    htrkseff[khis]=new TH1F(hisnam,"",100,0.,200.);
+    sprintf(hisnam,"haux%d",khis);
+    haux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"htrksaux%d",khis);
+    htrksaux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"htrks3Daux%d",khis);
+    htrks3Daux[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+  }
+  for(Int_t khis=0;khis<nbinsmult;khis++){
+    sprintf(hisnam,"hxm%d",khis);
+    hxm[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hym%d",khis);
+    hym[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hzm%d",khis);
+    hzm[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hxpm%d",khis);
+    hxpullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hypm%d",khis);
+    hypullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hzpm%d",khis);
+    hzpullm[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hmult%d",khis);
+    hmult[khis]=new TH1F(hisnam,"",100,0.,200.);
+  }
+
+  for(Int_t khis=0;khis<8;khis++){
+    sprintf(hisnam,"hxc%d",khis);
+    hxc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hyc%d",khis);
+    hyc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hzc%d",khis);
+    hzc[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hxpc%d",khis);
+    hxpullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hypc%d",khis);
+    hypullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hzpc%d",khis);
+    hzpullc[khis]=new TH1F(hisnam,"",100,-rangeHPull,rangeHPull);
+    sprintf(hisnam,"hcont%d",khis);
+    hcont[khis]=new TH1F(hisnam,"",100,0.,200.);
+  }
+
+  for(Int_t khis=0;khis<nbinz;khis++){
+    sprintf(hisnam,"hxz%d",khis);
+    hxz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hyz%d",khis);
+    hyz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hzz%d",khis);
+    hzz[khis]=new TH1F(hisnam,"",binsHRes,-rangeHRes,rangeHRes);
+    sprintf(hisnam,"hz%d",khis);
+    hz[khis]=new TH1F(hisnam,"",100,-20.,20.);   
+  }
+
+  Float_t totev=0,totevtriggered=0,nvtx3D=0,nvtxZ=0;
+  TFile *f=new TFile(fname.Data());
+  TList *cOutput = (TList*)f->Get("cOutput");
+  TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
+  Int_t NNNev=nt->GetEntries();
+  printf("Events = %d\n",NNNev);
+  Float_t xVtx,xdiffVtx,xerrVtx;
+  Float_t yVtx,ydiffVtx,yerrVtx;
+  Float_t zVtx,zdiffVtx,zerrVtx;
+  Float_t ntrklets,ncontrVtx,dndy,triggered,vtx3D;
+  Float_t ztrue,zref;
+  
+  TString sxx="x"; sxx.Append(vtxtype.Data());
+  nt->SetBranchAddress(sxx.Data(),&xVtx);
+  TString syy="y"; syy.Append(vtxtype.Data());
+  nt->SetBranchAddress(syy.Data(),&yVtx);
+  TString szz="z"; szz.Append(vtxtype.Data());
+  nt->SetBranchAddress(szz.Data(),&zVtx);
+  
+  TString xdiff="xdiff"; xdiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(xdiff.Data(),&xdiffVtx);
+  TString ydiff="ydiff"; ydiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(ydiff.Data(),&ydiffVtx);
+  TString zdiff="zdiff"; zdiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(zdiff.Data(),&zdiffVtx);
+    
+  TString xerr="xerr"; xerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(xerr.Data(),&xerrVtx);
+  TString yerr="yerr"; yerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(yerr.Data(),&yerrVtx);
+  TString zerr="zerr"; zerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(zerr.Data(),&zerrVtx);
+
+  TString trkstitle;
+  if(vtxtype.Contains("TPC")) {
+    nt->SetBranchAddress("nTPCin",&ntrklets);
+    trkstitle="TPC tracks pointing to beam pipe";
+  } else {
+    nt->SetBranchAddress("ntrklets",&ntrklets);
+    trkstitle="SPD tracklets";
+  }
+  TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
+  nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
+  nt->SetBranchAddress("dndygen",&dndy);
+  
+  nt->SetBranchAddress("ztrue",&ztrue);
+
+  nt->SetBranchAddress("triggered",&triggered);
+
+  nt->SetBranchAddress("SPD3D",&vtx3D);
+
+
+  // loop on events
+  for(Int_t iev=0;iev<NNNev;iev++){
+    nt->GetEvent(iev);
+    zref = (useztrue ? ztrue : zVtx);
+    if(!vtxtype.Contains("SPD")) vtx3D=1.;
+    if(triggered<0.5) continue; // not triggered
+    totevtriggered += 1.;
+    if(ncontrVtx>0) {
+      htot->Fill(zdiffVtx);
+      if(vtx3D>0.5) {
+       nvtx3D += 1.;
+      } else {
+       nvtxZ += 1.;
+      }
+    }
+
+    if(ncontrVtx>0 && vtx3D>0.5) {
+      gbeamxz->SetPoint(gbeamxz->GetN(),zVtx,xVtx);
+      gbeamyz->SetPoint(gbeamxz->GetN(),zVtx,yVtx);
+      gbeamxy->SetPoint(gbeamxz->GetN(),xVtx,yVtx);
+      
+      hbeamx->Fill(xVtx);
+      hbeamy->Fill(yVtx);
+      hbeamxz->Fill(zVtx,xVtx);
+      hbeamyz->Fill(zVtx,yVtx);
+      hbeamxy->Fill(xVtx,yVtx);
+    }
+
+    for(Int_t khis=0;khis<nbinseff;khis++){
+      if(dndy>=limeff[khis] && dndy<limeff[khis+1]){
+       hmulteff[khis]->Fill(dndy);
+       if(ncontrVtx>0) haux[khis]->Fill(zdiffVtx);
+      }
+      if(ntrklets>=limeff[khis] && ntrklets<limeff[khis+1]){
+       htrkseff[khis]->Fill(ntrklets);
+       if(ncontrVtx>0) htrksaux[khis]->Fill(zdiffVtx);
+       if(ncontrVtx>0 && vtx3D>0.5) htrks3Daux[khis]->Fill(zdiffVtx);
+      }
+    }
+    for(Int_t khis=0;khis<nbinsmult;khis++){
+      if(ntrklets>=limmult[khis] && ntrklets<limmult[khis+1]){
+       hmult[khis]->Fill(ntrklets);
+       if(ncontrVtx>0){
+         if(vtx3D>0.5) hxm[khis]->Fill(xdiffVtx);
+         if(vtx3D>0.5) hym[khis]->Fill(ydiffVtx);
+         hzm[khis]->Fill(zdiffVtx);
+         if(vtx3D>0.5) hxpullm[khis]->Fill(xdiffVtx/10000./xerrVtx);
+         if(vtx3D>0.5) hypullm[khis]->Fill(ydiffVtx/10000./yerrVtx);
+         hzpullm[khis]->Fill(zdiffVtx/10000./zerrVtx);
+       }
+      }
+    }
+    for(Int_t khis=0;khis<8;khis++){
+      if(ncontrVtx>=limcont[khis]&&ncontrVtx<limcont[khis+1]){
+       hcont[khis]->Fill(ncontrVtx);
+       if(ncontrVtx>0){        
+         if(vtx3D>0.5)hxc[khis]->Fill(xdiffVtx);
+         if(vtx3D>0.5)hyc[khis]->Fill(ydiffVtx);
+         hzc[khis]->Fill(zdiffVtx);
+         if(vtx3D>0.5)hxpullc[khis]->Fill(xdiffVtx/10000./xerrVtx);
+         if(vtx3D>0.5)hypullc[khis]->Fill(ydiffVtx/10000./yerrVtx);
+         hzpullc[khis]->Fill(zdiffVtx/10000./zerrVtx);
+       }
+      }
+    }
+    for(Int_t khis=0;khis<nbinz;khis++){
+      if(zref>=limitzt[khis]&&zref<limitzt[khis+1]){
+       hz[khis]->Fill(zref);
+       if(ncontrVtx>0){
+         if(vtx3D>0.5)hxz[khis]->Fill(xdiffVtx);
+         if(vtx3D>0.5)hyz[khis]->Fill(ydiffVtx);
+         hzz[khis]->Fill(zdiffVtx);
+       }
+      }
+    }
+  }
+  totev+=NNNev;
+
+  if(totev==0){
+    printf("Total number of events = 0\n");
+    return;
+  }
+  for(Int_t khis=0;khis<nbinseff;khis++){
+    Double_t x=hmulteff[khis]->GetMean();
+    Double_t ex=hmulteff[khis]->GetRMS();;
+    Float_t nEv=(Float_t)(hmulteff[khis]->GetEntries());
+    Float_t trkeff=-1.;
+    cout<<"Eff. dNch/dy bin "<<khis<<" ("<<limeff[khis]<<"-"<<limeff[khis+1]<<")  # Events ="<<nEv<<" with vertex ="<<haux[khis]->GetEntries()<<endl;
+    if(nEv>0) trkeff=(Float_t)(haux[khis]->GetEntries())/nEv;
+    geffm->SetPoint(khis,x,trkeff);
+    geffm->SetPointError(khis,ex,0.);
+  }
+  for(Int_t khis=0;khis<nbinseff;khis++){
+    Double_t x=htrkseff[khis]->GetMean();
+    Double_t ex=htrkseff[khis]->GetRMS();;
+    Float_t nEv=(Float_t)(htrkseff[khis]->GetEntries());
+    Float_t trkeff=-1.;
+    cout<<"Eff. trks bin "<<khis<<" ("<<limeff[khis]<<"-"<<limeff[khis+1]<<")  # Events ="<<nEv<<" with vertex ="<<htrksaux[khis]->GetEntries()<<endl;
+    if(nEv>0) trkeff=(Float_t)(htrksaux[khis]->GetEntries())/nEv;
+    gefftrks->SetPoint(khis,x,trkeff);
+    gefftrks->SetPointError(khis,ex,0.);
+    if(nEv>0) trkeff=(Float_t)(htrks3Daux[khis]->GetEntries())/nEv;
+    geff3Dtrks->SetPoint(khis,x,trkeff);
+    geff3Dtrks->SetPointError(khis,ex,0.);
+  }
+
+  for(Int_t khis=0;khis<nbinsmult;khis++){ 
+
+    c1->cd(khis+1);
+    if(hxm[khis]->GetEntries()<10) continue;
+    hxm[khis]->Draw();
+    hxm[khis]->Fit("gaus","Q0");
+    fitf= hxm[khis]->GetFunction("gaus");
+    Double_t avexg=fitf->GetParameter(1);
+    Double_t eavexg=fitf->GetParError(1);
+    Double_t rmsxg=fitf->GetParameter(2);
+    Double_t ermsxg=fitf->GetParError(2);
+    c1->cd(nbinsmult+khis+1);
+    hym[khis]->Draw();
+    hym[khis]->Fit("gaus","Q0");
+    fitf= hym[khis]->GetFunction("gaus");
+    Double_t aveyg=fitf->GetParameter(1);
+    Double_t eaveyg=fitf->GetParError(1);
+    Double_t rmsyg=fitf->GetParameter(2);
+    Double_t ermsyg=fitf->GetParError(2);
+    c1->cd(2*nbinsmult+khis+1);
+    hzm[khis]->Draw();
+    hzm[khis]->Fit("gaus","Q0");
+    fitf= hzm[khis]->GetFunction("gaus");
+    Double_t avezg=fitf->GetParameter(1);
+    Double_t eavezg=fitf->GetParError(1);
+    Double_t rmszg=fitf->GetParameter(2);
+    Double_t ermszg=fitf->GetParError(2);
+
+    cp->cd(khis+1);
+    hxpullm[khis]->Draw();
+    hxpullm[khis]->Fit("gaus","Q0");
+    fitf= hxpullm[khis]->GetFunction("gaus");
+    Double_t pullxg=fitf->GetParameter(2);
+    Double_t epullxg=fitf->GetParError(2);
+    cp->cd(nbinsmult+khis+1);
+    hypullm[khis]->Draw();
+    hypullm[khis]->Fit("gaus","Q0");
+    fitf= hypullm[khis]->GetFunction("gaus");
+    Double_t pullyg=fitf->GetParameter(2);
+    Double_t epullyg=fitf->GetParError(2);
+    cp->cd(2*nbinsmult+khis+1);
+    hzpullm[khis]->Draw();
+    hzpullm[khis]->Fit("gaus","Q0");
+    fitf= hzpullm[khis]->GetFunction("gaus");
+    Double_t pullzg=fitf->GetParameter(2);
+    Double_t epullzg=fitf->GetParError(2);
+
+    Double_t rmsxt=hxm[khis]->GetRMS();
+    Double_t rmsyt=hym[khis]->GetRMS();
+    Double_t rmszt=hzm[khis]->GetRMS();
+    Double_t ermsxt=hxm[khis]->GetRMSError();
+    Double_t ermsyt=hym[khis]->GetRMSError();
+    Double_t ermszt=hzm[khis]->GetRMSError();
+    Double_t avext=hxm[khis]->GetMean();
+    Double_t aveyt=hym[khis]->GetMean();
+    Double_t avezt=hzm[khis]->GetMean();
+    Double_t eavext=hxm[khis]->GetMeanError();
+    Double_t eaveyt=hym[khis]->GetMeanError();
+    Double_t eavezt=hzm[khis]->GetMeanError();
+    Double_t pullxt=hxpullm[khis]->GetRMS();
+    Double_t pullyt=hypullm[khis]->GetRMS();
+    Double_t pullzt=hzpullm[khis]->GetRMS();
+    Double_t epullxt=hxpullm[khis]->GetRMSError();
+    Double_t epullyt=hypullm[khis]->GetRMSError();
+    Double_t epullzt=hzpullm[khis]->GetRMSError();
+
+    Double_t x=hmult[khis]->GetMean();
+    if(hmult[khis]->GetEntries()==0) x=-1;
+    Double_t ex=hmult[khis]->GetRMS();;
+
+    gavexm->SetPoint(khis,x,avext);
+    gavexm->SetPointError(khis,ex,eavext);
+    gaveym->SetPoint(khis,x,aveyt);
+    gaveym->SetPointError(khis,ex,eaveyt);
+    gavezm->SetPoint(khis,x,avezt);
+    gavezm->SetPointError(khis,ex,eavezt);
+    grmsxm->SetPoint(khis,x,rmsxt);
+    grmsxm->SetPointError(khis,ex,ermsxt);
+    grmsym->SetPoint(khis,x,rmsyt);
+    grmsym->SetPointError(khis,ex,ermsyt);
+    grmszm->SetPoint(khis,x,rmszt);
+    grmszm->SetPointError(khis,ex,ermszt);
+
+    gavexmg->SetPoint(khis,x,avexg);
+    gavexmg->SetPointError(khis,ex,eavexg);
+    gaveymg->SetPoint(khis,x,aveyg);
+    gaveymg->SetPointError(khis,ex,eaveyg);
+    gavezmg->SetPoint(khis,x,avezg);
+    gavezmg->SetPointError(khis,ex,eavezg);
+    grmsxmg->SetPoint(khis,x,rmsxg);
+    grmsxmg->SetPointError(khis,ex,ermsxg);
+    grmsymg->SetPoint(khis,x,rmsyg);
+    grmsymg->SetPointError(khis,ex,ermsyg);
+    grmszmg->SetPoint(khis,x,rmszg);
+    grmszmg->SetPointError(khis,ex,ermszg);
+
+    gpullxm->SetPoint(khis,x,pullxt);
+    gpullxm->SetPointError(khis,ex,epullxt);
+    gpullym->SetPoint(khis,x,pullyt);
+    gpullym->SetPointError(khis,ex,epullyt);
+    gpullzm->SetPoint(khis,x,pullzt);
+    gpullzm->SetPointError(khis,ex,epullzt);
+
+    gpullxmg->SetPoint(khis,x,pullxg);
+    gpullxmg->SetPointError(khis,ex,epullxg);
+    gpullymg->SetPoint(khis,x,pullyg);
+    gpullymg->SetPointError(khis,ex,epullyg);
+    gpullzmg->SetPoint(khis,x,pullzg);
+    gpullzmg->SetPointError(khis,ex,epullzg);
+
+    Float_t nEv=hmult[khis]->GetEntries();
+    cout<<"Mult. bin "<<khis<<"  # Events ="<<nEv<<endl;
+  }
+
+  for(Int_t khis=0;khis<8;khis++){
+
+    Double_t rmsxt=hxc[khis]->GetRMS();
+    Double_t rmsyt=hyc[khis]->GetRMS();
+    Double_t rmszt=hzc[khis]->GetRMS();
+    Double_t ermsxt=hxc[khis]->GetRMSError();
+    Double_t ermsyt=hyc[khis]->GetRMSError();
+    Double_t ermszt=hzc[khis]->GetRMSError();
+    Double_t avext=hxc[khis]->GetMean();
+    Double_t aveyt=hyc[khis]->GetMean();
+    Double_t avezt=hzc[khis]->GetMean();
+    Double_t eavext=hxc[khis]->GetMeanError();
+    Double_t eaveyt=hyc[khis]->GetMeanError();
+    Double_t eavezt=hzc[khis]->GetMeanError();
+    Double_t pullxt=hxpullc[khis]->GetRMS();
+    Double_t pullyt=hypullc[khis]->GetRMS();
+    Double_t pullzt=hzpullc[khis]->GetRMS();
+    Double_t epullxt=hxpullc[khis]->GetRMSError();
+    Double_t epullyt=hypullc[khis]->GetRMSError();
+    Double_t epullzt=hzpullc[khis]->GetRMSError();
+
+    Double_t x=hcont[khis]->GetMean();
+    Double_t ex=hcont[khis]->GetRMS();;
+
+    gavexc->SetPoint(khis,x,avext);
+    gavexc->SetPointError(khis,ex,eavext);
+    gaveyc->SetPoint(khis,x,aveyt);
+    gaveyc->SetPointError(khis,ex,eaveyt);
+    gavezc->SetPoint(khis,x,avezt);
+    gavezc->SetPointError(khis,ex,eavezt);
+    grmsxc->SetPoint(khis,x,rmsxt);
+    grmsxc->SetPointError(khis,ex,ermsxt);
+    grmsyc->SetPoint(khis,x,rmsyt);
+    grmsyc->SetPointError(khis,ex,ermsyt);
+    grmszc->SetPoint(khis,x,rmszt);
+    grmszc->SetPointError(khis,ex,ermszt);
+
+    gpullxc->SetPoint(khis,x,pullxt);
+    gpullxc->SetPointError(khis,ex,epullxt);
+    gpullyc->SetPoint(khis,x,pullyt);
+    gpullyc->SetPointError(khis,ex,epullyt);
+    gpullzc->SetPoint(khis,x,pullzt);
+    gpullzc->SetPointError(khis,ex,epullzt);
+
+    Float_t nEv=hcont[khis]->GetEntries();
+    cout<<"Contrib. bin "<<khis<<"  # Events ="<<nEv<<endl;
+  }
+
+  for(Int_t khis=0; khis<nbinz; khis++){
+
+    Double_t rmsxt=hxz[khis]->GetRMS();
+    Double_t rmsyt=hyz[khis]->GetRMS();
+    Double_t rmszt=hzz[khis]->GetRMS();
+    Double_t ermsxt=hxz[khis]->GetRMSError();
+    Double_t ermsyt=hyz[khis]->GetRMSError();
+    Double_t ermszt=hzz[khis]->GetRMSError();
+    Double_t avext=hxz[khis]->GetMean();
+    Double_t aveyt=hyz[khis]->GetMean();
+    Double_t avezt=hzz[khis]->GetMean();
+    Double_t eavext=hxz[khis]->GetMeanError();
+    Double_t eaveyt=hyz[khis]->GetMeanError();
+    Double_t eavezt=hzz[khis]->GetMeanError();
+
+    Float_t nEv=hz[khis]->GetEntries();
+    Double_t x=-999.;
+    if(nEv>0) x=hz[khis]->GetMean();
+    Double_t ex=hz[khis]->GetRMS();;
+    
+    gavexz->SetPoint(khis,x,avext);
+    gavexz->SetPointError(khis,ex,eavext);
+    gaveyz->SetPoint(khis,x,aveyt);
+    gaveyz->SetPointError(khis,ex,eaveyt);
+    gavezz->SetPoint(khis,x,avezt);
+    gavezz->SetPointError(khis,ex,eavezt);
+    grmsxz->SetPoint(khis,x,rmsxt);
+    grmsxz->SetPointError(khis,ex,ermsxt);
+    grmsyz->SetPoint(khis,x,rmsyt);
+    grmsyz->SetPointError(khis,ex,ermsyt);
+    grmszz->SetPoint(khis,x,rmszt);
+    grmszz->SetPointError(khis,ex,ermszt);
+
+    c1z->cd(khis+1);
+    hxz[khis]->Draw();
+    hxz[khis]->Fit("gaus","Q0");
+    fitf= hxz[khis]->GetFunction("gaus");
+    Double_t avexg=fitf->GetParameter(1);
+    Double_t eavexg=fitf->GetParError(1);
+    Double_t rmsxg=fitf->GetParameter(2);
+    Double_t ermsxg=fitf->GetParError(2);
+    c1z->cd(1*nbinz+khis+1);
+    hyz[khis]->Draw();
+    hyz[khis]->Fit("gaus","Q0");
+    fitf= hyz[khis]->GetFunction("gaus");
+    Double_t aveyg=fitf->GetParameter(1);
+    Double_t eaveyg=fitf->GetParError(1);
+    Double_t rmsyg=fitf->GetParameter(2);
+    Double_t ermsyg=fitf->GetParError(2);
+    c1z->cd(2*nbinz+khis+1);
+    hzz[khis]->Draw();
+    hzz[khis]->Fit("gaus","Q0");
+    fitf= hzz[khis]->GetFunction("gaus");
+    Double_t avezg=fitf->GetParameter(1);
+    Double_t eavezg=fitf->GetParError(1);
+    Double_t rmszg=fitf->GetParameter(2);
+    Double_t ermszg=fitf->GetParError(2);
+
+    gavexzg->SetPoint(khis,x,avexg);
+    gavexzg->SetPointError(khis,ex,eavexg);
+    gaveyzg->SetPoint(khis,x,aveyg);
+    gaveyzg->SetPointError(khis,ex,eaveyg);
+    gavezzg->SetPoint(khis,x,avezg);
+    gavezzg->SetPointError(khis,ex,eavezg);
+    grmsxzg->SetPoint(khis,x,rmsxg);
+    grmsxzg->SetPointError(khis,ex,ermsxg);
+    grmsyzg->SetPoint(khis,x,rmsyg);
+    grmsyzg->SetPointError(khis,ex,ermsyg);
+    grmszzg->SetPoint(khis,x,rmszg);
+    grmszzg->SetPointError(khis,ex,ermszg);
+
+    Float_t zeff=-999.;
+    if(nEv>0) zeff=hzz[khis]->GetEntries()/nEv;
+    geffz->SetPoint(khis,x,zeff);
+    geffz->SetPointError(khis,ex,0.);
+    
+    cout<<"Z bin "<<khis<<"  # Events ="<<nEv<<endl;
+  }
+  Double_t efftrk=htot->GetEntries()/totevtriggered;
+
+  printf("EVENTS STATISTICS:\n Total: %d\n Triggered (MB1): %d\n Triggered and with vertex %d\n",totev,totevtriggered,htot->GetEntries());
+  if(vtxtype.Contains("SPD")) printf("  %d with Vertexer3D, %d with VertexerZ\n",nvtx3D,nvtxZ);
+  printf("Overall efficiency (for triggered evts) Vertexer%s = %f / %f = %f\n",vtxtype.Data(),htot->GetEntries(),totevtriggered,efftrk);
+
+  TFile* in = new TFile("vert-graphs.root","recreate");
+  gbeamxz->Write("gbeamxz");
+  gbeamyz->Write("gbeamyz");
+  grmsxm->Write();
+  grmsxmg->Write();
+  gpullxm->Write();
+  gpullxmg->Write();
+  gavexm->Write();
+  gavexmg->Write();
+  grmsym->Write();
+  grmsymg->Write();
+  gpullym->Write();
+  gpullymg->Write();
+  gaveym->Write();
+  gaveymg->Write();
+  grmszm->Write();
+  grmszmg->Write();
+  gpullzm->Write();
+  gpullzmg->Write();
+  gavezm->Write();
+  gavezmg->Write();
+  geffm->Write();
+  gefftrks->Write();
+  geff3Dtrks->Write();
+  grmsxc->Write();
+  gpullxc->Write();
+  gavexc->Write();
+  grmsyc->Write();
+  gpullyc->Write();
+  gaveyc->Write();
+  grmszc->Write();
+  gpullzc->Write();
+  gavezc->Write();
+  gavexz->Write();
+  gaveyz->Write();
+  gavezz->Write();
+  grmsxz->Write();
+  grmsyz->Write();
+  grmszz->Write();
+  gavexzg->Write();
+  gaveyzg->Write();
+  gavezzg->Write();
+  grmsxzg->Write();
+  grmsyzg->Write();
+  grmszzg->Write();
+
+  in->Close();
+
+  gStyle->SetOptTitle(0);
+
+  Char_t outgif[100];
+  TLine *lin0=new TLine(0,0,60,0);
+  lin0->SetLineStyle(2);
+
+  TCanvas *cg1=new TCanvas("cg1","Histo mean");
+  cg1->SetBottomMargin(0.14);
+  cg1->SetTopMargin(0.08);
+  cg1->SetLeftMargin(0.14);
+  cg1->SetRightMargin(0.08);
+  gavexm->SetMarkerStyle(20);
+  gavexm->SetMarkerColor(1);
+  gaveym->SetMarkerStyle(21);
+  gaveym->SetMarkerColor(2);
+  gavezm->SetMarkerStyle(22);
+  gavezm->SetMarkerColor(4);
+  TLegend *leg=new TLegend(0.18,0.70,0.25,0.90);
+  leg->SetBorderSize(0);
+  leg->SetFillStyle(0);
+  TLegendEntry *ent=leg->AddEntry(gavexm,"x","P");
+  ent->SetTextColor(1);
+  ent=leg->AddEntry(gaveym,"y","P");
+  ent->SetTextColor(2);
+  ent=leg->AddEntry(gavezm,"z","P");
+  ent->SetTextColor(4);
+  gavezm->GetXaxis()->SetLimits(0.,60.);
+  gavezm->SetMinimum(-rangeGrAve);
+  gavezm->SetMaximum(rangeGrAve);
+  gavezm->Draw("AP");
+  gavezm->GetXaxis()->SetTitle(trkstitle.Data());
+  gavezm->GetXaxis()->SetTitleSize(0.05);
+  gavezm->GetYaxis()->SetTitle("<Pos_{found}-Pos_{true}> [#mum]");
+  gavezm->GetYaxis()->SetTitleSize(0.05);
+  gavexm->Draw("PSAME");
+  gaveym->Draw("PSAME");
+  leg->Draw();
+  lin0->Draw();
+  sprintf(outgif,"vert%s-ave-mult.gif",vtxtype.Data());
+  if(optgif) cg1->SaveAs(outgif);
+
+
+  TCanvas *cg2=new TCanvas("cg2","Histo RMS");
+  cg2->SetBottomMargin(0.14);
+  cg2->SetTopMargin(0.08);
+  cg2->SetLeftMargin(0.14);
+  cg2->SetRightMargin(0.08);
+  grmsxm->SetMarkerStyle(20);
+  grmsxm->SetMarkerColor(1);
+  grmsym->SetMarkerStyle(21);
+  grmsym->SetMarkerColor(2);
+  grmszm->SetMarkerStyle(22);
+  grmszm->SetMarkerColor(4);
+  grmszm->SetMinimum(0);
+  grmszm->SetMaximum(rangeGrRms);
+  grmszm->GetXaxis()->SetLimits(0,60);
+  grmszm->Draw("AP");
+  grmszm->GetXaxis()->SetTitle(trkstitle.Data());
+  grmszm->GetXaxis()->SetTitleSize(0.05);
+  grmszm->GetYaxis()->SetTitle("Resolution [#mum]");
+  grmszm->GetYaxis()->SetTitleSize(0.05);
+  grmsym->Draw("PSAME");
+  grmsxm->Draw("PSAME");
+  grmsxmg->SetMarkerStyle(24);
+  grmsxmg->SetMarkerColor(1);
+  grmsymg->SetMarkerStyle(25);
+  grmsymg->SetMarkerColor(2);
+  grmszmg->SetMarkerStyle(26);
+  grmszmg->SetMarkerColor(4);
+  grmszmg->SetMinimum(0);
+  grmszmg->SetMaximum(rangeGrRms);
+  grmszmg->GetXaxis()->SetLimits(0,60);
+  grmszmg->Draw("PSAME");
+  grmszmg->GetXaxis()->SetTitle(trkstitle.Data());
+  grmszmg->GetXaxis()->SetTitleSize(0.05);
+  grmszmg->GetYaxis()->SetTitle("Resolution [#mum]");
+  grmszmg->GetYaxis()->SetTitleSize(0.05);
+  grmsymg->Draw("PSAME");
+  grmsxmg->Draw("PSAME");
+  leg->Draw();
+  sprintf(outgif,"vert%s-rms-mult.gif",vtxtype.Data());
+  if(optgif) cg2->SaveAs(outgif);
+
+
+
+
+  TCanvas *cg3=new TCanvas("cg3","Efficiency vs dNch/dy");
+  cg3->SetBottomMargin(0.14);
+  cg3->SetTopMargin(0.08);
+  cg3->SetLeftMargin(0.14);
+  cg3->SetRightMargin(0.08);
+  geffm->SetMarkerStyle(22);
+  geffm->SetMarkerColor(1);
+  geffm->GetXaxis()->SetLimits(0.,40.);
+  geffm->SetMinimum(0.);
+  geffm->SetMaximum(1.2);
+  geffm->Draw("AP");
+  geffm->GetXaxis()->SetTitle("MC dN_{ch}/dy in |y|<1");
+  geffm->GetXaxis()->SetTitleSize(0.05);
+  geffm->GetYaxis()->SetTitle("efficiency");
+  geffm->GetYaxis()->SetTitleSize(0.05);
+  sprintf(outgif,"vert%s-eff-mult.gif",vtxtype.Data());
+  if(optgif) cg3->SaveAs(outgif);
+
+
+  TCanvas *cg3b=new TCanvas("cg3b","Efficiency vs tracks");
+  cg3b->SetBottomMargin(0.14);
+  cg3b->SetTopMargin(0.08);
+  cg3b->SetLeftMargin(0.14);
+  cg3b->SetRightMargin(0.08);
+  gefftrks->SetMarkerStyle(22);
+  gefftrks->SetMarkerColor(1);
+  gefftrks->GetXaxis()->SetLimits(0.,40.);
+  gefftrks->SetMinimum(0.);
+  gefftrks->SetMaximum(1.2);
+  gefftrks->Draw("AP");
+  gefftrks->GetXaxis()->SetTitle(trkstitle.Data());
+  gefftrks->GetXaxis()->SetTitleSize(0.05);
+  gefftrks->GetYaxis()->SetTitle("efficiency");
+  gefftrks->GetYaxis()->SetTitleSize(0.05);
+  if(vtxtype.Contains("SPD")) {
+    geff3Dtrks->SetMarkerStyle(26);
+    geff3Dtrks->SetMarkerColor(1);
+    geff3Dtrks->Draw("P");
+  }
+  sprintf(outgif,"vert%s-eff-mult.gif",vtxtype.Data());
+  if(optgif) cg3b->SaveAs(outgif);
+
+
+  TCanvas *cg4=new TCanvas("cg4","Pulls");
+  cg4->SetBottomMargin(0.14);
+  cg4->SetTopMargin(0.08);
+  cg4->SetLeftMargin(0.14);
+  cg4->SetRightMargin(0.08);
+  gpullxm->SetMarkerStyle(20);
+  gpullxm->SetMarkerColor(1);
+  gpullym->SetMarkerStyle(21);
+  gpullym->SetMarkerColor(2);
+  gpullzm->SetMarkerStyle(22);
+  gpullzm->SetMarkerColor(4);
+  gpullzm->GetXaxis()->SetLimits(0,60);
+  gpullzm->SetMinimum(0);
+  gpullzm->SetMaximum(rangeGrPull);
+  gpullzm->Draw("AP");
+  gpullzm->GetXaxis()->SetTitle(trkstitle.Data());
+  gpullzm->GetXaxis()->SetTitleSize(0.05);
+  gpullzm->GetYaxis()->SetTitle("PULL");
+  gpullzm->GetYaxis()->SetTitleSize(0.05);
+  gpullxm->Draw("PSAME");
+  gpullym->Draw("PSAME");
+  gpullxmg->SetMarkerStyle(24);
+  gpullxmg->SetMarkerColor(1);
+  gpullymg->SetMarkerStyle(25);
+  gpullymg->SetMarkerColor(2);
+  gpullzmg->SetMarkerStyle(26);
+  gpullzmg->SetMarkerColor(4);
+  gpullzmg->GetXaxis()->SetLimits(0,60);
+  gpullzmg->SetMinimum(0);
+  gpullzmg->SetMaximum(rangeGrPull);
+  gpullzmg->Draw("PSAME");
+  gpullzmg->GetXaxis()->SetTitle(trkstitle.Data());
+  gpullzmg->GetXaxis()->SetTitleSize(0.05);
+  gpullzmg->GetYaxis()->SetTitle("PULL");
+  gpullzmg->GetYaxis()->SetTitleSize(0.05);
+  gpullxmg->Draw("PSAME");
+  gpullymg->Draw("PSAME");
+  TLine *lin=new TLine(0,1,60,1);
+  lin->SetLineStyle(2);
+  lin->Draw();
+  leg->Draw();
+  sprintf(outgif,"vert%s-pull-mult.gif",vtxtype.Data());
+  if(optgif) cg4->SaveAs(outgif);
+
+
+
+
+  TCanvas *cz1=new TCanvas("cz1","Efficiency vs. Z");
+  cz1->SetBottomMargin(0.14);
+  cz1->SetTopMargin(0.08);
+  cz1->SetLeftMargin(0.14);
+  cz1->SetRightMargin(0.08);
+  geffz->SetMarkerStyle(22);
+  geffz->SetMarkerColor(1);
+  geffz->GetXaxis()->SetLimits(-20,20.);
+  geffz->SetMinimum(0.);
+  geffz->SetMaximum(1.);
+  geffz->Draw("AP");
+  geffz->GetXaxis()->SetTitle("Z [cm]");
+  geffz->GetXaxis()->SetTitleSize(0.05);
+  geffz->GetYaxis()->SetTitle("efficiency");
+  geffz->GetYaxis()->SetTitleSize(0.05);
+  sprintf(outgif,"vert%s-eff-z.gif",vtxtype.Data());
+  if(optgif) cz1->SaveAs(outgif);
+
+
+  TLine *lin0z=new TLine(-20,0,20,0);
+  lin0z->SetLineStyle(2);
+  TCanvas *cz2=new TCanvas("cz2","Mean vs. Z");
+  cz2->SetBottomMargin(0.14);
+  cz2->SetTopMargin(0.08);
+  cz2->SetLeftMargin(0.14);
+  cz2->SetRightMargin(0.08);
+  gavexz->SetMarkerStyle(20);
+  gavexz->SetMarkerColor(1);
+  gaveyz->SetMarkerStyle(21);
+  gaveyz->SetMarkerColor(2);
+  gavezz->SetMarkerStyle(22);
+  gavezz->SetMarkerColor(4);
+  gavezz->GetXaxis()->SetLimits(-20,20.);
+  gavezz->SetMinimum(-rangeGrAve);
+  gavezz->SetMaximum(rangeGrAve);
+  gavezz->Draw("AP");
+  gavezz->GetXaxis()->SetTitle("Z [cm]");
+  gavezz->GetXaxis()->SetTitleSize(0.05);
+  gavezz->GetYaxis()->SetTitle("<Pos_{found}-Pos_{true}> [#mum]");
+  gavezz->GetYaxis()->SetTitleSize(0.05);
+  gavexz->Draw("PSAME");
+  gaveyz->Draw("PSAME");
+  lin0z->Draw();
+  gavexzg->SetMarkerStyle(24);
+  gavexzg->SetMarkerColor(1);
+  gavexzg->Draw("P");
+  gaveyzg->SetMarkerStyle(25);
+  gaveyzg->SetMarkerColor(2);
+  gaveyzg->Draw("P");
+  gavezzg->SetMarkerStyle(26);
+  gavezzg->SetMarkerColor(4);
+  gavezzg->Draw("P");
+  leg->Draw();
+  sprintf(outgif,"vert%s-ave-z.gif",vtxtype.Data());
+  if(optgif) cz2->SaveAs(outgif);
+
+
+  TCanvas *cz3=new TCanvas("cz3","Resolution vs. Z");
+  cz3->SetBottomMargin(0.14);
+  cz3->SetTopMargin(0.08);
+  cz3->SetLeftMargin(0.14);
+  cz3->SetRightMargin(0.08);
+  grmsxz->SetMarkerStyle(20);
+  grmsxz->SetMarkerColor(1);
+  grmsyz->SetMarkerStyle(21);
+  grmsyz->SetMarkerColor(2);
+  grmszz->SetMarkerStyle(22);
+  grmszz->SetMarkerColor(4);
+  grmszz->SetMinimum(0);
+  grmszz->SetMaximum(rangeGrRms);
+  grmszz->GetXaxis()->SetLimits(-20,20);
+  grmszz->Draw("AP");
+  grmszz->GetXaxis()->SetTitle("Z [cm]");
+  grmszz->GetXaxis()->SetTitleSize(0.05);
+  grmszz->GetYaxis()->SetTitle("Resolution [#mum]");
+  grmszz->GetYaxis()->SetTitleSize(0.05);
+  grmsxz->Draw("PSAME");
+  grmsyz->Draw("PSAME");
+  leg->Draw();
+  sprintf(outgif,"vert%s-rms-z.gif",vtxtype.Data());
+  if(optgif) cz3->SaveAs(outgif);
+
+  gStyle->SetPalette(1);
+  TCanvas *cbeam = new TCanvas("cbeam","Beam Long",800,800);
+  cbeam->Divide(2,2);
+  cbeam->cd(1);
+  hbeamx->GetYaxis()->SetTitle("X [cm]");
+  hbeamx->Draw();
+  cbeam->cd(3);
+  hbeamx->GetYaxis()->SetTitle("Y [cm]");
+  hbeamy->Draw();
+  cbeam->cd(2);
+  cbeam_2->SetLogz();
+  //gbeamxz->SetMarkerStyle(7);
+  hbeamxz->Draw("colz");
+  hbeamxz->GetXaxis()->SetTitle("Z [cm]");
+  hbeamxz->GetYaxis()->SetTitle("X [cm]");
+  cbeam_1->Update();
+  TPaveStats *st1=(TPaveStats*)hbeamxz->GetListOfFunctions()->FindObject("stats");
+  st1->SetX1NDC(0.13);
+  st1->SetX2NDC(0.33);
+  cbeam_2->Modified();
+  cbeam_2->Update();
+  cbeam->cd(4);
+  cbeam_4->SetLogz();
+  //gbeamyz->SetMarkerStyle(7);
+  hbeamyz->Draw("colz");
+  hbeamyz->GetXaxis()->SetTitle("Z [cm]");
+  hbeamyz->GetYaxis()->SetTitle("Y [cm]");
+  cbeam_4->Update();
+  TPaveStats *st2=(TPaveStats*)hbeamyz->GetListOfFunctions()->FindObject("stats");
+  st2->SetX1NDC(0.13);
+  st2->SetX2NDC(0.33);
+  cbeam_4->Modified();
+  cbeam_4->Update();
+  cbeam->Update();
+
+  TCanvas *cbeam2 = new TCanvas("cbeam2","Beam Transv",500,500);
+  cbeam2->SetLogz();
+  cbeam2->SetRightMargin(0.14);
+  //gbeamxy->SetMarkerStyle(7);
+  hbeamxy->Draw("colz");
+  hbeamxy->GetXaxis()->SetTitle("X [cm]");
+  hbeamxy->GetYaxis()->SetTitle("Y [cm]");
+  cbeam2->Update();
+  TPaveStats *st3=(TPaveStats*)hbeamxy->GetListOfFunctions()->FindObject("stats");
+  st3->SetX1NDC(0.13);
+  st3->SetX2NDC(0.33);
+  cbeam2->Modified();
+  cbeam2->Update();
+
+  return;
+}
+//----------------------------------------------------------------------------
+void ComputeVtxMean(TString vtxtype="TPC",
+                   TString fname="VertexESDwithMC.root",
+                   TString ntname="fNtupleVertexESD",
+                   Int_t nEventsToUse=10000,
+                   Int_t mincontr=10) {
+  //-----------------------------------------------------------------------
+  // Compute weighted mean and cov. matrix from the ntuple
+  //-----------------------------------------------------------------------
+  gStyle->SetOptStat(0);
+  //gStyle->SetOptFit(0);
+
+  Double_t diamondx=0.0200.,diamondy=0.0200.,diamondz=7.5.;
+
+  Double_t avx=0.;
+  Double_t avy=0.;
+  Double_t avz=0.;
+  Double_t wgtavx=0.;
+  Double_t wgtavy=0.;
+  Double_t wgtavz=0.;
+  Double_t sum=0.;
+  Double_t sumwgtx=0.;
+  Double_t sumwgty=0.;
+  Double_t sumwgtz=0.;
+  Double_t rmsx=0;
+  Double_t rmsy=0;
+  Double_t rmsz=0;
+  Double_t varx=0.;
+  Double_t vary=0.;
+  Double_t varz=0.;
+  Double_t covxy=0.;
+  Double_t covxz=0.;
+  Double_t covyz=0.;
+  Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
+
+  TH1F* hx = new TH1F("hx","",200,-1,1);
+  hx->SetXTitle("vertex x [#mu m]");
+  hx->SetYTitle("events");
+  TH1F* hy = new TH1F("hy","",200,-1,1);
+  hy->SetXTitle("vertex y [#mu m]");
+  hy->SetYTitle("events");
+  TH1F* hz = new TH1F("hz","",200,-20,20);
+  hz->SetXTitle("vertex z [cm]");
+  hz->SetYTitle("events");
+
+
+  TFile *f=new TFile(fname.Data());
+  TList *cOutput = (TList*)f->Get("cOutput");
+  TNtuple *nt=(TNtuple*)cOutput->FindObject(ntname.Data());
+  Int_t NNNev=nt->GetEntries();
+  printf("Events = %d\n",NNNev);
+  Float_t xVtx,xdiffVtx,xerrVtx;
+  Float_t yVtx,ydiffVtx,yerrVtx;
+  Float_t zVtx,zdiffVtx,zerrVtx;
+  Float_t ntrklets,ncontrVtx,dndy,triggered,vtx3D;
+  Float_t ztrue,zref;
+  
+  TString sxx="x"; sxx.Append(vtxtype.Data());
+  nt->SetBranchAddress(sxx.Data(),&xVtx);
+  TString syy="y"; syy.Append(vtxtype.Data());
+  nt->SetBranchAddress(syy.Data(),&yVtx);
+  TString szz="z"; szz.Append(vtxtype.Data());
+  nt->SetBranchAddress(szz.Data(),&zVtx);
+  
+  TString xdiff="xdiff"; xdiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(xdiff.Data(),&xdiffVtx);
+  TString ydiff="ydiff"; ydiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(ydiff.Data(),&ydiffVtx);
+  TString zdiff="zdiff"; zdiff.Append(vtxtype.Data());
+  nt->SetBranchAddress(zdiff.Data(),&zdiffVtx);
+    
+  TString xerr="xerr"; xerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(xerr.Data(),&xerrVtx);
+  TString yerr="yerr"; yerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(yerr.Data(),&yerrVtx);
+  TString zerr="zerr"; zerr.Append(vtxtype.Data());
+  nt->SetBranchAddress(zerr.Data(),&zerrVtx);
+
+  TString trkstitle;
+  if(vtxtype.Contains("TPC")) {
+    nt->SetBranchAddress("nTPCin",&ntrklets);
+    trkstitle="TPC tracks pointing to beam pipe";
+  } else {
+    nt->SetBranchAddress("ntrklets",&ntrklets);
+    trkstitle="SPD tracklets";
+  }
+  TString ntrks="ntrks"; ntrks.Append(vtxtype.Data());
+  nt->SetBranchAddress(ntrks.Data(),&ncontrVtx);
+  nt->SetBranchAddress("dndygen",&dndy);
+  
+  nt->SetBranchAddress("ztrue",&ztrue);
+
+  nt->SetBranchAddress("triggered",&triggered);
+
+  nt->SetBranchAddress("SPD3D",&vtx3D);
+
+  Int_t total=0;
+
+  Int_t divider=(Int_t)(NNNev/nEventsToUse);
+  // first loop on events
+  for(Int_t iev=0;iev<NNNev;iev++) {
+    if(iev%divider!=0) continue;
+    total++;
+    nt->GetEvent(iev);
+    if(!vtxtype.Contains("SPD")) vtx3D=1.;
+    if(vtx3D<0.5) continue;
+    if(triggered<0.5) continue; // not triggered
+    if(ncontrVtx<=0) continue; // no vertex
+
+    if(ncontrVtx<mincontr) continue;
+
+    avx += xVtx;
+    avy += yVtx;
+    avz += zVtx;
+    sum += 1.;
+    wgtavx += xVtx/xerrVtx/xerrVtx;
+    wgtavy += yVtx/yerrVtx/yerrVtx;
+    wgtavz += zVtx/zerrVtx/zerrVtx;
+    sumwgtx += 1./xerrVtx/xerrVtx;
+    sumwgty += 1./yerrVtx/yerrVtx;
+    sumwgtz += 1./zerrVtx/zerrVtx;
+     
+    hx->Fill(xVtx);
+    hy->Fill(yVtx);
+    hz->Fill(zVtx);
+  }
+  
+  avx /= sum;
+  avy /= sum;
+  avz /= sum;
+  wgtavx /= sumwgtx;
+  wgtavy /= sumwgty;
+  wgtavz /= sumwgtz;
+  ewgtavx = 1./TMath::Sqrt(sumwgtx);
+  ewgtavy = 1./TMath::Sqrt(sumwgty);
+  ewgtavz = 1./TMath::Sqrt(sumwgtz);
+  
+
+  // second loop on events
+  for(Int_t iev=0;iev<NNNev;iev++){
+    if(iev%divider!=0) continue;
+    nt->GetEvent(iev);
+    if(!vtxtype.Contains("SPD")) vtx3D=1.;
+    if(vtx3D<0.5) continue;
+    if(triggered<0.5) continue; // not triggered
+    if(ncontrVtx<=0) continue; // no vertex
+
+    if(ncontrVtx<mincontr) continue;
+  
+    varx += (xVtx-avx)*(xVtx-avx);
+    vary += (yVtx-avy)*(yVtx-avy);
+    varz += (zVtx-avz)*(zVtx-avz);
+    covxy += (xVtx-avx)*(yVtx-avy);
+    covxz += (xVtx-avx)*(zVtx-avz);
+    covyz += (yVtx-avy)*(zVtx-avz);
+  }
+  
+  varx /= sum;
+  vary /= sum;
+  varz /= sum;
+  covxy /= sum;
+  covxz /= sum;
+  covyz /= sum;
+  rmsx = TMath::Sqrt(varx);
+  rmsy = TMath::Sqrt(vary);
+  rmsz = TMath::Sqrt(varz);
+  eavx = rmsx/TMath::Sqrt(sum);
+  eavy = rmsy/TMath::Sqrt(sum);
+  eavz = rmsz/TMath::Sqrt(sum);
+  
+
+  printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
+  printf("Minimum number of contributors: %d\n",mincontr);
+  printf("Average:\n x = (%f +- %f) cm\n y = (%f +- %f) cm\n z = (%f +- %f) cm\n",avx,eavx,avy,eavy,avz,eavz);
+  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);
+  printf("RMS:\n x = %f cm\n y = %f cm\n z = %f cm\n",rmsx,rmsy,rmsz);
+
+  TCanvas *c = new TCanvas("c","c",0,0,1000,500);
+  c->Divide(3,1);
+  c->cd(1);
+  hx->Draw();
+  TF1 *gx = new TF1("gx","gaus",-1000,1000);
+  gx->SetLineColor(2);
+  hx->Fit(gx,"Q");
+  TF1 *gxx = (TF1*)gx->Clone("gxx");
+  gxx->FixParameter(2,diamondx);
+  gxx->SetLineStyle(2);
+  gxx->Draw("same");
+  c->cd(2);
+  hy->Draw();
+  TF1 *gy = new TF1("gy","gaus",-1000,1000);
+  gy->SetLineColor(2);
+  hy->Fit(gy,"Q");
+  TF1 *gyy = (TF1*)gy->Clone("gyy");
+  gyy->FixParameter(2,diamondy);
+  gyy->SetLineStyle(2);
+  gyy->Draw("same");
+  c->cd(3);
+  hz->Draw();
+  TF1 *gz = new TF1("gz","gaus",-10,10);
+  gz->SetLineColor(2);
+  hz->Fit(gz,"Q");
+  TF1 *gzz = (TF1*)gz->Clone("gzz");
+  gzz->FixParameter(2,diamondz);
+  gzz->SetLineStyle(2);
+  gzz->Draw("same");
+
+
+  return;
+}