Index: examples/ExampleClCount.C
===================================================================
--- examples/ExampleClCount.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
+++ examples/ExampleClCount.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
@@ -0,0 +1,173 @@
+#include <iostream>
+#include <TMath.h>
+#include <TVectorD.h>
+#include <TVector3.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TGraph.h>
+#include <TLine.h>
+#include <TRandom.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TLegend.h>
+#include "external/TrackCovariance/TrkUtil.h"
+
+
+void ExampleClCount(Int_t Nev = 1000, Double_t ang = 90., Int_t Opt=0)
+{
+	//
+	// Constants
+	Double_t angRad = ang*TMath::Pi() / 180.;	// To radians
+	Double_t KpMass = 0.493677;		// Charged Kaon mass
+	Double_t PiMass = 0.13957018;	// Charged Pion mass
+	Double_t Bz = 2.0;				// Field
+	TVector3 x(0.0, 0.0, 0.0);		// Track origin
+	//
+	Double_t pmin = 0.1;				// Momentum range
+	Double_t pmax = 50.;
+	//
+	// Setup chamber
+	//
+	TrkUtil *TU = new TrkUtil(Bz);
+	//cout << "TU done" << endl;
+	Double_t Rmin = 0.35;
+	Double_t Rmax = 2.0;
+	Double_t Zmin = -2.0;
+	Double_t Zmax = 2.0;
+	TU->SetDchBoundaries(Rmin, Rmax, Zmin, Zmax);
+	TU->SetGasMix(Opt);	// gas selection
+	char title[50];
+	Double_t y_min = 0.0;
+	Double_t y_max = 0.0;
+	switch (Opt)
+	{
+	case 0:
+	{
+		sprintf(title, "Helium 90 - Isobutane 10");	// He-Isobutane
+		y_min = 1500.; y_max = 3500;
+		break;
+	}
+	case 1:
+	{
+		sprintf(title, "Helium 100");				// pure He
+		y_min = 400.; y_max = 1000;
+		break;
+	}
+	case 2:
+	{
+		sprintf(title, "Argon 50 - Ethane 50");		// Argon - Ethane
+		y_min = 6000.; y_max = 8500.;
+		break;
+	}
+	case 3:
+	{
+		sprintf(title, "Argon 100");				// pure Argon
+		y_min = 4000.; y_max = 6000.;
+		break;
+	}
+	}
+	//cout << "DchBoundaries set" << endl;
+	//
+	// Histograms
+	//
+	TH2F *h_kaon = new TH2F("h_kaon", "Nclusters distribution", 200, 0., 50., 50, y_min, y_max);
+	TH2F *h_pion = new TH2F("h_pion", "Nclusters distribution", 200, 0., 50., 50, y_min, y_max);
+	//cout << "histograms booked" << endl;
+	//
+	// Fill plots
+	//
+	for (Int_t n = 0; n < Nev; n++)
+	{
+		Double_t R = gRandom->Rndm();
+		Double_t pval = pmin + R*(pmax - pmin);
+		Double_t px = pval*TMath::Sin(angRad);
+		Double_t pz = pval*TMath::Cos(angRad);
+		TVector3 p(px, 0.0, pz);
+		Double_t Q = 1.0;
+		TVectorD Par = TrkUtil::XPtoPar(x, p, Q, Bz);
+		//cout << "Parameters done. pval = "<< pval << endl;
+		//
+		Double_t Ncl = 0.0;
+		if (TU->IonClusters(Ncl, KpMass, Par))
+		{
+			//cout << "Kaon extracted. pval = "<< pval << endl;
+			h_kaon->Fill(pval, Ncl);
+			//cout << "Kaon filled. Ncl = "<<Ncl << endl;
+		}
+		if (TU->IonClusters(Ncl, PiMass, Par))
+		{
+			//cout << "Pion extracted. pval = "<<pval << endl;
+			h_pion->Fill(pval, Ncl);
+			//cout << "Pion filled. Ncl = "<<Ncl << endl;
+		}
+	}
+	//
+	// graph K/pi separation
+	//
+	const Int_t Nint = 501;
+	Double_t pmom[Nint];
+	Double_t pmn = 0.2; Double_t pmx = 100.;
+	Double_t stp = (TMath::Log(pmx) - TMath::Log(pmn)) / (Double_t)(Nint - 1);
+	for (Int_t i = 0; i < Nint; i++)pmom[i] = TMath::Exp(i * stp + TMath::Log(pmn));
+	Double_t SigDiff[Nint];
+	for (Int_t i = 0; i < Nint; i++)
+	{
+		Double_t bgK = pmom[i] / KpMass;
+		Double_t bgP = pmom[i] / PiMass;
+		Double_t px = pmom[i]*TMath::Sin(angRad);
+		Double_t pz = pmom[i]*TMath::Cos(angRad);
+		TVector3 p(px, 0.0, pz);
+		Double_t Q = 1.0;
+		TVectorD Par = TrkUtil::XPtoPar(x, p, Q, Bz);
+		Double_t tLen = TU->TrkLen(Par);
+		Double_t CluKp = TrkUtil::Nclusters(bgK,Opt)*tLen;
+		Double_t CluKpErr = TMath::Sqrt(CluKp);
+		Double_t CluPi = TrkUtil::Nclusters(bgP,Opt)*tLen;
+		Double_t CluPiErr = TMath::Sqrt(CluPi);
+		//cout << "Momentum = " << pmom[i] << ", length= " << tLen
+		//	<< ", CluKp = " << CluKp << ", CluPi = " << CluPi << endl;
+		//if (CluKp + CluPi > 0.0)SigDiff[i] = TMath::Abs(CluPi - CluKp) / TMath::Sqrt(CluKp + CluPi);
+		if (CluKp + CluPi > 0.0)SigDiff[i] = TMath::Abs(CluPi - CluKp) / (0.5*(CluKpErr + CluPiErr));
+	}
+	//
+	// Plots
+	//
+	TCanvas *cnv = new TCanvas("cnv", title, 50, 50, 800, 500);
+	cnv->Divide(2, 1);
+	cnv->cd(1);
+	gStyle->SetOptStat(0);
+	h_kaon->SetMarkerColor(kRed);
+	h_kaon->SetLineColor(kRed);
+	h_kaon->GetXaxis()->SetTitle("Momentum (GeV)");
+	h_kaon->Draw();
+	h_pion->SetMarkerColor(kBlack);
+	h_pion->SetLineColor(kBlack);
+	h_pion->Draw("SAME");
+	//
+	TLegend* lg = new TLegend(0.1, 0.9, 0.3, 0.70);
+	TString LgTitle = "Particle type:";
+	lg->SetHeader(LgTitle);
+	lg->AddEntry(h_pion, "#pi", "L");
+	lg->AddEntry(h_kaon, "#color[2]{K}", "L");
+	lg->Draw();
+	cnv->cd(2);
+	gPad->SetLogx();
+	TGraph *gr = new TGraph(Nint, pmom, SigDiff);
+	gr->SetMinimum(0.0);
+	gr->SetMaximum(8.0);
+	gr->GetXaxis()->SetLimits(0.,100.);
+	gr->SetTitle("K/#pi separation in nr. of #sigma");
+	gr->GetXaxis()->SetTitle("Momentum (GeV)");
+	gr->SetLineColor(kBlue);
+	gr->Draw("APC");
+	TLine *line = new TLine(0.0, 3.0, 100., 3.0);
+	line->Draw("SAME");
+	//
+	TCanvas* cnvf = new TCanvas("cnvf", "Interpolating function", 100, 100, 800, 500);
+	cnvf->Divide(1, 1);
+	cnvf->cd(1);
+	gPad->SetLogx();
+	TF1* f_ncl = new TF1("f_ncl", TU, &TrkUtil::funcNcl, 0.5, 1000., 1);
+	f_ncl->SetNpx(500);
+	f_ncl->Draw();
+}
Index: examples/ExamplePVtxFind.C
===================================================================
--- examples/ExamplePVtxFind.C	(revision a7bd0513d5985e067963efaabe4931d7311ed48c)
+++ examples/ExamplePVtxFind.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
@@ -51,7 +51,6 @@
 		// Load selected branches with data from specified event
 		treeReader->ReadEntry(entry);
