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

/** ExtractEfficiencies.cpp
 * 
 * Workhorse for determining sensor efficiency once tracks have been processed
 * and you're happy with the alignment.
 * 
 * Exports plots to supplied outpufile name.
 * 
 * This is a good place to start to understand the general code structure, and
 * as a general example for how you could write your own tool.
 * 
 * Usage: ExtractEfficiencies input file outputfile
 * 
 * Jamie Ballin, Imperial College London
 * 		February 2008
 * 
 */

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.018, 0.018);
	//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();
	std::vector<MapsSensor*>& sensors = mtm2.getSensors();
	MapsSensor* s1;
	MapsSensor* s2;
	for (std::vector<MapsSensor*>::const_iterator it = sensors.begin(); it
			!= sensors.end(); ++it) {
		if ((*it)->id() == 8) {
			s1 = *it;
		} else if ((*it)->id() == 6) {
			s2 = *it;
		}
	}

	//count efficient and inefficient frequency
	unsigned ineff(0), eff(0), deadArs(0);
	//counts specific to shapers and samplers
	unsigned shapEff(0), sampEff(0);
	unsigned shapIneff(0), sampIneff(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;
		//construct the three hit track
		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++;
				if (t->getHits().at(s4).first < 84)
					shapEff++;
				else
					sampEff++;
			} else {
				//fourth sensor wasn't efficient!
				//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++;
					//Where would the hit have been i.t.o. a pixel coordinate?
					MapsSensor::coord predHit(0, 0);
					s4->convertPhysicalToHit(
							t->findXYPrediction(s4->zPosition()), predHit);
					if (predHit.first < 84)
						shapIneff++;
					else
						sampIneff++;
				} else {
					//hit would have intersected with a dead area
					++deadArs;
				}
			}

		}
	}

	//All that follows is housekeeping and plot extraction...

	char output[100];
	strcpy(output, argv[2]);
	TFile f(output, "recreate");

	//Find the projected dead area defined by the 4 sensors
	TH2F p;
	MapsSensor::mutualSensorMask(sensors, p);
	
	f.mkdir("effCurves");
	
	//iterate over sensors to extract efficiency plots
	for (std::vector<MapsSensor*>::iterator it = sensors.begin(); it
			!= sensors.end(); it++) {
		MapsSensor* s = *it;
		TH1F h, r, q;
		TH1F g, k;
		TH2F j, l, m, n;
		f.cd("effCurves");
		//Get global sensor efficiency
		s->getEfficiencyCurve(h, 13, 80, 210, true, true);
		h.Write();
		//Just shapers please
		s->getEfficiencyCurve(r, 13, 80, 210, true, false);
		r.Write();
		//Now just samplers
		s->getEfficiencyCurve(q, 13, 80, 210, false, true);
		q.Write();
		f.cd("/");

		s->getEfficiencyTimestamps(g);
		s->getEfficiencyXYPlot(j);
		s->getInefficiencyXYPlot(m, true);
		s->getInefficiencyTimestamps(k);
		s->getFourthHitResidualPlot(l);
		s->getEfficiencyByGroup(n);

		g.Write();
		j.Write();
		k.Write();
		l.Write();
		m.Write();
		n.Write();
		p.Write();
		std::cout.width(4);
		std::cout << "\n" << *s << ": Efficiency: \n";
		//		for (unsigned tMult(7); tMult < 20; tMult++) {
		//			double efficiency = s->getEfficiency(tMult*10);
		//			double shapEff = s->getEfficiency(tMult*10, true, false);
		//			double sampEff = s->getEfficiency(tMult*10, false, true);
		//			if (efficiency > -0.9) {
		//				std::cout.precision(3);
		//
		//				std::cout << "\t Thresh: " << tMult*10 << "\t: " << efficiency
		//						* 100 << ",\t shap:\t" << shapEff * 100
		//						<< ",\t samp:\t" << sampEff * 100 << "\n";
		//			}
		//		}
		
		//Commented code above now obseleted by...
		s->printEfficiencies(std::cout);
	}

	//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 << "\t\tShaper efficient hits:\t" << shapEff << "\n";
	std::cout << "\t\tShaper inefficient hits:\t" << shapIneff << "\n";
	std::cout << "\t\tSampler efficient hits:\t" << sampEff << "\n";
	std::cout << "\t\tSampler inefficient hits:\t" << sampIneff << "\n";
	std::cout << "\tDead area intersections:\t" << deadArs << "\n";
	std::cout << "Global efficiency:\t" << 100 * static_cast<double>(eff)
			/ (eff+ineff) << "\n";
	std::cout << "\tShapers:\t" << 100 * static_cast<double>(shapEff)
			/ (shapEff+shapIneff) << "\n";
	std::cout << "\tSamplers:\t" << 100 * static_cast<double>(sampEff)
			/ (sampEff+sampIneff) << "\n";

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

}