// LSR BFEM Noise Analysis   September 9, 2001
//

// Include files
// Gaudi system includes
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/Algorithm.h"


// ntuple interface
#include "ntupleWriterSvc/INTupleWriterSvc.h"

// if use the gui
#include "GuiSvc/IGuiSvc.h"
#include "gui/GuiMgr.h"
#include "gui/DisplayControl.h"

// TDS class declarations, raw and recon data
#include "digiRootData/DigiEvent.h"
#include "GlastEvent/Hits/SiLayers.h"
#include "GlastEvent/TopLevel/Event.h"
#include "CalRecon/CalADCLogs.h"

#include "GlastEvent/Digi/AcdDigi.h"

#include "TkrRecon/SiRecObjs.h"
#include "TkrRecon/SiClusters.h"
#include "TkrRecon/GFgamma.h"
#include "CalRecon/CalRecLogs.h"
#include "CalRecon/CsIClusters.h"

#include "xml/IFile.h"

#include <cassert>
#include <algorithm>
#include <map>
#include <string>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>

#include "GaudiKernel/RndmGenerators.h"

//For the Histograms

#include "GaudiKernel/IHistogram1D.h"
#include "Gaudikernel/IHistogramSvc.h"
static IHistogram1D  *h_summary, *h_pos, *h_distNear, *h_dist1Hit, *h_distNHit, *h_totHits,
                     *h_totClusts, *h_clustSize, *h_distRndm, *h_posRndm, *h_clustsSel,
                     *h_logChi;

//keep track of strips hit
static int hits[41600];
static int summary[100];



// Define the class here instead of in a header file: not needed anywhere but here!
//------------------------------------------------------------------------------
/** 
A simple algorithm.

  
*/
class UserAlg : public Algorithm {
public:
    UserAlg(const std::string& name, ISvcLocator* pSvcLocator);
    StatusCode initialize();
    StatusCode execute();
    StatusCode finalize();
    
private: 
    /// Read in the event list provided in the job options file
    void processEventList();
    //! number of times called
    int m_count; 
    //! constants from the "instrument.xml" file.
    // xml::IFile * m_ini; 
    //! access to the Gui Service for display of 3-d objects
    IGuiSvc*    m_guiSvc;
    //! access the ntupleWriter service to write out to ROOT ntuples
    INTupleWriterSvc *m_ntupleWriteSvc;
    //! parameter to store the logical name of the ROOT file to write to
    std::string m_tupleName;
    /// List of Event Ids to Pause on
    StringArrayProperty m_eventList;
    /// keeps track of run numbers and event ids specified in job options file
    std::map< int, std::vector<int> > m_idMap;
    /// Flag which denotes if the user provided run numbers in the eventList
    bool m_runNumFlag;
    ///hit strip list
    int m_processed;

    // for ntuple
    float m_chi1, m_chi2;
};
//------------------------------------------------------------------------

// necessary to define a Factory for this algorithm
// expect that the xxx_load.cxx file contains a call     
//     DLL_DECL_ALGORITHM( UserAlg );

static const AlgFactory<UserAlg>  Factory;
const IAlgFactory& UserAlgFactory = Factory;

//------------------------------------------------------------------------
//! ctor
UserAlg::UserAlg(const std::string& name, ISvcLocator* pSvcLocator)
:Algorithm(name, pSvcLocator)
,m_count(0), m_guiSvc(0), m_runNumFlag(true)
{
    // declare properties with setProperties calls
    declareProperty("tupleName",  m_tupleName="");
    /// List of strings, of the form:  {"runNum_eventId", "runNum_eventId"...}
    declareProperty("eventList", m_eventList);
    
}


