/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//                                                                          //
// This macro calculates photon fluxes as a function of                     //
//    - energy and time                                                     //
//    - energy and Theta                                                    //
//                                                                          //
// It is assumed that the data have been taken in the wobble mode.          //
// This means that there is only one set of data, from which both           //
// 'on' and 'off' data are constructed.                                     //
//                                                                          //
// Necessary input from MC :                                                //
//    - migration matrix (E_est, E_true) as a functioin of Theta            //
//    - effective collection areas as a function of energy and Theta        //
//                                                                          //
//                                                                          // 
// The input from MC has to refer to the wobble mode too.                   // 
//                                                                          //
// The macro includes :                                                     // 
//    - the calculation of Hillas parameters                                //
//    - the calculation of the effective on time                            //
//    - the unfolding of the distributions in the estimated energy E_est    //
//                                                                          //
// For testing purposes (as long as no real data, no energy estimator,      //
// no migration matrices and no effective collection areas are available)   //
//    - event times are generated (according to a Poissonian with dead time)//
//    - a dummy energy estimator is provided                                //
//    - migration matrices are generated                                    //
//    - a kind of effective collection areas are calculated                 //
//                                                                          // 
//////////////////////////////////////////////////////////////////////////////
void flux()
{ 
    //--------------------------------------------------------------
    // empty lists of parameter containers and tasks
    //
    MParList  parlist;
    MTaskList tasklist;

    //--------------------------------------------------------------
    // Geometry information of the MAGIC camera
    //
    MGeomCamMagic geomcam;

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

    //--------------------------------------------------------------
    // Setup binning for the histograms
    //
    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(100, 0, 1);   // 100 bins from 0 to 1 deg

    MBinning binslength("BinningLength");
    binslength.SetEdges(100, 0, 1);  // 100 bins from 0 to 1 deg

    MBinning binsdist("BinningDist");
    binsdist.SetEdges(100, 0, 2);    // 100 bins from 0 to 2 deg

    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, 100);

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


    //--------------------------------------------------------------
    // Fill list of parameter containers
    //

    parlist.AddToList(&tasklist);
    parlist.AddToList(&geomcam);

    parlist.AddToList(&source);
    parlist.AddToList(&antisrc);

    parlist.AddToList(&binswidth);
    parlist.AddToList(&binslength);
    parlist.AddToList(&binsdist);

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

    //--------------------------------------------------------------
    // Setup the tasks
    //

    //.....................................
    // input file
    //
    // next statement commented out
    //    MReadMarsFile reader("Events", "~/data/Gamma*.root");

  //MReadMarsFile reader("Events", "/hd31/rkb/Camera/mpi/Gamma*.root");
  //reader.AddFile("/hd31/rkb/Camera/mpi/prot_N_*.root");

    // next 2 statements added
    MReadMarsFile reader("Events");
    reader.AddFile("protondata/gamma_10_cn*.root");

    reader.EnableBranch("MMcEvt.fTheta");

    //.....................................
    // generation of event time
    MMcTimeGenerate   rand;

    //.....................................
    // collection area
    MMcCollectionAreaCalc acalc;

    MMcPedestalCopy   pcopy;
    MMcPedestalNSBAdd pnsb;
    MCerPhotCalc      ncalc;
    MImgCleanStd      clean;
    MBlindPixelCalc   blind;

    //.....................................
    // source independent image parameters
    MHillasCalc       hcalc;

    //.....................................
    // source dependent image parameters
    MHillasSrcCalc    hsrc1("Source",     "HillasSrc");
    MHillasSrcCalc    hsrc2("AntiSource", "HillasAntiSrc");

    //.....................................
    // energy estimation
    MEnergyEstimate   estim;

    //.....................................
    // migration matrix (MC)
    MFillH eestetrue("MHMcEnergyMigration",  "HillasSrc");


    //.....................................
    // plots of width and length
    MFillH hfill1h("MHHillas",  "MHillas");

    //.....................................
    // star map
    MFillH hfill1m("MHStarMap", "MHillas");

    //.....................................
    // plots of alpha and dist
    MFillH hfill2s("HSource     [MHHillasSrc]", "HillasSrc");
    MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
 
    //.....................................
    // filters
    MFTriggerLvl1 lvl1;

    const Float_t alpha0 = 15; // [deg]
    MFAlpha fsrc ("HillasSrc",     '>', alpha0);
    MFAlpha fasrc("HillasAntiSrc", '>', alpha0);

    //.....................................
    // 1-D profile plots     (<Theta>,  time)
    //                       (<Theta>, Theta)
    MFillH fillthetabartime ("ThetabarTime  [MHThetabarTime]",  "MMcEvt");
    // fillthetabartime.SetFilter(&lvl1);
    MFillH fillthetabartheta("ThetabarTheta [MHThetabarTheta]", "MMcEvt");
    // fillthetabartheta.SetFilter(&lvl1);

    //.....................................
    // 2-D plots     (Delta(time),  time)
    //               (Delta(time), Theta)
    MFillH filldtimetime ("EffOnTime  [MHTimeDiffTime]",  "MMcEvt");
    // filldtime.SetFilter(&lvl1);
    MFillH filldtimetheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");
    // fillontheta.SetFilter(&lvl1);

    //.....................................
    // 3-D plots     (alpha, E_est, time) 
    MFillH fillsptime  ("FluxSrcTime   [MHAlphaEnergyTime]",  "HillasSrc");
    fillsptime.SetFilter(&fasrc);
    MFillH fillasptime ("FluxASrcTime  [MHAlphaEnergyTime]",  "HillasAntiSrc");
    fillasptime.SetFilter(&fsrc);

    //.....................................
    // 3-D plots     (alpha, E_est, Theta) 
    MFillH fillsptheta ("FluxSrcTheta  [MHAlphaEnergyTheta]", "HillasSrc");
    fillsptheta.SetFilter(&fasrc);
    MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");
    fillasptheta.SetFilter(&fsrc);

    //.....................................
    // MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt");
    // MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt");
    // fethetasel.SetFilter(&lvl1);

    // MFillH fetimeall ("AllTime  [MHEnergyTime]",  "MMcEvt");
    // MFillH fetimesel ("SelTime  [MHEnergyTime]",  "MMcEvt");
    // fetimesel.SetFilter(&lvl1);


    //---------------------------------------------------------------------
    // Setup the task list
    //
    tasklist.AddToList(&reader);
    tasklist.AddToList(&rand);

    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(&eestetrue);

    tasklist.AddToList(&hfill1h);
    tasklist.AddToList(&hfill1m);
    tasklist.AddToList(&hfill2s);
    tasklist.AddToList(&hfill2a);

    tasklist.AddToList(&fillthetabartime);
    tasklist.AddToList(&fillthetabartheta);

    tasklist.AddToList(&filldtimetime);
    tasklist.AddToList(&filldtimetheta);

    // tasklist.AddToList(&lvl1);

    // tasklist.AddToList(&fethetaall);
    // tasklist.AddToList(&fethetasel);
    // tasklist.AddToList(&fetimeall);
    // tasklist.AddToList(&fetimesel);

    tasklist.AddToList(&fsrc);
    tasklist.AddToList(&fasrc);

    tasklist.AddToList(&fillsptime);
    tasklist.AddToList(&fillasptime);

    tasklist.AddToList(&fillsptheta);
    tasklist.AddToList(&fillasptheta);

 
    //----------------------------------------------------------------------
    // Event loop
    //
    MEvtLoop magic;
    magic.SetParList(&parlist);

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


    tasklist.PrintStatistics(10);

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

    //----------------------------------------------------------------------
    // average Theta versus time
    //           and versus Theta
    //
    MHThetabarTime  *bartime  = (MHThetabarTime*)parlist.FindObject("ThetabarTime");
    MHThetabarTheta *bartheta = (MHThetabarTheta*)parlist.FindObject("ThetabarTheta");
    
    /*
    bartime->DrawClone();
    bartheta->DrawClone();
    */


    //----------------------------------------------------------------------
    // this is something like the collection area
    //         as a function of time
    //     and as a function of Theta
     /*
     MHEnergyTime  &alltime  = *(MHEnergyTime*)parlist.FindObject("AllTime");
     MHEnergyTime  &seltime  = *(MHEnergyTime*)parlist.FindObject("SelTime");
     MHEnergyTime collareatime;
     collareatime.Divide(&seltime, &alltime);

     MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta");
     MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta");
     MHEnergyTheta collareatheta;
     collareatheta.Divide(&seltheta, &alltheta);
     */

     //----------------------------------------------------------------------
     // Effective on time

     //....................................
     // get plots of time differences for different intervals in time
     // and plots of time differences for different intervals in Theta 

     MHTimeDiffTime  *dtimetime  = (MHTimeDiffTime*)parlist.FindObject("EffOnTime");
     MHTimeDiffTheta *dtimetheta = (MHTimeDiffTheta*)parlist.FindObject("EffOnTheta");
     /*
     dtimetime->DrawClone();
     dtimetheta->DrawClone();
     */


     //....................................
     // calculate effective on time for different intervals in time
     //                         and for different intervals in Theta      
     MHEffOnTimeTime  effontime;
     MHEffOnTimeTheta effontheta;
     effontime.SetupFill(&parlist);
     effontheta.SetupFill(&parlist);

     effontime.Calc(dtimetime->GetHist());
     effontheta.Calc(dtimetheta->GetHist());

     // plot effective on time versus time
     //                    and versus Theta
     
     effontime.DrawClone();
     effontheta.DrawClone();
     


     //----------------------------------------------------------------------
     // 3D-plots of alpha, E_est and time
     //      and of alpha, E_est and Theta 
     // with alpha calculated relative to the source position
     //                   and relative to the anti-source position

     MHAlphaEnergyTime  *evtsptime   = (MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime");
     MHAlphaEnergyTime  *evtasptime  = (MHAlphaEnergyTime*)parlist.FindObject("FluxASrcTime");
     MHAlphaEnergyTheta *evtsptheta  = (MHAlphaEnergyTheta*)parlist.FindObject("FluxSrcTheta");
     MHAlphaEnergyTheta *evtasptheta = (MHAlphaEnergyTheta*)parlist.FindObject("FluxASrcTheta");

     /* this test shows that the addresses
             'evtsptime' and  '&evtspobj' are identical
     MHAlphaEnergyTime  &evtspobj     = *(MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime");
     cout << "evtsptime = " << evtsptime << endl;
     cout << "&evtspobj = " << &evtspobj << endl;
     */

     /*
      evtsptime->SetTitle("3D-plot Al-En-time, alpha\_{asp} .gt. alpha\_0");
      evtsptime->DrawClone();
      evtasptime->SetTitle("3D-plot Al-En-time, alpha\_{sp} .gt. alpha\_0");
      evtasptime->DrawClone();

      evtsptheta->SetTitle("3D-plot Al-En-Theta, alpha\_{asp} .gt. alpha\_0");
      evtsptheta->DrawClone();
      evtasptheta->SetTitle("3D-plot Al-En-Theta, alpha\_{sp} .gt. alpha\_0");
      evtasptheta->DrawClone();
     */                


     //----------------------------------------------------------------------
     // Difference between source and antisource (= 'gamma' sample) 
     // for             3D-plots of alpha, E_est and time
     //                      and of alpha, E_est and Theta 
      
     MHAlphaEnergyTime  gammatime("Al-En-time","3D-plot Al-En-time");
     MHAlphaEnergyTheta gammatheta("Al-En-Theta","3D-plot Al-En-Theta");

     gammatime.SetupFill(&parlist);
     gammatime.Subtract(evtsptime, evtasptime);
     gammatheta.SetupFill(&parlist);
     gammatheta.Subtract(evtsptheta, evtasptheta);


     /*
      gammatime.SetTitle("3D-plot Al-En-time for \'gamma\' sample");
      gammatime.DrawClone();

      gammatheta.SetTitle("3D-plot Al-En-Theta for \'gamma\' sample");
      gammatheta.DrawClone();
     */                

     //----------------------------------------------------------------------
     // get number of 'gammas' in the alpha range (lo, up) as a function of
     //                 E_est and time
     //        or       E_est and Theta 

     Axis_t lo = -10;
     Axis_t up =  10;

     //     TH2D &evttime  = *gammatime.GetAlphaProjection (lo, up);
     //     TH2D &evttheta = *gammatheta.GetAlphaProjection(lo, up);

     TH2D *evttime  = gammatime.DrawAlphaProjection (lo, up);
     TH2D *evttheta = gammatheta.DrawAlphaProjection(lo, up);


     //----------------------------------------------------------------------
     // plot migration matrix:  E-true E-est Theta

     parlist.FindObject("MHMcEnergyMigration")->DrawClone();

     return;

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

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

         /* UNFOLDING */

         //hist->Divide(collareatime);
         hist.Scale(1./effontime.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 &evttime;
     delete &evttheta;

     return;

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

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