-		Int_t Ntr = 0;	// # of starting tracks used in primary vertex
 		Int_t NtrG = branchTrack->GetEntries();
-		//std::cout << "Event opened containing " << NtrG << " tracks" << std::endl;
+		if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl;
 		TVectorD** pr = new TVectorD * [NtrG];
 		TMatrixDSym** cv = new TMatrixDSym * [NtrG];
@@ -63,4 +62,7 @@
 		//
 		Double_t Nprim = 0.0;
+		Double_t xpv = 0.0;		// Init true primary vertex position
+		Double_t ypv = 0.0;
+		Double_t zpv = 0.0;
 		if (branchTrack->GetEntries() > 0)
 		{
@@ -80,16 +82,12 @@
 				//std::cout << "Got track parameters for track " << it << std::endl;
 				//
-				// Load tracks for vertex fit if impact parameters is not ridiculous
-				Double_t Dmax = 1.0;	// max is 1 mm
-				if (TMath::Abs(obsD0) < Dmax) {
-					Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
-					TVectorD obsPar(5, oPar);	// Fill observed parameters
-					pr[Ntr] = new TVectorD(obsPar);
-					cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
-					Ntr++;
-					//std::cout << "Track loaded Ntr= " << Ntr << std::endl;
-				}
-				//
-				// Get associated generated particle
+				// Load all tracks for vertex fit 
+				Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
+				TVectorD obsPar(5, oPar);	// Fill observed parameters
+				pr[it] = new TVectorD(obsPar);
+				cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
+				//std::cout << "Track loaded it= " << it << std::endl;
+			//
+			// Find true primary vertex
 				GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
 				//std::cout << "GenParticle pointer "<<gp << std::endl;
@@ -99,25 +97,83 @@
 				Double_t y = gp->Y;
 				Double_t z = gp->Z;
-				//std::cout << "Got position of origin " << std::endl;
-				//
-				// Count tracks originating from the primary vertex
-				if (x == 0.0 && y == 0.0) Nprim++;
+				Bool_t prim = kTRUE;	// Is primary?
+				Int_t mp = gp->M1;	// Mother
+				while (mp > 0) {
+					GenParticle* gm =
+						(GenParticle*)branchGenPart->At(mp);
+					Double_t xm = gm->X;
+					Double_t ym = gm->Y;
+					Double_t zm = gm->Z;
+					if (x != xm || y != ym || z != zm) {
+						prim = kFALSE;
+						break;
+					}
+					else mp = gm->M1;
+				}
+				if (prim) {		// It's a primary track
+					Nprim++;
+					xpv = x;	// Store true primary
+					ypv = y;
+					zpv = z;
+				}
 
 			}		// End loop on tracks
 		}
 		//
