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

/////////////////////////////////////////////////////////////////////////////
//
//  MJSimulation
//
//
// Force reading a corsika file even if the footer (RUNE-section) is missing
// by setting fForceMode to kTRUE or from the resource file by
//
//    ForceMode: Yes
//
//
// In case of a pedestal or calibration run the artificial trigger can
// be "switched off" and the cosmics trrigger "switched on" by setting
// fForceTrigger to kTRUE or from the resource file by
//
//    ForceTrigger: Yes
//
//
/////////////////////////////////////////////////////////////////////////////
#include "MJSimulation.h"

#include <TEnv.h>

// Core
#include "MLog.h"
#include "MLogManip.h"

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

#include "MStatusDisplay.h"

// Tasks
#include "MCorsikaRead.h"
#include "MContinue.h"
#include "MFillH.h"
#include "MGeomApply.h"
#include "MHillasCalc.h"
#include "MImgCleanStd.h"
#include "MWriteRootFile.h"

#include "MSimAbsorption.h"
#include "MSimReflector.h"
#include "MSimPointingPos.h"
#include "MSimPSF.h"
#include "MSimGeomCam.h"
#include "MSimSignalCam.h"
#include "MSimAPD.h"
#include "MSimExcessNoise.h"
#include "MSimCamera.h"
#include "MSimTrigger.h"
#include "MSimReadout.h"
#include "MSimRandomPhotons.h"
#include "MSimBundlePhotons.h"
#include "MSimCalibrationSignal.h"

// Histograms
#include "MBinning.h"

#include "MHn.h"
#include "MHCamEvent.h"
#include "MHPhotonEvent.h"

// Container
#include "MRawRunHeader.h"
#include "MParameters.h"
#include "MReflector.h"
#include "MParEnv.h"
#include "MPulseShape.h"

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

ClassImp(MJSimulation);

using namespace std;

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

Bool_t MJSimulation::CheckEnvLocal()
{
    fForceMode    = GetEnv("ForceMode",    fForceMode);
    fForceTrigger = GetEnv("ForceTrigger", fForceTrigger);

    return kTRUE;
}

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

    TObjArray cont;
    cont.Add(const_cast<TEnv*>(GetEnv()));
    //cont.Add(const_cast<MSequence*>(&fSequence));

    if (fDisplay)
    {
//        TString title = "--  Reflector: ";
//        title += fSequence.GetSequence();
//        title += "  --";
//        fDisplay->SetTitle(title, kFALSE);
        fDisplay->SetTitle("Ceres", kFALSE);

        cont.Add(fDisplay);
    }

//    const TString oname = Form("reflector%08d.root", fSequence.GetSequence());
    const TString oname = "ceres.root";
    return WriteContainer(cont, oname, "RECREATE");
}

void MJSimulation::SetupHist(MHn &hist) const
{
    hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
    hist.InitName("Energy");
    hist.InitTitle("Energy;E [GeV]");
    hist.SetLog(kTRUE, kTRUE, kFALSE);

    hist.AddHist("MPhotonEvent.GetNumExternal");
    hist.InitName("Size");
    hist.InitTitle("Size;S [#]");
    hist.SetLog(kTRUE, kTRUE, kFALSE);

    hist.AddHist("MCorsikaEvtHeader.fX/100","MCorsikaEvtHeader.fY/100");
    hist.SetDrawOption("colz");
    hist.InitName("Impact;Impact;Impact");
    hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");

    hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
    hist.InitName("Height");
    hist.InitTitle("FirstInteractionHeight;h [km]");

    hist.AddHist("MCorsikaEvtHeader.fAz*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()");
    hist.InitName("SkyOrigin;Az;Zd");
    hist.InitTitle("Sky Origin;Az [deg];Zd [deg]");
    hist.SetDrawOption("colz");

    TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
    TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
    TString cos  = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";

    TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()";

    hist.AddHist(form);
    hist.InitName("ViewCone");
    hist.InitTitle("Incident Angle;\\alpha [\\deg]");
}

