/** @file check.cpp functions for checking track reconstruction quality */

#include <algorithm>
#include <map>
#include <math.h>
#include <set>
#include <stdio.h>
#include <utility>

#include "check.h"

void copy_tracks(vector<Track> &tracks, map<Track, int> &track_counts) {
	typedef map<Track, int>::iterator TrackIter;
	for(TrackIter it = track_counts.begin(); it != track_counts.end(); it++) {
		const Track& t = it->first;
		int cur_count = it->second;
		tracks.push_back(Track(t.x0, t.y0, t.r, cur_count, t.evtid));
	}	
}  // copy_tracks

void extract_tracks
(vector<vector<Track> > &tracks, const vector<EventPoint> &eps) {
	tracks.clear();
	map<Track, int> track_counts;
	int prev_ievt = -1;
	for(int i = 0; i < eps.size(); i++) {
		EventPoint ep = eps[i];
		int ievt = ep.evtid;
		if(ievt != prev_ievt) {
			if(prev_ievt >= 0) {
				tracks.push_back(vector<Track>());
				copy_tracks(tracks[prev_ievt], track_counts);
				track_counts.clear();
			}
			prev_ievt = ievt;
		}
		float x0 = ep.z, r = ep.r, y0 = sqrt(r * r - x0 * x0);
		Track track = Track(x0, y0, r, 0, ievt);
		track_counts[track]++;
	}  // for each event point
	// copy the tracks for the last event
	if(!track_counts.empty()) {
		tracks.push_back(vector<Track>());
		copy_tracks(*tracks.rbegin(), track_counts);
	}
}  // extract_tracks

void sort_by_event
(vector<vector<Track> > &evt_tracks, const vector<Track> &tracks) {
	map<int, vector<Track> > evt_m;
	for(vector<Track>::const_iterator it = tracks.begin(); 
			it != tracks.end();	it++)
		evt_m[it->evtid].push_back(*it);

	// find the maximum event number
	int max_ievt = -1;
	for(map<int, vector<Track> >::iterator it = evt_m.begin(); it != evt_m.end();
			it++) {
		max_ievt = max(max_ievt, it->first);
	}
	int nevt = max_ievt + 1;
	
	// push the tracks to event's vectors
	evt_tracks.resize(nevt);
	for(map<int, vector<Track> >::iterator it = evt_m.begin(); it != evt_m.end();
			it++)
		evt_tracks[it->first] = it->second;
	
}  // sort_by_event

void compute_rates
(const vector< vector<Track> > &orig_tracks, 
const vector< vector<Track > > &rec_tracks,
float eps) {

	// print some tracks
	// for(int ievt = 10; ievt < 20; ievt++) {
	// 	printf("event %d, original tracks:\n", ievt);
	// 	for(int i = 0; i < orig_tracks[ievt].size(); i++) {
	// 		const Track &t = orig_tracks[ievt][i];
	// 		printf("(%lf %lf %lf), %d\n", t.x0, t.y0, t.r, t.npoints);
	// 	}
	// 	printf("event %d, reconstructed tracks:\n", ievt);
	// 	for(int i = 0; i < rec_tracks[ievt].size(); i++) {
	// 		const Track &t = rec_tracks[ievt][i];
	// 		printf("(%lf %lf %lf), %d\n", t.x0, t.y0, t.r, t.npoints);
	// 	}
	// }

	// count of tracks in the arrays
	int norig = 0, nrec = 0;
	// count for correct tracks, failed tracks, false positives
	int nsucc = 0, nfail = 0, nfalsepos = 0;
	for(int ievt = 0; ievt < orig_tracks.size(); ievt++) {
		const vector<Track> &orig_evt_tracks = orig_tracks[ievt];
		norig += orig_evt_tracks.size();
		if(ievt >= rec_tracks.size()) {
			// event has not been reconstructed at all
			nfail += orig_evt_tracks.size();
			continue;
		}
		const vector<Track> &rec_evt_tracks = rec_tracks[ievt];
		nrec += rec_evt_tracks.size();
		// event has been reconstructed, compare tracks
		// for each original track, find whether it has been reconstructed
		for(int iorig = 0; iorig < orig_evt_tracks.size(); iorig++) {
			bool reconstructed = false;
			for(int irec = 0; irec < rec_evt_tracks.size(); irec++) {
				if(equal_eps(orig_evt_tracks[iorig], rec_evt_tracks[irec], eps)) {
					reconstructed = true;
					break;
				}
			}  // for each reconstructed track
			if(reconstructed) 
				nsucc++;
			else
				nfail++;
		}  // for each original track

		// for each reconstructed track, find whether it corresponds to the original
		// track 
		for(int irec = 0; irec < rec_evt_tracks.size(); irec++) {
			bool corresponds = false;
			for(int iorig = 0; iorig < orig_evt_tracks.size(); iorig++) {
				if(equal_eps(orig_evt_tracks[iorig], rec_evt_tracks[irec], eps)) {
					corresponds = true;
					break;
				}		
			}
			if(!corresponds)
				nfalsepos++;
		}  // for each reconstructed track
	}  // for each event

	// count rates
	int ndup = nrec - nsucc - nfalsepos;
	double succ_rate = nsucc / (double)norig;
	double fail_rate = nfail / (double)norig;
	double falsepos_rate = nfalsepos / (double)norig;
	double dup_rate = ndup / (double)norig;
	printf("%d tracks original, %d tracks reconstructed\n", norig, nrec);
	printf("%d success, %d failed, %d false positives, %d duplicates\n", 
				 nsucc, nfail, nfalsepos, ndup);
	printf("success rate %lf, fail rate %lf, " 
				 "false positive rate %lf, duplicate rate %lf\n",
				 succ_rate, fail_rate, falsepos_rate, dup_rate);
}  // compute_rates