#include "MapsTrackManager.hh"
#include "MapsTrack.hh"
#include "MapsSensor.hh"
#include "MapsException.hh"
#include <cmath>
#include <algorithm>
#include <iostream>
#include "TRandom2.h"
#include "TFile.h"
#include "TH1F.h"
#include "Operators.h"

int main(int argc, const char **argv) {

	if (argc != 3) {
		std::cout << "Usage: ExtractEfficiencies <input file> <outputfile>\n";
		return 0;
	}
	std::cout << "Welcome to ExtractEfficiencies...\n";
	
	//Define chi squared probability cuts
	double cut(0.05);
	//Define errors for chi squared calculation
	std::pair<double, double> error(0.1, 0.1);
	//std::pair<double, double> error(0.026, 0.02);
	std::cout << "\tUsing cuts for probability of: \t" << cut << "\n";
	std::cout << "\tError parameters: \t" << error << "\n";
	
	std::cout << "\tRecreating from root file...\n";
	char input[100];
	strcpy(input, argv[1]);
	
	//Create a maps track manager to recreate the tracks
	MapsTrackManager mtm2;
	//Reinitialise all the tracks and sensors in mtm2
	mtm2.recreateFromRootFile(input);

	std::vector<MapsTrack*>& tracks = mtm2.getTracks();

	//count efficient and inefficient frequency
	unsigned ineff(0), eff(0), deadArs(0);
	unsigned candidateTracks(0);

	for (std::vector<MapsTrack*>::iterator it = tracks.begin(); it
			!= tracks.end(); it++) {
		MapsTrack* t = *it;

		//this is not obligatory, but highly advised!
		t->setTrackError(error);
		
		//see if we have a good three hit track first: for tracks which are just three hits anyway
		//t3 will be a copy of t
		MapsTrack t3;
		t->make3HitTrack(t3);
		t3.setTrackError(error);

		//evaluate quality of the 3 hit track
		if (t3.chiSqProb(0) > cut && t3.chiSqProb(1) > cut) {
			//it's a real track then
			candidateTracks++;

			//get the original track's fourth sensor
			MapsSensor* s4 = t->fourthSensor();

			//pose the ultimate question...
			std::pair<double, double> resid;
			unsigned fourthConfirm = t->getFourthHitResidual(t3, resid);

			if (fourthConfirm == 0 && t->chiSqProb(0) > cut && t->chiSqProb(1)
					> cut) {
				//fourth sensor was efficient
				//one could also cut on the residual value, if you wanted
				s4->registerTrackConfirm(t->fourthSensorThresh(), t->getHits().at(s4), t->timeStamp(), resid);
				eff++;
			} else {
				//fourth sensor wasn't efficien!
				//but only if it could have been
				if(!s4->isDeadArea(t->findXYPrediction(s4->zPosition()))) {
					s4->registerTrackMiss(t->fourthSensorThresh(), t->findXYPrediction(s4->zPosition()), t->timeStamp());
					ineff++;
				} else {
					++deadArs;
				}
			}

		}
	}
	
	//All that follows is housekeeping and plot extraction...
	
	char output[100];
	strcpy(output, argv[2]);
	TFile f(output, "recreate");

	//iterate over sensors to extract efficiency plots	
	std::vector<MapsSensor*>& sensors = mtm2.getSensors();
	for (std::vector<MapsSensor*>::iterator it = sensors.begin(); it
			!= sensors.end(); it++) {
		MapsSensor* s = *it;
		TH1F h;
		TH1F g, k;
		TH2F j, l, m, n;
		s->getEfficiencyCurve(h, 20, 0, 200);
		s->getEfficiencyTimestamps(g);
		s->getEfficiencyXYPlot(j);
		s->getInefficiencyXYPlot(m, true);
		s->getInefficiencyTimestamps(k);
		s->getFourthHitResidualPlot(l);
		s->getEfficiencyByGroup(n);
		h.Write();
		g.Write();
		j.Write();
		k.Write();
		l.Write();
		m.Write();
		n.Write();
		std::cout.precision(4);
		std::cout << "\n" << *s << ": Efficiency: \n";
		for (unsigned tMult(7); tMult < 20; tMult++) {
			double efficiency = s->getEfficiency(tMult*10);
			if (efficiency > -0.9)
				std::cout << "\t Thresh: " << tMult*10 << "\t: " << efficiency
						* 100 << "\n";
		}
	}

	//save plots to disk
	std::cout << "\tWriting root file: " << output << "\n";
	f.Write();
	f.Close();

	//close up shop
	std::cout << "ExtractEfficiencies: summary:\n";
	std::cout << "\tTotal candidate tracks:\t" << candidateTracks << "\n";
	std::cout << "\tEfficient hits:\t" << eff << "\n";
	std::cout << "\tInefficient hits:\t" << ineff << "\n";
	std::cout << "\tDead area intersections:\t" << deadArs << "\n";

	std::cout << "Have a nice day.\n";

}