Bool_t MJSimulation::Process(const MArgs &args)
{
    /*
    if (!fSequence.IsValid())
    {
        *fLog << err << "ERROR - Sequence invalid!" << endl;
        return kFALSE;
    }
    */

    //if (!HasWritePermission(GetPathOut()))
    //    return kFALSE;

    *fLog << inf;
    fLog->Separator(GetDescriptor());

    if (!CheckEnv())
        return kFALSE;

    *fLog << warn << "FIXME: Monte Carlo simulation: Sequences not supported yet.";
    //*fLog << fSequence.GetFileName() << endl;
    *fLog << endl;

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

    //MDirIter iter;
    //if (fSequence.GetRuns(iter, MSequence::kCalibrated)<=0)
    //    return kFALSE;

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

    // setup TaskList
    MTaskList tasks;
    plist.AddToList(&tasks);

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

    // ---------- FIXME FIXME FIXME -----------
    // FIXME: Find a better name (What if both cams are identical?)
    MParEnv env1("GeomAPDs");  // Inheritance!
    MParEnv env2("MGeomCam");  // Inheritance!
    env1.SetClassName("MGeomCam");
    env2.SetClassName("MGeomCam");
    plist.AddToList(&env1);
    plist.AddToList(&env2);

    // FIXME: Allow the user to create other reflectors
    MReflector reflector;
    plist.AddToList(&reflector);

    plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");

    MPulseShape shape;
    plist.AddToList(&shape);

    MRawRunHeader header;
    header.SetValidMagicNumber();
    //header.InitFadcType(3);

    header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTData,        1, 1);
    if (args.GetNumArguments()==1)
    {
        if (!args.GetArgumentStr(0).CompareTo("pedestal", TString::kIgnoreCase))
        {
            header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTPedestal,    1, 1);
            header.SetSourceInfo("Pedestal");
        }
        if (!args.GetArgumentStr(0).CompareTo("calibration", TString::kIgnoreCase))
        {
            header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTCalibration,    1, 1);
            header.SetSourceInfo("Calibration");
        }
        if (!args.GetArgumentStr(0).CompareTo("pointrun", TString::kIgnoreCase))
            header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTPointRun,    1, 1);
    }

    // FIXME: Move to MSimPointingPos, MSimCalibrationSignal
    header.SetObservation("On", "MonteCarlo");
    plist.AddToList(&header);
    // ++++++++ FIXME FIXME FIXME +++++++++++++

    /*
    MPedestalCam pedcam;
    pedcam.Init(geomcam.GetNumPixels());
    for (UInt_t i=0; i<geomcam.GetNumPixels(); i++)
        pedcam[i].Set(128./header.GetScale(), 22.5/header.GetScale());
    plist.AddToList(&pedcam);
    */

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

    MCorsikaRead read;
    read.SetForceMode(fForceMode);

    for (int i=0; i<args.GetNumArguments();i ++)
        read.AddFile(args.GetArgumentStr(i));

    MSimAbsorption absapd("AbsorptionAPDs");
    MSimAbsorption absmir("AbsorptionMirrors");
    MSimAbsorption cones("AbsorptionCones");
    cones.SetUseTheta();

    MSimPointingPos pointing;

    MSimReflector reflect;
    reflect.SetNameGeomCam("GeomAPDs");
