/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MCT1ReadPreProc
//
// Reads a output file of the CT1 preproc.
//
//  Input Containers:
//   -/-
//
//  Output Containers:
//    MCerPhotEvt     the data container for all data.
//    MPedestalCam    ct1 pedestals
//    MMcEvt          monte carlo data container for MC files
//    MMcTrig         mc data container for trigger information
//    MSrcPosCam      source position in the camera
//    MBlindPixels    Array holding blind pixels
//
/////////////////////////////////////////////////////////////////////////////

#include "MCT1ReadPreProc.h"

#include <fstream.h>

#include <TList.h>
#include <TSystem.h>

#define LINUX
#define HISTO void
#define HBOOK_FILE int
#include "defines.h"
#include "structures.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"
#include "MCerPhotEvt.h"

#include "MPedestalPix.h"
#include "MPedestalCam.h"

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MBlindPixels.h"

#include "MMcEvt.hxx"
#include "MMcTrig.hxx"

ClassImp(MCT1ReadPreProc);

// --------------------------------------------------------------------------
//
// Default constructor. Creates an array which stores the file names of
// the files which should be read. If a filename is given it is added
// to the list.
//
MCT1ReadPreProc::MCT1ReadPreProc(const char *fname, const char *name,
                                 const char *title) : fIn(NULL)
{
    fName  = name  ? name  : "MCT1ReadPreProc";
    fTitle = title ? title : "Task to loop over events in CT1 ascii file";

    //
    // remember file name for opening the file in the preprocessor
    //
    fFileNames = new TList;
    fFileNames->SetOwner();

    if (fname)
        AddFile(fname);
}

// --------------------------------------------------------------------------
//
// Delete the filename list and the input stream if one exists.
//
MCT1ReadPreProc::~MCT1ReadPreProc()
{
    delete fFileNames;
    if (fIn)
        delete fIn;
}

// --------------------------------------------------------------------------
//
// Add this file as the last entry in the chain
//
void MCT1ReadPreProc::AddFile(const char *txt)
{
    TNamed *name = new TNamed(txt, "");
    fFileNames->AddLast(name);
}

