Index: examples/ExamplePVtxFind.C
===================================================================
--- examples/ExamplePVtxFind.C	(revision ebf40fd3f3ff25430347fd821d73d468acce4df6)
+++ examples/ExamplePVtxFind.C	(revision 46b3e017290a3d0e25f3c45d70c9fc3d53b01af8)
@@ -1,5 +1,5 @@
 /*
 Example of using vertex fitter class to fit primary vertex
-assumed to be generated in (0,0,0)
+assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
 */
 
@@ -44,4 +44,5 @@
 	TH1D* hTrFoundG = new TH1D("hTrFoundG", "Found GOOD primary tracks", 41, -0.5, 40.5);
 	TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
+	TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
 	//
 	// Loop over all events
@@ -87,4 +88,19 @@
 				pr[it] = new TVectorD(obsPar);
 				cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
+				// Check covariance
+				/*
+				TMatrixDSymEigen Eign(*cv[it]);
+				TVectorD lambda = Eign.GetEigenValues();
+				Bool_t good = kTRUE;
+				for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;
+				if(!good){
+					std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "
+					<<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "
+					<<lambda(3)<<", "<<lambda(4)<<std::endl;
+					TMatrixDSym Ctmp = *cv[it];
+					std::cout<<"Error D (dir - cov) "<<trk->ErrorD0
+					<<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;
+				}
+				*/
 				//std::cout << "Track loaded it= " << it << std::endl;
 			//
@@ -120,4 +136,5 @@
 			}		// End loop on tracks
 		}
+		std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
 		//
 		// Find primary vertex
@@ -132,4 +149,9 @@
 		covpvc(1, 1) = 2.55e-05 * 2.55e-05;
 		covpvc(2, 2) = 0.64 * 0.64;
+		if(Nprim == 0){
+			xpv = xpvc(0);
+			ypv = xpvc(1);
+			zpv = xpvc(2);
+		}
 		//
 		//
@@ -137,9 +159,12 @@
 		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.);
+		TVectorD** PrSk = new TVectorD * [1];
+		TMatrixDSym** CvSk = new TMatrixDSym * [1];
+		Double_t MaxChi2 = 2.0+3*2.;
 		for (Int_t n = 0; n < NtrG; n++) {
-			Vskim->AddTrk(pr[n], cv[n]);
+			PrSk[0] = new TVectorD(*pr[n]);
+			CvSk[0] = new TMatrixDSym(*cv[n]);
+			VertexFit* Vskim = new VertexFit(1,PrSk, CvSk);
+			Vskim->AddVtxConstraint(xpvc, covpvc);
 			Double_t Chi2One = Vskim->GetVtxChi2();
 			//std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
@@ -149,23 +174,25 @@
 				nSkim++;
 			}
-			Vskim->RemoveTrk(0);
-		}
+			// Cleanup
+			delete Vskim;
+		}
+		delete PrSk[0];
+		delete CvSk[0];
+		delete[] PrSk;
+		delete[] CvSk;
+		//
 		// 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) {
+			TVectorD** PrFit = new TVectorD * [nSkim];
+			TMatrixDSym** CvFit = new TMatrixDSym * [nSkim];
 			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) {
+			}
+			delete[] nSkimmed;
+			Int_t Nfound = nSkim;
+			if(entry%500 ==0)std::cout << "Found tracks "<<Nfound << std::endl;
+			const Int_t MaxFound = 100; Double_t Chi2LL[MaxFound]; Double_t *Chi2L = Chi2LL;
 			VertexFit* Vtx = new VertexFit(nSkim, PrFit, CvFit);
 			//std::cout << "Vertex fit created " << std::endl;
@@ -173,5 +200,5 @@
 			//
 			// Remove tracks with large chi2
-			Double_t MaxChi2Fit = 9.;
+			Double_t MaxChi2Fit = 4.5;
 			Bool_t Done = kFALSE;
 			while (!Done) {
@@ -180,5 +207,5 @@
 				TVectorD Chi2List = Vtx->GetVtxChi2List();	// Get contributions to Chi2
 				//std::cout << "After Chi2List.  " << std::endl; Chi2List.Print();
-				Double_t* Chi2L = new Double_t[Nfound];
+				//Double_t* Chi2L = new Double_t[Nfound];
 				Chi2L = Chi2List.GetMatrixArray();
 				Int_t iMax = TMath::LocMax(Nfound, Chi2L);
@@ -186,5 +213,5 @@
 				Double_t Chi2Mx = Chi2L[iMax];
 				//std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
-				if (Chi2Mx > MaxChi2Fit && Nfound > 2) {
+				if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
 					//std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
 					Vtx->RemoveTrk(iMax);
@@ -195,7 +222,4 @@
 					Done = kTRUE;
 				}
-				//std::cout << "Done =  " << Done << std::endl;
-				//delete[] Chi2L;
-				//std::cout << "Array Chi2L removed " << std::endl;
 			}
 			//
@@ -203,8 +227,8 @@
 			//
 			// Require minimum number of tracks in vertex
-			Int_t Nmin = 2;
+			Int_t Nmin = 1;
 			if (Nfound >= Nmin) {
 				TVectorD xvtx = Vtx->GetVtx();
-				//std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
+				std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
 				TMatrixDSym covX = Vtx->GetVtxCov();
 				Double_t Chi2 = Vtx->GetVtxChi2();
@@ -222,8 +246,14 @@
 				hTrPrim->Fill(Nprim);
 				hTrFound->Fill((Double_t)Nfound);
+				hTrDiff->Fill((Double_t)Nfound-Nprim);
 				//std::cout << "Histograms filled " << std::endl;
 			}
 			//