-		// Fit primary vertex
-		//
-		Int_t Nfound = Ntr;
-		//std::cout << "Found tracks "<<Nfound << std::endl;
-		Int_t MinTrk = 2;	// Minumum # tracks for vertex fit
-		Double_t MaxChi2 = 6.;
-		if (Ntr >= MinTrk) {
-			VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
+		// Find primary vertex
+		//
+		//Beam constraint
+		TVectorD xpvc(3);
+		xpvc(0) = 1.0;
+		xpvc(1) = -2.0;
+		xpvc(2) = 10.0;
+		TMatrixDSym covpvc(3); covpvc.Zero();
+		covpvc(0, 0) = 0.0097 * 0.0097;
+		covpvc(1, 1) = 2.55e-05 * 2.55e-05;
+		covpvc(2, 2) = 0.64 * 0.64;
+		//
+		//
+		// Skim tracks
+		Int_t nSkim = 0;
+		Int_t* nSkimmed = new Int_t[NtrG];
+		VertexFit* Vskim = new VertexFit();
+		Vskim->AddVtxConstraint(xpvc, covpvc);
+		Double_t MaxChi2 = 2.0+3*sqrt(2.);
+		for (Int_t n = 0; n < NtrG; n++) {
+			Vskim->AddTrk(pr[n], cv[n]);
+			Double_t Chi2One = Vskim->GetVtxChi2();
+			//std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
+			if (Chi2One < MaxChi2) {
+				nSkimmed[nSkim] = n;
+				//std::cout << "nSkimmed[" << nSkim << "] = " << n << std::endl;
+				nSkim++;
+			}
+			Vskim->RemoveTrk(0);
+		}
+		// Load tracks for primary fit
+		Int_t MinTrk = 1;	// Minumum # tracks for vertex fit
+		TVectorD** PrFit = new TVectorD * [nSkim];
+		TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
+		if (nSkim >= MinTrk) {
+			for (Int_t n = 0; n < nSkim; n++) {
+				PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
+				CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
+				TVectorD test = *PrFit[n];
+				//std::cout << "Test *PrFit[" << n << "](0) = " << test(0) << std::endl;
+			}
+		}
+		delete[] nSkimmed;
+		delete Vskim;
+		Int_t Nfound = nSkim;
+		if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
+		if (nSkim >= MinTrk) {
+			VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
 			//std::cout << "Vertex fit created " << std::endl;
+			Vtx->AddVtxConstraint(xpvc, covpvc);
 			//
 			// Remove tracks with large chi2
+			Double_t MaxChi2Fit = 9.;
 			Bool_t Done = kFALSE;
-			while (!Done){
+			while (!Done) {
 				//std::cout << "After while " << std::endl;
 				// Find largest Chi2 contribution
@@ -130,6 +186,8 @@
 				Double_t Chi2Mx = Chi2L[iMax];
 				//std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
-				if (Chi2Mx > MaxChi2 && Nfound > 2) {
-					Vtx->RemoveTrk(iMax); 
+				if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
+					//std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
+					Vtx->RemoveTrk(iMax);
+					//std::cout << "After remove." << std::endl;
 					Nfound--;
 				}
@@ -145,5 +203,5 @@
 			//
 			// Require minimum number of tracks in vertex
-			Int_t Nmin = 4;
+			Int_t Nmin = 2;
 			if (Nfound >= Nmin) {
 				TVectorD xvtx = Vtx->GetVtx();
@@ -151,8 +209,8 @@
 				TMatrixDSym covX = Vtx->GetVtxCov();
 				Double_t Chi2 = Vtx->GetVtxChi2();
-				Double_t Ndof = 2 * (Double_t)Nfound - 3;
-				Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
-				Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
-				Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
+				Double_t Ndof = 2 * (Double_t)Nfound;
+				Double_t PullX = (xvtx(0) - xpv) / TMath::Sqrt(covX(0, 0));
+				Double_t PullY = (xvtx(1) - ypv) / TMath::Sqrt(covX(1, 1));
+				Double_t PullZ = (xvtx(2) - zpv) / TMath::Sqrt(covX(2, 2));
 				//
 				// Fill histograms
@@ -173,8 +231,12 @@
 		//
 		// Cleanup
-		for (Int_t i = 0; i < Ntr; i++) delete pr[i];
-		for (Int_t i = 0; i < Ntr; i++) delete cv[i];
+		for (Int_t i = 0; i < NtrG; i++) delete pr[i];
+		for (Int_t i = 0; i < NtrG; i++) delete cv[i];
+		for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
+		for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
 		delete[] pr;
 		delete[] cv;
+		delete[] PrFit;
+		delete[] CvFit;
 	}
 	//
@@ -195,5 +257,5 @@
 	CnvN->cd(1);
 	hTrPrim->Draw();
-	CnvN-> cd(2);
+	CnvN->cd(2);
 	hTrFound->SetLineColor(kRed);
 	hTrFound->Draw();
Index: examples/ExamplePVtxFit.C
===================================================================
--- examples/ExamplePVtxFit.C	(revision a7bd0513d5985e067963efaabe4931d7311ed48c)
+++ examples/ExamplePVtxFit.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
@@ -14,4 +14,5 @@
 
 #endif
+
 
 
@@ -45,8 +46,12 @@
 		// Load selected branches with data from specified event
 		treeReader->ReadEntry(entry);
+		// 
 		Int_t Ntr = 0;	// # of tracks from primary vertex
 		Int_t NtrG = branchTrack->GetEntries();
 		TVectorD** pr = new TVectorD * [NtrG];
 		TMatrixDSym** cv = new TMatrixDSym * [NtrG];
+		//
+		// True vertex
+		Double_t xpv, ypv, zpv;
 		// If event contains at least 1 track
 		//
@@ -60,4 +65,10 @@
 				// Get associated generated particle
 				GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
+				TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
+				TVector3 pg(gp->Px,gp->Py,gp->Pz);
+				Double_t Q = (Double_t)gp->Charge;
+				Double_t Bz = 2.0;
+				TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
+				TVectorD genPar = TrkUtil::ParToMm(genParM);
 				//
 				// Position of origin in mm
@@ -65,8 +76,26 @@
 				Double_t y = gp->Y;
 				Double_t z = gp->Z;
+				Bool_t prim = kTRUE;	// Is primary?
+				Int_t mp = gp->M1;	// Mother
+				while(mp>0){
+					GenParticle* gm = 
+					(GenParticle*)branchGenPart->At(mp);
+					Double_t xm = gm->X;
+                                	Double_t ym = gm->Y;
+                                	Double_t zm = gm->Z;
+					if(x!=xm || y!=ym || z!=zm){
+						prim = kFALSE;
+						break;
+					}else mp = gm->M1;
+				}
 				//
 				// group tracks originating from the primary vertex
-				if (x == 0.0 && y == 0.0)
+				if (prim)
 				{
+					//std::cout<<"Event: "<<entry<<", track: "<<it;
+					//std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 
+					xpv = x;
+					ypv = y;
+					zpv = z; 
 					//
 					// Reconstructed track parameters
@@ -75,4 +104,7 @@
 					Double_t obsC = trk->C;
 					Double_t obsZ0 = trk->DZ;
+					//std::cout<<"Z0 track = "<< obsZ0
+					//<<", gen Z0 = "<<genPar(3)
+					//<<", gen cot = "<<genPar(4)<<std::endl;
 					Double_t obsCtg = trk->CtgTheta;
 					Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
@@ -96,7 +128,11 @@
 			Double_t Chi2 = Vtx->GetVtxChi2();
 			Double_t Ndof = 2 * (Double_t)Ntr - 3;
-			Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
-			Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
-			Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
+			Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
+			Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
+			Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
+			//std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
+			//<<ypv<<", "<<zpv<<std::endl;
+			//std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
+                        //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
 			//
 			// Fill histograms
Index: examples/ExamplePVtxFitEC.C
===================================================================
--- examples/ExamplePVtxFitEC.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
+++ examples/ExamplePVtxFitEC.C	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
@@ -0,0 +1,176 @@
+/*
+Example of using vertex fitter class to fit primary vertex
+assumed to be generated in (0,0,0)
+*/
+
+#ifdef __CLING__
+R__LOAD_LIBRARY(libDelphes)
+
+#include "classes/DelphesClasses.h"
+#include "external/ExRootAnalysis/ExRootTreeReader.h"
+#include "modules/TrackCovariance.h"
+#include "external/TrackCovariance/TrkUtil.h"
+#include "external/TrackCovariance/VertexFit.h"
+
+#endif
+
+
+
+//------------------------------------------------------------------------------
+
+void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
+{
+	// Create chain of root trees
+	TChain chain("Delphes");
+	chain.Add(inputFile);
+
+	// Create object of class ExRootTreeReader
+	ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
+	Long64_t numberOfEntries = treeReader->GetEntries();
+
+	// Get pointers to branches used in this analysis
+	TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
+	TClonesArray* branchTrack = treeReader->UseBranch("Track");
+
+	// Book histograms
+	Int_t Nbin = 100;
+	TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
+	TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
+	TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
+	TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
+	//
+	// Loop over all events
+	Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
+	for (Int_t entry = 0; entry < Nev; ++entry)
+	{
+		// Load selected branches with data from specified event
+		treeReader->ReadEntry(entry);
+		// 
+		Int_t Ntr = 0;	// # of tracks from primary vertex
+		Int_t NtrG = branchTrack->GetEntries();
+		TVectorD** pr = new TVectorD * [NtrG];
+		TMatrixDSym** cv = new TMatrixDSym * [NtrG];
+		//
+		// True vertex
+		Double_t xpv, ypv, zpv;
+		// If event contains at least 1 track
+		//
+		if (branchTrack->GetEntries() > 0)
+		{
+			// Loop on tracks
+			for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
+			{
+				Track* trk = (Track*)branchTrack->At(it);
+				//
+				// Get associated generated particle
+				GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
+				TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
+				TVector3 pg(gp->Px,gp->Py,gp->Pz);
+				Double_t Q = (Double_t)gp->Charge;
+				Double_t Bz = 2.0;
+				TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
+				TVectorD genPar = TrkUtil::ParToMm(genParM);
+				//
+				// Position of origin in mm
+				Double_t x = gp->X;
+				Double_t y = gp->Y;
+				Double_t z = gp->Z;
+				Bool_t prim = kTRUE;	// Is primary?
+				Int_t mp = gp->M1;	// Mother
+				while(mp>0){
+					GenParticle* gm = 
+					(GenParticle*)branchGenPart->At(mp);
+					Double_t xm = gm->X;
+                                	Double_t ym = gm->Y;
+                                	Double_t zm = gm->Z;
+					if(x!=xm || y!=ym || z!=zm){
+						prim = kFALSE;
+						break;
+					}else mp = gm->M1;
+				}
+				//
+				// group tracks originating from the primary vertex
+				if (prim)
+				{
+					//std::cout<<"Event: "<<entry<<", track: "<<it;
+					//std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl; 
+					xpv = x;
+					ypv = y;
+					zpv = z; 
+					//
+					// Reconstructed track parameters
+					Double_t obsD0 = trk->D0;
+					Double_t obsPhi = trk->Phi;
+					Double_t obsC = trk->C;
+					Double_t obsZ0 = trk->DZ;
+					//std::cout<<"Z0 track = "<< obsZ0
+					//<<", gen Z0 = "<<genPar(3)
+					//<<", gen cot = "<<genPar(4)<<std::endl;
+					Double_t obsCtg = trk->CtgTheta;
+					Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
+					TVectorD obsPar(5, oPar);	// Fill observed parameters
+					//
+					pr[Ntr] = new TVectorD(obsPar);
+					//std::cout<<"Cov Matrix:"<<std::endl;
+					//trk->CovarianceMatrix().Print();
+					cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
+					Ntr++;
+				}
+			}		// End loop on tracks
+		}
+		//
+		// Fit primary vertex with beam constraint
+		//
+		//Beam constraint
+		TVectorD xpvc(3);
+		xpvc(0) = 1.0;
+		xpvc(1) = -2.0;
+		xpvc(2) = 10.0;
+		TMatrixDSym covpvc(3); covpvc.Zero();
+		covpvc(0, 0) = 0.0097 * 0.0097;
+		covpvc(1, 1) = 2.55e-05 * 2.55e-05;
+		covpvc(2, 2) = 0.64 * 0.64;
+		Int_t MinTrk = 2;	// Minumum # tracks for vertex fit
+		if (Ntr >= MinTrk) {
+			VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
+			Vtx->AddVtxConstraint(xpvc, covpvc);
+			TVectorD xvtx = Vtx->GetVtx();
+			TMatrixDSym covX = Vtx->GetVtxCov();
+			Double_t Chi2 = Vtx->GetVtxChi2();
+			Double_t Ndof = 2 * (Double_t)Ntr;
+			Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
+			Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
+			Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
+			//std::cout<<"**** True  vertex (x, y, z) "<<xpv<<", "
+			//<<ypv<<", "<<zpv<<std::endl;
+			//std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
+                        //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
+			//
+			// Fill histograms
+			hXpull->Fill(PullX);
+			hYpull->Fill(PullY);
+			hZpull->Fill(PullZ);
+			hChi2->Fill(Chi2 / Ndof);
+		}
+
+		//std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
+		//
+		// Cleanup
+		for (Int_t i = 0; i < Ntr; i++) delete pr[i];
+		for (Int_t i = 0; i < Ntr; i++) delete cv[i];
+		delete[] pr;
+		delete[] cv;
+	}
+	//
+	// Show resulting histograms
+	//
+	TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
+	Cnv->Divide(2, 2);
+	Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hXpull->Fit("gaus"); hXpull->Draw();
+	Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hYpull->Fit("gaus"); hYpull->Draw();
+	Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hZpull->Fit("gaus"); hZpull->Draw();
+	Cnv->cd(4); hChi2->Draw();
+}
Index: examples/Pythia8/ee_zh.cmnd
===================================================================
--- examples/Pythia8/ee_zh.cmnd	(revision a7bd0513d5985e067963efaabe4931d7311ed48c)
+++ examples/Pythia8/ee_zh.cmnd	(revision 4acf2fd03ea62013a5ce2d50c6c501206c4c8906)
@@ -1,3 +1,4 @@
 
+! number of events to generate
 Main:numberOfEvents = 1000         ! number of events to generate
 
@@ -7,11 +8,10 @@
 
 ! Vertex smearing :
-!Beams:allowVertexSpread = on
-!Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
-!Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
-!Beams:sigmaVertexZ = 0.64      !  0.64 mm
-!Beams:sigmaTime    = 0.64      !  0.64 mm
-ï¿Œ
 
+Beams:allowVertexSpread = on
+Beams:sigmaVertexX = 9.70e-3   !  13.7 mum / sqrt2
+Beams:sigmaVertexY = 25.5E-6   !  36.1 nm / sqrt2
+Beams:sigmaVertexZ = 0.64      !  0.64 mm
+Beams:sigmaTime    = 0.64      !  0.64 mm
 
 ! Higgsstrahlung process
@@ -20,3 +20,3 @@
 ! 5) Force the Z decays to muons
 23:onMode = off
-23:onIfAny = 11 -11
+23:onIfAny = 13 -13
