#include <cstdlib>
#include <iostream>
#include <vector>
#include <map>
#include <string>
using namespace std;

#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TStopwatch.h"
#include "TLatex.h"
#include "TH1.h"
#include "TH2.h"
#include "TChain.h"

#if not defined(__CINT__) || defined(__MAKECINT__)
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TMVA/MethodCuts.h"
#endif

using namespace TMVA;

void TMVA_L3BDTApplication(float myBDTcut = -0.05, TString path = "./", TString myFileName = "hltTree.root") 
{   
   
	//
	// create the Reader object
	//
	TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );

	// create a set of variables and declare them to the reader
	// - the variable names must corresponds in name and type to 
	// those given in the weight file(s) that you use
	float Ephoton_truth;     // From DBeamPhoton for MC data only (GeV)
	int numParticlesGeneratedByType[50]; // Number of particles of different types 
	
	float Nstart_counter;      // Number of start counter hits
	float Ntof;                // Number of TOF hits
	float Ntof_point;          // Number of TOF points
	
	float Nbcal_showers;      // Number of BCAL showers
	float EbcalShowers;       // Total energy in BCAL (Showers)
	
	float Nfcal_showers;      // Number of FCAL showers
	float EfcalShowers;       // Total energy in FCAL (Showers)
	
	float Ntrack_candidates;   // Number of track candidates
	float Ptot_tracks_cand;         // Scaler sum of total momentum from tracks
	
	bool L1_fired;
	
	float NtagRF;
	float Ntagm;
	float Ntagh;
	
	reader->AddVariable("Nstart_counter",&Nstart_counter);
	reader->AddVariable("Nfcal_showers",&Nfcal_showers);
	reader->AddVariable("EfcalShowers",&EfcalShowers);
	
	reader->AddVariable("Ntof_point",&Ntof_point);
	reader->AddVariable("Nbcal_showers",&Nbcal_showers);
	reader->AddVariable("EbcalShowers",&EbcalShowers);
	
	reader->AddVariable("Ntrack_candidates",&Ntrack_candidates);
	reader->AddVariable("Ptot_tracks_cand",&Ptot_tracks_cand);

	//
	// book the MVA methods
	//
	TString dir    = "weights/";
	TString prefix = "L3_BDT";
	TString weightfile = dir + prefix + ".weights.xml";
	reader->BookMVA( "BDT method", weightfile ); 
	
	//
	// book histograms
	//
	TH1F* hEphotonTruthAll = new TH1F("hEphotonTruthAll","",120,0.,12.);
	TH1F* hEphotonTruthAllNoAccept = new TH1F("hEphotonTruthAllNoAccept","",120,0.,12.);
	TH1F* hEphotonTruthAllProton = new TH1F("hEphotonTruthAllProton","",120,0.,12.);
	TH1F* hEphotonTruthAllNeutron = new TH1F("hEphotonTruthAllNeutron","",120,0.,12.);
	TH1F* hEphotonTruthL1 = new TH1F("hEphotonTruthL1","",120,0.,12.);
	TH1F* hEphotonTruthL1NoAccept = new TH1F("hEphotonTruthL1NoAccept","",120,0.,12.);
	TH1F* hEphotonTruthL1_tagged = new TH1F("hEphotonTruthL1_tagged","",120,0.,12.);
	TH1F* hEphotonTruthL1Proton = new TH1F("hEphotonTruthL1Proton","",120,0.,12.);
	TH1F* hEphotonTruthL1Neutron = new TH1F("hEphotonTruthL1Neutron","",120,0.,12.);
	TH1F* hEphotonTruthL3 = new TH1F("hEphotonTruthL3","",120,0.,12.);
	TH1F* hEphotonTruthL3NoAccept = new TH1F("hEphotonTruthL3NoAccept","",120,0.,12.);
	TH1F* hEphotonTruthL3Proton = new TH1F("hEphotonTruthL3Proton","",120,0.,12.);
	TH1F* hEphotonTruthL3Neutron = new TH1F("hEphotonTruthL3Neutron","",120,0.,12.);
	
	TH2F* hEphotonTruthAll_2D = new TH2F("hEphotonTruthAll_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthAllNoAccept_2D = new TH2F("hEphotonTruthAllNoAccept_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthAllProton_2D = new TH2F("hEphotonTruthAllProton_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthAllNeutron_2D = new TH2F("hEphotonTruthAllNeutron_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL1_2D = new TH2F("hEphotonTruthL1_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL1NoAccept_2D = new TH2F("hEphotonTruthL1NoAccept_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL1Proton_2D = new TH2F("hEphotonTruthL1Proton_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL1Neutron_2D = new TH2F("hEphotonTruthL1Neutron_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL3_2D = new TH2F("hEphotonTruthL3_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL3NoAccept_2D = new TH2F("hEphotonTruthL3NoAccept_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL3Proton_2D = new TH2F("hEphotonTruthL3Proton_2D","",120,0.,12.,1000,-1.0,1.0);
	TH2F* hEphotonTruthL3Neutron_2D = new TH2F("hEphotonTruthL3Neutron_2D","",120,0.,12.,1000,-1.0,1.0);
	
	TH2F* hTopologyTruthHighAll = new TH2F("hTopologyTruthHighAll",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentAll = new TH2F("hTopologyTruthCoherentAll",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighAllProton = new TH2F("hTopologyTruthHighAllProton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentAllProton = new TH2F("hTopologyTruthCoherentAllProton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighAllNeutron = new TH2F("hTopologyTruthHighAllNeutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentAllNeutron = new TH2F("hTopologyTruthCoherentAllNeutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	
	TH2F* hTopologyTruthHighL1 = new TH2F("hTopologyTruthHighL1",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL1 = new TH2F("hTopologyTruthCoherentL1",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighL1Proton = new TH2F("hTopologyTruthHighL1Proton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL1Proton = new TH2F("hTopologyTruthCoherentL1Proton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighL1Neutron = new TH2F("hTopologyTruthHighL1Neutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL1Neutron = new TH2F("hTopologyTruthCoherentL1Neutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	
	TH2F* hTopologyTruthHighL3 = new TH2F("hTopologyTruthHighL3",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL3 = new TH2F("hTopologyTruthCoherentL3",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighL3Proton = new TH2F("hTopologyTruthHighL3Proton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL3Proton = new TH2F("hTopologyTruthCoherentL3Proton",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthHighL3Neutron = new TH2F("hTopologyTruthHighL3Neutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	TH2F* hTopologyTruthCoherentL3Neutron = new TH2F("hTopologyTruthCoherentL3Neutron",";# Thrown Photons; # Thrown Charged",15,0.,15.,10,0.,10.);
	
	TH1F* hNphotonsAll = new TH1F("hNphotonsAll","",50,0,50);
	
	// monitoring for L3 trigger
	TH1F* hBdtClass = new TH1F("hBdtClass","Boosted Decision Tree classifier for L3 decision; BDT Classifier",100,-1.,1.);
	TH1F* hTrackMomKeep = new TH1F("hTrackMomKeep","Track momentum sum for events kept by L3; Track momentum sum (GeV/c)",120,0.,12.);
	TH1F* hTrackMomDiscard = new TH1F("hTrackMomDiscard","Track momentum sum for events discarded by L3; Track momentum sum (GeV/c)",120,0.,12.);
	TH2F* hEbcal_EfcalKeep = new TH2F("hEbcal_EfcalKeep","BCAL energy sum vs FCAL energy sum for events kept by L3; FCAL energy sum; BCAL energy sum",120,0.,6.,120,0.,6.);
	TH2F* hEbcal_EfcalDiscard = new TH2F("hEbcal_EfcalDiscard","BCAL energy sum vs FCAL energy sum for events discarded by L3; FCAL energy sum; BCAL energy sum",120,0.,6.,120,0.,6.);
	
	// correlations of thrown vs reco for right and wrong trigger decisions
	TH2F* hNTrack2DHighKeep = new TH2F("hNTrack2DHighKeep","hNTrack2DHighKeep;Thrown # charged tracks; Reconstructed # tracks",10,0,10,10,0,10);
	TH2F* hTrack2DHighKeep = new TH2F("hTrack2DHighKeep","hTrack2DHighKeep;Thrown Momentum; Reconstructed Momentum",100,0.,10.,100,0.,10.);
	TH2F* hBCAL2DHighKeep = new TH2F("hBCAL2DHighKeep","hBCAL2DHighKeep;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hFCAL2DHighKeep = new TH2F("hFCAL2DHighKeep","hFCAL2DHighKeep;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hNTrack2DHighDiscard = new TH2F("hNTrack2DHighDiscard","hNTrack2DHighDiscard;Thrown # charged tracks; Reconstructed # tracks",10,0,10,10,0,10);
	TH2F* hTrack2DHighDiscard = new TH2F("hTrack2DHighDiscard","hTrack2DHighDiscard;Thrown Momentum; Reconstructed Momentum",100,0.,10.,100,0.,10.);
	TH2F* hBCAL2DHighDiscard = new TH2F("hBCAL2DHighDiscard","hBCAL2DHighDiscard;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hFCAL2DHighDiscard = new TH2F("hFCAL2DHighDiscard","hFCAL2DHighDiscard;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	
	TH2F* hNTrack2DLowKeep = new TH2F("hNTrack2DLowKeep","hNTrack2DLowKeep;Thrown # charged tracks; Reconstructed # tracks",10,0,10,10,0,10);
	TH2F* hTrack2DLowKeep = new TH2F("hTrack2DLowKeep","hTrack2DLowKeep;Thrown Momentum; Reconstructed Momentum",100,0.,10.,100,0.,10.);
	TH2F* hBCAL2DLowKeep = new TH2F("hBCAL2DLowKeep","hBCAL2DLowKeep;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hFCAL2DLowKeep = new TH2F("hFCAL2DLowKeep","hFCAL2DLowKeep;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hNTrack2DLowDiscard = new TH2F("hNTrack2DLowDiscard","hNTrack2DLowDiscard;Thrown # charged tracks; Reconstructed # tracks",10,0,10,10,0,10);
	TH2F* hTrack2DLowDiscard = new TH2F("hTrack2DLowDiscard","hTrack2DLowDiscard;Thrown Momentum; Reconstructed Momentum",100,0.,10.,100,0.,10.);
	TH2F* hBCAL2DLowDiscard = new TH2F("hBCAL2DLowDiscard","hBCAL2DLowDiscard;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	TH2F* hFCAL2DLowDiscard = new TH2F("hFCAL2DLowDiscard","hFCAL2DLowDiscard;Thrown Energy; Reconstructed Energy",100,0.,10.,100,0.,10.);
	
	// more on events lost due to tracking
	TH1F* hTrackResidHighKeep = new TH1F("hTrackResidHighKeep","hTrackResidHighKeep; Charged track momentum residual (reco - thrown)",100,-10.,10.);
	TH1F* hTrackResidHighDiscard = new TH1F("hTrackResidHighDiscard","hTrackResidHighDiscard; Charged track momentum residual (reco - thrown)",100,-10.,10.);
	TH2F *hMissingTrackHighDiscard = new TH2F("hMissingTrackHighDiscard","hMissingTrackHighDiscard; thrown #theta (degrees); thrown momentum (GeV)",130,0,130,100,0,10);
	
	TChain *theTree = new TChain("theTree");
	theTree->AddFile(path+"hltTree.root",-1,"TrainingL3_Tree");
	
	theTree->SetBranchAddress("Ephoton_truth",&Ephoton_truth);
	theTree->SetBranchAddress("numParticlesGeneratedByType",numParticlesGeneratedByType);
	
	theTree->SetBranchAddress("Nstart_counter",&Nstart_counter);
	theTree->SetBranchAddress("Ntof",&Ntof);
	theTree->SetBranchAddress("Ntof_point",&Ntof_point); 
	
	theTree->SetBranchAddress("Nbcal_showers",&Nbcal_showers);	
	theTree->SetBranchAddress("EbcalShowers",&EbcalShowers);
   
	theTree->SetBranchAddress("Nfcal_showers",&Nfcal_showers);
	theTree->SetBranchAddress("EfcalShowers",&EfcalShowers); 
	
	theTree->SetBranchAddress("Ntrack_candidates",&Ntrack_candidates);
	theTree->SetBranchAddress("Ptot_tracks_cand",&Ptot_tracks_cand);
	
	theTree->SetBranchAddress("L1_fired",&L1_fired);
	
	Int_t Accept = 1;
	//theTree->SetBranchAddress("Accept",&Accept);
	
	double thrownTrackP = 0;
	double thrownTrackNeutronP = 0;
	double thrownNeutronP = 0;
	double thrownBCALE = 0;
	double thrownFCALE = 0;

	std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
	TStopwatch sw;
	sw.Start();
	int offset = -10000;
	for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
		
		if (ievt%10000 == 0){
			std::cout << "--- ... Processing event: " << ievt << std::endl;
			offset+=10000;
		}
		
		theTree->GetEntry(ievt);     
		thrownTrackNeutronP = thrownTrackP+thrownNeutronP;
	
		// deal with ints where required
		Nfcal_showers = (int)Nfcal_showers;
		Nbcal_showers = (int)Nbcal_showers;
		Ntrack_candidates = (int)Ntrack_candidates;
		
		//if(EfcalShowers!=EfcalShowers) EfcalShowers = 0;
		//if(EbcalShowers!=EbcalShowers) EbcalShowers = 0;
		//if(Ptot_tracks_cand!=Ptot_tracks_cand) Ptot_tracks_cand = 0;

		// cut used in training
		if(Ptot_tracks_cand>1e2) {
			//cout<<Ptot_tracks_cand<<endl;
			continue;
		}
		
		bool coherentPeak = false;
		bool highEnergy = false;
		if(Ephoton_truth > 7.0) highEnergy=true;
		if(Ephoton_truth >= 8.4 && Ephoton_truth <= 9.0) coherentPeak=true;
		
		int nPhotons = numParticlesGeneratedByType[1];
		hNphotonsAll->Fill(nPhotons);
		
		double BDT = reader->EvaluateMVA("BDT method");
		double minBDT = myBDTcut;
		
		// Egamma distribution with NO acceptance cuts
		hEphotonTruthAllNoAccept->Fill(Ephoton_truth);
		hEphotonTruthAllNoAccept_2D->Fill(Ephoton_truth,BDT);
		if(L1_fired) {
			hEphotonTruthL1NoAccept->Fill(Ephoton_truth);
			hEphotonTruthL1NoAccept_2D->Fill(Ephoton_truth,BDT);
			hEphotonTruthL3NoAccept_2D->Fill(Ephoton_truth,BDT);
			if(BDT > minBDT) hEphotonTruthL3NoAccept->Fill(Ephoton_truth);
		}
		
		// proton, anti-proton, pi+, pi-, K+, K-
		Int_t nCharged = numParticlesGeneratedByType[14] + numParticlesGeneratedByType[8] + numParticlesGeneratedByType[11] + numParticlesGeneratedByType[15] + numParticlesGeneratedByType[9] + numParticlesGeneratedByType[12];
		if(nCharged%2 == 0) continue; // skip even charged tracks for now, is bug in Geant?
		
		// acceptance criteria for thrown tracks in HLT.C
		if(!Accept) continue;
		
		hEphotonTruthAll->Fill(Ephoton_truth);
		hEphotonTruthAll_2D->Fill(Ephoton_truth,BDT);
		if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0) {
			hEphotonTruthAllProton->Fill(Ephoton_truth);
			hEphotonTruthAllProton_2D->Fill(Ephoton_truth,BDT);
			if(highEnergy) hTopologyTruthHighAllProton->Fill(nPhotons,nCharged);
			if(coherentPeak) hTopologyTruthCoherentAllProton->Fill(nPhotons,nCharged);
		}
		else {
			hEphotonTruthAllNeutron->Fill(Ephoton_truth);
			hEphotonTruthAllNeutron_2D->Fill(Ephoton_truth,BDT);
			if(highEnergy) hTopologyTruthHighAllNeutron->Fill(nPhotons,nCharged);
			if(coherentPeak) hTopologyTruthCoherentAllNeutron->Fill(nPhotons,nCharged);
		}
		if(highEnergy) hTopologyTruthHighAll->Fill(nPhotons,nCharged);
		if(coherentPeak) hTopologyTruthCoherentAll->Fill(nPhotons,nCharged);
		hBdtClass->Fill(BDT);
		
		// start cuts for L1 and L3 trigger here
		if(L1_fired) {
			//bool L1 = true; //cout<<"passed L1"<<endl;
		}
		else continue;
		
		// trigger monitoring plots
		if(BDT > minBDT){
			hTrackMomKeep->Fill(Ptot_tracks_cand);
			if(!highEnergy) hEbcal_EfcalKeep->Fill(EfcalShowers,EbcalShowers);
		}
		else {
			hTrackMomDiscard->Fill(Ptot_tracks_cand);
			if(!highEnergy) hEbcal_EfcalDiscard->Fill(EfcalShowers,EbcalShowers);
		}
		
		hEphotonTruthL1->Fill(Ephoton_truth);
		hEphotonTruthL1_2D->Fill(Ephoton_truth,BDT);
		
		if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0){
			hEphotonTruthL1Proton->Fill(Ephoton_truth);
			hEphotonTruthL1Proton_2D->Fill(Ephoton_truth,BDT);
			if(highEnergy) hTopologyTruthHighL1Proton->Fill(nPhotons,nCharged);
			if(coherentPeak) hTopologyTruthCoherentL1Proton->Fill(nPhotons,nCharged);
		}
		else {
			hEphotonTruthL1Neutron->Fill(Ephoton_truth);
			hEphotonTruthL1Neutron_2D->Fill(Ephoton_truth,BDT);
			if(highEnergy) hTopologyTruthHighL1Neutron->Fill(nPhotons,nCharged);
			if(coherentPeak) hTopologyTruthCoherentL1Neutron->Fill(nPhotons,nCharged);
		}
		if(highEnergy) hTopologyTruthHighL1->Fill(nPhotons,nCharged);
		if(coherentPeak) hTopologyTruthCoherentL1->Fill(nPhotons,nCharged);
		
		hEphotonTruthL3_2D->Fill(Ephoton_truth,BDT); 
		if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0) hEphotonTruthL3Proton_2D->Fill(Ephoton_truth,BDT);
		else hEphotonTruthL3Neutron_2D->Fill(Ephoton_truth,BDT);
		//fill some plots
		if(BDT > minBDT) {
			hEphotonTruthL3->Fill(Ephoton_truth);
			if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0) {
				hEphotonTruthL3Proton->Fill(Ephoton_truth);
				if(highEnergy) hTopologyTruthHighL3Proton->Fill(nPhotons,nCharged);
				if(coherentPeak) hTopologyTruthCoherentL3Proton->Fill(nPhotons,nCharged);
			}
			else {
				hEphotonTruthL3Neutron->Fill(Ephoton_truth);
				if(highEnergy) hTopologyTruthHighL3Neutron->Fill(nPhotons,nCharged);
				if(coherentPeak) hTopologyTruthCoherentL3Neutron->Fill(nPhotons,nCharged);
			}
			if(coherentPeak) hTopologyTruthCoherentL3->Fill(nPhotons,nCharged);
			if(highEnergy) hTopologyTruthHighL3->Fill(nPhotons,nCharged);
			
			if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0 && nPhotons==0) {
				if(coherentPeak) {
					if(thrownTrackP>=8.4 && thrownTrackP<=9.0) {
						hTrackResidHighKeep->Fill(thrownTrackP-Ptot_tracks_cand);
					}
					hNTrack2DHighKeep->Fill(nCharged,Ntrack_candidates);
					hTrack2DHighKeep->Fill(thrownTrackP,Ptot_tracks_cand);
					hBCAL2DHighKeep->Fill(thrownBCALE,EbcalShowers);
					hFCAL2DHighKeep->Fill(thrownFCALE,EfcalShowers);
				}
				if(!highEnergy) {
					hNTrack2DLowKeep->Fill(nCharged,Ntrack_candidates);
					hTrack2DLowKeep->Fill(thrownTrackP,Ptot_tracks_cand);
					hBCAL2DLowKeep->Fill(thrownBCALE,EbcalShowers);
					hFCAL2DLowKeep->Fill(thrownFCALE,EfcalShowers);
				}
			}
		}
		else {
			if(numParticlesGeneratedByType[14]>0 && numParticlesGeneratedByType[13] == 0 && nPhotons==0) {
				if(coherentPeak) {
					if(thrownTrackP>=8.4 && thrownTrackP<=9.0) {
						cout<<"Event: "<<ievt-offset<<" E_gamma="<<Ephoton_truth<<endl;
						hTrackResidHighDiscard->Fill(Ptot_tracks_cand-thrownTrackP);
						//hMissingTrackHighDiscard->Fill(thrownMaxMissingTrackTheta*180./3.14,thrownMaxMissingTrackP);
					}
					hNTrack2DHighDiscard->Fill(nCharged,Ntrack_candidates);
					hTrack2DHighDiscard->Fill(thrownTrackP,Ptot_tracks_cand);
					hBCAL2DHighDiscard->Fill(thrownBCALE,EbcalShowers);
					hFCAL2DHighDiscard->Fill(thrownFCALE,EfcalShowers);
				}
				if(!highEnergy) {
					hNTrack2DLowDiscard->Fill(nCharged,Ntrack_candidates);
					hTrack2DLowDiscard->Fill(thrownTrackP,Ptot_tracks_cand);
					hBCAL2DLowDiscard->Fill(thrownBCALE,EbcalShowers);
					hFCAL2DLowDiscard->Fill(thrownFCALE,EfcalShowers);
				}
			}
		}
	}

	// get elapsed time
	sw.Stop();
	std::cout << "--- End of event loop: "; sw.Print();
	
	//
	// write histograms
	//
	TString outFile = "TMVApp_"; outFile+=myFileName; 
	TFile *target  = new TFile( outFile,"RECREATE" );
	
	hEphotonTruthAll->Write();
	hEphotonTruthAllNoAccept->Write();
	hEphotonTruthAllProton->Write();
	hEphotonTruthAllNeutron->Write();
	hEphotonTruthL1->Write();
	hEphotonTruthL1NoAccept->Write();
	hEphotonTruthL1_tagged->Write();
	hEphotonTruthL1Proton->Write();
	hEphotonTruthL1Neutron->Write();
	hEphotonTruthL3->Write();
	hEphotonTruthL3NoAccept->Write();
	hEphotonTruthL3Proton->Write();
	hEphotonTruthL3Neutron->Write();

	hEphotonTruthAll_2D->Write();
	hEphotonTruthAllNoAccept_2D->Write();
	hEphotonTruthAllProton_2D->Write();
	hEphotonTruthAllNeutron_2D->Write();
	hEphotonTruthL1_2D->Write();
	hEphotonTruthL1NoAccept_2D->Write();
	hEphotonTruthL1Proton_2D->Write();
	hEphotonTruthL1Neutron_2D->Write();
	hEphotonTruthL3_2D->Write();
	hEphotonTruthL3NoAccept_2D->Write();
	hEphotonTruthL3Proton_2D->Write();
	hEphotonTruthL3Neutron_2D->Write();
	
	hTopologyTruthHighAll->Write();
	hTopologyTruthCoherentAll->Write();
	hTopologyTruthHighL1->Write();
	hTopologyTruthCoherentL1->Write();
	hTopologyTruthHighL3->Write();
	hTopologyTruthCoherentL3->Write();
	
	hTopologyTruthHighAllProton->Write();
	hTopologyTruthCoherentAllProton->Write();
	hTopologyTruthHighL1Proton->Write();
	hTopologyTruthCoherentL1Proton->Write();
	hTopologyTruthHighL3Proton->Write();
	hTopologyTruthCoherentL3Proton->Write();
	
	hTopologyTruthHighAllNeutron->Write();
	hTopologyTruthCoherentAllNeutron->Write();
	hTopologyTruthHighL1Neutron->Write();
	hTopologyTruthCoherentL1Neutron->Write();
	hTopologyTruthHighL3Neutron->Write();
	hTopologyTruthCoherentL3Neutron->Write();
	
	hNphotonsAll->Write();
	
	hBdtClass->Write();
	hTrackMomKeep->Write();
	hTrackMomDiscard->Write();
	hEbcal_EfcalKeep->Write();
	hEbcal_EfcalDiscard->Write();
	
	hNTrack2DHighKeep->Write();
	hTrack2DHighKeep->Write();
	hBCAL2DHighKeep->Write();
	hFCAL2DHighKeep->Write();
	hNTrack2DLowKeep->Write();
	hTrack2DLowKeep->Write();
	hBCAL2DLowKeep->Write();
	hFCAL2DLowKeep->Write();
	
	hNTrack2DHighDiscard->Write();
	hTrack2DHighDiscard->Write();
	hBCAL2DHighDiscard->Write();
	hFCAL2DHighDiscard->Write();
	hNTrack2DLowDiscard->Write();
	hTrack2DLowDiscard->Write();
	hBCAL2DLowDiscard->Write();
	hFCAL2DLowDiscard->Write();
	
	hTrackResidHighKeep->Write();
	hTrackResidHighDiscard->Write();
	hMissingTrackHighDiscard->Write();
	
	target->Close();
	std::cout << "--- Created root file: \"TMVApp.root\" containing the MVA output histograms" << std::endl;
	
	delete reader;
	std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
	
}