/* ======================================================================== *\
!
! *
! * 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, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MJStar
//
// Resource file entries are case sensitive!
//
/////////////////////////////////////////////////////////////////////////////
#include "MJStar.h"

#include <TEnv.h>
#include <TFile.h>

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

#include "MDirIter.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MEvtLoop.h"

#include "MStatusDisplay.h"

#include "MH3.h"
#include "MHVsTime.h"
#include "MHCamEvent.h"
#include "MHCamEventRot.h"
#include "MBinning.h"

#include "MReadReports.h"
#include "MReadMarsFile.h"
#include "MF.h"
#include "MFDataMember.h"
#include "MFDeltaT.h"
#include "MFSoftwareTrigger.h"
#include "MContinue.h"
#include "MGeomApply.h"
#include "MEventRateCalc.h"
#include "MImgCleanStd.h"
#include "MSrcPosCalc.h"
#include "MHillasCalc.h"
#include "MMuonSearchParCalc.h"
#include "MMuonCalibParCalc.h"
#include "MFillH.h"
#include "MWriteRootFile.h"

#include "MMuonSetup.h"
#include "MObservatory.h"
#include "MPointingPosCalc.h"

ClassImp(MJStar);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
//
MJStar::MJStar(const char *name, const char *title) : fMuonAnalysis(kTRUE)
{
    fName  = name  ? name  : "MJStar";
    fTitle = title ? title : "Standard analysis and reconstruction";
}

Bool_t MJStar::CheckEnvLocal()
{
    DisableMuonAnalysis(!GetEnv("MuonAnalysis", fMuonAnalysis));
    return kTRUE;
}

Bool_t MJStar::WriteResult()
{
    if (fPathOut.IsNull())
    {
        *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
        return kTRUE;
    }

    const TString oname = Form("%s/star%08d.root", (const char*)fPathOut, fSequence.GetSequence());

    *fLog << inf << "Writing to file: " << oname << endl;

    TFile file(oname, "RECREATE");

    *fLog << inf << " - MStatusDisplay..." << flush;
    if (fDisplay && fDisplay->Write()<=0)
    {
        *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
        return kFALSE;
    }
    *fLog << inf << "ok." << endl;

    return kTRUE;
}

Bool_t MJStar::Process(Bool_t ismc)
{
    if (!fSequence.IsValid())
    {
        *fLog << err << "ERROR - Sequence invalid!" << endl;
        return kFALSE;
    }

    //if (!CheckEnv())
    //    return kFALSE;

    CheckEnv();

    // --------------------------------------------------------------------------------

    *fLog << inf;
    fLog->Separator(GetDescriptor());
    *fLog << "Calculate image parameters of sequence ";
    *fLog << fSequence.GetName() << endl;
    *fLog << endl;

    // --------------------------------------------------------------------------------

    MDirIter iter;
    const Int_t n0 = fSequence.SetupDatRuns(iter, MSequence::kCalibrated, fPathData);
    const Int_t n1 = fSequence.GetNumDatRuns();
    if (n0==0)
    {
        *fLog << err << "ERROR - No input files of sequence found!" << endl;
        return kFALSE;
    }
    if (n0!=n1)
    {
        *fLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
        if (fLog->GetDebugLevel()>4)
        {
            *fLog << dbg << "Files which are searched:" << endl;
            iter.Print();
        }
        return kFALSE;
    }

    // Setup Parlist
    MParList plist;
    plist.AddToList(this); // take care of fDisplay!

    MObservatory obs;
    plist.AddToList(&obs);

    MMuonSetup muonsetup;
    plist.AddToList(&muonsetup);

    // Setup binnings for muon analysis
    MBinning bins1("BinningRadius");
    MBinning bins2("BinningArcWidth");
    MBinning bins3("BinningRingBroadening");
    MBinning bins4("BinningSizeVsArcRadius");
    MBinning bins5("BinningMuonWidth");
    MBinning bins6("BinningArcPhi");
    plist.AddToList(&bins1);
    plist.AddToList(&bins2);
    plist.AddToList(&bins3);
    plist.AddToList(&bins4);
    plist.AddToList(&bins5);
    plist.AddToList(&bins6);


    // Setup Tasklist
    MTaskList tlist;
    plist.AddToList(&tlist);

    MReadReports readreal;
    readreal.AddTree("Events", "MTime.", MReadReports::kMaster);
    readreal.AddTree("Drive", MReadReports::kRequired);
    readreal.AddTree("CC");
    //read.AddTree("Trigger");
    //read.AddTree("Camera");
    //read.AddTree("Currents");
    readreal.AddFiles(iter);

    MReadMarsFile readmc("Events");
    readmc.DisableAutoScheme();
    readmc.AddFiles(iter);

    // ------------------ Setup general tasks ----------------

    MFDeltaT               fdeltat;
    MContinue              cont(&fdeltat, "FilterDeltaT", "Filter events with wrong timing");
    cont.SetInverted();

    MGeomApply             apply; // Only necessary to craete geometry
    MEventRateCalc         rate;
    //MEventRateCalc         rate1; // 5min
    rate.SetNumEvents(1200);
    //rate1.SetNumEvents(60000);
    //rate1.SetNameEventRate("MEventRate2");
    //rate1.SetNameTimeRate("MTimeRate2");

    /*
    MEventRateCalc         rate10000;
    rate10000.SetNameEventRate("MEventRate10000");
    rate10000.SetNumEvents(10000);
 */
    //MBadPixelsMerge        merge(&badpix);

    MFSoftwareTrigger swtrig;
    MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
    contsw.SetInverted();


    MImgCleanStd           clean;
    clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");

    MSrcPosCalc            poscalc;
    MHillasCalc            hcalc;
    hcalc.Disable(MHillasCalc::kCalcConc);

    // ------------------ Setup histograms and fill tasks ----------------
    MHCamEvent evt0a(0, "Cleaned", "Signal after Cleaning;;S [\\gamma]");
    MHCamEvent evt0b(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
    evt0b.SetThreshold(0);

    //MHCamEventRot evt0r("UsedRot", "Pixels marked Used (derotated)");
    //evt0r.SetThreshold(0);

    MH3 h1("MEventRate.fRate");
    h1.SetName("MHEventRate");
    h1.SetTitle("Event Rate distribution;R [Hz];Counts");
    h1.SetLogy();
/*
    MH3 h12("MEventRate10000.fRate");
    h12.SetName("MHEventRate");
    h12.SetLogy();
 */
    MBinning b1("BinningMHEventRate");
    b1.SetEdges(150, 0, 1500);
    plist.AddToList(&b1);

    MHVsTime hvs("MEventRate.fRate");
    hvs.SetTitle("Rate per 500 events;;R [Hz]");
    hvs.SetNumEvents(500);
    hvs.SetMinimum(0);

    //MContinue cont1("MEventRate2.fRate/MEventRate.fRate>1.1");
    //MContinue cont2("MEventRate.fRate/MEventRate2.fRate>1.1");

    MFillH fillvs(&hvs,           "MTime",        "FillEventRate10s");

    MFillH fill0a(&evt0a,         "MSignalCam",   "FillSignalCam");
    MFillH fill0b(&evt0b,         "MSignalCam",   "FillCntUsedPixels");
    //MFillH fill0r(&evt0r,         "MCerPhotEvt",  "FillCntUsedRotated");
    MFillH fill1("MHHillas",      "MHillas",      "FillHillas");
    MFillH fill2("MHHillasExt",   "",             "FillHillasExt");
    MFillH fill3("MHHillasSrc",   "MHillasSrc",   "FillHillasSrc");
    MFillH fill4("MHImagePar",    "MImagePar",    "FillImagePar");
    MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
    //MFillH fill6("MHImageParTime","MImageParTime","FillImageParTime");
    //MFillH fill7("MHNewImagePar2","MNewImagePar2","FillNewImagePar2");
    MFillH fill8(&h1,             "",             "FillEventRate");
    MFillH fill9("MHEffectiveOnTime", "MTime",    "FillEffOnTime");
    //MFillH fillb(&h12, "", "FillEvtRate2");
    //MFillH fill9("MHCerPhot");

    //fill0r.SetDrawOption("colz");
    fill8.SetNameTab("EvtRate");
    fill9.SetNameTab("EffOnTime");

    // ------------------ Setup write task ----------------

    const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data()));
    MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW");
    MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW");
    // Data
    write.AddContainer("MHillas",       "Events");
    write.AddContainer("MHillasExt",    "Events");
    write.AddContainer("MHillasSrc",    "Events");
    write.AddContainer("MImagePar",     "Events");
    write.AddContainer("MNewImagePar",  "Events");
    //write.AddContainer("MNewImagePar2", "Events");
    //write.AddContainer("MImageParTime", "Events");
    write.AddContainer("MRawEvtHeader", "Events");
    write.AddContainer("MPointingPos",  "Events");

    // Run Header
    write.AddContainer("MRawRunHeader", "RunHeaders");
    write.AddContainer("MBadPixelsCam", "RunHeaders");
    write.AddContainer("MGeomCam",      "RunHeaders");
    write.AddContainer("MObservatory",  "RunHeaders");

    // Muon Setup
    write.AddContainer("BinningRadius",            "RunHeaders");
    write.AddContainer("BinningArcWidth",          "RunHeaders");
    write.AddContainer("BinningRingBroadening",    "RunHeaders");
    write.AddContainer("BinningSizeVsArcRadius",   "RunHeaders");
    write.AddContainer("MMuonSetup",               "RunHeaders");

    if (ismc)
    {
        // Monte Carlo Data
        write.AddContainer("MMcEvt",               "Events");
        write.AddContainer("MMcTrig",              "Events");
        // Monte Carlo Run Headers
        write.AddContainer("MMcRunHeader",         "RunHeaders");
        write.AddContainer("MMcTrigHeader",        "RunHeaders");
        write.AddContainer("MMcFadcHeader",        "RunHeaders");
        write.AddContainer("MMcConfigRunHeader",   "RunHeaders");
        write.AddContainer("MMcCorsikaRunHeader",  "RunHeaders");
    }
    else
    {
        write.AddContainer("MTime",                "Events");
        // Drive
        write.AddContainer("MReportDrive",         "Drive");
        write.AddContainer("MTimeDrive",           "Drive");
        // Effective On Time
        write.AddContainer("MEffectiveOnTime",     "EffectiveOnTime");
        write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime");
    }

    // What to write in muon tree
    writem.AddContainer("MMuonSearchPar", "Muons");
    writem.AddContainer("MMuonCalibPar",  "Muons");
    writem.AddContainer("MHillas",        "Muons");
    writem.AddContainer("MHillasExt",     "Muons");
    writem.AddContainer("MHillasSrc",     "Muons");
    writem.AddContainer("MImagePar",      "Muons");
    writem.AddContainer("MNewImagePar",   "Muons");
    writem.AddContainer("MRawEvtHeader",  "Muons");
    writem.AddContainer("MPointingPos",   "Muons");
    if (ismc)
    {
        // Monte Carlo Data
        writem.AddContainer("MMcEvt",     "Muons");
        writem.AddContainer("MMcTrig",    "Muons");
    }

    if (ismc)
        if (fMuonAnalysis)
            writem.AddCopySource("OriginalMC");
        else
            write.AddCopySource("OriginalMC");

    MTaskList tlist2("Events");
    tlist2.AddToList(&apply);
    if (!ismc)
        tlist2.AddToList(&cont);
    tlist2.AddToList(&contsw);
    if (!ismc)
    {
        // Calibration events don't enter star at all.
        tlist2.AddToList(&rate);
        tlist2.AddToList(&fillvs);
        tlist2.AddToList(&fill8);
        tlist2.AddToList(&fill9);
    }
    //tlist2.AddToList(&fillb);
    tlist2.AddToList(&clean);
    tlist2.AddToList(&fill0a);
    tlist2.AddToList(&fill0b);
    //tlist2.AddToList(&fill0r);
    tlist2.AddToList(&poscalc);
    tlist2.AddToList(&hcalc);
    tlist2.AddToList(&fill1);
    tlist2.AddToList(&fill2);
    tlist2.AddToList(&fill3);
    tlist2.AddToList(&fill4);
    tlist2.AddToList(&fill5);

    // ----------------------- Muon Analysis ----------------------
    MF fmuon1("MHillas.fSize>150", "MuonPreCut");
    MF fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&"
              "(MMuonSearchPar.fDeviation<45)", "MuonSearchCut");
    MF fmuon3("(MMuonCalibPar.fArcPhi>190)  && (MMuonSearchPar.fDeviation<35) &&"
              "(MMuonCalibPar.fArcWidth<0.20) && (MMuonCalibPar.fArcWidth>0.04)",
              "MuonFinalCut");
    MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
    writem.SetFilter(&fmuon4);

    MMuonSearchParCalc muscalc;
    muscalc.SetFilter(&fmuon1);

    MMuonCalibParCalc mcalc;
    mcalc.SetFilter(&fmuon2);

    MFillH fillmuon("MHSingleMuon", "", "FillMuon");
    MFillH fillmpar("MHMuonPar",    "", "FillMuonPar");
    fillmuon.SetFilter(&fmuon2);
    fillmpar.SetFilter(&fmuon3);

    if (fMuonAnalysis)
    {
        tlist2.AddToList(&fmuon1);
        tlist2.AddToList(&muscalc);
        tlist2.AddToList(&fmuon2);
        tlist2.AddToList(&fillmuon);
        tlist2.AddToList(&mcalc);
        tlist2.AddToList(&fmuon3);
        tlist2.AddToList(&fillmpar);
        tlist2.AddToList(&fmuon4);
        tlist2.AddToList(&writem);
    }
    // ------------------------------------------------------------

    MFillH fillw("MHWeather", "MTimeCC", "FillWeather");
    //tlist2.AddToList(&fill6);
    //tlist2.AddToList(&fill7);
    //tlist2.AddToList(&fill9);

    MPointingPosCalc pcalc;

    tlist.AddToList(ismc ? (MTask*)&readmc : (MTask*)&readreal);
    tlist.AddToList(&pcalc,  ismc ? "Events" : "Drive");
    tlist.AddToList(&tlist2, "Events");
    if (!ismc)
        tlist.AddToList(&fillw,  "CC");
    tlist.AddToList(&write);

    // Create and setup the eventloop
    MEvtLoop evtloop(fName);
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(fDisplay);
    evtloop.SetLogStream(fLog);
    if (!SetupEnv(evtloop))
        return kFALSE;

    // Execute first analysis
    if (!evtloop.Eventloop(fMaxEvents))
    {
        *fLog << err << GetDescriptor() << ": Failed." << endl;
        return kFALSE;
    }

    if (!WriteResult())
        return kFALSE;

    *fLog << all << GetDescriptor() << ": Done." << endl;
    *fLog << endl << endl;

    return kTRUE;
}