//  MSimStarField stars;

    MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
    MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");

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

    MBinning binse( 100,     1,   100000, "BinningEnergy", "log");
    MBinning binss( 100,     1, 10000000, "BinningSize",   "log");
    MBinning binsi(  50,  -250,      250, "BinningImpact");
    MBinning binsh(  50,     0,       12, "BinningHeight");
    MBinning binsaz(360,  -360,      360, "BinningAz");
    MBinning binszd( 70,     0,       70, "BinningZd");
    MBinning binsvc(155,     0,       31, "BinningViewCone");
    MBinning binsew(150,     0,       25, "BinningEvtWidth");
    MBinning binstr("BinningTrigPos");

    plist.AddToList(&binse);
    plist.AddToList(&binss);
    plist.AddToList(&binsi);
    plist.AddToList(&binsh);
    plist.AddToList(&binszd);
    plist.AddToList(&binsaz);
    plist.AddToList(&binsvc);
    plist.AddToList(&binsew);
    plist.AddToList(&binstr);

    MHn mhn1, mhn2, mhn3;
    SetupHist(mhn1);
    SetupHist(mhn2);
    SetupHist(mhn3);

    MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-MPulseShape.GetPulseWidth");
    mhtp.SetName("TrigPos");
    mhtp.SetTitle("Trigger position w.r.t. the first photon hitting an APD");

    MH3 mhew("MPhotonStatistics.fTimeLast-MPhotonStatistics.fTimeFirst - MRawRunHeader.GetFreqSampling/1000.*MRawRunHeader.fNumSamplesHiGain - 2*MPulseShape.GetPulseWidth");
    mhew.SetName("EvtWidth");
    mhew.SetTitle("Time between first and last photon hitting an APD");

    MFillH fillh1(&mhn1, "", "FillH1");
    MFillH fillh2(&mhn2, "", "FillH2");
    MFillH fillh3(&mhn3, "", "FillH3");
    MFillH filltp(&mhtp, "", "FillTriggerPos");
    MFillH fillew(&mhew, "", "FillEvtWidth");
    fillh1.SetNameTab("H1", "Distribution of Muons as simulated");
    fillh2.SetNameTab("H2", "Distribution of Muons as available after all");
    fillh3.SetNameTab("H3", "Distribution after trigger");
    filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
    fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting an APD");

    MHPhotonEvent planeG(1);     // Get from MaxImpact
    MHPhotonEvent plane0(2);     // Get from MReflector
    MHPhotonEvent plane1(2);
    MHPhotonEvent plane2(2);
    MHPhotonEvent plane3(2);
    MHPhotonEvent plane4(2);
    MHPhotonEvent planeF1(3);   // Get from MGeomCam
    MHPhotonEvent planeF2(header.IsPointRun()?4:3); // Get from MGeomCam

    //MHPSF psf;

    MFillH fillG(&planeG,   "MPhotonEvent", "FillGround");
    MFillH fill0(&plane0,   "MirrorPlane0", "FillReflector");
    MFillH fill1(&plane1,   "MirrorPlane1", "FillCamShadow");
    MFillH fill2(&plane2,   "MirrorPlane2", "FillCandidates");
    MFillH fill3(&plane3,   "MirrorPlane3", "FillReflected");
    MFillH fill4(&plane4,   "MirrorPlane4", "FillFocal");
    MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
    MFillH fillF2(&planeF2, "MPhotonEvent", "FillAPDs");
    //MFillH fillP(&psf,      "MPhotonEvent", "FillPSF");

    fillG.SetNameTab("Ground",     "Photon distribution at ground");
    fill0.SetNameTab("Reflector",  "Photon distribution at reflector plane w/o camera shadow");
    fill1.SetNameTab("CamShadow",  "Photon distribution at reflector plane w/ camera shadow");
    fill2.SetNameTab("Candidates", "'Can hit' photon distribution at reflector plane w/ camera shadow");
    fill3.SetNameTab("Reflected",  "Photon distribution after reflector projected to reflector plane");
    fill4.SetNameTab("Focal",      "Photon distribution at focal plane");
    fillF1.SetNameTab("Camera",    "Photon distribution which hit the detector");
    fillF2.SetNameTab("APDs",      "Photon distribution after cones");
    //fillP.SetNameTab("PSF",        "Photon distribution after cones for all mirrors");

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

    const TString rule1(Form("s/cer([0-9]+)/%s\\/ref$1.root/", Esc(fPathOut).Data()));
    const TString rule2(Form("s/cer([0-9]+)/%s\\/sig$1.root/", Esc(fPathOut).Data()));
    const TString rule3(Form("s/cer([0-9]+)/%s\\/cam$1.root/", Esc(fPathOut).Data()));

    MWriteRootFile write1( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
    MWriteRootFile write2( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
    MWriteRootFile write3( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");

    write1.SetName("WriteRef");
    write2.SetName("WriteSig");
    write3.SetName("WriteCam");

    write1.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
    write1.AddContainer("MPhotonEvent",      "Events");

    write2.AddContainer("MRawRunHeader",             "RunHeaders");
    write2.AddContainer("MGeomCam",                  "RunHeaders");
    write2.AddContainer("MSignalCam",                "Events");
    write2.AddContainer("MRawEvtHeader",             "Events");
    write2.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders");
    /* Monte Carlo Headers
    write.AddContainer("MMcRunHeader",              "RunHeaders", kFALSE);
    write.AddContainer("MMcFadcHeader",             "RunHeaders", kFALSE);
    write.AddContainer("MMcTrigHeader",             "RunHeaders", kFALSE);
    write.AddContainer("MMcConfigRunHeader",        "RunHeaders", kFALSE);
    write.AddContainer("MMcCorsikaRunHeader",       "RunHeaders", kFALSE);
    // Monte Carlo
    write.AddContainer("MMcEvt",                    "Events",     kFALSE);
    write.AddContainer("MMcTrig",                   "Events",     kFALSE);
    */// Data tree
    //write.AddContainer("MPedPhotFromExtractor",     "Events");
    //write.AddContainer("MPedPhotFromExtractorRndm", "Events");
    //write.AddContainer("MTime",                     "Events",     kFALSE);
    //write.AddContainer("MRawEvtHeader",             "Events");
    //write.AddContainer("MTriggerPattern",           "Events");

    write3.AddContainer("MRawEvtData",   "Events");
    write3.AddContainer("MRawRunHeader", "RunHeaders");
    write3.AddContainer("MGeomCam",      "RunHeaders");
    write3.AddContainer("MRawEvtHeader", "Events");
    write3.AddContainer("MPedestalCam",  "RunHeaders", kFALSE);

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

    MGeomApply apply;

    MSimPSF simpsf;

    MSimGeomCam simgeom;
    simgeom.SetNameGeomCam("GeomAPDs");

    MSimCalibrationSignal simcal;
    simcal.SetNameGeomCam("GeomAPDs");

    //  Dark Current: 4MHz = 0.004 GHz
    //  Light of night sky at La Palma: ~ 0.2 ph / cm^2 / sr / ns
    //  5deg camera: 6e-3 sr
    //  NSB = 0.2*6e-3

    MSimRandomPhotons  simnsb;
//    simnsb.SetFreq(0.005, 0.004);    // 5MHz/cm^2, 4MHz
    simnsb.SetFreq(0, 0.04);           // 40MHz fixed for each APD
    simnsb.SetNameGeomCam("GeomAPDs");

    // FIXME: How do we handle star-light? We need the rate but we also
    //        need to process it through the mirror!

    MSimAPD simapd;
    simapd.SetFreq(0.04); // Must be identical to MSimRanodmPhotons!!
    simapd.SetNameGeomCam("GeomAPDs");

    MSimExcessNoise   simexcnoise;
    MSimBundlePhotons simsum;
    MSimSignalCam     simsignal;
    MSimCamera        simcam;
    MSimTrigger       simtrig;
    MSimReadout       simdaq;

    MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");

    MParameterD pulpos("IntendedPulsePos");
//    pulpos.SetVal(40);  // [ns]
    plist.AddToList(&pulpos);

    MParameterD trigger("TriggerPos");
    trigger.SetVal(0);
    plist.AddToList(&trigger);

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

    MImgCleanStd clean(7, 4.5);
    clean.SetMethod(MImgCleanStd::kAbsolute);

    //clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");

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

    MHCamEvent evt0a(/*10*/0, "Signal",    "Average signal;;S [ph]");
    MHCamEvent evt0c(/*10*/0, "MaxSignal", "Maximum signal;;S [ph]");
    MHCamEvent evt0d(/*11*/8, "ArrTm",     "Time after first photon;;T [ns]");
    evt0a.SetErrorSpread(kFALSE);
    evt0c.SetCollectMax();

    MFillH fillx0a(&evt0a,             "MSignalCam",      "FillAvgSignal");
    MFillH fillx0c(&evt0c,             "MSignalCam",      "FillMaxSignal");
    MFillH fillx0d(&evt0d,             "MSignalCam",      "FillArrTm");
    MFillH fillx1("MHHillas",          "MHillas",         "FillHillas");
    //MFillH fillx2("MHHillasExt",       "",                "FillHillasExt");
    MFillH fillx3("MHHillasSrc",       "MHillasSrc",      "FillHillasSrc");
    //MFillH fillx4("MHImagePar",        "MImagePar",       "FillImagePar");
    //MFillH fillx5("MHNewImagePar",     "MNewImagePar",    "FillNewImagePar");

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

    if (header.IsDataRun())
    {
        tasks.AddToList(&read);
        tasks.AddToList(&pointing);
        //    if (header.IsPointRun())
        //        tasks.AddToList(&stars);
        //tasks.AddToList(&print);
        tasks.AddToList(&fillh1);
        tasks.AddToList(&fillG);
        if (!header.IsPointRun())
        {
            tasks.AddToList(&absapd);
            tasks.AddToList(&absmir);
        }
        tasks.AddToList(&reflect);
        if (!header.IsPointRun())
        {
            tasks.AddToList(&fill0);
            tasks.AddToList(&fill1);
            tasks.AddToList(&fill2);
            tasks.AddToList(&fill3);
            tasks.AddToList(&fill4);
            tasks.AddToList(&fillF1);
        }
        tasks.AddToList(&cones);
        tasks.AddToList(&fillF2);
        //if (header.IsPointRun())
        //    tasks.AddToList(&fillP);
        tasks.AddToList(&cont1);
        if (!header.IsPointRun())
            tasks.AddToList(&fillh2);
    }
    // -------------------------------
    tasks.AddToList(&apply);
    if (header.IsDataRun())
    {
        tasks.AddToList(&simpsf);
        tasks.AddToList(&simgeom);
        tasks.AddToList(&cont2);
    }
    if (header.IsPedestalRun() || header.IsCalibrationRun())
        tasks.AddToList(&simcal);
    tasks.AddToList(&simnsb);
    tasks.AddToList(&simapd);
    tasks.AddToList(&simexcnoise);
    tasks.AddToList(&simsum);
    tasks.AddToList(&simcam);
    if (header.IsDataRun() || fForceTrigger)
        tasks.AddToList(&simtrig);
    tasks.AddToList(&conttrig);
    tasks.AddToList(&simsignal);  // What do we do if signal<0?
    tasks.AddToList(&simdaq);
    if (!fPathOut.IsNull() && !HasNullOut())
    {
        tasks.AddToList(&write1);
        if (!header.IsPedestalRun())
            tasks.AddToList(&write2);
        tasks.AddToList(&write3);
    }
    // -------------------------------
    if (header.IsDataRun())
        tasks.AddToList(&fillh3);
    tasks.AddToList(&filltp);
    if (header.IsDataRun())
        tasks.AddToList(&fillew);
    if (!header.IsPedestalRun())
    {
        tasks.AddToList(&fillx0a);
        tasks.AddToList(&fillx0c);
        //tasks.AddToList(&clean);
        tasks.AddToList(&hcalc);
        tasks.AddToList(&fillx0d);
        tasks.AddToList(&fillx1);
        //tasks.AddToList(&fillx2);
        tasks.AddToList(&fillx3);
        //tasks.AddToList(&fillx4);
        //tasks.AddToList(&fillx5);
    }

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

    tasks.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

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

    if (binstr.IsDefault())
        binstr.SetEdgesLin(150, -shape.GetPulseWidth(),
                           header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetPulseWidth());

    header.Print();

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

    //-------------------------------------------
    // FIXME: Display Spline     in tab
    // FIXME: Display Reflector  in tab
    // FIXME: Display Routing(?) in tab
    // FIXME: Display Camera(s)  in tab
    //-------------------------------------------

    if (!WriteResult())
        return kFALSE;

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

    return kTRUE;
}