//------------------------------------------------------------------------
//! set parameters and attach to various perhaps useful services.
StatusCode UserAlg::initialize(){
    StatusCode  sc = StatusCode::SUCCESS;
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "initialize" << endreq;
    
    // Use the Job options service to set the Algorithm's parameters
    setProperties();
    
    processEventList();
    
    if( m_tupleName.empty()) {log << MSG::ERROR << "tupleName property not set!"<<endreq;
    return StatusCode::FAILURE;}
    // now try to find the GlastDevSvc service
    
    // get the Gui service (not required)
    if (service("GuiSvc", m_guiSvc).isFailure ()){
        log << MSG::WARNING << "No GuiSvc, so no display" << endreq;
    }
    
    // get a pointer to our ntupleWriterSvc
    if (service("ntupleWriterSvc", m_ntupleWriteSvc).isFailure()) {
        log << MSG::ERROR << "writeJunkAlg failed to get the ntupleWriterSvc" << endreq;
        return StatusCode::FAILURE;
    }

    // define histograms
    h_summary  = histoSvc()->book("/stat","Summary", "Total Hits/Strip",
        100, -.5, 99.5);
    h_pos  = histoSvc()->book("/stat","HitPosition", "Transverse Position of Hit",
        100, -16.,16.);
    h_distNear  = histoSvc()->book("/stat","NearDist", "Distance clusters from track",
        100, 0, 1.);
    h_dist1Hit  = histoSvc()->book("/stat","Dist1Hit", "Distance of 1-hit clusters",
        100, 0, 50.);
    h_distNHit  = histoSvc()->book("/stat","DistNHit", "Distance of multihit clusters",
        100, 0, 50.);
    h_distRndm  = histoSvc()->book("/stat","DistRndm", "Distance to a random Hit",
        100, 0, 50.);
    h_posRndm  = histoSvc()->book("/stat","PosRndm", "Position of random Hit",
        100, -16.,16.);
    h_totHits  = histoSvc()->book("/stat","TotHits", "Total # Hits/Event",
        100, -.5, 99.5);
    h_clustsSel  = histoSvc()->book("/stat","SelClusts", "Clusters/Event for all selected events",
        100, -.5, 99.5);
    h_totClusts  = histoSvc()->book("/stat","ToTClusts", "Total # Clusters/Event",
        100, -.5, 99.5);
    h_clustSize  = histoSvc()->book("/stat","ClustSize", "Cluster Size",
        100, -.5, 99.5);
    h_logChi = histoSvc()->book("/stat", "LogChiSq", "Log of Max Chisquared", 49, -5., 2.);

    return sc;
}

