#ifndef OAROOTRACKER_IROOTRACKERHISTOGRAM_HXX
#define OAROOTRACKER_IROOTRACKERHISTOGRAM_HXX

#include <map>
#include <set>
#include <string>

#include <IDatum.hxx>

#include <TLorentzVector.h>

class TFile;
class THnBase;
class TH1;

namespace COMET{
    class IRooTrackerHistogram;
    class IOutputRooTrackerHistogram;

    OA_EXCEPTION(EIRooTrackerHistogram,COMET::EoaCore);
    OA_EXCEPTION(ECannotGeneratePrimary,EIRooTrackerHistogram);
    OA_EXCEPTION(ENegativeBinContent,EIRooTrackerHistogram);
    OA_EXCEPTION(ECannotAddHistograms,EIRooTrackerHistogram);
    OA_EXCEPTION(ECannotDrawHistogram,EIRooTrackerHistogram);
}
 
class COMET::IRooTrackerHistogram: public IDatum{
    public:
        struct PrimaryParticle_t{
            int pid;
            double time, posX, posY, posZ;
            double energy, momX, momY, momZ;
            PrimaryParticle_t():
                pid(0),
                time(-9999), posX(-9999), posY(-9999), posZ(-9999),
                energy(-9999),momX(-9999), momY(-9999), momZ(-9999){
            }
            PrimaryParticle_t(int _pid):
                pid(_pid),
                time(-9999), posX(-9999), posY(-9999), posZ(-9999),
                energy(-9999),momX(-9999), momY(-9999), momZ(-9999){
            }
            PrimaryParticle_t(int _pid,double _time,double _x, double _y, double _z,
                    double _energy,double _px,double _py,double _pz):
                pid(_pid),
                time(_time), posX(_x), posY(_y), posZ(_z),
                energy(_energy),momX(_px), momY(_py), momZ(_pz){
            }
            void Print(const std::string& prefix="")const;
        };
        typedef std::vector<PrimaryParticle_t> PrimaryList;

        enum Dimensions_t{
            kPosX=0, kPosY, kPosZ,
            kTime,
            kMomX, kMomY, kMomZ,
            kEnergy,
            kNumDimensions
        };

        // SpecificTrack stores sets of requirement to a track
        struct SpecificTrack_t{
            static const char* GetParameterNames() {
                const char* names = "count time posX posY posZ energy momX momY momZ";
                return names;
            }
            struct Range_t {
                double min, max;
                Range_t(): min(0), max(-1) {}
                Range_t(double _min, double _max): min(_min), max(_max) {}
                void Set(const double _min, const double _max){ min = _min; max = _max; } 
                inline bool IsSet() const { return min < max; }
                inline bool IsIn(const double val) const { return ((IsSet())? ( min < val && val < max ) : true); }
            };
            int pid;
            Range_t count;
            Range_t posX, posY, posZ, momX, momY, momZ; // For vector components
            Range_t energy, time;                       // For scalar values
            SpecificTrack_t(int _pid = 0) : pid(_pid) {}
            inline bool IsSpecificPID() const {return (pid != 0);}
            inline bool CheckNTracks(int _count) const { 
                return count.IsSet()? (count.min <= (int)_count && _count <= (int)count.max) : true; }
            inline bool CheckTrackProperties(double* properties) const {
                if(time.IsIn(properties[kTime]) && energy.IsIn(properties[kEnergy]) &&
                        posX.IsIn(properties[kPosX]) && posY.IsIn(properties[kPosY]) && posZ.IsIn(properties[kPosZ]) &&
                        momX.IsIn(properties[kMomX]) && momY.IsIn(properties[kMomY]) && momZ.IsIn(properties[kMomZ]) )
                    return true;
                return false;
            }
            bool SetRangeWithStr(std::string name, double min, double max){
                if     (name == "posX")   posX.  Set(min, max);
                else if(name == "posY")   posY.  Set(min, max);
                else if(name == "posZ")   posZ.  Set(min, max);
                else if(name == "momX")   momX.  Set(min, max);
                else if(name == "momY")   momY.  Set(min, max);
                else if(name == "momZ")   momZ.  Set(min, max);
                else if(name == "energy") energy.Set(min, max);
                else if(name == "time")   time.  Set(min, max);
                else if(name == "count")  count. Set(min, max);
                else return false;
                return true;
            }
        }; // End of SpecificTrack_t
        typedef std::vector<SpecificTrack_t> SpecificTracks;

        typedef int PDG_t;
        typedef std::map<PDG_t,THnBase*> ParticleDistributionHistograms;
        typedef std::map<PDG_t,TH1*> ParticleFrequencyHistograms;
        typedef std::map<PDG_t,int> ParticleTally;

    public:

        IRooTrackerHistogram():IDatum(GetNameOnFile().c_str(),"Histogrammed version of RooTracker data"){
        }
        virtual ~IRooTrackerHistogram(){}

        static IRooTrackerHistogram* Open(const std::string& filename);

        static std::string DimensionStr(Dimensions_t dim);
        static Dimensions_t DimensionStr(const std::string& dim);

        const PrimaryList& GeneratePrimaries();
        PrimaryParticle_t GenerateOnePrimary();

        inline const PrimaryList& GetPrimaryList() const {return fGeneratedPrimaries;}
        inline Int_t GetNGeneratedPrimaries(){return fGeneratedPrimaries.size(); }
  