// --------------------------------------------------------------------------
//
// Print data from the header to the screen and analyse the header data,
// means store and copy the needed data into Mars structures and
// data members
//
void MCT1ReadPreProc::ProcessHeader(const struct outputpars &outpars)
{
    fNumPixels = outpars.inumpixels; // number of pixels in the camera

    if (fNumPixels>iMAXNUMPIX || fNumPixels==0)
        fNumPixels=iMAXNUMPIX;

    //
    // ------------------- Output some stuff -----------------------
    //

    // int     itelescope;       // number of the CT which took the data
    *fLog << inf << "Telescope: CT" << outpars.itelescope;

    // float   flongitude_deg;   // longitude (counted positive towards West) of CT position */
    // float   flatitude_deg;    // latitude (counted positive towards North) of CT position */
    *fLog << " located @ Longitude=" << outpars.flongitude_deg;
    *fLog << "deg  Latitude=" << outpars.flatitude_deg << "deg" << endl;

    // int     irunnum;          // run number (from parameters file)
    // enum    onoroff {NEITHER_ON_NOR_OFF, OFF_SOURCE, ON_SOURCE} eruntype; // runtype
    *fLog << "Run:       #" << outpars.irunnum << "  (";
    switch (outpars.eruntype)
    {
    case NEITHER_ON_NOR_OFF: *fLog << "unknown";    break;
    case OFF_SOURCE:         *fLog << "off-source"; break;
    case ON_SOURCE:          *fLog << "on-source";  break;
    default:                 *fLog << (int)outpars.eruntype; break;
    }
    *fLog << ", ";
    switch (outpars.etrackmode)
    {
    case NORMAL:  *fLog << "normal tracking";  break;
    case REVERSE: *fLog << "reverse tracking"; break;
    case DUNNO:   *fLog << "unknown tracking"; break;
    default:      *fLog << (int)outpars.etrackmode; break;
    }
    *fLog << ")" << endl;

    //double  dsourcera_hours;  // right ascension of observed source in hours
    //double  dsourcedec_deg;   // declination of observed source in degrees
    *fLog << "Source:    RA=" << outpars.dsourcera_hours << "h  DEC=";
    *fLog << outpars.dsourcedec_deg << "deg" << endl;

    *fLog << "Pixels:    " << fNumPixels << endl;

    //int     inummuonpixels;   // number of pixels in the muon shield
    //int     inumcointdcs;     // number of coincidence tdcs recorded in the runfile
    //float   fpixdiameter_deg; // smallest pixel diameter (degrees) (from parameters file) */

    // enum    axes {RA, DEC, ALT, AZ} ese1_is; // name of the axis to which shaft encoder 1 is attached (implies the type of mount)
    *fLog << "Shaftencoder 1 @ ";
    switch (outpars.ese1_is)
    {
    case RA:  *fLog << "RA";  break;
    case DEC: *fLog << "DEC"; break;
    case ALT: *fLog << "ALT"; break;
    case AZ:  *fLog << "AZ";  break;
    default:  *fLog << (int)outpars.ese1_is; break;
    }
    *fLog << endl;

    // int     isezeropos[2];       // zero position of shaftencoders 1 and 2 (from parameters file)
    *fLog << "SE Zero:   SE(1)=" << outpars.isezeropos[0] << "  ";
    *fLog << "SE(2)=" << outpars.isezeropos[1] << endl;

    // int     iaz_rev_track_corr;  // correction for the azimuth shaft encoder (ALT/AZ mount only) in reverse tracking mode
    // int     ialt_rev_track_corr; // correction for the altitude shaft encoder (ALT/AZ mount only) in reverse tracking mode
    *fLog << "Reverse tracking corrections: SE(az)=" << outpars.iaz_rev_track_corr;
    *fLog << "  SE(alt)=" << outpars.ialt_rev_track_corr << endl;

    // float   fbendingcorr;        // bending correction factor (ALT/AZ mount only)
    // float   fextinction;         // atmospheric extinction (typically taken from the Carlsberg Meridian Circle data)
    *fLog << "Bending: Correction factor=" << outpars.fbendingcorr << "  ";
    *fLog << "Extinction=" << outpars.fextinction << endl;

    // Boolean bdontusepix[iMAXNUMPIX]; // bdontusepix is set true if the pixel should not be used in image analysis, otherwise it is true;
    fBlinds->Clear();
    *fLog << "Don't use pixels: ";
    for (int i=0; i<iMAXNUMPIX; i++)
        if (outpars.bdontusepix[i])
        {
            *fLog << i << " ";
            fBlinds->SetPixelBlind(i);
        }
    *fLog << endl;

    *fLog << "Exclude: ";
    // Boolean bexcludepix[iMAXNUMPIX];
    for (int i=0; i<iMAXNUMPIX; i++)
        if (outpars.bdontusepix[i])
            *fLog << i << " ";
    *fLog << endl;
    /* bexcludepix[] is set TRUE (== exclude from pedestal, Laser
     * calibration and the further analysis) when the Mean value
     * of a pixel inside a pedestal Run is larger than 50 or ( || )
     * if the pedestal RMS value of this pixel is larger than 5.0
     * This is reflected in the (new for versions >= 0.4)
     * variable "pixonoff" in the ntuple written by preproc:
     * preproc.nt.hbook
     *
     * When the pixel is excluded by the user it will get a -2 otherwise
     * pixonoff = 0.
     * Additive to this a -1 is added when preproc excludes the pixel
     * for a given Run. So the actual value tells you whether you caught
     * it already by hand or not.
     *
     * A plot of pixonoff may also be handy to tell you the status of your
     * ADC equipment. */

    // float   fphotoel_per_adccnt[iMAXNUMPIX]; // conversion factors for the pixel signals */
    /*
    float padc = outpars.fphotoel_per_adccnt[0];
    *fLog << "Phe/ADC (pixel 0): " << padc << endl;
    for (int i=0; i<iMAXNUMPIX; i++)
        *fLog << outpars.fphotoel_per_adccnt[i] << " ";
    *fLog << endl;
    */
    /*
     --- USEFULL? NEEDED? ---
     int     irubminusutc_usecs;              // difference between rubidium clock and UTC in microseconds
     int     isum_thresh_phot;                // threshold for the total sum of photoelectrons filter
     int     i2out_thresh_phot;               // threshold for the two-pixels-out-of-all software
     int     imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; // thresholds for the muon cut
     Boolean bmuon_suppression;               // "Dolby" noise reduction flag
     float   ftolerated_pointerror_deg;       // maximum tolerated pointing error in the position
     */

    // float fxpointcorr_deg; // pointing correction (to be added along the camera x axis) e.g. from point run */
    // float fypointcorr_deg; // pointing correction (to be added along the camera y axis) e.g. from point run */
    *fLog << "Pointing correction: dx=" << outpars.fxpointcorr_deg << "deg  ";
    *fLog << "dy=" << outpars.fxpointcorr_deg << "deg" << endl;

    /*
     fSrcPos->SetXY( outpars.fxpointcorr_deg/fGeom->GetConvMm2Deg(),
     outpars.fypointcorr_deg/fGeom->GetConvMm2Deg());
     */

    /*
     --- USEFULL? NEEDED? ---
     float   fcamera_align_angle_deg;         // the angle between the camera y-axis and the meridian when a culminating object is observed (defined counter-clockwise looking at the sky)
     int     iratecalc_numevents_odd;         // number of events used in the rate calculation (must be odd)
     int     inumpedfile;                     // number of the pedestal file used
     int     inumpedrun;                      // number of the pedestal run used in the file (starting at 0)
     int     inumcalfile;                     // number of the calibration file used
     int     inumlaserrun;                    // number of the laserrun used in the file (starting at 0)
     int     inumtellogfile;                  // number of the TelLog file to be used
     int     inumtellogrun;                   // number of the tellog entry (Runnumber) used from the log file
     int     imcparticle;                     // CORSIKA-coded Monte Carlo particle type.
    */

    // ----- preprocessing results -----

    // int     istart_mjdate_day;                 // MJD of run start (first event) */
    // int     iend_mjdate_day;                   // MJD of run end (last event) */
    // int     irunduration_secs;                 // difference between start and end time (secs) */
    *fLog << "Run Time: From " << outpars.istart_mjdate_day << " to ";
    *fLog << outpars.iend_mjdate_day << " (MJD),  Duration=";
    *fLog << outpars.irunduration_secs/3600 << "h";
    *fLog << (outpars.irunduration_secs/60)%60 << "m";
    *fLog << outpars.irunduration_secs%60 << "s" << endl;

    /*
     --- USEFULL? NEEDED? ---
     int     iproc_mjdate;                      // MJD of data processing (i.e. creation of this file)
     */

    // int     iproc_evts;                        // number of events processed */
    *fLog << "Number of processed events: " << outpars.iproc_evts << endl;

    // --- USEFULL? NEEDED? ---
    // double  dactual_sourcera_hours;            // for off runs: the false source (that should have been) observed */

    // float   frms_pedsig_phot[iMAXNUMPIX];      // standard deviation of the calibrated signals from the pedestal run */
    fPedest->InitSize(iMAXNUMPIX);
    for (Int_t i=0; i<iMAXNUMPIX; i++)
        (*fPedest)[i].SetMeanRms(outpars.frms_pedsig_phot[i]);

    // Used to communicate the mean over all pixels
    // pedestal RMS into the Runs NTuple, as it might
    // be used for e.g. quality cuts.
    // float   fpedrms_mean;
    *fLog << "Pedestal RMS: " << outpars.fpedrms_mean << endl;

    // The average current over the active pixels
    // for this run. This is done separately for
    // ON and OF runs.
    //float   fcurrent_mean;

    // enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK} eerrortolerance;
    /* 0 == "cautious", exits on any reason (but tells in
     * the .err file,
     * 1 == "goodphysics", exits when physics could be affected
     * by the error,
     * 2 == "tank", never exits except on coredumps and when
     * all files have been processed. Do not use such files for
     * physics analysis!
     *
     * NOTE: the capital letter words are the enums, the small letter
     * words must be used inside the parameter file. */

    // enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG} emctriggerflag;
    /* all: all events which survive the filter are written to the
     *      events NTuple.
     * flag: When Dorota's triggerflag is set to 1 for a particular
     *       event, it shall be written to the output. All others shall
     *       just be disregarded. (Default)
     * noflag: Opposite of 'flag': only events with triggerflag = 0 shall
     *         be treated further. */

    if (outpars.bmontecarlo)
        *fLog << "File is a Monte Carlo file." << endl;

    fIsMcFile = outpars.bmontecarlo==TRUE;
}

