#include "Riostream.h" 
#include "TFile.h" 
#include <stdio.h>
#include "TNtuple.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"

void myCompareNTuples3kPureD0sHelices(){
	//cout << "Hello World" << endl;
	ifstream geantData;
	ifstream muDstData;

//Array keys
	Int_t i = 0;
	Int_t j=0;

//Variable keys from Geant
	Int_t  FileIdGeant = 0;
 	Int_t  EventIdGeant = 0;
	Int_t  RunIdGeant = 0; 
	Float_t  VertexXGeant = 0.0; 
	Float_t  VertexYGeant = 0.0; 
	Float_t  VertexZGeant = 0.0; 
	Float_t  ParentIdGeant = 0.0; 
	Float_t  ParentPtGeant = 0.0; 
	Float_t  ParentPxGeant = 0.0; 
	Float_t  ParentPyGeant = 0.0; 
	Float_t  ParentPzGeant = 0.0; 
	Float_t  ParentEtaGeant = 0.0;	 
	Float_t CosThetaStarGeant = 0.0;
	Float_t VrtxDXGeant = 0.0;
	Float_t VrtxDYGeant = 0.0;
	Float_t VrtxDZGeant = 0.0;
	Int_t  KeyKGeant = 0; 
	Float_t  PtKGeant = 0.0; 
	Float_t  PxKGeant = 0.0; 
	Float_t  PyKGeant = 0.0; 
	Float_t  PzKGeant = 0.0; 
	Float_t  EtaKGeant = 0.0; 
	Int_t  KeyPiGeant = 0; 
	Float_t  PtPiGeant = 0.0; 
	Float_t  PxPiGeant = 0.0; 
	Float_t  PyPiGeant = 0.0; 
	Float_t  PzPiGeant = 0.0; 
	Float_t  EtaPiGeant = 0.0;
	Float_t  dcaXYKGeant = 0.0;
	Float_t  dcaZKGeant = 0.0;
	Float_t  dcaXYPiGeant = 0.0;
	Float_t  dcaZPiGeant = 0.0;
	Float_t decayLengthDGeant = 0.0;
			
	//Variables Stored in muDstData.txt
	Int_t FileIdMuDst = 0;
	Int_t EventIdMuDst = 0;
	Int_t RunIdMuDst = 0;
	Float_t PrimVertexXMuDst = 0.0;
	Float_t PrimVertexYMuDst = 0.0;
	Float_t PrimVertexZMuDst = 0.0;
	Float_t PtDMuDst = 0.0;
	Float_t PxDMuDst = 0.0;
	Float_t PyDMuDst = 0.0;
	Float_t PzDMuDst = 0.0;
	Float_t EtaDMuDst = 0.0;
	Float_t MDMuDst = 0.0;
	Float_t CosThetaStarMuDst = 0.0;
	Float_t DvtxXMuDst = 0.0;
	Float_t DvtxYMuDst = 0.0;
	Float_t DvtxZMuDst = 0.0;
	Int_t TrackKMuDst = 0;
	Float_t PtKMuDst = 0.0;
	Float_t PzKMuDst = 0.0;
	Float_t EtaKMuDst = 0.0;
	Float_t dcaXYKMuDst = 0.0;
	Float_t dcaZKMuDst = 0.0;
	Int_t TrackPiMuDst = 0;
	Float_t PtPiMuDst = 0.0;
	Float_t PzPiMuDst = 0.0;
	Float_t EtaPiMuDst = 0.0;
	Float_t dcaXYPiMuDst = 0.0;
	Float_t dcaZPiMuDst = 0.0;
	Int_t noTracksGlMuDst = 0;
	Int_t noTracksPrMuDst = 0;
	Float_t decayLengthDMuDst = 0.0;
	Float_t sigmaXYK = 0.0;
	Float_t sigmaXYPi = 0.0;
	Float_t slength = 0.0;
	Float_t PxKMuDst = 0.0;
	Float_t PyKMuDst = 0.0;
	Float_t PxPiMuDst = 0.0;
	Float_t PyPiMuDst = 0.0;
	Float_t StopRDMuKpi = 0.0;
	Float_t secVtxXOrigMuKpi = 0.0;
	Float_t secVtxYOrigMuKpi = 0.0;
	Float_t secVtxZOrigMuKpi = 0.0;
	Float_t secVtxXHel1MuKpi = 0.0;
	Float_t secVtxYHel1MuKpi = 0.0;
	Float_t secVtxZHel1MuKpi = 0.0;
	Float_t decayLengthHel1MuKpi = 0.0;
	Float_t stopRDHel1MuKpi = 0.0;
	Float_t secVtxXHel2MuKpi = 0.0;
	Float_t secVtxYHel2MuKpi = 0.0;
	Float_t secVtxZHel2MuKpi = 0.0;
	Float_t decayLengthHel2MuKpi = 0.0;
	Float_t stopRDHel2MuKpi = 0.0;
	
	Int_t nEntriesGeant = 30;
	Int_t nEntriesMuDst = 46;
	Int_t nEntriesCompare = nEntriesGeant + nEntriesMuDst;
	
	//Create Array to Store Geant
	Float_t tupleGeant[50000][100];	//Array to store the data that will go into the NTuple
	for(Int_t line=0;line<50000;line++){
		for(Int_t entry=0;entry<100;entry++){
			tupleGeant[line][entry]=0.0;
		}
	}
	Float_t tupleGe[100];	//Array to store the data that will go into the NTuple
	for(Int_t entry=0;entry<100;entry++){
		tupleGe[entry]=0.0;
	}

	//Create Array to store MuDst
	Float_t tupleMuDst[100000][100];	//Array to store the data that will go into the NTuple
	for(Int_t line=0;line<100000;line++){		
		for(Int_t entry=0;entry<100;entry++){
			tupleMuDst[line][entry]=0.0;
		}
	}

	Float_t tupleMu[100];	//Array to store the data that will go into the NTuple for MuDst full
	for(Int_t entry=0;entry<100;entry++){
		tupleMu[entry]=0.0;
	}
	
	//Create Array to strore Comparison Data
	Float_t tupleCompare[100];	//Array to store the data that will go into the NTuple for MuDst full
	for(Int_t entry=0;entry<100;entry++){
		tupleMu[entry]=0.0;
	}
	
	//Create Array to strore Unmatched Data
	Float_t tupleNoMatchGeant[100];	//Array to store the data that will go into the NTuple for MuDst full
	for(Int_t entry=0;entry<100;entry++){
		tupleMu[entry]=0.0;
	}

	//Open Files for reading and writing
	TFile *fMuD = new TFile("myCompare3kPureD0sHelicesLinear.root","recreate");

	//Create Ntuples
	TNtuple *nTupleGeant = new TNtuple("GeantFull","GeantFull","FileIdGeant:EventIdGeant:RunIdGeant:VertexXGeant:VertexYGeant:VertexZGeant:ParentIdGeant:ParentPxGeant:ParentPyGeant:ParentPzGeant:ParentEtaGeant:CosThetaStarGeant:VrtxDXGeant:VrtxDYGeant:VrtxDZGeant:KeyKGeant:PxKGeant:PyKGeant:PzKGeant:EtaKGeant:dcaXYKGeant:dcaZKGeant:KeyPiGeant:PxPiGeant:PyPiGeant:PzPiGeant:EtaPiGeant:dcaXYPiGeant:dcaZPiGeant:decayLengthDGeant");
	TNtuple *nTupleCompare = new TNtuple("compare","compare","FileIdGeant:EventIdGeant:RunIdGeant:VertexXGeant:VertexYGeant:VertexZGeant:ParentIdGeant:ParentPxGeant:ParentPyGeant:ParentPzGeant:ParentEtaGeant:CosThetaStarGeant:VrtxDXGeant:VrtxDYGeant:VrtxDZGeant:KeyKGeant:PxKGeant:PyKGeant:PzKGeant:EtaKGeant:dcaXYKGeant:dcaZKGeant:KeyPiGeant:PxPiGeant:PyPiGeant:PzPiGeant:EtaPiGeant:dcaXYPiGeant:dcaZPiGeant:decayLengthDGeant:FileIdMuDst:EventIdMuDst:RunIdMuDst:PrimVertexXMuDst:PrimVertexYMuDst:PrimVertexZMuDst:PxDMuDst:PyDMuDst:PzDMuDst:EtaDMuDst:MDMuDst:CosThetaStarMuDst:TrackKMuDst:PzKMuDst:EtaKMuDst:dcaXYKMuDst:dcaZKMuDst:TrackPiMuDst:PzPiMuDst:EtaPiMuDst:dcaXYPiMuDst:dcaZPiMuDst:noTracksGlMuDst:noTracksPrMuDst:decayLengthDMuDst:sigmaXYK:sigmaXYPi:slength:PxKMuDst:PyKMuDst:PxPiMuDst:PyPiMuDst:StopRDMuKpi:secVtxXOrigMuKpi:secVtxYOrigMuKpi:secVtxZOrigMuKpi:secVtxXHel1MuKpi:secVtxYHel1MuKpi:secVtxZHel1MuKpi:decayLengthHel1MuKpi:stopRDHel1MuKpi:secVtxXHel2MuKpi:secVtxYHel2MuKpi:secVtxZHel2MuKpi:decayLengthHel2MuKpi:stopRDHel2MuKpi");
	TNtuple *nTupleMuDst = new TNtuple("MuDstFull","MuDstFull","FileIdMuDst:EventIdMuDst:RunIdMuDst:PrimVertexXMuDst:PrimVertexYMuDst:PrimVertexZMuDst:PxDMuDst:PyDMuDst:PzDMuDst:EtaDMuDst:MDMuDst:CosThetaStarMuDst:TrackKMuDst:PzKMuDst:EtaKMuDst:dcaXYKMuDst:dcaZKMuDst:TrackPiMuDst:PzPiMuDst:EtaPiMuDst:dcaXYPiMuDst:dcaZPiMuDst:noTracksGlMuDst:noTracksPrMuDst:decayLengthDMuDst:sigmaXYK:sigmaXYPi:slength:PxKMuDst:PyKMuDst:PxPiMuDst:PyPiMuDst:StopRDMuKpi:secVtxXOrigMuKpi:secVtxYOrigMuKpi:secVtxZOrigMuKpi:secVtxXHel1MuKpi:secVtxYHel1MuKpi:secVtxZHel1MuKpi:decayLengthHel1MuKpi:stopRDHel1MuKpi:secVtxXHel2MuKpi:secVtxYHel2MuKpi:secVtxZHel2MuKpi:decayLengthHel2MuKpi:stopRDHel2MuKpi");
	TNtuple *nTupleNoMatchGeant = new TNtuple("unmatchedGeant","unmatchedGeant","FileIdGeant:EventIdGeant:RunIdGeant:VertexXGeant:VertexYGeant:VertexZGeant:ParentIdGeant:ParentPxGeant:ParentPyGeant:ParentPzGeant:ParentEtaGeant:CosThetaStarGeant:VrtxDXGeant:VrtxDYGeant:VrtxDZGeant:KeyKGeant:PxKGeant:PyKGeant:PzKGeant:EtaKGeant:dcaXYKGeant:dcaZKGeant:KeyPiGeant:PxPiGeant:PyPiGeant:PzPiGeant:EtaPiGeant:dcaXYPiGeant:dcaZPiGeant:decayLengthDGeant");

	//Book some histos
	//TH1D *dVtxXMuDst  = new TH1D("decayVtxXMuDst","Decay Vertex X - Primary Vertex X MuDst",100,-0.2,0.2);
	//TH1D *dVtxYMuDst  = new TH1D("decayVtxYMuDst","Decay Vertex Y - Primary Vertex Y MuDst",100,-0.2,0.2);
	//TH1D *dVtxZMuDst  = new TH1D("decayVtxZMuDst","Decay Vertex Z - Primary Vertex Z MuDst",100,-0.2,0.2);
	//TH1D *dVtxXGeant  = new TH1D("decayVtxXGeant","Decay Vertex X - Primary Vertex X Geant",100,-0.2,0.2);
	//TH1D *dVtxYGeant  = new TH1D("decayVtxYGeant","Decay Vertex Y - Primary Vertex Y Geant",100,-0.2,0.2);
	//TH1D *dVtxZGeant  = new TH1D("decayVtxZGeant","Decay Vertex Z - Primary Vertex Z Geant",100,-0.2,0.2);
	//TH1D *dLMuDst = new TH1D("decayLengthMuDst","Decay Length MuDst",100,0.0,0.12);
	//TH1D *dLGeant = new TH1D("decayLengthGeant","Decay Length Geant",100,0.0,0.12);
	
	//TH1D *dVtxXRes  = new TH1D("decayVtxXRes","Decay Vertex X Geant - Decay Vertex X MuDst",100,-0.2,0.2);
	//TH1D *dVtxYRes  = new TH1D("decayVtxYRes","Decay Vertex Y Geant - Decay Vertex Y MuDst",100,-0.2,0.2);
	//TH1D *dVtxZRes  = new TH1D("decayVtxZRes","Decay Vertex Z Geant - Decay Vertex Z MuDst",100,-0.2,0.2);
	//TH1D *dLRes = new TH1D("decayLengthRes","Decay Length Geant - Decay Length MuDst",100,-0.12,0.12);
	
	//program start
	geantData.open("geantData3kPureD0sGlobal.txt");
	muDstData.open("muDstData3kPureD0sHelicesLinearGlobal2.txt");

	while (!geantData.eof()) {

		//File 3 is geant
		geantData >> FileIdGeant 
				  >> EventIdGeant 
				  >> RunIdGeant 
				  >> VertexXGeant 
				  >> VertexYGeant 
				  >> VertexZGeant 
				  >> ParentIdGeant 
				  >> ParentPtGeant 
				  >> ParentPxGeant 
				  >> ParentPyGeant 
				  >> ParentPzGeant 
				  >> ParentEtaGeant 
				  >> CosThetaStarGeant 
				  >> VrtxDXGeant 
				  >> VrtxDYGeant 
				  >> VrtxDZGeant 
				  >> KeyKGeant 
				  >> PtKGeant 
				  >> PxKGeant 
				  >> PyKGeant 
				  >> PzKGeant 
				  >> EtaKGeant 
				  >> dcaXYKGeant 
				  >> dcaZKGeant 
				  >> KeyPiGeant 
				  >> PtPiGeant 
				  >> PxPiGeant 
				  >> PyPiGeant 
				  >> PzPiGeant 
				  >> EtaPiGeant 
				  >> dcaXYPiGeant 
				  >> dcaZPiGeant 
				  >> decayLengthDGeant;
				
		if (i%4000==0)cout << "FileId: " << FileIdGeant << " Event Id: " << EventIdGeant << " i is : " << i << endl;
		//fill nTuple and save to array for Geant
			
		tupleGeant[i][0] = (Float_t)FileIdGeant;
		tupleGeant[i][1] = (Float_t)EventIdGeant;
		tupleGeant[i][2] = (Float_t)RunIdGeant;
		tupleGeant[i][3] = VertexXGeant;
		tupleGeant[i][4] = VertexYGeant;
		tupleGeant[i][5] = VertexZGeant;
		tupleGeant[i][6] = ParentIdGeant;
		tupleGeant[i][7] = ParentPxGeant;
		tupleGeant[i][8] = ParentPyGeant;
		tupleGeant[i][9] = ParentPzGeant;
		tupleGeant[i][10] = ParentEtaGeant;
		tupleGeant[i][11] = CosThetaStarGeant;
		tupleGeant[i][12] = VrtxDXGeant;
		tupleGeant[i][13] = VrtxDYGeant;
		tupleGeant[i][14] = VrtxDZGeant;
		tupleGeant[i][15] = (Float_t)KeyKGeant;
		tupleGeant[i][16] = PxKGeant;
		tupleGeant[i][17] = PyKGeant;
		tupleGeant[i][18] = PzKGeant;
		tupleGeant[i][19] = EtaKGeant;
		tupleGeant[i][20] = dcaXYKGeant;
		tupleGeant[i][21] = dcaZKGeant;
		tupleGeant[i][22] = (Float_t)KeyPiGeant;
		tupleGeant[i][23] = PxPiGeant;
		tupleGeant[i][24] = PyPiGeant;
		tupleGeant[i][25] = PzPiGeant;
		tupleGeant[i][26] = EtaPiGeant;
		tupleGeant[i][27] = dcaXYPiGeant;
		tupleGeant[i][28] = dcaZPiGeant;
		tupleGeant[i][29] = decayLengthDGeant;

		//Fill all geant Data into nTupleGeant
		for(Int_t k=0;k<nEntriesGeant;k++){
			tupleGe[k]=tupleGeant[i][k];
		}
		
		nTupleGeant->Fill(tupleGe);  //fill ntuple with values from array
		
		i++;
			
		//if (geantData.eof())break;
	}
	Int_t maxI = i+1;
	cout << "maxI = " << maxI << endl;
			
	//save muDst data to array

	while (!muDstData.eof()) {
		//Read MuDst Data
		muDstData >> FileIdMuDst
				  >> EventIdMuDst
				  >> RunIdMuDst
				  >> PrimVertexXMuDst
				  >> PrimVertexYMuDst
				  >> PrimVertexZMuDst
				  >> PtDMuDst
				  >> PxDMuDst
				  >> PyDMuDst
				  >> PzDMuDst
				  >> EtaDMuDst
				  >> MDMuDst
				  >> DvtxXMuDst
				  >> DvtxYMuDst
				  >> DvtxZMuDst
				  >> CosThetaStarMuDst
				  >> TrackKMuDst
				  >> PtKMuDst
				  >> PzKMuDst
				  >> EtaKMuDst 
				  >> dcaXYKMuDst 
				  >> dcaZKMuDst
				  >> TrackPiMuDst
				  >> PtPiMuDst
				  >> PzPiMuDst
				  >> EtaPiMuDst 
				  >> dcaXYPiMuDst 
				  >> dcaZPiMuDst
				  >> noTracksGlMuDst
				  >> noTracksPrMuDst
				  >> decayLengthDMuDst
				  >> sigmaXYK
				  >> sigmaXYPi
				  >> slength
				  >> PxKMuDst
				  >> PyKMuDst
		          >> PxPiMuDst
				  >> PyPiMuDst
				  >> StopRDMuKpi
				  >> secVtxXOrigMuKpi
				  >> secVtxYOrigMuKpi
				  >> secVtxZOrigMuKpi
				  >> secVtxXHel1MuKpi
				  >> secVtxYHel1MuKpi
				  >> secVtxZHel1MuKpi
				  >> decayLengthHel1MuKpi
				  >> stopRDHel1MuKpi
				  >> secVtxXHel2MuKpi
				  >> secVtxYHel2MuKpi
				  >> secVtxZHel2MuKpi
				  >> decayLengthHel2MuKpi
				  >> stopRDHel2MuKpi;

		if (j%4000==0) cout << "FileId: " << FileIdMuDst << " Event Id: " << EventIdMuDst <<  " j is : " << j <<  endl;

		//fill array
		tupleMuDst[j][0] = (Float_t)FileIdMuDst;  //turn Int_t type into Float_t
		tupleMuDst[j][1] = (Float_t)EventIdMuDst;
		tupleMuDst[j][2] = (Float_t)RunIdMuDst;
		tupleMuDst[j][3] = PrimVertexXMuDst;
		tupleMuDst[j][4] = PrimVertexYMuDst;
		tupleMuDst[j][5] = PrimVertexZMuDst;
		tupleMuDst[j][6] = PxDMuDst;
		tupleMuDst[j][7] = PyDMuDst;
		tupleMuDst[j][8] = PzDMuDst;
		tupleMuDst[j][9] = EtaDMuDst;
		tupleMuDst[j][10] = MDMuDst;
		tupleMuDst[j][11] = CosThetaStarMuDst;
		tupleMuDst[j][12] = (Float_t)TrackKMuDst;
		tupleMuDst[j][13] = PzKMuDst;
		tupleMuDst[j][14] = EtaKMuDst;
		tupleMuDst[j][15] = dcaXYKMuDst;
		tupleMuDst[j][16] = dcaZKMuDst;
		tupleMuDst[j][17] = (Float_t)TrackPiMuDst;
		tupleMuDst[j][18] = PzPiMuDst;
		tupleMuDst[j][19] = EtaPiMuDst;
		tupleMuDst[j][20] = dcaXYPiMuDst;
		tupleMuDst[j][21] = dcaZPiMuDst;
		tupleMuDst[j][22] = (Float_t)noTracksGlMuDst;
		tupleMuDst[j][23] = (Float_t)noTracksPrMuDst;
		tupleMuDst[j][24] = decayLengthDMuDst;
		tupleMuDst[j][25] = sigmaXYK;
		tupleMuDst[j][26] = sigmaXYPi;
		tupleMuDst[j][27] = slength;
		tupleMuDst[j][28] = PxKMuDst;
		tupleMuDst[j][29] = PyKMuDst;
		tupleMuDst[j][30] = PxPiMuDst;
		tupleMuDst[j][31] = PyPiMuDst;
		tupleMuDst[j][32] = StopRDMuKpi;
		tupleMuDst[j][33] = secVtxXOrigMuKpi;
		tupleMuDst[j][34] = secVtxYOrigMuKpi;
		tupleMuDst[j][35] = secVtxZOrigMuKpi;
		tupleMuDst[j][36] = secVtxXHel1MuKpi;
		tupleMuDst[j][37] = secVtxYHel1MuKpi;
		tupleMuDst[j][38] = secVtxZHel1MuKpi;
		tupleMuDst[j][39] = decayLengthHel1MuKpi;
		tupleMuDst[j][40] = stopRDHel1MuKpi;
		tupleMuDst[j][41] = secVtxXHel2MuKpi;
		tupleMuDst[j][42] = secVtxYHel2MuKpi;
		tupleMuDst[j][43] = secVtxZHel2MuKpi;
		tupleMuDst[j][44] = decayLengthHel2MuKpi;
		tupleMuDst[j][45] = stopRDHel2MuKpi;

		//Fill all MuDst Data into nTupleMuDst
		for(Int_t l=0;l<nEntriesMuDst;l++){
			tupleMu[l]=tupleMuDst[j][l];
		}

		nTupleMuDst->Fill(tupleMu); //fill ntuple with values from array
		
		j++;
		
		//if (muDstData.eof()) break;
	}
	
	Int_t maxJ = j+1;
	cout << "maxJ = " << maxJ << endl;
//End of reading into Arrays
 		
//Now for comparisons and matching Geant and MuDst events
	Bool_t FoundMatch = 0;
	Int_t nLineGeant = 0;
	Int_t nMatches = 0;
	Int_t nLineMuDst = 0;
	Int_t nLineSearch = 0;
	while(nLineGeant<maxI){
		//cout << "FileIdGeant = " << tupleGeant[nLineGeant][0] << " EventIdGeant = " << tupleGeant[nLineGeant][1] << " RunIdGeant = " << tupleGeant[nLineGeant][2] << endl;
		while(nLineMuDst<maxJ){
			//cout << "nLineGeant = " << nLineGeant << " nLineMuDst = " << nLineMuDst << " nMatches = " << nMatches << endl;
			//Compare FileId, EventID, and RunID
			if(tupleGeant[nLineGeant][0]==tupleMuDst[nLineMuDst][0] && 
			   tupleGeant[nLineGeant][1]==tupleMuDst[nLineMuDst][1] &&
			   tupleGeant[nLineGeant][2]==tupleMuDst[nLineMuDst][2]){
				FoundMatch = 1;
				if(nMatches%100==0) cout << "found match!  nLineGeant = " << nLineGeant << " nLineMuDst = " << nLineMuDst << endl;
				nMatches++;
				
				for(Int_t m=0;m<nEntriesCompare;m++){
					if(m<nEntriesGeant){
						tupleCompare[m] = tupleGeant[nLineGeant][m];
					}
					else{
						tupleCompare[m] = tupleMuDst[nLineMuDst][m-nEntriesGeant];
					}
				}
				nTupleCompare->Fill(tupleCompare);
				
				//dVtxXMuDst->Fill(tupleMuDst[nLineMuDst][12]-tupleMuDst[nLineMuDst][3]);
				//dVtxYMuDst->Fill(tupleMuDst[nLineMuDst][13]-tupleMuDst[nLineMuDst][4]);
				//dVtxZMuDst->Fill(tupleMuDst[nLineMuDst][14]-tupleMuDst[nLineMuDst][5]);
				//dVtxXGeant->Fill(tupleGeant[nLineGeant][13]-tupleGeant[nLineGeant][3]);
				//dVtxYGeant->Fill(tupleGeant[nLineGeant][14]-tupleGeant[nLineGeant][4]);
				//dVtxZGeant->Fill(tupleGeant[nLineGeant][15]-tupleGeant[nLineGeant][5]);
				//dLMuDst->Fill(tupleMuDst[nLineMuDst][30]);
				//dLGeant->Fill(tupleGeant[nLineGeant][32]);
				
				//dVtxXRes->Fill((tupleGeant[nLineGeant][13]-tupleGeant[nLineGeant][3])-(tupleMuDst[nLineMuDst][12]-tupleMuDst[nLineMuDst][3]));
				//dVtxYRes->Fill((tupleGeant[nLineGeant][14]-tupleGeant[nLineGeant][4])-(tupleMuDst[nLineMuDst][13]-tupleMuDst[nLineMuDst][4]));
				//dVtxZRes->Fill((tupleGeant[nLineGeant][15]-tupleGeant[nLineGeant][5])-(tupleMuDst[nLineMuDst][14]-tupleMuDst[nLineMuDst][5]));
				//dLRes->Fill(tupleGeant[nLineGeant][32]-tupleMuDst[nLineMuDst][30]);
				
				nLineSearch = nLineMuDst++;
				break;
			}
			//if(nLineMuDst%4000==0){
				//cout << "nLineGeant = " << nLineGeant << " nLineMuDst = " << nLineMuDst << endl;
			//}
			if(nLineMuDst == maxJ){
				for(Int_t n=0;n<nEntriesGeant;n++){
						tupleNoMatchGeant[n] = tupleGeant[nLineGeant][n];
				}
				nTupleNoMatchGeant->Fill(tupleNoMatchGeant);
			}
			nLineMuDst++;
			//if(tupleGeant[nLineGeant][0]!=tupleMuDst[nLineMuDst][0]) break;
			/*if(FoundMatch){
				nLineMuDst = 0;
				break;
			}*/
		}
		nLineMuDst = nLineSearch;
		nLineGeant++;
	}
	cout << "Geant has " << maxI << " lines.  MuDst has " << maxJ << " lines." << endl;
	cout << "Comparison complete.  Found " << nMatches << " matches." << endl;
	
	//Make Some Plots
	/*TCanvas *secVtxPlots = new TCanvas("secVtxPlots","Secondary Vertex Location wrt Primary Vertex",900,600);
	secVtxPlots->Divide(3,3);
	secVtxPlots->cd(1);
	dVtxXMuDst->Draw();
	secVtxPlots->cd(2);
	dVtxYMuDst->Draw();
	secVtxPlots->cd(3);
	dVtxZMuDst->Draw();
	secVtxPlots->cd(4);
	dVtxXGeant->Draw();
	secVtxPlots->cd(5);
	dVtxYGeant->Draw();
	secVtxPlots->cd(6);
	dVtxZGeant->Draw();
	secVtxPlots->cd(7);
	dVtxXRes->Draw();
	secVtxPlots->cd(8);
	dVtxYRes->Draw();
	secVtxPlots->cd(9);
	dVtxZRes->Draw();
	
	TCanvas *decayLengthPlots = new TCanvas("decayLengthPlots","Decay Length Plots and Resolution",900,600);
	decayLengthPlots->Divide(2,2);
	decayLengthPlots->cd(1);
	dLMuDst->Draw();
	decayLengthPlots->cd(2);
	dLGeant->Draw();
	decayLengthPlots->cd(3);
	dLRes->Draw();*/
	
	fMuD->Write();
	geantData.close();
	muDstData.close();
			
}