        void SetSpecificTracks(const SpecificTracks& tracks) { fSpecificTracks = tracks; }
        void ClearSpecificTracks() { fSpecificTracks.clear(); }
        inline SpecificTracks& GetSpecificTracks() {return fSpecificTracks; }
        inline Int_t GetNSpecificTracks(){ return fSpecificTracks.size(); }
        const PrimaryList& GenerateSpecificPrimaries();
  
        std::set<int> GetDescribedParticles()const;
        const TH1* GetFrequencyHist(int pid)const;
        const THnBase* GetDistributionHist(int pid)const;
        void PrintParticles()const;
        /// Get or make a new directory under the optional parent directory
        TDirectory* GetElseMakeDirectory(const std::string& name, 
                                         TDirectory* parentDir=NULL);
        /// Draw projections of distribution histograms and frequency histograms 
        /// to a file.  By default, it writes to the current file, but another 
        /// file can be specified.
        void Draw(TFile* aFile, bool only1DHistos=false);

        /// Unhide the other, virtual Draw from TObject.
        using TObject::Draw;

        static std::string GetNameOnFile(){return "RooTrackHist";}
        template <typename iterator>
        void BlackListParticles(iterator begin, iterator end){
            for(iterator it=begin; it!=end; ++it){
                fParticleDistributions.erase(*it);
                fParticleFrequencies.erase(*it);
            }
        }
        template <typename iterator>
        void WhiteListParticles(const iterator& begin, const iterator& end){
            ParticleFrequencyHistograms copy=fParticleFrequencies;
            for(ParticleFrequencyHistograms::const_iterator i_particle=copy.begin(); i_particle!=copy.end(); ++i_particle){
                iterator found=std::find(begin,end,i_particle->first);
                if(found!=end){
                    fParticleDistributions.erase(i_particle->first);
                    fParticleFrequencies.erase(i_particle->first);
                }
            }
        }
    private:
        void InitializeParticleWeights();
        void DrawByPID(int pid, bool only1DHistos);

    protected:
        ParticleDistributionHistograms fParticleDistributions;
        ParticleFrequencyHistograms fParticleFrequencies;
        ParticleTally fParticlesInThisEvent;  // ought to be in IOutputRooTrackerHistogram

        PrimaryList fGeneratedPrimaries; // ! 
        typedef std::map<double,int> ParticleWeights;
        ParticleWeights fParticleWeights;

        SpecificTracks fSpecificTracks; // !

        ClassDef(IRooTrackerHistogram,1);
};

class COMET::IOutputRooTrackerHistogram:public IRooTrackerHistogram{
    private:
        struct AxisProperties{
            int bins;
            double low, high;
            std::string title;
            AxisProperties():
                bins(-1),low(-1),high(-1),title("blah"){}
            AxisProperties(int _b,double _l,double _h, const char* _title):
                bins(_b),low(_l),high(_h),title(_title){}
        };
        typedef std::map<Dimensions_t,AxisProperties> AxisList;

    public:
        IOutputRooTrackerHistogram(const std::string& filename):
            fFile(NULL),
            fEventsIncluded(0),
            fTotalEvents(0),
            fFilename(filename) {
                BuildDefaultAxes();
            }

        bool SetupAxis(const std::string& axis_specification);

        bool InitialiseForWriting(const char* mode="CREATE");

        void SetOutputFileName(const std::string& filename){
            fFilename=filename;
        }
        std::string GetOutputFileName() const{
            return fFilename;
        }
        /// Fills in the frequency histograms with events where particles were 
        /// not seen by comparing the integral of the frequency histogram to the 
        /// total number of events read
        void SetFreqTotalEvents(int totalEvents);
        /// Sets the POT equivalent to the number of events read
        void Finalize(const int eventsRead);
        virtual void Close();
        /// Write the underlying distribution histograms (THnSparse) and 
        /// frequency histograms (TH1) to a new file
        void WriteHistogramsToFile(TFile* aFile);
        int GetNumEventsIncluded()const{return fEventsIncluded;}
        int GetNumTotalEvents()const{return fTotalEvents;}

        bool Add(const IRooTrackerHistogram& hist, 
                 double scale=1.0, 
                 bool skipCheck=false);
        /// Checks if hB can be added to hA via hA.Add(hB)
        /// Check that all the binning for the frequency and distribution
        /// histograms are consistent.  This runs the CanAdd(TH1*,TH1*) for the 
        /// frequency histograms and a projection of each dimension of the 
        /// THnSparse histogram
        static bool CanAdd(const IRooTrackerHistogram& hA,
                           const IRooTrackerHistogram& hB);
        /// Checks that the binning for a 1D histogram are consistent.
        /// This means the same number of bins, the same axis limits, and the 
        /// same bin limits.
        /// NOTE: Copies ROOT's TH1::CheckConsistency, only simiplified for the 
        /// one dimensional case
        static bool CanAdd(const TH1* hA, const TH1* hB);

    public:

        void BeginEvent(){fParticlesInThisEvent.clear();}
        void EndEvent();

        void AddParticle(int pid,
                   const TLorentzVector& position,
                   const TLorentzVector& momentum,
                   double weight=1);

    private:
        void BuildDefaultAxes();

        THnBase* GetElseMakeDistribution(int pid);
        TH1* GetElseMakeFrequency(int pid);

    private:
        TFile* fFile;
        int fEventsIncluded;
        int fTotalEvents;
        std::string fFilename;
        AxisList fAxisProperties;
};

#endif // OAROOTRACKER_IROOTRACKERHISTOGRAM_HXX