// read muDst & geant.root in sync
//#define DEBUG
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "Riostream.h"
#include <stdio.h>
#include "TSystem.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TStyle.h"
#include "TF1.h"
#include "TTree.h"
#include "TChain.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TCanvas.h"
#include "TMinuit.h"
#include "TSpectrum.h"
#include "TString.h"
#include "TLine.h"
#include "TText.h"
#include "TROOT.h"
#include "TList.h"
#include "TPolyMarker.h"
#include "StBichsel/Bichsel.h"
#include "BetheBloch.h"
#include "TDirIter.h"
#include "TTreeIter.h"
#include "TLorentzVector.h"
#include "TVector3.h"
//#include "StThreeVectorF.hh"
//#include "StChain.h"
#else
class TSystem;
class TMath;
class TH1;
class TH2;
class TH3;
class TProfile;
class TStyle;
class TF1;
class TTree;
class TChain;
class TFile;
class TNtuple;
class TCanvas;
class TMinuit;
class TSpectrum;
class TString;
class TLine;
class TText;
class TROOT;
class TList;
class TPolyMarker;
class StChain;
#endif
class StChain;
StChain *chain=0;
//#include "tables/St_g2t_ssd_hit_Table.h"

void mygeant_out_allLD0s(const char* listName="1kPureD0geant.list",
					  const char* outName="myresults_geant.root"){
	//cout << "In mygeant_out_all2.C" << endl;
	ifstream inFile; 
	//Char_t *t = new Char_t[200]; 
	//Int_t NFiles = 0; 
	//inFile.open(listName); 
	/*while(NFiles < 10){
		inFile >> t;  
		cout << "NFiles = "<< NFiles << "  t=" <<t << endl;  
		if (!inFile.good()) break; 
		//Char_t *file = t;
		NFiles++; 
	}*/
	//const Int_t MAXFILES = NFiles; 
	//inFile.close();  
	//Int_t max = 400*(NFiles+1);
	//printf("maxfiles =%d maxEvent = %d\n",MAXFILES,max);
	//cout << "\n fileList :" <<listName<<endl;  
  
	TFile *myFile = new TFile(outName,"RECREATE");
	Float_t tuple[40];
	for(int counter=0;counter<40;counter++){
		tuple[counter]=0.0;
	}
	TNtuple *nTuple = new TNtuple("GeantTree","GeantTree","fileId:eventId:runId:primVrtxX:primVrtxY:primVrtxZ:parentGeantId:ptD0:pXD0:pYD0:pZD0:etaD0:cosThetaStar:secVrtxX:secVrtxY:secVrtxZ:KeyK:ptK:pXK:pYK:pZK:etaK:dcaXYK:dcaZK:KeyPi:ptPi:pXPi:pYPi:pZPi:etaPi:dcaXYPi:dcaZPi");
  
	FILE *readMe = fopen("READMEgeant.txt","w");{		// create file with the layout for data stored in text file
		fprintf(readMe,"The data stored in each row corresponds to data for the parent and daughter particles for one event. \n");
		fprintf(readMe,"Each variable will be separated by a tab space.  The list is as follows: \n");
		fprintf(readMe,"FileId \t EventId \t RunId \t VertexX \t VertexY \t VertexZ \t ParentId \t ParentPt \t ParentPx \t ParentPy \t ParentPz \t ParentEta \t cosThetaStar \t VrtxDX \t VrtxDY \t VrtxDZ \t KeyK \t PtK \t PxK \t PyK \t PzK \t EtaK \t dcaXYK \t dcaZK \t KeyPi \t PtPi \t PxPi \t PyPi \t PzPi \t EtaPi dcaXYPi \t dcaZPi \t decayLengthD \n");
	}
	fclose(readMe);
	cout << "************** README Written ************" << endl;
	//book histos here
  
	//TH1F *h1      = new TH1F("h1","D0 p_{T}",200,0,10);
	//TH1F *h2      = new TH1F("h2","D0 #eta",160,-4,4);
	//TH1F *h3      = new TH1F("h3","Particles Id",500,0,50);
	//TH1F *h4      = new TH1F("h4","distribution of the angle between pt of kaon and pion",360,0,360);
	//TH1F *h6      = new TH1F("h6","distribution of decaylength",200,-0.1,0.1);
	//TH1F *h7      = new TH1F("h7","event no. distribution",400,0,400);
	//TH1F *h8      = new TH1F("h8","Pz of the parent D0 from geant",100,-5,5);
	//TH1F *h9      = new TH1F("h9","Phi distibution of the parent D0",200,-10,10);
  	//TH1F *h10     = new TH1F("h10","distribution of primary vertex position_R",100,-5,5);
	//TH1F *h11     = new TH1F("h11","distribution of decay vertex position_R",100,-5,5);
	//TH1F *h12     = new TH1F("h12","geant ID",100,0,100);
	//TH1F *h122    = new TH1F("h122","geant ID for first track =37",100,0,100);
	//TH1F *hpTD0GEANT_[MAXFILES];
	//TH1F *hDecay_[MAXFILES];
	//TString title;
	//Char_t *name   = new Char_t[20];
	//Char_t *name2  = new Char_t[20];
	//for(Int_t ii=0;ii<MAXFILES;ii++){
		//title = Form("p_{T} of D0 per event for the file : %i",ii);
		//sprintf(name,"%s%d","hpTD0GEANT_",ii);
		//hpTD0GEANT_[ii] =  new TH1F(name,title,500,0,500);
		//title = Form("decay per event for the file : %i",ii);
		//sprintf(name2,"%s%d","hDecay_",ii);
		//hDecay_[ii] =  new TH1F(name2,title,500,0,500);
	//}  
  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  cout << " loading done " << endl;
	Int_t NumFiles = 0;
	Int_t nFilesUsed = -1;
	Char_t *ti = new Char_t[200];
	//cout << listName << endl; 
	inFile.open(listName);
	FILE *geantData = fopen("geantData.txt","w");{		// create file with the layout for data stored in text file
 
	Int_t D0found = 0;
	//while(NumFiles < MAXFILES){
		
	ifstream keysFile;
	keysFile.open("miniMcKeys.txt");
		
	while(NumFiles < 1){
		TString geantFile;
		Int_t fileIdIntGeant = 0;
		cout << NumFiles << endl;
		inFile >> ti;  
		cout << "ti  = " <<ti << endl;
		//cout << inFile.good() << endl;
		TString fileIdGeant(ti);
		fileIdGeant.ReplaceAll("/star/institutions/ksu/fisyak/D0B/rcf1272_","");
		fileIdGeant.ReplaceAll("_400evts_37.geant.root","");
		fileIdIntGeant = atoi(fileIdGeant);
		//cout << "fileIdIntGeant = " << fileIdIntGeant << endl;
		if (!inFile.good()) break;
		//cout << "made it past break" << endl; 
		//Char_t *file = ti;
		geantFile = ti;
		// create chain    
		chain = new StChain("StChain");   
		//chain->SetDebug();
		//chain->PrintInfo();
    
		//set up maker in read mode
		StIOMaker* IOMk = new StIOMaker("IO","r",geantFile,"bfcTree");
		//StIOMaker* IOMk = new StIOMaker("IO","r",file,"bfcTree");
		IOMk->SetDebug();
		IOMk->SetIOMode("r");
		IOMk->SetBranch("*",0,"0");
		IOMk->SetBranch("geantBranch",0,"r");
		printf("I'm opening file # %d\n",NumFiles);

		chain->Init();
		//chain->ls(3);
	
		Int_t FileId = 0;
		Int_t mEventId = 0;
		Int_t mRunId = 0;
		Float_t mVertexX = 0.0;
		Float_t mVertexY = 0.0;
		Float_t mVertexZ = 0.0;
		Float_t StopRDMc = 0.0;
		Float_t PtDMc = 0.0;
		Float_t PxDMc = 0.0;
		Float_t PyDMc = 0.0;
		Float_t PzDMc = 0.0;
		Float_t EtaDMc = 0.0;
		Float_t PhiDMc = 0.0;
		Int_t GeantK = 0;
		Float_t StopRK = 0;
		Int_t KeyK = 0;
		Int_t RecoKeyK = 0;
		Float_t PtKPr = 0.0;
		Float_t PxKPr = 0.0;
		Float_t PyKPr = 0.0;
		Float_t PzKPr = 0.0;
		Float_t EtaKPr = 0.0;
		Float_t PhiKPr = 0.0;
		Float_t DCAKPr = 0.0;
		Float_t DCAxyKPr = 0.0;
		Float_t DCAzKPr = 0.0;
		Float_t PtKGl = 0.0;
		Float_t PxKGl = 0.0;
		Float_t PyKGl = 0.0;
		Float_t PzKGl = 0.0;
		Float_t EtaKGl = 0.0;
		Float_t PhiKGl = 0.0;
		Float_t DCAKGl = 0.0;
		Float_t DCAxyKGl = 0.0;
		Float_t DCAzKGl = 0.0;
		Int_t GeantPi = 0;
		Float_t StopRPi = 0.0;
		Int_t KeyPi = 0;
		Int_t RecoKeyPi = 0;
		Float_t PtPiPr = 0.0;
		Float_t PxPiPr = 0.0;
		Float_t PyPiPr = 0.0;
		Float_t PzPiPr = 0.0;
		Float_t EtaPiPr = 0.0;
		Float_t PhiPiPr = 0.0;
		Float_t DCAPiPr = 0.0;
		Float_t DCAxyPiPr = 0.0;
		Float_t DCAzPiPr = 0.0;
		Float_t PtPiGl = 0.0;
		Float_t PxPiGl = 0.0;
		Float_t PyPiGl = 0.0;
		Float_t PzPiGl = 0.0;
		Float_t EtaPiGl = 0.0;
		Float_t PhiPiGl = 0.0;
		Float_t DCAPiGl = 0.0;
		Float_t DCAxyPiGl = 0.0;
		Float_t DCAzPiGl = 0.0;
	
		Int_t trackK = 0;
		Int_t trackPi = 0;
	
		keysFile >> FileId;
		keysFile >> mEventId;
		//cout << "Find mEventId = " << mEventId << endl;
		keysFile >> mRunId;
		keysFile >> mVertexX;
		keysFile >> mVertexY;
		keysFile >> mVertexZ;
		keysFile >> StopRDMc;
		keysFile >> PtDMc;
		keysFile >> PxDMc;
		keysFile >> PyDMc;
		keysFile >> PzDMc;
		keysFile >> EtaDMc;
		keysFile >> PhiDMc;
		keysFile >> GeantK;
		keysFile >> StopRK;
		keysFile >> KeyK;
		keysFile >> RecoKeyK;
		//cout << "Find Track RecoKeyK-1 = " << RecoKeyK-1 << endl;
		keysFile >> PtKPr;
		keysFile >> PxKPr;
		keysFile >> PyKPr;
		keysFile >> PzKPr;
		keysFile >> EtaKPr;
		keysFile >> PhiKPr;
		keysFile >> DCAKPr;
		keysFile >> DCAxyKPr;
		keysFile >> DCAzKPr;
		keysFile >> PtKGl;
		keysFile >> PxKGl;
		keysFile >> PyKGl;
		keysFile >> PzKGl;
		keysFile >> EtaKGl;
		keysFile >> PhiKGl;
		keysFile >> DCAKGl;
		keysFile >> DCAxyKGl;
		keysFile >> DCAzKGl;
		keysFile >> GeantPi;
		keysFile >> StopRPi;
		keysFile >> KeyPi;
		keysFile >> RecoKeyPi;
		//cout << "Find Track RecoKeyPi-1 = " <<RecoKeyPi-1 << endl;
		keysFile >> PtPiPr;
		keysFile >> PxPiPr;
		keysFile >> PyPiPr;
		keysFile >> PzPiPr;
		keysFile >> EtaPiPr;
		keysFile >> PhiPiPr;
		keysFile >> DCAPiPr;
		keysFile >> DCAxyPiPr;
		keysFile >> DCAzPiPr;
		keysFile >> PtPiGl;
		keysFile >> PxPiGl;
		keysFile >> PyPiGl;
		keysFile >> PzPiGl;
		keysFile >> EtaPiGl;
		keysFile >> PhiPiGl;
		keysFile >> DCAPiGl;
		keysFile >> DCAxyPiGl;
		keysFile >> DCAzPiGl;
	
		Bool_t foundEvent = 0;
	
		//---------------------------------------------------
		while ( 1) {// loop over events
		
			chain->Clear();
			int stat = chain->Make();
			if(stat) break;
	  
			// Access to geant-tables .......................
      
			St_DataSet* Event = chain->GetDataSet("geant");
			//Event->ls(0);
			St_DataSetIter geant(Event);
			St_g2t_vertex  *g2t_vertexTablePointer  =  (St_g2t_vertex *)  geant("g2t_vertex");
			St_g2t_event   *g2t_eventTablePointer   =  (St_g2t_event *)  geant("g2t_event"); 
      
			g2t_vertex_st  *vertexTable = 0;
			g2t_event_st *eventTable = 0;
			if(g2t_eventTablePointer){
				eventTable = g2t_eventTablePointer->GetTable();
				Int_t Nevents = g2t_eventTablePointer->GetNRows();
				//cout << "Nevents = " << Nevents << endl;
			}
			else
				cout << "Table g2t_event Not found in Dataset " << geant.Pwd()->GetName() << endl;
			if (!g2t_vertexTablePointer)continue;
			if (g2t_vertexTablePointer){
				vertexTable  = g2t_vertexTablePointer->GetTable();
				Int_t Nvertex = g2t_vertexTablePointer->GetNRows();
				//printf("Nvertex  = %d\n",Nvertex);
			}
			else
				cout << "Table g2t_vertex Not found in Dataset " << geant.Pwd()->GetName() << endl;
			St_g2t_track   *g2t_trackTablePointer   =  (St_g2t_track   *) geant("g2t_track");
			if(g2t_trackTablePointer){
				g2t_track_st *trackTable = g2t_trackTablePointer->GetTable();    
				Int_t NTracks            = g2t_trackTablePointer->GetNRows();
				//printf("NTracks  = %d\n",NTracks);
			}
			else                       
				cout << "Table g2t_tracks Not found in Dataset " << geant.Pwd()->GetName() << endl;
			if (!g2t_trackTablePointer)continue;
			cout << "eventTable[0].n_event = " << eventTable[0].n_event << " mEventId = " << mEventId << " fileIdIntGeant = " << fileIdIntGeant << " FileId = " << FileId << " eventTable[0].n_run+1 = " << eventTable[0].n_run+1 << " Run = " << mRunId << endl;
			if (eventTable[0].n_event == mEventId && fileIdIntGeant == FileId && eventTable[0].n_run+1 == mRunId) foundEvent = 1;
			//cout << "Found Event" << endl;
			//if(mEventId > 5) continue;
	  
			St_g2t_vertex *vertex = (St_g2t_vertex *) chain->FindObject("g2t_vertex");
			g2t_vertex_st *vrtx = vertex->GetTable();
	  
			//cout << "VertexX = " << vrtx[0].ge_x[0] << " VertexY = " << vrtx[0].ge_x[1] << " VertexZ = " << vrtx[0].ge_x[2] << endl;
	  
			if(!foundEvent) continue;
			
			//Int_t VrtxDId = 0;
			//VrtxDId = trackTable[trackK].start_vertex_p;
			//if(VrtxDId != trackTable[trackPi].start_vertex_p) continue;
			//cout << "VrtxDId = " << VrtxDId << endl;
	  
			if (foundEvent){  
				trackK = (KeyK-1);
				trackPi = (KeyPi-1);
				
				Float_t dcaXYK = 0.0;
				Float_t dcaXYPi = 0.0;
				Float_t rK = 0.0;
				Float_t rPi = 0.0;
				Float_t dcaZK = 0.0;
				Float_t dcaZPi = 0.0;
				Float_t x0 = 0.0;
				Float_t y0 = 0.0;
				Float_t x1K = 0.0;
				Float_t y1K = 0.0;
				Float_t x1Pi = 0.0;
				Float_t y1Pi = 0.0;
				Float_t x2 = 0.0;
				Float_t y2 = 0.0;
				Float_t ptK = 0.0;
				Float_t ptPi = 0.0;
				Float_t distVtxK = 0.0;
				Float_t distVtxPi = 0.0;
			
				// dcaZ calculation
				// Consider particle moving in circle.  Pt points in a tangent to the circle.  Radius of circle r = Pt/(0.3*B). B = mag field strength
				// Draw line from secondary vertex (x,y) to center of circle perpendicular to Pt of course.  Find (x,y) of center of circle.
				// dcaZ = r - (distance between prim vtx and center of circle)
				// (x0,y0) = prim vtx	(x1,y1) = center of circle	(x2,y2) = secondary vtx
				x0 = vrtx[0].ge_x[0];
				y0 = vrtx[0].ge_x[1];
				x2 = vrtx[1].ge_x[0];
				y2 = vrtx[1].ge_x[1];
								
				ptK = TMath::Sqrt(trackTable[trackK].p[0]*trackTable[trackK].p[0]+trackTable[trackK].p[1]*trackTable[trackK].p[1]);
				ptPi = TMath::Sqrt(trackTable[trackPi].p[0]*trackTable[trackPi].p[0]+trackTable[trackPi].p[1]*trackTable[trackPi].p[1]);
				rK = ptK/(0.3*4.9798);
				rPi = ptPi/(0.3*4.9798);
				
				x1K = x2 - TMath::Sqrt(rK*rK/((trackTable[trackK].p[0]*trackTable[trackK].p[0]/trackTable[trackK].p[1]*trackTable[trackK].p[1])+1));
				y1K = (trackTable[trackK].p[0]/trackTable[trackK].p[1])*(x2 - x1K) + y2;
				x1Pi = x2 + TMath::Sqrt(rPi*rPi/((trackTable[trackPi].p[0]*trackTable[trackPi].p[0]/trackTable[trackPi].p[1]*trackTable[trackPi].p[1])+1));
				y1Pi = (trackTable[trackPi].p[0]/trackTable[trackPi].p[1])*(x2 - x1Pi) + y2;
				
				distVtxK = TMath::Sqrt((x0 - x1K)*(x0 - x1K) + (y0 - y1K)*(y0 - y1K));
				distVtxPi = TMath::Sqrt((x0 - x1Pi)*(x0 - x1Pi) + (y0 - y1Pi)*(y0 - y1Pi));
				
				dcaZK = rK - distVtxK;
				dcaZPi = rPi - distVtxPi;
				
				dcaXYK = ((x2-x0)*trackTable[trackK].p[1]-(y2-y0)*trackTable[trackK].p[0])/TMath::Abs(ptK);
				dcaXYPi = ((x2-x0)*trackTable[trackPi].p[1]-(y2-y0)*trackTable[trackPi].p[0])/TMath::Abs(ptPi);
				
				//cout << "trackTable[trackK].id = " << trackTable[trackK].id << " KeyK = " << KeyK << " trackTable[trackPi].id = " << trackTable[trackPi].id << " KeyPi = " << KeyPi << endl;
				//cout << "eventTable[0].n_run = " << eventTable[0].n_run << " mRunId = " << mRunId << endl;
				//cout << "SecondaryVertexX = " << vrtx[1].ge_x[0] << " SecondaryVertexY = " << vrtx[1].ge_x[1] << " SecondaryVertexZ = " << vrtx[1].ge_x[2] << endl;
				//cout << "DCAxyK = " << dcaXYK << endl;
				//cout << "DCAxyPi = " << dcaXYPi << endl;
				//cout << "Rk = " << rK << endl;
				//cout << "DCAzK = " << dcaZK << endl;
				//cout << "Rpi = " << rPi << endl;
				//cout << "DCAzPi = " << dcaZPi << endl << endl;
				
				// Cos(ThetaStar) calculation.  ThetaStar is Theta for the Kaon in the D0 center of mass frame.
				TLorentzVector pD0;
				TLorentzVector pK;
				Float_t phiD0 = 0.0;
				Float_t phiK = 0.0;
				//Float_t ptK = 0.0;
				Float_t etaK = 0.0;
				//Float_t ptPi = 0.0;
				Float_t etaPi = 0.0;
				Float_t cosThetaStar;
				TVector3 boostD0;
				phiD0 = TMath::ATan2(trackTable[0].p[1],trackTable[0].p[0]);
				pD0.SetPtEtaPhiE(trackTable[0].pt, trackTable[0].eta, phiD0, trackTable[0].e);
				//cout << "pD0pt = " << trackTable[0].pt << endl;
				//cout << "pD0eta = " << trackTable[0].eta << endl;
				//cout << "pD0e = " << trackTable[0].e << endl;
				boostD0 = pD0.BoostVector();
				cout << "boostD0x = " << boostD0.X() << " boostD0y = " << boostD0.Y() << " boostD0z = " << boostD0.Z() << endl;
				TVector3 b(-boostD0.X(),-boostD0.Y(),-boostD0.Z());
				cout << "bx = " << b.X() << " by = " << b.Y() << " bz = " << b.Z() << endl;
				//ptK = TMath::Sqrt(trackTable[trackK].p[0]*trackTable[trackK].p[0]+trackTable[trackK].p[1]*trackTable[trackK].p[1]);
				phiK = TMath::ATan2(trackTable[trackK].p[1],trackTable[trackK].p[0]);
				etaK = -TMath::Log(TMath::Tan(TMath::ATan2(ptK,trackTable[trackK].p[2])/2));
				pK.SetPtEtaPhiE(ptK, etaK, phiK, trackTable[trackK].e);
				cout << "pKpt = " << ptK << endl;
				cout << "pKeta = " << etaK << endl;
				cout << "pKe = " << trackTable[trackK].e << endl;
				pK.Boost(b);
				cout << "pKx = " << pK.X() << endl;
				TVector3 daughter(pK.Vect());
				TVector3 mother(pD0.Vect());
				cosThetaStar = daughter.Dot(mother)/(daughter.Mag()*mother.Mag());
				
				//ptPi = TMath::Sqrt(trackTable[trackPi].p[0]*trackTable[trackPi].p[0]+trackTable[trackPi].p[1]*trackTable[trackPi].p[1]);
				etaPi = -TMath::Log(TMath::Tan(TMath::ATan2(ptPi,trackTable[trackPi].p[2])/2));
				
				cout << "cosThetaStar = " << cosThetaStar << endl << endl;
				
				if(KeyK == trackTable[trackK].id && KeyPi == trackTable[trackPi].id && mRunId == eventTable[0].n_run+1){
					//cout << "GeantEventId = " << eventTable[0].n_event << endl;
					//cout << "mEventId = " << mEventId <<" EventId = " << currentEvent << " KeyK = " << KeyK << " trackK geantId = " << trackTable[trackK].ge_pid << " KeyPi = " << KeyPi << " trackPi geantId = " << trackTable[trackPi].ge_pid << endl;
					//cout << "x = " << trackTable[trackK].p[0] << " PxPi = " << trackTable[trackPi].p[0] << endl;
					//cout << "PtD = " << trackTable[0].pt << " PzD = " << trackTable[0].p[2] << " EtaD = " << trackTable[0].eta << endl;
					//cout << "CosThetaStar = " << cosThetaStar << endl;
					//fprintf(geantData,"FileId \t EventId \t ParentId \t ParentPt \t ParentPz \t ParentEta \t GeantId1 \t Eta1 \t Pt1 \t Pz1 \t Key1 \t GeantId2 \t Eta2 \t Pt2 \t Pz2 \t Key2 \t");
					
					Float_t decayLengthD = TMath::Sqrt((vrtx[0].ge_x[0]-vrtx[1].ge_x[0])*(vrtx[0].ge_x[0]-vrtx[1].ge_x[0])+(vrtx[0].ge_x[1]-vrtx[1].ge_x[1])*(vrtx[0].ge_x[1]-vrtx[1].ge_x[1])+(vrtx[0].ge_x[2]-vrtx[1].ge_x[2])*(vrtx[0].ge_x[2]-vrtx[1].ge_x[2]));
					
					fprintf(geantData,"%i \t %i \t %i \t %.3f \t %.3f \t %.3f \t %i \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %i \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %i \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t %.3f \t  %.4f \n",
						fileIdIntGeant,eventTable[0].n_event,eventTable[0].n_run+1,vrtx[0].ge_x[0],vrtx[0].ge_x[1],vrtx[0].ge_x[2],trackTable[0].ge_pid,trackTable[0].pt,trackTable[0].p[0],trackTable[0].p[1],trackTable[0].p[2],trackTable[0].eta,cosThetaStar,vrtx[1].ge_x[0],vrtx[1].ge_x[1],vrtx[1].ge_x[2],
						KeyK,TMath::Sqrt(trackTable[trackK].p[0]*trackTable[trackK].p[0]+trackTable[trackK].p[1]*trackTable[trackK].p[1]),trackTable[trackK].p[0],trackTable[trackK].p[1],trackTable[trackK].p[2],etaK,dcaXYK,dcaZK,
						KeyPi,ptPi,trackTable[trackPi].p[0],trackTable[trackPi].p[1],trackTable[trackPi].p[2],etaPi,dcaXYPi,dcaZPi,decayLengthD);
					
					tuple[0] = fileIdIntGeant;
					tuple[1] = eventTable[0].n_event;
					tuple[2] = eventTable[0].n_run+1;
					tuple[3] = vrtx[0].ge_x[0];
					tuple[4] = vrtx[0].ge_x[1];
					tuple[5] = vrtx[0].ge_x[2];
					cout << "vrtxX = " << vrtx[0].ge_x[0] << " vrtxY = " << vrtx[0].ge_x[1] << " vrtxZ = " << vrtx[0].ge_x[2] << endl;
					cout << "vrtxXVrtxDId = " << vrtx[1].ge_x[0] << " vrtxYVrtxDId = " << vrtx[1].ge_x[1] << " vrtxZVrtxDId = " << vrtx[1].ge_x[2] << endl;
					cout << "vrtxXVrtxDId-1 = " << vrtx[1-1].ge_x[0] << " vrtxYVrtxDId-1 = " << vrtx[1-1].ge_x[1] << " vrtxZVrtxDId-1 = " << vrtx[1-1].ge_x[2] << endl;
					tuple[6] = trackTable[0].ge_pid;
					tuple[7] = trackTable[0].pt;
					tuple[8] = trackTable[0].p[0];
					tuple[9] = trackTable[0].p[1];
					tuple[10] = trackTable[0].p[2];
					tuple[11] = trackTable[0].eta;
					tuple[12] = cosThetaStar;
					tuple[13] = vrtx[1].ge_x[0];
					tuple[14] = vrtx[1].ge_x[1];
					tuple[15] = vrtx[1].ge_x[2];
					tuple[16] = KeyK;
					tuple[17] = ptK;
					tuple[18] = trackTable[trackK].p[0];
					tuple[19] = trackTable[trackK].p[1];
					tuple[20] = trackTable[trackK].p[2];
					tuple[21] = etaK;
					tuple[22] = dcaXYK;
					tuple[23] = dcaZK;
					tuple[24] = KeyPi;
					tuple[25] = ptPi;
					tuple[26] = trackTable[trackPi].p[0];
					tuple[27] = trackTable[trackPi].p[1];
					tuple[28] = trackTable[trackPi].p[2];
					tuple[29] = etaPi;
					tuple[30] = dcaXYPi;
					tuple[31] = dcaZPi;
					
					nTuple->Fill(tuple);
		
					keysFile >> FileId;
					keysFile >> mEventId;
					//cout << "Find mEventId = " << mEventId << endl;
					keysFile >> mRunId;
					keysFile >> mVertexX;
					keysFile >> mVertexY;
					keysFile >> mVertexZ;
					keysFile >> StopRDMc;
					keysFile >> PtDMc;
					keysFile >> PxDMc;
					keysFile >> PyDMc;
					keysFile >> PzDMc;
					keysFile >> EtaDMc;
					keysFile >> PhiDMc;
					keysFile >> GeantK;
					keysFile >> StopRK;
					keysFile >> KeyK;
					keysFile >> RecoKeyK;
					//cout << "Find Track RecoKeyK-1 = " << RecoKeyK-1 << endl;
					keysFile >> PtKPr;
					keysFile >> PxKPr;
					keysFile >> PyKPr;
					keysFile >> PzKPr;
					keysFile >> EtaKPr;
					keysFile >> PhiKPr;
					keysFile >> DCAKPr;
					keysFile >> DCAxyKPr;
					keysFile >> DCAzKPr;
					keysFile >> PtKGl;
					keysFile >> PxKGl;
					keysFile >> PyKGl;
					keysFile >> PzKGl;
					keysFile >> EtaKGl;
					keysFile >> PhiKGl;
					keysFile >> DCAKGl;
					keysFile >> DCAxyKGl;
					keysFile >> DCAzKGl;
					keysFile >> GeantPi;
					keysFile >> StopRPi;
					keysFile >> KeyPi;
					keysFile >> RecoKeyPi;
					//cout << "Find Track RecoKeyPi-1 = " <<RecoKeyPi-1 << endl;
					keysFile >> PtPiPr;
					keysFile >> PxPiPr;
					keysFile >> PyPiPr;
					keysFile >> PzPiPr;
					keysFile >> EtaPiPr;
					keysFile >> PhiPiPr;
					keysFile >> DCAPiPr;
					keysFile >> DCAxyPiPr;
					keysFile >> DCAzPiPr;
					keysFile >> PtPiGl;
					keysFile >> PxPiGl;
					keysFile >> PyPiGl;
					keysFile >> PzPiGl;
					keysFile >> EtaPiGl;
					keysFile >> PhiPiGl;
					keysFile >> DCAPiGl;
					keysFile >> DCAxyPiGl;
					keysFile >> DCAzPiGl;
	  
					foundEvent = 0;
					D0found++;
				}
			}  
		}
		//if(NumFiles > 0) break;
		NumFiles++;
		delete IOMk;
		delete chain;
		//if(eventTable[0].n_event == 1){
			//nFilesUsed++;
		//}
		//if(nFilesUsed > 2) break;
	}
		cout << "*********** Found " << D0found << " D0s. **************" << endl;
	}
	fclose(geantData);
	cout << "Geant Data Closed " << endl;
	inFile.close();
	myFile->Write();
	myFile->Close();
	cout << "Finished" << endl;
}