// --------------------------------------------------------------------------
//
// Read CT1 PreProc File Header:
//
void MCT1ReadPreProc::ReadHeader()
{
    char cheadertitle[iHEADERTITLELENGTH];
    fIn->read(cheadertitle, iHEADERTITLELENGTH);

    *fLog << cheadertitle << flush;

    // cTITLE_TEMPLATE "PREPROC V%f/S%f CT %d RUN %d %d PROCMJD %d\n"
    struct outputpars outpars;

    Float_t fpreprocversion, structversion;
    sscanf(cheadertitle, cTITLE_TEMPLATE,
           &fpreprocversion,    &structversion,
           &outpars.itelescope, &outpars.irunnum,
           &outpars.eruntype,   &outpars.iproc_mjdate);

    if (STRUCT_VERSION != structversion)
    {
        *fLog << warn << "WARNING: Version of C-structures of file (V";
        *fLog << structversion << ") not identical with current structures (V";
        *fLog << STRUCT_VERSION << ")" << endl;
    }

    fIn->read((Byte_t*)&outpars, sizeof(struct outputpars));

    ProcessHeader(outpars);
}

void MCT1ReadPreProc::ReadFooter()
{
    char cheadertitle[iHEADERTITLELENGTH];
    fIn->read(cheadertitle, iHEADERTITLELENGTH);
    /*
     ssscanf(cheadertitle, cEND_EVENTS_TEMPLATE,
             &filterres.ifilter_passed_evts);
     */

    struct filterresults filterres;
    fIn->read((Byte_t*)&filterres, sizeof(struct filterresults));
    /*
     int   imax_alt_arcs;            // maximum altitude reached during the run
     int   iaz_at_max_alt_arcs;      // azimuth at the time the max. alt. was reached
     int   itimeaverage_alt_arcs;    // altitude averaged over the runtime
     int   icoord_inconsist_evts;    // number of events with time-coordinate inconsistency in this run
     int   ifilter_passed_evts;      // number of events which passed the filter
     int   imuon_fail_evts;          // number of events rejected as muons (other filters passed)
     int   i2out_fail_evts;          // number of events which failed in the two out of all pixels software trigger
     int   imuon_and_2out_fail_evts; // number of events failing in both muon and 2out filter
     int   isum_fail_evts;           // number of events which failed the sum-of-all-calibrated ADC counts filter
     int   isum_and_muon_fail_evts;  // number of events which failed in both the sum and the muon filter
     int   isum_and_2out_fail_evts;  // number of events which failed in both the sum and the 2out filter
     int   iall_filters_fail_evts;   // number of events which failed in all filters
     float favg_event_rate_hz;       // average rate before filtering
     float fstddev_event_rate_hz;    // standard deviation of the rate before filtering
     */
}