+			// Clean
 			delete Vtx;
+			for (Int_t i = 0; i < nSkim; i++) delete PrFit[i];
+			for (Int_t i = 0; i < nSkim; i++) delete CvFit[i];
+			delete[] PrFit;
+			delete[] CvFit;
 		}
 
@@ -233,10 +263,6 @@
 		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;
 	}
 	//
@@ -254,5 +280,5 @@
 	//
 	TCanvas* CnvN = new TCanvas("CnvN", "Primary tracks found", 100, 100, 900, 500);
-	CnvN->Divide(2, 1);
+	CnvN->Divide(2, 2);
 	CnvN->cd(1);
 	hTrPrim->Draw();
@@ -260,3 +286,5 @@
 	hTrFound->SetLineColor(kRed);
 	hTrFound->Draw();
+	CnvN->cd(3);
+	hTrDiff->Draw();
 }
Index: examples/ExamplePVtxFitEC.C
===================================================================
--- examples/ExamplePVtxFitEC.C	(revision ebf40fd3f3ff25430347fd821d73d468acce4df6)
+++ examples/ExamplePVtxFitEC.C	(revision 46b3e017290a3d0e25f3c45d70c9fc3d53b01af8)
@@ -1,5 +1,5 @@
 /*
 Example of using vertex fitter class to fit primary vertex
-assumed to be generated in (0,0,0)
+assumed to be generated with Pythia8/ee_zh_smr-shf.cmn
 */
 
@@ -21,4 +21,20 @@
 void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
 {
+	//
+	// Beam constraint 
+	// (consistent with generation with ee_zh_smr-shf.cmnd)
+	//
+	// Mean beam position
+	TVectorD xpvc(3);
+	xpvc(0) = 1.0;
+	xpvc(1) = -2.0;
+	xpvc(2) = 10.0;
+	// Interaction region covariance
+	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;
+
+	//
 	// Create chain of root trees
 	TChain chain("Delphes");
@@ -35,8 +51,14 @@
 	// Book histograms
 	Int_t Nbin = 100;
+	// Vertex fit pulls
 	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.);
+	// Generation check
+	TH1D* hXvpull = new TH1D("hXvpull", "Pull X generated vertex", Nbin, -10., 10.);
+	TH1D* hYvpull = new TH1D("hYvpull", "Pull Y generated vertex", Nbin, -10., 10.);
+	TH1D* hZvpull = new TH1D("hZvpull", "Pull Z generated vertex", Nbin, -10., 10.);
+
 	//
 	// Loop over all events
@@ -65,10 +87,13 @@
 				// 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 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z); 	// mm -> meters
 				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);
+				TVectorD genPar = TrkUtil::ParToMm(genParM); 		// -> back to mm
+
+				// 
+				// Check if original track is from primary vertex
 				//
 				// Position of origin in mm
@@ -79,20 +104,20 @@
 				Int_t mp = gp->M1;	// Mother
 				while(mp>0){
-					GenParticle* gm = 
-					(GenParticle*)branchGenPart->At(mp);
+					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;
+						prim = kFALSE;	// Not primary
 						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; 
+					//
+					// Store true primary vertex for this event
 					xpv = x;
 					ypv = y;
@@ -104,7 +129,4 @@
 					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 };
@@ -112,6 +134,4 @@
 					//
 					pr[Ntr] = new TVectorD(obsPar);
-					//std::cout<<"Cov Matrix:"<<std::endl;
-					//trk->CovarianceMatrix().Print();
 					cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
 					Ntr++;
@@ -122,14 +142,5 @@
 		// 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
+		Int_t MinTrk = 3;	// Minumum # tracks for vertex fit
 		if (Ntr >= MinTrk) {
 			VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
@@ -139,11 +150,14 @@
 			Double_t Chi2 = Vtx->GetVtxChi2();
 			Double_t Ndof = 2 * (Double_t)Ntr;
+			delete Vtx;
+			//
 			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;
+			//
+			Double_t PullXv = (xpvc(0)-xpv) / TMath::Sqrt(covpvc(0, 0));
+			Double_t PullYv = (xpvc(1)-ypv) / TMath::Sqrt(covpvc(1, 1));
+			Double_t PullZv = (xpvc(2)-zpv) / TMath::Sqrt(covpvc(2, 2));
+			//
 			//
 			// Fill histograms
@@ -152,7 +166,10 @@
 			hZpull->Fill(PullZ);
 			hChi2->Fill(Chi2 / Ndof);
+			//
+			hXvpull->Fill(PullXv);
+			hYvpull->Fill(PullYv);
+			hZvpull->Fill(PullZv);
 		}
 
-		//std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
 		//
 		// Cleanup
@@ -162,4 +179,5 @@
 		delete[] cv;
 	}
+
 	//
 	// Show resulting histograms
@@ -174,3 +192,12 @@
 	hZpull->Fit("gaus"); hZpull->Draw();
 	Cnv->cd(4); hChi2->Draw();
+	//
+	TCanvas* Cnv1 = new TCanvas("Cnv1", "Generated primary vertex pulls", 100, 100, 900, 500);
+	Cnv1->Divide(3, 1);
+	Cnv1->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hXvpull->Fit("gaus"); hXvpull->Draw();
+	Cnv1->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hYvpull->Fit("gaus"); hYvpull->Draw();
+	Cnv1->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
+	hZvpull->Fit("gaus"); hZvpull->Draw();
 }