//------------------------------------------------------------------------
//! process an event
StatusCode UserAlg::execute()
{
    StatusCode  sc = StatusCode::SUCCESS;
    MsgStream   log( msgSvc(), name() );
    log << MSG::INFO << "executing " << ++m_count << " time" << endreq;
    
    // Here we are adding a new column to our ROOT ntuple
    sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "NumCalls", m_count);
    sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "chi1", m_chi1);
    sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "chi2", m_chi2);
    
    // Check the event Id for this event, if this id and run number combination
    // is in our eventList, we will pause the GUI.
    SmartDataPtr<Event> event(eventSvc(), "/Event/Header");
    unsigned int eventId = event->event();
    unsigned int runNum = 0;
    if (m_runNumFlag) runNum = event->run();
    std::vector<int> &eventList = m_idMap[runNum];
    // Search to see if this event id is among the list of ids we want to pause on
    int *loc = std::find(eventList.begin(), eventList.end(), eventId);
    // If we found it, we set the GUI to PAUSE
    if (loc != eventList.end()) {
        m_guiSvc->guiMgr()->pause();
    }

    SmartDataPtr<SiRecObjs> tkrReconData(eventSvc(), "/Event/TkrRecon/SiRecObjs");
    if (tkrReconData != 0) {

        // find "gammas" with 2 and only 2 tracks
        int numGammas = tkrReconData->numGammas();
        if (numGammas!=1) return sc;
        int numParticles = tkrReconData->numParticles();
        if (numParticles) return sc;
        GFgamma* gamma = tkrReconData->Gamma(0);
        if (gamma->empty()) return sc;

        bool only_best =  !gamma->getBest(SiCluster::X)->empty() &
                          !gamma->getBest(SiCluster::Y)->empty() &
                           gamma->getPair(SiCluster::X)->empty() &
                           gamma->getPair(SiCluster::Y)->empty() ;

        // if you want 2-particle gammas
        //
        //bool four_track = !gamma->getBest(SiCluster::X)->empty() &
        //                  !gamma->getBest(SiCluster::Y)->empty() &
        //                  !gamma->getPair(SiCluster::X)->empty() &
        //                  !gamma->getPair(SiCluster::Y)->empty() ;

        m_chi1 = -9999.;
        m_chi2 = -9999.;

        if(!only_best) return sc;

        // check for track chiSquared and # planes
        int nplanes_min = 1000;
        float chiSq_max = -1.;
        float chiSq[2];

        for (int view = 0; view<2; view++) {
            GFtrack* xbest = gamma->getBest((SiCluster::view) view);
            int nplanes = xbest->kplanelist.size();
            if (nplanes<nplanes_min) nplanes_min = nplanes;
            chiSq[view] = xbest->chiSquare();
            if (chiSq[view]>chiSq_max) chiSq_max = chiSq[view];
        }
        m_chi1 = chiSq[0];
        m_chi2 = chiSq[1];
        

        h_logChi->fill(log10(chiSq_max));
        if (chiSq_max>0.01) return sc;
        if (nplanes_min<5)  return sc;

        std::cout << "chiSq_max = " << chiSq_max << std::endl;
        //m_guiSvc->guiMgr()->pause();

        m_processed++;

        // get parameters of line between first and last hits
        float x1[2], x2[2], z1[2], z2[2], dxdz[2];

        for (view = 0; view<2; view++) {
            GFtrack* xbest = gamma->getBest((SiCluster::view) view);               
            int nplanes = xbest->kplanelist.size();               
            float chiSq = xbest->chiSquare();
            //std::cout << "nplanes " << nplanes << " chisq " << chiSq << std::endl;
            for (int iplane = 0; iplane<nplanes;iplane++) {
                KalPlane kplane = xbest->kplanelist[iplane];
                Point hitPoint = kplane.getPoint(KalHit::MEAS);
                //std::cout <<  "     " << hitPoint << std::endl;

                if (iplane==0) {
                    x1[view] = hitPoint.x();
                    z1[view] = hitPoint.z();
                } else if (iplane==nplanes-1) {
                    x2[view] = hitPoint.x();
                    z2[view] = hitPoint.z();
                }

            }
            dxdz[view] = (x2[view]-x1[view])/(z2[view]-z1[view]);
        }
        // go through the  clusters and throw out the ones 
        //    within a certain distance of the track

        SmartDataPtr<SiClusters> clusterData(eventSvc(), "/Event/TkrRecon/SiClusters");

        int totalHits = 0;
        int totalClusters = 0;
        int maxNStrips = -1;
        float maxDist = -1.;
        int nclust = clusterData->nHits();
        h_clustsSel->fill(nclust, 1.);

        // count the hits and clusters first
        for (int ihit = 0; ihit<nclust; ihit++){
            SiCluster* cluster = clusterData->getHit(ihit);
            Point hitPoint = cluster->position();
            int view = cluster->v();
            int plane = cluster->plane();
            float transverse; 
            if (view==0) {transverse = hitPoint.x();}
            else {transverse = hitPoint.y();}
            float projected_pt = dxdz[view]*(hitPoint.z()-z1[view]) + x1[view];
            float distance = projected_pt - transverse;
            if (fabs(distance)> maxDist) maxDist = fabs(distance);

            Rndm::Numbers uniform(randSvc(), Rndm::Flat(-16.,16.));
            float position = uniform();
            h_distRndm->fill(fabs(position-projected_pt),1.);
            h_posRndm->fill(position,1.);
            if (fabs(distance)>15.0) {
                totalClusters++;
                int strip0 = cluster->firstStrip();
                int stripf = cluster->lastStrip();
                int nstrips = stripf-strip0+1;
                if (nstrips>maxNStrips) maxNStrips = nstrips;
                totalHits += nstrips;
            }
        }

        //make cuts on total #clusters, #hits if you want
        //if (totalClusters>3) return sc;

        // now score the survivors

        totalHits = 0;
        totalClusters = 0;

        for (ihit = 0; ihit<nclust; ihit++){
            SiCluster* cluster = clusterData->getHit(ihit);
            Point hitPoint = cluster->position();
            int view = cluster->v();
            int plane = cluster->plane();
            float transverse;
            if (view==0) {transverse = hitPoint.x();}
            else {transverse = hitPoint.y();}
            float distance = dxdz[view]*(hitPoint.z()-z1[view]) + x1[view] - transverse;
            h_pos->fill(transverse, 1.);
            h_distNear->fill(fabs(distance),1.);

            int strip0 = cluster->firstStrip();
            int stripf = cluster->lastStrip();
            int nstrips = stripf-strip0+1;

            if (nstrips==1) h_dist1Hit->fill(fabs(distance),1.);
            if (nstrips>1)  h_distNHit->fill(fabs(distance),1.);
            if ((fabs(distance)>15.0)) {
                totalClusters++;
            // add to the cluster list and write it out
                std::cout << " Run" << runNum << " ev " << eventId << " cluster " << cluster->plane() << " " << view 
                    << " " << cluster->strip() << " dist " << distance 
                    << " nclusts " << cluster->lastStrip()-cluster->firstStrip()+1 << std::endl;
                int base = 1600*(view + 2*(12-plane));
                if (nstrips>maxNStrips) maxNStrips = nstrips;
                h_clustSize->fill(nstrips,1.);
                if (nstrips>0) {
                    for (int istrip= strip0; istrip<stripf+1; istrip++) {
                        hits[istrip+base]++;
                        totalHits++;
                    }
                }
            }
        }

        h_totHits->fill(totalHits,1.);
        h_totClusts->fill(totalClusters,1.);
   
        // can use the lines below to pause the display on some condition
        //if (totalHits>5) {
        //    std::cout << "total hits" << totalHits << std::endl;
        //    m_guiSvc->guiMgr()->pause();
        //}
    }

    // Here's how to retrieve TKR raw data from the TDS
