///////////////////////////////////////////////////////////////////////////////////////////////
//Create Histograms with cuts applied //
////////////////////////////////////////////////////////////////////////////////////////////////

#include <TChain.h>
#include <TTree.h>
#include <TH1F.h>
#include <TFile.h>

#include <TH1D.h>
#include <TH2D.h>
#include <map>
#include <string>
#include <iostream>
#include <fstream>

#include "TRandom3.h"
#include <TH1D.h>
#include <TCanvas.h>
#include "TMath.h"


using namespace std;


void LoadNtuple(TChain *chain, TFile *fL, string Ecut){

    //Declaring histograms of the quantities I am interested in
    //----------------------------------------------------------------------//

    TH1D *henergy = new TH1D("henergy", "henergy", 150,0,15);
    TH2D *hrenergy = new TH2D("hrenergy", "hrenergy",60, 0, 6000,100,0,10);   
    
    TH1D *hdx = new TH1D("hdx", "hdx", 200, -4000, 4000);
    TH1D *hdy = new TH1D("hdy", "hdy", 200, -4000, 4000);
    TH1D *hdz = new TH1D("hdz", "hdz", 200, -4000, 4000);
    TH1D *hdr = new TH1D("hdr", "hdr", 200, -4000, 4000);

///////////////////////////////////////////////////////////////////////////////////////  

    //Switch Back to file to write:
    fL->cd();
    
    //Create New Tree//
    TTree* outtree = new TTree("tree","tree");
    outtree->SetName("tree");

    double energy;
    double posx, posy, posz;
    bool fitValid;
    //Setting addresses to branches we are interested in using
    chain->SetBranchAddress("energy", &energy);
    chain->SetBranchAddress("posx", &posx);
    chain->SetBranchAddress("posy", &posy);
    chain->SetBranchAddress("posz", &posz);
    chain->SetBranchAddress("fitValid", &fitValid);
    
    //Define new branches
    double posr;
    outtree->Branch("energy", &energy);
    outtree->Branch("posx", &posx);
    outtree->Branch("posy", &posy);
    outtree->Branch("posz", &posz);
    outtree->Branch("posr", &posr);
    
    //Get number of entries
    ulong entries = chain->GetEntries();
    
  // Create header to store cut information
  TTree *header = new TTree("header", "cut info");
  int accepted = 0,accepted_valid = 0;

  double efficiency = 0,efficiency_valid = 0;

  int total = 0;
  header->Branch("efficiency", &efficiency);
  header->Branch("efficiency_valid", &efficiency_valid);
  header->Branch("accepted", &accepted);
  header->Branch("accepted_valid", &accepted_valid);
  header->Branch("total", &total);
   
    cerr<<"Begin processing: "<<entries<<" events"<<endl;
    for(ulong i=0; i<entries; i++) //for each event
    {
        total++;
        if(i%10000==0) //just a counter
        cerr<<"Proc: "<<i/double(entries)*100<<"%\r"<<flush;

        chain->GetEvent(i);
            
        //Fitter produces valid reconstruction result (All Values do now)
        if (!fitValid) continue; 
            accepted_valid+=1;
        
        //Energy Cut exactly at 3MeV (if specified --can make this more flexible)
        if (Ecut=="3"){
            double Ecut_double = atof(Ecut.c_str());
            if (energy<Ecut_double)continue;
        }
        
        accepted+=1; //get total number of events that passes all cuts
        
       
    posr = TMath::Sqrt(posx*posx + posy*posy + posz*posz);
        
    //Fill Histograms//    

    henergy->Fill(energy);
    henergy->GetXaxis()->SetTitle("E (RSP) [MeV]"); 

    hrenergy->Fill(posr,energy);
    hrenergy->GetYaxis()->SetTitle("E (RSP) [MeV]");
        
    hdx->Fill(posx);
    hdy->Fill(posy);
    hdz->Fill(posz);    
    hdr->Fill(posr);  
      
    outtree->Fill();
        
    }
    
    efficiency_valid = double(accepted_valid)/total;
    efficiency = double(accepted)/total;
    cout << "accepted      : " << accepted  << "  total: " << total << endl;
    cout << "efficiency    : " << efficiency<< "  efficiency_valid " << efficiency_valid<< endl;

    header->Fill();

    fL->Write(); //write file
    fL->Close();
}


void ChainMyEsim(string Ecut){
 
    /////////////////////////////////////////////////////////////////////////////////////
    string outputfilename_mc = "simout.root";
    const char *coutputfilename_mc = outputfilename_mc.c_str();
    TFile *f_outputfilename_mc = new TFile(coutputfilename_mc, "RECREATE");
    TChain *chainMC = new TChain("output");
    string Sim1 = "Sim1.root";
    string Sim2 = "Sim2.root";
    chainMC->Add(Sim1.c_str());
    chainMC->Add(Sim2.c_str());
    /////////////////////////////////////////////////////////////////////////////////////
    LoadNtuple(chainMC,f_outputfilename_mc, Ecut);
        
}