#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <time.h>
#include <cstdio>
#include <cstdlib>
#include <unistd.h>
#include <cstring>
#include <unistd.h>

#include "TTreeIndex.h"
#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"

#include "Event.h"
#include "Detector.h"
#include "Report.h"
#include "Vector.h"
#include "counting.hh"

using namespace std;

//Get effective volumes with error bars

class EarthModel;

int main(int argc, char **argv)
{
	gStyle->SetOptStat(111111);
	gStyle->SetOptDate(1);
	
	if(argc<4){
		cout<<"Not enough arguments! Abort run."<<endl;
		cout<<"Run like: ./veff volradius voldepth AraOut1.root AraOut2.rot ..."<<endl;
		return -1;
	}

	double volradius = atof(argv[1]);
	double voldepth = atof(argv[2]);

	TChain *AraTree = new TChain("AraTree");
	TChain *AraTree2 = new TChain("AraTree2");

	for(int i=3; i<argc;i++){
		AraTree->Add(string(argv[i]).c_str());
		AraTree2->Add(string(argv[i]).c_str());
	}

	Report *report = 0;
	Event *event=0;
	AraTree2->SetBranchAddress("report",&report);
	AraTree2->SetBranchAddress("event",&event);

	int totnthrown = AraTree2->GetEntries();
	cout << "Total number of events: " << totnthrown << endl;

	int NBINS=10;
	double eventsfound_binned[NBINS];
	for(int i=0; i<NBINS; i++) eventsfound_binned[i]=0.;

	double totweight=0;
	for(int iEvt2=0; iEvt2<totnthrown; iEvt2++){
		AraTree2->GetEntry(iEvt2);
		if(report->stations[0].Global_Pass<=0) continue;
		double weight = event->Nu_Interaction[0].weight;
		if(weight > 1.0) continue;
		totweight += weight;

		int index_weights = Counting::findWeightBin(log10(weight));
		if(index_weights<NBINS) eventsfound_binned[index_weights]++;
	}

	double error_minus=0.;
	double error_plus=0.;
	Counting::findErrorOnSumWeights(eventsfound_binned,error_plus,error_minus);

	double vtot = TMath::Pi() * double(volradius) * double(volradius) * double(voldepth) / 1.e9; //answer in km^3
	double veff = vtot * totweight / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
	double veff_p  = vtot * (error_plus) / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
	double veff_m  = vtot * (error_minus) / double(totnthrown) * 4. * TMath::Pi(); //answer in km^3 sr
                                                                                                                          
	printf("volthrown: %.6f \n  totweight: %.6f + %.6f - %.6f \n  Veff: %.6f + %.6f - %.6f \n",
			vtot,
			totweight, error_plus, error_minus,
			veff, veff_p, veff_m
			);
	return 0;
}