/*
    SmartDataPtr<SiLayers> tkrRawLayerData(eventSvc(),"/Event/TkrRecon/SiLayers");
    
    if (tkrRawLayerData != 0) {
        int nlayers = tkrRawLayerData->num();
        
        for (int ilayer = 0; ilayer < nlayers; ilayer++) {
            SiLayer* layer  = tkrRawLayerData->Layer(ilayer);
            
            int      klayer = layer->layer();
            int      iview  = layer->view();
            int      nHits  = layer->nstrips();
             for (int i = 0; i<nHits; i++) {
                int strip = layer->idstrip();
            }               
        }
    }
*/
    
    
    SmartDataPtr<CalADCLogs> calRawData(eventSvc(), "/Event/CalRecon/CalADCLogs");
    if (calRawData != 0) {
        
    }
    
    SmartDataPtr<CalRecLogs> calReconData(eventSvc(), "/Event/CalRecon/CalRecLogs");
    if (calReconData != 0) {
        
    }
    
    SmartDataPtr<AcdDigiVector> acdRawData(eventSvc(), "/Event/Digi/AcdDigis");
    if (acdRawData != 0) {
        
    }
    
    
    return sc;
}
/

//------------------------------------------------------------------------
//! clean up, summarize
StatusCode UserAlg::finalize(){
    StatusCode  sc = StatusCode::SUCCESS;
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "finalize after " << m_count << " calls." << endreq;

    //Write out some summary info
    std::ofstream res("results.txt");

/
    res << m_count << " events, " << m_processed << " analyzed." << std::endl;

    for (int i = 0; i<41600; i++) {
        int nhits = hits[i];
        if (nhits>99) nhits = 99;
        summary[nhits]++;
    }

    res << "Summary of hits:" << std::endl;
    for(i=0;i<100;i++) {
        if (summary[i]) {
            res << i << " " << summary[i] << std::endl;
        }
    }

    //This one doesn't work, for some reason... too late maybe?
    res << "Filling: ";
    for (i=0;i<100;i++) {
        h_summary->fill(i*1.,summary[i]*1.);
    }
    res << std::endl;

    //List of potentially hot strips
    res << "hits" << std::endl;
    for (i = 0; i<41600; i++) {
        int nhits = hits[i];
        if (nhits>99) nhits = 99;
        if (nhits>1) {
            int temp = i/1600;
            int layer = temp/2;
            int view = temp%2;
            int strip = i%1600;
            float occupancy = hits[i]*1./m_processed;
            res << "hit layer, view, strip: " 
                << layer << " " << view << " " << strip << " nhits " << nhits << " Occupancy " << occupancy 
                << std::endl;
        }
    }

    res.close();
 
    return sc;
}



void UserAlg::processEventList() {
    const std::vector<std::string>& theIds = m_eventList.value( );
    std::vector<std::string>::const_iterator it;
    std::vector<std::string>::const_iterator itend = theIds.end( );
    for (it = theIds.begin(); it != itend; it++) {
        int len = (*it).size();
        int delimPos = (*it).find_first_of('_');
        int runNum;
        if (delimPos >= 0) {
            runNum = atoi((*it).substr(0, delimPos-1).c_str());
        } else {
            int runNum = 0;
            m_runNumFlag = false;
        }
        
        int eventId = atoi((*it).substr(delimPos+1, len).c_str());
        std::vector<int>& curList = m_idMap[runNum];
        curList.push_back(eventId);
        
    }
}