#ifndef MAPSTRACK_HH_
#define MAPSTRACK_HH_
#include <map>
#include <iostream>
#include <vector>
#include "ToString.h"
#include "MapsSensor.hh"
#include "MapsException.hh"
#include "Operators.h"

/** MapsTrack.hh
 * 
 * Jamie Ballin, Imperial College London
 * 		February 2008
 * 
 * 
 * A track is defined by at least 3 hits in 3 sensors. Only one hit per sensor
 * is allowed. At most it is defined by 4 hits in 4 sensors (though only slight
 * revisions would be required to change this).
 * 
 * Note that MapsTracks are lazy: they do not tell each of the assosicated MapsSensors of
 * their residuals unless you explicitly tell them to do so. Telling sensors whether they
 * have been efficient or not is entirely up to the user; see examples in tests/, especially
 * ExtractEfficiencies.cpp code.
 * 
 */
class MapsTrack {
public:
	typedef std::pair<int, int> coord;
	typedef std::pair<double, double> physXY;

	MapsTrack() :
		myT(0), myFourthThresh(0) {
	}
	;

	/** Constructor */
	MapsTrack(unsigned bx, std::map<MapsSensor*, MapsTrack::coord> hits,
			MapsSensor* fourthSensor = 0, unsigned fourthThreshold = 0) :
		myT(bx), myFourthThresh(fourthThreshold),
				myFourthSensor(fourthSensor), myHits(hits), myTrackError(std::pair<double, double>(0.1, 0.1)), myPixelPitch(0.05) {
	}
	;

	virtual ~MapsTrack();

	/** the timestamp of this track */
	inline unsigned timeStamp() const {
		return myT;
	}
	;
	/** returns a map of hits making the track */
	const std::map<MapsSensor*, coord>& getHits();

	/** The standard deviation of the hits in x */
	double sigmaX() const;

	/** The standard deviation of the hits in y */
	double sigmaY() const;

	/** Returns the arithmetic mean of the x hits */
	double meanX() const;

	/** Returns the arithmetic mean of the y hits */
	double meanY() const;

	inline double pixelPitch() const {
		return myPixelPitch;
	}
	;

	inline std::pair<double, double> trackError() const {
		return myTrackError;
	}
	;

	/**
	 * Use this to set the size of the pixel to anything other than 0.05 mm.
	 */
	void setPixelPitch(const double& pitch);

	/**
	 * Use this to set the sigma parameter used in the chi squared computations.
	 */
	void setTrackError(const std::pair<double, double>& error);

	/** 
	 * Getter to tell you the fourth sensor for this track 
	 */
	MapsSensor* fourthSensor() const;

	/**
	 * Getter for the fourth sensor threshold.
	 */
	inline unsigned fourthSensorThresh() const {
		return myFourthThresh;
	}
	;
	/**
	 * Sets/overwrites the fourth sensor's threshold.
	 */
	void setFourthSensorThresh(const unsigned& threshold);

	/**
	 * A track is defined by at least 3 hits; this returns the first sensor 
	 * not found in the supplied list, and zero if all were found.
	 */
	MapsSensor* missingSensor(std::vector<MapsSensor*> sensors) const;

	/**
	 * Determines the residual i.t.o. the difference between the real hit in sensor s 
	 * and the track defined by its left and right extrapolation points
	 * 
	 * Throws an exception if the track doesn't have a hit for the supplied sensor.
	 */
	std::pair<double, double> simpleResidual(MapsSensor* s) const
			throw(MapsException);

	/**
	 * Convenience method for calling simpleResidual() for all sensors :-| 
	 */
	std::map<MapsSensor*, std::pair<double, double> > allResiduals();

	/**
	 * Sets each sensor's residuals based on this track.
	 */
	void tellSensorsOfResiduals();

	/**
	 * Sets sensor residuals conditionally: this allows you to define a fixed coordinate
	 * system. Usually the left most and right most sensors are set as the fixed references.
	 */
	void tellSensorsOfResiduals(const MapsSensor* const requiredLeft,
			const MapsSensor* const requiredRight);

	/**
	 * Sets each sensor's efficiency hits etc. based on this track.
	 */
	void tellSensorsOfHits();

