#include <memory>
#include <cmath>
#include <iostream>
#include <string>
#include <set>

#include <tut.h>

#include <TVector3.h>
#include <HEPUnits.hxx>
#include <IGeomInfo.hxx>
#include <IHitSelection.hxx>
#include <IHit.hxx>
#include <IMCHit.hxx>
#include <IG4HitSegment.hxx>

#include "tutGenerateEvents.hxx"

namespace tut {
    struct baseP0DHits {
        baseP0DHits() {
            // Run before each test.
            SimG4::GenerateEvents();
        }
        ~baseP0DHits() {
            // Run after each test.
        }
    };

    // Declare the test
    typedef test_group<baseP0DHits>::object testP0DHits;
    test_group<baseP0DHits> groupP0DHits("P0DHits");

    // Check that the SimG4 events with SimDetectorResponse run exist.
    template<> template<>
    void testP0DHits::test<1> () {
        ensure("COMET MC Events exist.", SimG4::Events.size()>0);        
    }

    // Test that hits were generated.
    template<> template<>
    void testP0DHits::test<2> () {
        int p0dCount = 0;
        for (SimG4::EventVector::iterator e = SimG4::Events.begin();
             e != SimG4::Events.end(); 
             ++e) {
            COMET::IHandle<COMET::IHitSelection> p0d = (*e)->GetHitSelection("p0d");
            if (p0d && p0d->size()>0) ++p0dCount;
        }
        ensure_greaterthan("P0D hits exist",p0dCount,0);
    }

    // Check that all of the P0D hits have reasonable values.
    template<> template<>
    void testP0DHits::test<3>() {
        int p0dEvtCount = 0;
        for (SimG4::EventVector::iterator e = SimG4::Events.begin();
             e != SimG4::Events.end(); ++e) {
            COMET::IHandle<COMET::IHitSelection> p0d((*e)->GetHitSelection("p0d"));
            if (!p0d) continue;
            ++p0dEvtCount;
            int zeroTimes = 0;
            int totalTimes = 0;
            for (COMET::IHitSelection::iterator h = p0d->begin();
                 h != p0d->end(); ++h) {
                double charge = (*h)->GetCharge();
                double time = (*h)->GetTime();
                ensure_greaterthan("P0D hit charge positive",charge,0);
                ensure_lessthan("P0D hit charge valid",charge,667);
                ensure_lessthan("P0D hit time valid",time,1*unit::second);
                ensure_greaterthan("P0D hit time valid",time,-1*unit::second);
                if (std::abs(time) < 0.1*unit::nanosecond) ++zeroTimes;
                ++totalTimes;
            }
            ensure_lessthan("Few P0D zero times",zeroTimes,1+totalTimes/1000);
        }
        ensure_greaterthan("Events with P0D hits", p0dEvtCount, 0);
    }

    // Check that all of the P0D hits are in P0D bars.
    template<> template<>
    void testP0DHits::test<4>() {
        COMET::IGeomInfo& info = COMET::IGeomInfo::Get();
        for (SimG4::EventVector::iterator e = SimG4::Events.begin();
             e != SimG4::Events.end(); ++e) {
            COMET::IHandle<COMET::IHitSelection> p0d((*e)->GetHitSelection("p0d"));
            if (!p0d) continue;
            for (COMET::IHitSelection::iterator h = p0d->begin();
                 h != p0d->end(); ++h) {
                std::string name = info.NodeName((*h)->GetGeomId());
                ensure("Hit in P0D",
                       name.find("/P0D_") != std::string::npos
                       && name.find("/Bar_") != std::string::npos);
            }
        }
    }

    // Check that the PD hits are close to the deposited energy from tracks,
    // and that hits from noise exist.
    template<> template<>
    void testP0DHits::test<5>() {
        int withContributors = 0;
        int withoutContributors = 0;
        for (SimG4::EventVector::iterator e = SimG4::Events.begin();
             e != SimG4::Events.end(); ++e) {
            COMET::IHandle<COMET::IHitSelection> p0d
                = (*e)->GetHitSelection("p0d");
            if (!p0d) continue;
            
            for (COMET::IHitSelection::iterator h = p0d->begin();
                 h != p0d->end(); ++h) {
                COMET::IHandle<COMET::IMCHit> mcHit(*h);
                ensure("Hit is a COMET::IMCHit",mcHit);
                if (mcHit->GetContributors().size() < 1) {
                    ++withoutContributors;
                    continue;
                }
                ++withContributors;
                TVector3 pos(mcHit->GetPosition());
                for (COMET::IMCHit::ContributorContainer::const_iterator c
                         = mcHit->GetContributors().begin();
                     c != mcHit->GetContributors().end();
                     ++c) {
                    COMET::IG4HitSegment* seg = dynamic_cast<COMET::IG4HitSegment*>(*c);
                    ensure("COMET::IMCHit contributor is a COMET::IG4HitSegment*",seg);
                    if (mcHit->IsXHit()) {
                        double x = 0.5*(seg->GetStartX()+seg->GetStopX());
                        ensure_distance("Hit close in X",pos.X(),x,2*unit::cm);
                    }
                    if (mcHit->IsYHit()) {
                        double y = 0.5*(seg->GetStartY()+seg->GetStopY());
                        ensure_distance("Hit close in Y",pos.Y(),y,2*unit::cm);
                    }
                    if (mcHit->IsZHit()) {
                        double z = 0.5*(seg->GetStartZ()+seg->GetStopZ());
                        ensure_distance("Hit close in Z",pos.Z(),z,2*unit::cm);
                    }
                }
            }
        }
        ensure_greaterthan("P0D Hits from particle trajectories exist",
                           withContributors,0);
    }

    // Check that muons moving along z axis create hits in consecutive
    // scintillator planes.
    template<> template<>
    void testP0DHits::test<6>() {
        COMET::ICOMETEvent* event = SimG4::Events.front();
            COMET::IHandle<COMET::IHitSelection> p0d(event->GetHitSelection("p0d"));
        std::set<unsigned int> layerIndex;
        for (COMET::IHitSelection::iterator h = p0d->begin();
             h != p0d->end(); ++h) {
            const TVector3& pos = (*h)->GetPosition();
            unsigned int index = COMET::IGeomInfo::P0D().ActivePlane(pos.Z());
            layerIndex.insert(index);
        }
        
        int gaps = 0;
        for (std::set<unsigned int>::iterator iter = layerIndex.begin();
             iter != layerIndex.end(); ++iter) {
            std::set<unsigned int>::iterator next = iter;
            ++next;
            if (next == layerIndex.end()) break;
            unsigned int delta = *next - *iter;
            if (delta>1) ++gaps;
            ensure_lessthan("Hits in consecutive layers", delta, 3);
        }
        ensure_lessthan("Few gaps in track",gaps,2);
    }

};