/* ======================================================================== *\
!
! *
! * 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/2002 <mailto:tbretz@uni-sw.gwdg.de>
!   Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

void flux()
{ 
    //
    // first we have to create our empty lists
    //
    MParList  parlist;
    MTaskList tasklist;

    parlist.AddToList(&tasklist);

    //
    // Define the two source positions
    //
    MSrcPosCam source("Source");
    MSrcPosCam antisrc("AntiSource");
    source.SetXY(0, 0);
    antisrc.SetXY(+240, 0);
    parlist.AddToList(&source);
    parlist.AddToList(&antisrc);

    //
    // Setup binning for the histograms
    //
    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(45, -90, 90);

    MBinning binse("BinningE");
    binse.SetEdgesLog(10, 10, 1e3);

    MBinning binstheta("BinningTheta");
    binstheta.SetEdges(7, -2.5, 32.5);

    MBinning binstime("BinningTime");
    binstime.SetEdges(5, 0, 50);

    MBinning binsdifftime("BinningTimeDiff");
    binsdifftime.SetEdges(50, 0, 0.1);

    parlist.AddToList(&binsalpha);
    parlist.AddToList(&binse);
    parlist.AddToList(&binstheta);
    parlist.AddToList(&binstime);
    parlist.AddToList(&binsdifftime);

    //
    // Setup our tasks
    //
    MReadMarsFile reader("Events", "~/data/Gamma*.root");
    reader.EnableBranch("MMcEvt.fTheta");

    MGeomCamMagic geomcam;
    parlist.AddToList(&geomcam);

    MMcTimeGenerate   rand;
    MMcPedestalCopy   pcopy;
    MMcPedestalNSBAdd pnsb;
    MCerPhotCalc      ncalc;
    MImgCleanStd      clean;
    MBlindPixelCalc   blind;
    MHillasCalc       hcalc;
    MHillasSrcCalc    hsrc1("Source",     "HillasSrc");
    MHillasSrcCalc    hsrc2("AntiSource", "HillasAntiSrc");
    MEnergyEstimate   estim;

    MFillH hfill1h("MHHillas",  "MHillas");
    MFillH hfill1m("MHStarMap", "MHillas");
    MFillH hfill2s("HSource     [MHHillasSrc]", "HillasSrc");
    MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
 
    MMcCollectionAreaCalc acalc;

    const Float_t alpha0 = 15; // [deg]

    MFAlpha fsrc ("HillasSrc",     '>', alpha0);
    MFAlpha fasrc("HillasAntiSrc", '>', alpha0);

    MFillH fillsp      ("FluxSrcTime   [MHAlphaEnergyTime]",  "HillasSrc");
    MFillH fillasp     ("FluxASrcTime  [MHAlphaEnergyTime]",  "HillasAntiSrc");
    MFillH fillsptheta ("FluxSrcTheta  [MHAlphaEnergyTheta]", "HillasSrc");
    MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");

    MFTriggerLvl1 lvl1;
    MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt");
    MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt");
    MFillH fetimeall ("AllTime  [MHEnergyTime]",  "MMcEvt");
    MFillH fetimesel ("SelTime  [MHEnergyTime]",  "MMcEvt");

    fetimesel.SetFilter(&lvl1);
    fethetasel.SetFilter(&lvl1);

    fillsp.SetFilter(&fasrc);
    fillasp.SetFilter(&fsrc);
    fillsptheta.SetFilter(&fsrc);
    fillasptheta.SetFilter(&fasrc);

    MFillH fillontime ("EffOnTime  [MHTimeDiffTime]",  "MMcEvt");
    MFillH fillontheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");

    //
    // Setup Task list
    //
    tasklist.AddToList(&reader);
    tasklist.AddToList(&rand);
    tasklist.AddToList(&fillontime);
    tasklist.AddToList(&fillontheta);
    tasklist.AddToList(&acalc);
    tasklist.AddToList(&pcopy);
    tasklist.AddToList(&pnsb);
    tasklist.AddToList(&ncalc);
    tasklist.AddToList(&clean);
    tasklist.AddToList(&blind);
    tasklist.AddToList(&hcalc);
    tasklist.AddToList(&hsrc1);
    tasklist.AddToList(&hsrc2);
    tasklist.AddToList(&estim);
    tasklist.AddToList(&hfill1h);
    tasklist.AddToList(&hfill1m);
    tasklist.AddToList(&hfill2s);
    tasklist.AddToList(&hfill2a);
    tasklist.AddToList(&fsrc);
    tasklist.AddToList(&fasrc);
    tasklist.AddToList(&fillsp);
    tasklist.AddToList(&fillasptheta);
    tasklist.AddToList(&fillasp);
    tasklist.AddToList(&fillsptheta);
    tasklist.AddToList(&lvl1);
    tasklist.AddToList(&fethetaall);
    tasklist.AddToList(&fethetasel);
    tasklist.AddToList(&fetimeall);
    tasklist.AddToList(&fetimesel);
 
    //
    // set up the loop for the processing
    //
    MEvtLoop magic;
    magic.SetParList(&parlist);

    //
    // Start to loop over all events
    //
    if (!magic.Eventloop())
        return;

    tasklist.PrintStatistics();

    /*
     parlist.FindObject("HSource")->DrawClone();;
     parlist.FindObject("MHHillas")->DrawClone();
     parlist.FindObject("MHStarMap")->DrawClone();
     parlist.FindObject("HAntiSource")->DrawClone();
     parlist.FindObject("MHMcCollectionArea")->DrawClone();
     */

    /*
     MHEnergyTime  &alltime  = *(MHEnergyTime*)parlist.FindObject("AllTime");
     MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta");
     MHEnergyTime  &seltime  = *(MHEnergyTime*)parlist.FindObject("SelTime");
     MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta");

     MHEnergyTime collareatime;
     MHEnergyTheta collareatheta;
     collareatime.Divide(&seltime, &alltime);
     collareatheta.Divide(&seltheta, &alltheta);
     */

    MHTimeDiffTime  &effontime  = *(MHTimeDiffTime*)parlist.FindObject("EffOnTime");
    MHTimeDiffTheta &effontheta = *(MHTimeDiffTheta*)parlist.FindObject("EffOnTheta");

     /*
      effontime.DrawClone();
      effontheta.DrawClone();
     */

     MHEffOnTimeTime  ontime;
     MHEffOnTimeTheta ontheta;
     ontime.SetupFill(&parlist);
     ontheta.SetupFill(&parlist);

     ontime.Calc(effontime.GetHist());
     ontheta.Calc(effontheta.GetHist());

     /*
      ontime.DrawClone();
      ontheta.DrawClone();
      */

     MHAlphaEnergyTime  &fluxsp       = *(MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime");
     MHAlphaEnergyTime  &fluxasp      = *(MHAlphaEnergyTime*)parlist.FindObject("FluxASrcTime");
     MHAlphaEnergyTheta &fluxsptheta  = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxSrcTheta");
     MHAlphaEnergyTheta &fluxasptheta = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxASrcTheta");

     /*
      fluxsp.DrawClone();
      fluxasp.DrawClone();
      fluxsptheta.DrawClone();
      fluxasptheta.DrawClone();
      */

     MHAlphaEnergyTime  resulttime;
     MHAlphaEnergyTheta resulttheta;
     resulttime.Substract(&fluxsp, &fluxasp);
     resulttheta.Substract(&fluxsptheta, &fluxasptheta);

     /*
      resulttime.DrawClone();
      resulttheta.DrawClone();
      */

     TH2D &projecttime  = *resulttime.GetAlphaProjection(-10, 10);
     TH2D &projecttheta = *resulttheta.GetAlphaProjection(-10, 10);

     projecttime.SetTitle("Number of Gammas vs. EnergyEst and Time (Alpha integrated between -10, 10deg)");
     projecttheta.SetTitle("Number of Gammas vs. EnergyEst and Theta (Alpha integrated between -10, 10deg)");

     /*
      TCanvas *c = new TCanvas("Unfold", "To be unfolded", 350, 500);
      c->Divide(1, 2);
      c->cd(1);
      projecttime.DrawCopy();
      c->cd(2);
      projecttheta.DrawCopy();
      */

     return;

     for (int i=1; i<=binstime.GetNumBins(); i++)
     {
         if (ontime.GetHist()->GetBinContent(i)==0)
             continue;

         TH1D &hist = *projecttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i);

         /* UNFOLDING */

         //hist->Divide(collareatime);
         hist.Scale(1./ontime.GetHist()->GetBinContent(i));

         for (int j=1; j<=binse.GetNumBins(); j++)
             hist.SetBinContent(j, hist.GetBinContent(j)/hist.GetBinWidth(j));

         hist.SetName("Flux");
         hist.SetTitle("Flux[Gammas/s/m^2/GeV] vs. EnergyTrue for a fixed Time");

         char n[100];
         sprintf(n, "Canv%d", j);
         c= new TCanvas(n, "Title");
         hist.DrawCopy();
     }

     delete &projecttime;
     delete &projecttheta;

     return;

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

     MHMcCollectionArea carea;
     TH1D *collareatime  = carea.GetHist();  // FIXME!
     TH1D *collareatheta = carea.GetHist();  // FIXME!
}