	/** 
	 * Define the theta to be the angle in radians relative to a 
	 * straight (non-skew) line through the sensors. This uses the complete track fit.
	 */
	double theta() const;

	/**
	 * Uses the Root API to evaluate whether the observed chi2 for a 
	 * "correct model" is more than the chi2 of the track
	 */
	double chiSqProb(const unsigned& dimension) const;

	/**
	 * Overwrites/sets this track's hits with those supplied.
	 */
	void setHits(std::map<MapsSensor*, coord> hits);

	/**
	 * One line description of the track.
	 */
	friend std::ostream& operator<<(std::ostream& s, const MapsTrack& mt);

	/**
	 * Prints a complete description of this track to the specfified output stream
	 */
	friend std::ostream& diagnose(std::ostream& s, const MapsTrack& mt);

	/**
	 * Prints the hits sensor by sensor.
	 */
	std::ostream& printCoords(std::ostream& s) const;
	

	/**
	 * Returns the sensor that should be used as the left sensor in a 
	 * track extrapolation to the supplied sensor
	 */
	MapsSensor* leftExtrapPoint(MapsSensor*) const;

	/**
	 * Returns the sensor that should be used as the right sensor in a 
	 * track extrapolation to the supplied sensor
	 */
	MapsSensor* rightExtrapPoint(MapsSensor*) const;

	/**
	 * Determines the 'p' coefficients for this track, as defined by
	 *      chi^2 = sum_i [x_i - (p_0 +z_i*p_1)]^2/sigma_i^2
	 * where i is the layer index and the assumption
	 *      sigma_i -> sigma
	 * is applied.
	 * 
	 * Dimension: 	= 0 for x
	 * 				= 1 for y
	 */
	std::pair<double, double> fitParameters(const unsigned& dimension) const;

	/**
	 * Evaluates the chi squared based on the full track fit.
	 * 
	 * Dimension: 	= 0 for x
	 * 				= 1 for y
	 */
	double chiSq(const unsigned& dimension) const;

	/**
	 * Makes a three hit track from tracks with four hits, ignoring the ''fourth sensor''.
	 * Makes a copy of this track if we only have 3 hits anyway.
	 */
	void make3HitTrack(MapsTrack& mt);

	/**
	 * Evaluates the difference in pixel coordinates between the extrapolated three hit track 
	 * and the fourth sensor's hit. You supply the answer :-) Returns something other than zero 
	 * if there was no fourth hit for comparison.
	 */
	unsigned getFourthHitResidual(const MapsTrack& threeHitTrack,
			std::pair<double, double>& answer);

	/**
	 * Uses the supplied track's fit to extrapolate to the supplied z position.
	 */
	std::pair<double, double> findXYPrediction(const double& zpos) const;
	
	const std::map<MapsSensor*, physXY >& getGlobalHits() const;
	
	void eraseGlobalHits();

protected:

	MapsTrack& operator=(MapsTrack& mt) {
		return *this;
	}
	

private:

	//returns the sensor with a hit that is more "left" than any other
	MapsSensor* leftSensor() const;
	//and more "right" than any other
	MapsSensor* rightSensor() const;

	MapsTrack(MapsTrack& mt);
	
	unsigned myT;
	unsigned myFourthThresh;
	MapsSensor* myFourthSensor;
	std::map<MapsSensor*, coord> myHits;
	
	//this is the same as myHits, but with each coord converted to the global system
	//access it with getPhysicalHits()
	mutable std::map<MapsSensor*, physXY > myGlobalHits;

	std::pair<double, double> myTrackError;
	double myPixelPitch;

	//we don't want to tell the sensors twice of these things despite what application code may do
	bool myToldSensorsResids;
	bool myToldSensorsHits;

};

//std::ostream& operator<<(std::ostream& s, const MapsTrack::coord& c);

//std::ostream& operator<<(std::ostream& s, const std::pair<double, double>& p);

std::ostream& operator<<(std::ostream& s, const MapsTrack& mt);

std::ostream& diagnose(std::ostream& s, const MapsTrack& mt);

//#include "MapsTrack.cc"
#endif /**MAPSTRACK_HH_*/