/* ======================================================================== *\
!
! *
! * 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                                    //
//    - dummy effective collection areas are used                           //
//                                                                          // 
//////////////////////////////////////////////////////////////////////////////
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 the binning for the histograms
    //

    //...............................................................
    // These are NOT the binnings for the flux determination

    MBinning binswidth("BinningWidth");
    binswidth.SetEdges(100, 0, 0.5);   // 100 bins from 0 to 0.5 deg

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

    MBinning binscamera("BinningCamera");
    binscamera.SetEdges(90, -2.25, 2.25);  // 90 bins from -2.25 to 2.25 deg

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

    MBinning binsalpha("BinningAlpha");
    binsalpha.SetEdges(100, -100, 100);

    MBinning binsasym("BinningAsym");
    binsasym.SetEdges(100, -0.5, 0.5);   // 100 bins from -0.5 to 0.5 deg

    MBinning binsasymn("BinningAsymn");
    binsasymn.SetEdges(100, -0.005, 0.005);   // 100 bins from -0.005 to 0.005 deg

    MBinning binsm3long("BinningM3Long");
    binsm3long.SetEdges(100, -0.5, 0.5);   // 100 bins from -0.5 to 0.5 deg

    MBinning binsm3trans("BinningM3Trans");
    binsm3trans.SetEdges(100, -0.5, 0.5);   // 100 bins from -0.5 to 0.5 deg

    //...............................................................
    // These ARE the binnings for the flux determination

    MBinning binsalphaflux("BinningAlphaFlux");
    binsalphaflux.SetEdges(100, 0, 100);

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

    parlist.AddToList(&binsdist);

    parlist.AddToList(&binsalpha);

    parlist.AddToList(&binsasym);
    parlist.AddToList(&binsasymn);
    parlist.AddToList(&binsm3long);
    parlist.AddToList(&binsm3trans);

    parlist.AddToList(&binsalphaflux);
    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");

    //reader.AddFile("/hd31/rkb/Camera/mpi/Gamma*.root");
    //reader.AddFile("/hd31/rkb/Camera/mpi/prot_N_*.root");
    //reader.AddFile("/hd31/rkb/Camera/mpi/prot_NS_*.root");
    //reader.AddFile("protondata/gamma_10_n*.root");
    reader.AddFile("protondata/gamma_10_cn*.root");

    reader.EnableBranch("MMcEvt.fTheta");

    //.....................................
    // filters
    MFTriggerLvl1 lvl1;

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

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

    //.....................................
    // source independent image parameters
    MHillasCalc       hcalc;
    MHillasExt        hext;
    parlist.AddToList(&hext);

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

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

    //.....................................
    // migration matrix (MC)
    MFillH eestetrue("MHMcEnergyMigration",  "HillasSrc");
    eestetrue.SetTitle("Task to determine the migration matrix for the energy");

    //.....................................
    // plots of source independent parameters (length, width, delta, size, center)
    MFillH hfill1h("MHHillas",  "MHillas");
    hfill1h.SetTitle("Task to plot the source independent image parameters");
    hfill1h.SetFilter(&lvl1);

    MFillH hfill1x("MHHillasExt", "MHillas");
    hfill1x.SetTitle("Task to plot the extended image parameters");
    hfill1x.SetFilter(&lvl1);

    //.....................................
    // plots of source dependent parameters (alpha, dist)
    MFillH hfill2s("HSource     [MHHillasSrc]", "HillasSrc");
    hfill2s.SetTitle("Task to plot the source dependent image parameters (Source)");
    hfill2s.SetFilter(&lvl1);
    MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
    hfill2a.SetTitle("Task to plot the source dependent image parameters (AntiSource)");
    hfill2a.SetFilter(&lvl1); 

    //.....................................
    // star map
    MFillH hfill1m("MHStarMap", "MHillas");
    hfill1m.SetTitle("Task to plot the Star map");

    //.....................................
    const Float_t alpha0 = 15; // [deg]
    MFAlpha fsrc ("HillasSrc",     '>', alpha0);
    fsrc.SetTitle("Task to check alphaSrc > alpha0");

    MFAlpha fasrc("HillasAntiSrc", '>', alpha0);
    fasrc.SetTitle("Task to check alphaAntiSrc > alpha0");

    //.....................................
    // 1-D profile plots     (<Theta>,  time)
    //                       (<Theta>, Theta)
    MFillH fillthetabartime ("ThetabarTime  [MHThetabarTime]",  "MMcEvt");
    fillthetabartime.SetTitle("Task to plot <Theta> vs time");
    // fillthetabartime.SetFilter(&lvl1);

    MFillH fillthetabartheta("ThetabarTheta [MHThetabarTheta]", "MMcEvt");
    fillthetabartheta.SetTitle("Task to plot <Theta> vs theta");
    // fillthetabartheta.SetFilter(&lvl1);

    //.....................................
    // 2-D plots     (Delta(time),  time)
    //               (Delta(time), Theta)
    MFillH filldtimetime ("EffOnTime  [MHTimeDiffTime]",  "MMcEvt");
    filldtimetime.SetTitle("Task to plot Delta(time) vs time");
    // filldtime.SetFilter(&lvl1);

    MFillH filldtimetheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");
    filldtimetheta.SetTitle("Task to plot Delta(time) vs theta");
    // fillontheta.SetFilter(&lvl1);

    //.....................................
    // 3-D plots     (alpha, E_est, time) 
    MFillH fillsptime  ("SrcTime   [MHAlphaEnergyTime]",  "HillasSrc");
    fillsptime.SetTitle("Task for 3D plots (alpha, E_est, time) (Source)");
    fillsptime.SetFilter(&fasrc);

    MFillH fillasptime ("ASrcTime  [MHAlphaEnergyTime]",  "HillasAntiSrc");
    fillasptime.SetTitle("Task for 3D plots (alpha, E_est, time) (AntiSource)");
    fillasptime.SetFilter(&fsrc);

    //.....................................
    // 3-D plots     (alpha, E_est, Theta) 
    MFillH fillsptheta ("SrcTheta  [MHAlphaEnergyTheta]", "HillasSrc");
    fillsptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (Source)");
    fillsptheta.SetFilter(&fasrc);

    MFillH fillasptheta("ASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");
    fillasptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (AntiSource)");
    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(&lvl1);
    tasklist.AddToList(&rand);

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

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

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

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



    // 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);

    
    MHHillas *mhhillas = (MHHillas*)parlist.FindObject("MHHillas");
    mhhillas->SetTitle("Source indep. image parameters");
    mhhillas->DrawClone();
    

    
    MHHillasExt *mhhillasext = (MHHillasExt*)parlist.FindObject("MHHillasExt");
    mhhillasext->SetTitle("Extended image parameters");
    mhhillasext->DrawClone();
    
    
    MHHillasSrc *hsource = (MHHillasSrc*)parlist.FindObject("HSource");
    hsource->SetTitle("Source dependent image parameters (Source)");
    hsource->DrawClone();

    MHHillasSrc *hantisource = (MHHillasSrc*)parlist.FindObject("HAntiSource");
    hantisource->SetTitle("Source dependent image parameters (AntiSource)");
    hantisource->DrawClone();
    

    /*    
    parlist.FindObject("MHStarMap")->DrawClone();
    */

    //----------------------------------------------------------------------
    // average Theta versus time
    //           and versus Theta
    // This is needed for selecting later the right collection area 
    // (at the right Theta) for each bin in time or Theta
    //
    MHThetabarTime  *bartime  = (MHThetabarTime*)parlist.FindObject("ThetabarTime");
    MHThetabarTheta *bartheta = (MHThetabarTheta*)parlist.FindObject("ThetabarTheta");
    
    /*    
    bartime->DrawClone();
    bartheta->DrawClone();
    */

     //----------------------------------------------------------------------
     // 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      

     char *varnametime; 
     varnametime = "Time";
     char *unittime;
     unittime    =  "[ s]";
     MHEffOnTime effontime(varnametime, unittime);

     char *varnametheta; 
     varnametheta = "Theta";
     char *unittheta;
     unittheta    =  "[  \\circ]";
     MHEffOnTime effontheta(varnametheta, unittheta);

     effontime.SetupFill(&parlist);
     effontheta.SetupFill(&parlist);


     // Draw == 0  don't draw the individual distributions of time differences
     //      != 0  draw them
     Bool_t Draw=kFALSE;
     effontime.Calc (dtimetime->GetHist(), Draw);
     effontheta.Calc(dtimetheta->GetHist(),Draw);


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


     //======================================================================
     //
     // A : fully differential plots (Alpha, E-est, Var)
     //
     //----------------------------------------------------------------------
     // 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("SrcTime");
     MHAlphaEnergyTime  *evtasptime  = (MHAlphaEnergyTime*)parlist.FindObject("ASrcTime");
     MHAlphaEnergyTheta *evtsptheta  = (MHAlphaEnergyTheta*)parlist.FindObject("SrcTheta");
     MHAlphaEnergyTheta *evtasptheta = (MHAlphaEnergyTheta*)parlist.FindObject("ASrcTheta");


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

      evtsptheta->SetTitle("Source : 3D-plot Al-En-Theta, alpha\_{asp} .gt. alpha\_0");
      evtsptheta->DrawClone();
      evtasptheta->SetTitle("AntiSource : 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 
     // 3rd and 4th argument are name and title of the resulting histogram

     MHGamma gamma;

     Bool_t Draw = kFALSE;
     TH3D *hsubtime  = gamma.Subtract( evtsptime->GetHist(), 
                                       evtasptime->GetHist(),
                        "Al-En-time", "3D-plot of Alpha,E-est,time", Draw);    

     TH3D *hsubtheta = gamma.Subtract( evtsptheta->GetHist(), 
                                       evtasptheta->GetHist(),
                        "Al-En-time", "3D-plot of Alpha,E-est,Theta", Draw);
    

     //----------------------------------------------------------------------
     // 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 =   0;    // [deg]
     Axis_t up =  10;    // [deg]
     Bool_t Draw = kFALSE;
     const TH2D &evttime  = *(gamma.GetAlphaProjection(hsubtime, 
                              lo, up, Draw));

     const TH2D &evttheta = *(gamma.GetAlphaProjection(hsubtheta, 
                              lo, up, Draw));
     

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


     /*
     MHMcEnergyMigration *migrationmatrix = 
         (MHMcEnergyMigration*)parlist.FindObject("MHMcEnergyMigration");
     migrationmatrix->SetTitle("Migration matrix");;
     migrationmatrix->DrawClone();
     */


     //----------------------------------------------------------------------
     // determine gamma flux from distributions of E-est :
     //   -  do unfolding to get the distributions in E-unfold
     //   -  divide by bin width in energy
     //                effective ontime and
     //                effective collection area
     //      to get the absolute gamma flux
    
     //.............................................
     // get flux spectrum for different bins in Time
     //
     // use dummy histogram *aeff which eventually will contain
     //     the effective collection area as a function of Etrue and Theta
     TH2D aeff;

     // arguments of MHFlux constructor are :
     //   - 2D histogram containing the no.of gammas as a function of
     //     E-est and Var
     //   - name  of variable Var ("Time" or "Theta")
     //   - units of variable Var ( "[s]" or "[\\circ]")

     // Draw == kTRUE  draw the no.of photons vs. E-est 
     //                for the individual bins of the variable Var
     Bool_t Draw=kTRUE;
     MHFlux fluxtime(evttime, Draw, "Time", "  [  s]");
     fluxtime.Unfold(Draw);
     fluxtime.CalcFlux(effontime.GetHist(), bartime.GetHist(),
                       &aeff, Draw);
     
     fluxtime.DrawClone();
     


     //..............................................
     // get flux spectrum for different bins in Theta
     Bool_t Draw=kTRUE;
     MHFlux fluxtheta(evttheta, Draw, "Theta", "  [  \\circ]");
     fluxtheta.Unfold(Draw);
     fluxtheta.CalcFlux(effontheta.GetHist(), bartheta.GetHist(),
                       &aeff, Draw);
     
     fluxtheta.DrawClone();
     

return;
     //----------------------------------------------------------------------


     //======================================================================
     //
     // B : plots integrated over E-est
     //
     //----------------------------------------------------------------------
     // integrate 3D-plot of (alpha, E_est and Time) over E-est
     //    to get 2D-plot of (alpha, Time) 
                                   
     Draw = kFALSE;                        
     TH2D *evttimesp  = evtsptime.IntegrateEest ("AlphaSP vs. Time", Draw);
     TH2D *evttimeasp = evtasptime.IntegrateEest("AlphaASP vs. Time",Draw);


     //======================================================================
     //
     // C : plots integrated over Time
     //
     //----------------------------------------------------------------------
     // integrate 3D-plot of (alpha, E_est and Time) over the Time
     //    to get 2D-plot of (alpha, E_est) 
                                   
     Draw = kFALSE;                        
     TH2D *evteestsp  = evtsptime.IntegrateTime ("AlphaSP vs. E-est", Draw);
     TH2D *evteestasp = evtasptime.IntegrateTime("AlphaASP vs. E-est",Draw);



     //======================================================================
     //
     // D : plots integrated over E-est and Time
     //

     //----------------------------------------------------------------------
     // integrate 3D-plot of (alpha, E_est and Time) over E-est and Time
     //    to get 1D-plot of (alpha) 
                                   
     Draw = kFALSE;                        
     TH1D *evtalphasp  = evtsptime.IntegrateEestTime ("AlphaSP", Draw);
     TH1D *evtalphaasp = evtasptime.IntegrateEestTime("AlphaASP",Draw);

     //----------------------------------------------------------------------
     // list of tasks for processing

     cout << " " << endl;
     cout << "List of tasks for processing :" << endl;
     cout << "==============================" << endl;
     cout << " " << endl;

     // was folgt geht nicht :

     // Error: Function Next() is not defined in current scope  FILE:macros/flux.C LINE:636

     //while ( (task=(MTask*)Next()) )
     //{
     //  cout << tasklist.GetName() << " :   " << tasklist.GetTitle() << endl;
     //}

     return;

     //======================================================================                
}
//===========================================================================


