Index: examples/ExamplePVtxFind.C
===================================================================
--- examples/ExamplePVtxFind.C	(revision 46d34429082d61a6b55a5d706d0f7041cd5a0609)
+++ examples/ExamplePVtxFind.C	(revision dd263e4676f6a1140877fe9775508bee89267127)
@@ -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();
 }