// --------------------------------------------------------------------------
//
// This opens the next file in the list and deletes its name from the list.
//
Bool_t MCT1ReadPreProc::OpenNextFile()
{
    //
    // open the input stream and check if it is really open (file exists?)
    //
    if (fIn)
        delete fIn;
    fIn = NULL;

    //
    // Check for the existance of a next file to read
    //
    TNamed *file = (TNamed*)fFileNames->First();
    if (!file)
        return kFALSE;

    //
    // open the file which is the first one in the chain
    //
    const char *name = file->GetName();

    fIn = new ifstream(gSystem->ExpandPathName(name));

    if (!(*fIn))
    {
        *fLog << dbginf << "Cannot open file '" << name << "'" << endl;
        return kFALSE;
    }

    *fLog << "Open file: '" << name << "'" << endl;

    //
    // Remove this file from the list of pending files
    //
    fFileNames->Remove(file);

    *fLog << inf << "-----------------------------------------------------------------------" << endl;
    const Int_t sizef = sizeof(struct filterresults)+iHEADERTITLELENGTH;
    const Int_t sizeh = sizeof(struct outputpars)   +iHEADERTITLELENGTH;

    fIn->seekg(0, ios::end);
    const Int_t filesize = fIn->tellg();

    fIn->seekg(filesize-sizef, ios::beg);
    ReadFooter();
    fIn->seekg(0, ios::beg);
    ReadHeader();

    struct eventrecord event;
    const int size1   = sizeof(event)-sizeof(event.spixsig_10thphot);
    const int size2   = sizeof(event.spixsig_10thphot[0])*fNumPixels;
    const int evtsize = size1+size2;

    *fLog << "Storage for events: " << filesize-sizef-sizeh << "b" << endl;
    *fLog << "Size of one Event:  " << evtsize << "b" << endl;

    fNumEvents = (filesize-sizef-sizeh)/evtsize-1;
    fNumEvent = 0;

    *fLog << "Calculated Number of Events: " << (float)(filesize-sizef-sizeh)/evtsize << " " <<  (filesize-sizef-sizeh)%evtsize << endl;

    *fLog << "-----------------------------------------------------------------------" << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Open the first file in the list. Check for the output containers or create
// them if they don't exist.
//
// Initialize the size of the MPedestalCam container to 127 pixels (CT1 camera)
//
Bool_t MCT1ReadPreProc::PreProcess(MParList *pList)
{
    //
    //  look for the MCerPhotEvt class in the plist
    //
    fNphot = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt");
    if (!fNphot)
        return kFALSE;

    //
    //  look for the pedestal class in the plist
    //
    fPedest = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
    if (!fPedest)
        return kFALSE;

    //
    //  look for the pedestal class in the plist
    //
    fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
    if (!fBlinds)
        return kFALSE;

    //
    //  look for the source position in the camera
    //
    fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam");
    if (!fSrcPos)
        return kFALSE;

    //
    //  look for the camera geometry
    //
    fGeom = (MGeomCam*)pList->FindCreateObj("MGeomCamCT1", "MGeomCam");
    if (!fGeom)
        return kFALSE;

    //
    //  look for the mc event class
    //
    fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt");
    if (!fMcEvt)
        return kFALSE;

    //
    //  look for the mc trigger class
    //
    fMcTrig = (MMcTrig*)pList->FindCreateObj("MMcTrig");
    if (!fMcTrig)
        return kFALSE;

    //
    // Try to open at least one (the first) file
    //
    if (!OpenNextFile())
        return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Analyse the event data, means store and copy the needed data into
// Mars structures and data members
//
void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
{
    /*
     --- USEFULL? NEEDED? ---
     int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
     int   isecfrac_200ns;     // fractional part of isecs_since_midday
     short snot_ok_flags;      // the bits in these two bytes are flags for additional information on the event: Everything OK =: all Bits = 0
     int   ialt_arcs;          // altitude (arcseconds)
     int   iaz_arcs;           // azimuth (arcseconds)
     int   ipreproc_alt_arcs;  // "should be" alt according to preproc (arcseconds)
     int   ipreproc_az_arcs;   // "should be" az according to preproc (arcseconds)

     // for ALT-AZ mount telescopes: rotation angle of the field of
     // view; this angle is defined mathematically positive looking
     // towards the sky as the angle between the hour circle through
     // the object being tracked and the line through pixel 1 and 2
     int   ifieldrot_arcs;

     // event rate in milli Hertz before filtering calculated by
     // iratecalc_numevents_odd/(time[i+iratecalc_numevents_odd/2] -
     // time[i-iratecalc_numevents_odd/2])
     // For the first and the last iratecalc_numevents_odd/2
     // events the rate is assumed to be constant
     unsigned short srate_millihz;

     // This is the angle between the observation of this event and the
     // culmination point. It is going to be written into the events NTuple.
     float fhourangle;
     */
    /*
     *fLog << event.isecs_since_midday << "s ";
     *fLog << event.ialt_arcs << "s ";
     *fLog << event.iaz_arcs << "s ";
     *fLog << event.ipreproc_alt_arcs << "s ";
     *fLog << event.ipreproc_az_arcs << "s ";
     *fLog << event.ifieldrot_arcs << "s ";
     *fLog << event.fhourangle << endl;
     */

    //
    // read in the number of cerenkov photons and add the 'new' pixel
    // too the list with it's id, number of photons and error
    //
    fNphot->InitSize(fNumPixels);

    // number of photoelectrons measured in each pixel only the
    // actual number of pixels (outputpars.inumpixels) is written out
    // short spixsig_10thphot[iMAXNUMPIX];
    for (Int_t i=0; i<fNumPixels; i++)
    {
        if (event.spixsig_10thphot[i]>0)
            fNphot->AddPixel(i, 0.1*event.spixsig_10thphot[i],
                             (*fPedest)[i].GetMeanRms());
    }

    if (!fIsMcFile)
        return;

    fMcEvt->SetPartId(event.imcparticle);        // corsika particle type
    fMcEvt->SetEnergy(event.fmcenergy_tev*1000); // simulated energy
    fMcEvt->SetImpact(event.imcimpact_m*100);    // simulated impact

    fMcTrig->SetFirstLevel(event.imctriggerflag); // MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */

    //float fmcsize_phel;   // Simulated SIZE
}

// --------------------------------------------------------------------------
//
// Check for the event number and depending on this number decide if
// pedestals or event data has to be read.
//
// If the end of the file is reached try to open the next in the list. If
// there is now next file stop the eventloop.
//
Bool_t MCT1ReadPreProc::Process()
{
    //
    // check if we are done. Try to open the next file in chain.
    // If it was possible start reading. If not break the event loop
    //
    if (fNumEvent++==fNumEvents-1)
        return OpenNextFile() ? kCONTINUE : kFALSE;

    // event data to be read from the file
    struct eventrecord event;

    const int size1 = sizeof(event)-sizeof(event.spixsig_10thphot);
    const int size2 = sizeof(event.spixsig_10thphot[0])*fNumPixels;

    // read the eventrecord with the recorded number of stored pixels
    fIn->read((Byte_t*)&event, size1+size2);

    ProcessEvent(event);

    return kTRUE;
}