/* ======================================================================== *\
!
! *
! * 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): Eva Domingo,     12/2004 <mailto:domingo@ifae.es>
!              Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MFindDisp                                                               //
//                                                                         //
// Class for otimizing the estimation of Disp                              //
//                                                                         //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MFindDisp.h"

#include <math.h>            // fabs 

#include <TArrayD.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TMinuit.h>
#include <TStopwatch.h>
#include <TVirtualFitter.h>

#include "MBinning.h"
#include "MContinue.h"
#include "MDispParameters.h"
#include "MDispCalc.h"
#include "MDataElement.h"
#include "MDataMember.h"

#include "MEvtLoop.h"
#include "MFSelFinal.h"
#include "MF.h"
#include "MFDisp.h"
#include "MFEventSelector.h"
#include "MFEventSelector2.h"
#include "MFillH.h"
#include "MFEventSelector.h"
#include "MGeomCamMagic.h"
#include "MH3.h"
#include "MHDisp.h"
#include "MHFindSignificance.h"
#include "MHMatrix.h"
#include "MHOnSubtraction.h"

#include "MLog.h"
#include "MLogManip.h"
#include "MMatrixLoop.h"
#include "MMinuitInterface.h"
#include "MParameters.h"
#include "MParList.h"
#include "MProgressBar.h"
#include "MReadMarsFile.h"
#include "MReadTree.h"
#include "MStatusDisplay.h"
#include "MTaskList.h"


ClassImp(MFindDisp);

using namespace std;


//------------------------------------------------------------------------
// fcnDisp 
//------------------------------------------------------------------------
//
// - calculates the quantity to be minimized (by TMinuit)
//
// - the quantity to be minimized is defined in 
//   MHDisp::Fill() and MHDisp::Finalize()
//
// - the parameters to be varied in the minimization are the parameters
//   appearing in the parametrization of Disp (MDispCalc::Calc())
//
//------------------------------------------------------------------------

static void fcnDisp(Int_t &npar, Double_t *gin, Double_t &f, 
                    Double_t *par, Int_t iflag)
{
    cout <<  "entry fcnDisp" << endl;

    // save pointer to the MINUIT object (for optimizing Disp) for the case
    // that it is overwritten by the pointer to another Minuit object
    //    TMinuit *savePointer = gMinuit;
    //-------------------------------------------------------------


    MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();

    MParList *plistfcn   = (MParList*)evtloopfcn->GetParList();
    MTaskList *tasklistfcn   = (MTaskList*)plistfcn->FindObject("MTaskList");

    // get needed Disp containers from the existing parList
    MDispParameters *dispparams = (MDispParameters*)plistfcn->FindObject("MDispParameters");
    if (!dispparams)
    {
        gLog << "fcnDisp : MDispParameters object '" << "MDispParameters"
            << "' not found... aborting" << endl;
        return;
    }

    MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp");
    if (!hdisp)
    {
        gLog << "fcnDisp : MHDisp object '" << "MHDisp"
            << "' not found... aborting" << endl;
        return;
    }

    MParameterD *minpar = (MParameterD*)plistfcn->FindObject("MinimizationParameter");

    //
    // transfer current parameter values to MDispParameters
    //
    // Attention : npar is the number of variable parameters
    //                  not the total number of parameters
    //
    Double_t fMin,   fEdm,   fErrdef;
    Int_t    fNpari, fNparx, fIstat;
    gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);

    dispparams->SetParameters(TArrayD(fNparx, par));

    // checking that parameters have been properly set
    TArrayD checkparameters = dispparams->GetParameters();
    gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ",  " 
         << fNparx  << endl;
    gLog << "fcnDisp : i, par, checkparameters =" << endl;
    for (Int_t i=0; i<fNparx; i++)
    {
      gLog << i << ",  " << par[i] << ",  " << checkparameters[i] << endl;
    }

    //
    // loop over the events in the training matrices to compute the parameter 
    // to minimize for the current Disp parameter values
    evtloopfcn->Eventloop();

    tasklistfcn->PrintStatistics(0, kTRUE);

    //-------------------------------------------
    // get the Minimization parameter value for the current values 
    // of the Disp parameters
    if (minpar)
      f = minpar->GetVal();
    else
    {
      gLog << "fcnDisp : MParameterD object '" << "MinimizationParameter"
	   << "' not found... aborting" << endl;
      return;
    }
}


// --------------------------------------------------------------------------
//
// Default constructor.
//
MFindDisp::MFindDisp(MFDisp *fdisp, const char *name, const char *title)
  : fHowManyTrain(10000), fHowManyTest(10000)
{
    fName  = name  ? name  : "MFindDisp";
    fTitle = title ? title : "Optimizer of Disp";

    //---------------------------
    // camera geometry is needed for conversion mm ==> degree
    fCam = new MGeomCamMagic; 

    // initialize MFDisp filter cuts to default (no cutting) values
    // (done at MFDisp constructor)
    fDispFilter = fdisp;

    // matrices to contain the training/test samples
    fMatrixTrainCalc  = new MHMatrix("MatrixTrainCalc");
    fMatrixTrainHists = new MHMatrix("MatrixTrainHists");
    fMatrixTestCalc   = new MHMatrix("MatrixTestCalc");
    fMatrixTestHists  = new MHMatrix("MatrixTestHists");

    // objects of MDispCalc where the first part of the matrices mapping is defined
    fDispCalcTrain = new MDispCalc("DispTrain","MDispParameters");
    fDispCalcTest  = new MDispCalc("DispTest","MDispParametersTest");

    // objects of MHDisp where the second part of the matrices mapping is defined
    fHDispTrain = new MHDisp("DispTrain");
    fHDispTest  = new MHDisp("DispTest");

    // define columns of matrices
    // Train matrix
    fDispCalcTrain->InitMapping(fMatrixTrainCalc);
    fHDispTrain->InitMapping(fMatrixTrainHists);
    // Test matrix
    fDispCalcTest->InitMapping(fMatrixTestCalc);
    fHDispTest->InitMapping(fMatrixTestHists);
}


// --------------------------------------------------------------------------
//
// Default destructor.
//
MFindDisp::~MFindDisp()
{
    delete fCam;
    delete fMatrixTrainCalc;
    delete fMatrixTrainHists;
    delete fMatrixTestCalc;
    delete fMatrixTestHists;
    delete fDispCalcTrain;
    delete fDispCalcTest;
    delete fHDispTrain;
    delete fHDispTest;
}


// --------------------------------------------------------------------------
//
// Define the matrices 'fMatrixTrainCalc' and 'fMatrixTrainHists' 
// for the TRAINING sample
//
// alltogether 'howmanytrain' events are read from file 'nametrain';
// the events are selected according to a target distribution 'hreftrain'
//
//
Bool_t MFindDisp::DefineTrainMatrix(
			  const TString &nametrain, MH3 &hreftrain, 
			  const Int_t howmanytrain, const TString &filetrain, 
			  Int_t iflag)
{
    if (nametrain.IsNull() || howmanytrain <= 0)
        return kFALSE;

    *fLog   << "=============================================" << endl;
    *fLog   << "Fill TRAINING Matrices from file '" << nametrain << endl;
    *fLog   << "     select " << howmanytrain << " events " << endl;
    if (!hreftrain.GetHist().GetEntries()==0)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << hreftrain.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }
    *fLog   << "=============================================" << endl;


    MParList  plist;
    MTaskList tlist;

    // initialize display to check the reference distribution
    // for the event selection (just if iflag==1)
    MStatusDisplay *display = NULL;
    if (iflag)
      display = new MStatusDisplay;

    MReadMarsFile read("Events", nametrain);
    read.DisableAutoScheme();

    // apply cuts for event selection
    MContinue contdisp(fDispFilter);
    contdisp.SetName("ContFilterSelector2");

    // choose a reference distribution
    MFEventSelector2 seltrain(hreftrain);
    seltrain.SetNumMax(howmanytrain);
    seltrain.SetName("selectTrain");

    MFillH filltraincalc(fMatrixTrainCalc);
    filltraincalc.SetFilter(&seltrain);
    filltraincalc.SetName("fillMatrixTrainCalc");

    MFillH filltrainhists(fMatrixTrainHists);
    filltrainhists.SetFilter(&seltrain);
    filltrainhists.SetName("fillMatrixTrainHists");

    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrainCalc);
    plist.AddToList(fMatrixTrainHists);

    //******************************
    // entries in MTaskList 

    tlist.AddToList(&read);
    if (fDispFilter != NULL)
      tlist.AddToList(&contdisp);
    tlist.AddToList(&seltrain);
    tlist.AddToList(&filltraincalc);
    tlist.AddToList(&filltrainhists);

    //******************************

    MProgressBar bar;
    MEvtLoop evtloop;
    if (display != NULL)
      evtloop.SetDisplay(display);
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrain");
    evtloop.SetProgressBar(&bar);

    if (!evtloop.Eventloop())
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);


    // print the filled Training Matrices
    fMatrixTrainCalc->Print("SizeCols");
    fMatrixTrainHists->Print("SizeCols");

    // check that number of generated events is compatible with the resquested
    Int_t howmanygeneratedcalc = fMatrixTrainCalc->GetM().GetNrows();
    if (TMath::Abs(howmanygeneratedcalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events ("
	    << howmanygeneratedcalc 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    Int_t howmanygeneratedhists = fMatrixTrainHists->GetM().GetNrows();
    if (TMath::Abs(howmanygeneratedhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events ("
	    << howmanygeneratedhists 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }


    *fLog << "TRAINING Matrices were filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out training matrices

    if (filetrain != "")
    {
      TFile filetr(filetrain, "RECREATE");
      fMatrixTrainCalc->Write();
      fMatrixTrainHists->Write();
      filetr.Close();

      *fLog << "MFindDisp::DefineTrainMatrix; Training matrices were written onto file '"
            << filetrain << "'" << endl;
    }


    if (display != NULL)
      delete display;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Define the matrices 'fMatrixTestCalc' and 'fMatrixTestHists' 
// for the TEST sample
//
// alltogether 'howmanytest' events are read from file 'nametest'
// the events are selected according to a target distribution 'hreftest'
//
//
Bool_t MFindDisp::DefineTestMatrix(
			  const TString &nametest, MH3 &hreftest,
	                  const Int_t howmanytest, const TString &filetest,
			  Int_t iflag)
{
    if (nametest.IsNull() || howmanytest<=0)
        return kFALSE;

    *fLog   << "=============================================" << endl;
    *fLog   << "Fill TEST Matrices from file '" << nametest << endl;
    *fLog   << "     select " << howmanytest << " events " << endl;
    if (!hreftest.GetHist().GetEntries()==0)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << hreftest.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }


    MParList  plist;
    MTaskList tlist;

    // initialize display to check the reference distribution
    // for the event selection (just if iflag==1)
    MStatusDisplay *display = NULL;
    if (iflag)
      display = new MStatusDisplay;

    MReadMarsFile read("Events", nametest);
    read.DisableAutoScheme();

    // apply cuts for event selection
    MContinue contdisp(fDispFilter);
    contdisp.SetName("ContFilterSelector2");
 
    // choose a reference distribution 
    MFEventSelector2 seltest(hreftest);
    seltest.SetNumMax(howmanytest);
    seltest.SetName("selectTest");
 
    MFillH filltestcalc(fMatrixTestCalc);
    filltestcalc.SetFilter(&seltest);
    filltestcalc.SetName("fillMatrixTestCalc");

    MFillH filltesthists(fMatrixTestHists);
    filltesthists.SetFilter(&seltest);
    filltesthists.SetName("fillMatrixTestHists");

    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTestCalc);
    plist.AddToList(fMatrixTestHists);

    //******************************
    // entries in MTaskList 

    tlist.AddToList(&read);
    if (fDispFilter != NULL)
      tlist.AddToList(&contdisp);
    tlist.AddToList(&seltest);
    tlist.AddToList(&filltestcalc);
    tlist.AddToList(&filltesthists);

    //******************************

    MProgressBar bar;
    MEvtLoop evtloop;
    if (display != NULL)
      evtloop.SetDisplay(display);
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTest");
    evtloop.SetProgressBar(&bar);

    if (!evtloop.Eventloop())
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);


    // print the filled Test Matrices
    fMatrixTestCalc->Print("SizeCols");
    fMatrixTestHists->Print("SizeCols");

    // check that number of generated events is compatible with the resquested
    const Int_t howmanygeneratedcalc = fMatrixTestCalc->GetM().GetNrows();
    if (TMath::Abs(howmanygeneratedcalc-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindDisp::DefineTestMatrix; no.of generated events ("
	    << howmanygeneratedcalc
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }

    const Int_t howmanygeneratedhists = fMatrixTestHists->GetM().GetNrows();
    if (TMath::Abs(howmanygeneratedhists-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindDisp::DefineTestMatrix; no.of generated events ("
	    << howmanygeneratedhists
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }

    *fLog << "TEST Matrices were filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out test matrices

    if (filetest != "")
    {
      TFile filete(filetest, "RECREATE");
      fMatrixTestCalc->Write();
      fMatrixTestHists->Write();
      filete.Close();

      *fLog << "MFindDisp::DefineTestMatrix; Test matrices were written onto file '"
            << filetest << "'" << endl;
    }


    if (display != NULL)
      delete display;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Define the matrices for the TRAINING and TEST sample respectively
//
//   - this is for the case that training and test samples are generated
//     from the same input file (which has the name 'name')
//
Bool_t MFindDisp::DefineTrainTestMatrix(
			  const TString &name, MH3 &href,
	                  const Int_t howmanytrain, const Int_t howmanytest,
                          const TString &filetrain, const TString &filetest,
			  Int_t iflag)
{
    *fLog   << "=============================================" << endl;
    *fLog   << "Fill TRAINING and TEST Matrices from file '" << name << endl;
    *fLog   << "     select "   << howmanytrain 
	    << " training and " << howmanytest << " test events " << endl;
    if (!href.GetHist().GetEntries()==0)
    {
      *fLog << "     according to a distribution given by the MH3 object '"
            << href.GetName() << "'" << endl;
    }
    else
    {
      *fLog << "     randomly" << endl;
    }


    MParList  plist;
    MTaskList tlist;

    // initialize display to check the reference distribution
    // for the event selection (just if iflag==1)
    MStatusDisplay *display = NULL;
    if (iflag)
      display = new MStatusDisplay;

    MReadMarsFile read("Events", name);
    read.DisableAutoScheme();

    // apply cuts for event selection
    MContinue contdisp(fDispFilter);
    contdisp.SetName("ContFilterSelector2");

    // choose a reference distribution 
    MFEventSelector2 selector(href);
    selector.SetNumMax(howmanytrain+howmanytest);
    selector.SetName("selectTrainTest");
    selector.SetInverted();
    MContinue cont(&selector);
    cont.SetName("ContTrainTest");

    // choose randomly the events to fill the Training Matrix
    Double_t prob =  ( (Double_t) howmanytrain )
                   / ( (Double_t)(howmanytrain+howmanytest) );
    MFEventSelector split;
    split.SetSelectionRatio(prob);

    MFillH filltraincalc(fMatrixTrainCalc);
    filltraincalc.SetFilter(&split);
    filltraincalc.SetName("fillMatrixTrainCalc");

    MFillH filltrainhists(fMatrixTrainHists);
    filltrainhists.SetFilter(&split);
    filltrainhists.SetName("fillMatrixTrainHists");

    // consider this event as candidate for a test event 
    // only if event was not accepted as a training event
    MContinue conttrain(&split);
    conttrain.SetName("ContTrain");

    MFillH filltestcalc(fMatrixTestCalc);
    filltestcalc.SetName("fillMatrixTestCalc");

    MFillH filltesthists(fMatrixTestHists);
    filltesthists.SetName("fillMatrixTestHists");


    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrainCalc);
    plist.AddToList(fMatrixTrainHists);
    plist.AddToList(fMatrixTestCalc);
    plist.AddToList(fMatrixTestHists);

    //******************************
    // entries in MTaskList 

    tlist.AddToList(&read);
    if (fDispFilter != NULL)
      tlist.AddToList(&contdisp);
    tlist.AddToList(&cont);
    tlist.AddToList(&filltraincalc);
    tlist.AddToList(&filltrainhists);
    tlist.AddToList(&conttrain);
    tlist.AddToList(&filltestcalc);
    tlist.AddToList(&filltesthists);

    //******************************

    MProgressBar bar;
    MEvtLoop evtloop;
    if (display != NULL)
      evtloop.SetDisplay(display);
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrainTest");
    evtloop.SetProgressBar(&bar);

    Int_t maxev = -1;
    if (!evtloop.Eventloop(maxev))
      return kFALSE;

    tlist.PrintStatistics(0, kTRUE);


    // print the filled Train Matrices
    fMatrixTrainCalc->Print("SizeCols");
    fMatrixTrainHists->Print("SizeCols");

    // check that number of generated events is compatible with the resquested
    const Int_t generatedtraincalc = fMatrixTrainCalc->GetM().GetNrows();
    if (TMath::Abs(generatedtraincalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtraincalc
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }
    const Int_t generatedtrainhists = fMatrixTrainHists->GetM().GetNrows();
    if (TMath::Abs(generatedtrainhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtrainhists
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }


    // print the filled Test Matrices
    fMatrixTestCalc->Print("SizeCols");
    fMatrixTestHists->Print("SizeCols");

    // check that number of generated events is compatible with the resquested
    const Int_t generatedtestcalc = fMatrixTestCalc->GetM().GetNrows();
    if (TMath::Abs(generatedtestcalc-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtestcalc
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }
    const Int_t generatedtesthists = fMatrixTestHists->GetM().GetNrows();
    if (TMath::Abs(generatedtesthists-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtesthists
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }


    *fLog << "TRAINING and TEST Matrices were filled" << endl;
    *fLog << "=============================================" << endl;


    //----------------------------
    // write out training matrices

    if (filetrain != "")
    {
      TFile filetr(filetrain, "RECREATE");
      fMatrixTrainCalc->Write();
      fMatrixTrainHists->Write();
      filetr.Close();

      *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrices were written onto file '"
            << filetrain << "'" << endl;
    }

    //--------------------------
    // write out test matrices

    if (filetest != "")
    {
      TFile filete(filetest, "RECREATE");
      fMatrixTestCalc->Write();
      fMatrixTestHists->Write();
      filete.Close();

      *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrices were written onto file '"
            << filetest << "'" << endl;
    }

    if (display != NULL)
      delete display;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Read training and test matrices from files
//
//
Bool_t MFindDisp::ReadMatrix(const TString &filetrain, const TString &filetest)
{
  //--------------------------
  // read in training matrices

  TFile filetr(filetrain);
  fMatrixTrainCalc->Read("MatrixTrainCalc");
  fMatrixTrainHists->Read("MatrixTrainHists");
  fMatrixTrainCalc->Print("SizeCols");
  fMatrixTrainHists->Print("SizeCols");

  *fLog << "MFindDisp::ReadMatrix; Training matrices were read in from file '"
        << filetrain << "'" << endl;
  filetr.Close();


  //--------------------------
  // read in test matrices

  TFile filete(filetest);
  fMatrixTestCalc->Read("MatrixTestCalc");
  fMatrixTestHists->Read("MatrixTestHists");
  fMatrixTestCalc->Print("SizeCols");
  fMatrixTestHists->Print("SizeCols");

  *fLog << "MFindDisp::ReadMatrix; Test matrices were read in from file '"
        << filetest << "'" << endl;
  filete.Close();

  return kTRUE;  
}


//------------------------------------------------------------------------
//
// Steering program for optimizing Disp
// ------------------------------------
//
//      the criterion for the 'optimum' is defined in 
//      MHDisp::Fill()  and  MHDisp::Finalize(); 
//      for example : smallest sum (over all events) of d^2, where d is the 
//                    distance between the estimated source position 
//                    (calculated using the current value of Disp) and
//                    the true source position     
//
// The various steps are :
//
// - setup the event loop to be executed for each call to fcnDisp 
//
// - call TMinuit to do the minimization :
//        the fcnDisp function calculates the parameter to minimize 
//                            for the current Disp parameter values;
//        for this - Disp is calculated in the event loop by calling 
//                   MDispCalc::Process() ==> MDispCalc::Calc() 
//                 - the Minimization parameter contributions are summed up 
//                   in the event loop by calling MHDisp::Fill()
//                 - after the event loop the final value of the Minimization
//                   parameter is calculated by calling MHDisp::Finalize()
//
// Needed as input : (to be set by the Set functions)
//
// - fFilenameParam      name of file to which optimum values of the 
//                       parameters are written
//
// - for the minimization, the starting values of the parameters are taken  
//     - from the file parDispInit (if it is != "")
//     - or from the arrays params and/or steps 
//     - or from the DispParameters constructor
//
//----------------------------------------------------------------------
Bool_t MFindDisp::FindParams(TString parDispInit,
                             TArrayD &params, TArrayD &steps)
{
    // Setup the event loop which will be executed in the 
    //                 fcnDisp function  of MINUIT
    //
    // parDispInit is the name of the file containing the initial values 
    // of the parameters; 
    // if parDispInit = ""   'params' and 'steps' are taken as initial values
    //

    *fLog << "=============================================" << endl;
    *fLog << "Setup event loop for fcnDisp" << endl;


    if (fMatrixTrainCalc == NULL || fMatrixTrainHists == NULL)
    {
      *fLog << "MFindDisp::FindParams; training matrices are not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTrainCalc->GetM().GetNrows() <= 0  ||  fMatrixTrainHists->GetM().GetNrows() <= 0)
    {
      *fLog << "MFindDisp::FindParams; training matrices have no entries"
            << endl;
      return kFALSE;
    }

    //---------------------------------------------------------
    MParList  parlistfcn;
    MTaskList tasklistfcn;

    // loop over rows of matrix
    MMatrixLoop loopcalc(fMatrixTrainCalc);
    MMatrixLoop loophists(fMatrixTrainHists);

    //--------------------------------
    // create container for the Disp parameters
    // and set them to their initial values
    MDispParameters *dispparams = fDispCalcTrain->GetDispParameters();

    // take initial values from file parDispInit
    if (parDispInit != "")
    {
      TFile inparam(parDispInit);
      dispparams->Read("MDispParameters");
      inparam.Close();
      *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file "
            << parDispInit << endl;
    }

    // take initial values from 'params' and/or 'steps'
    else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
    {
      if (params.GetSize()  != 0)
      {
        *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'"
              << endl;
        dispparams->SetParameters(params);
      }
      if (steps.GetSize()  != 0)
      {
        *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'"
              << endl;
        dispparams->SetStepsizes(steps);
      }
    }
    else
    {
        *fLog << "MFindDisp::FindParams; initial values and step sizes are taken "
	      << "from the MDispParameters constructor" << endl;
    }

    // fill the plots for Disp and sum up the Minimization parameter contributions
    MFillH filldispplots("MHDisp", "");

    //******************************
    // entries in MParList
    
    parlistfcn.AddToList(&tasklistfcn);
    parlistfcn.AddToList(dispparams);
    parlistfcn.AddToList(fHDispTrain);
    parlistfcn.AddToList(fCam);
    parlistfcn.AddToList(fMatrixTrainCalc);
    parlistfcn.AddToList(fMatrixTrainHists);

    //******************************
    // entries in MTaskList

    tasklistfcn.AddToList(&loopcalc);
    tasklistfcn.AddToList(&loophists);
    tasklistfcn.AddToList(fDispCalcTrain);
    tasklistfcn.AddToList(&filldispplots);

    //******************************

    MEvtLoop evtloopfcn("EvtLoopFCN");
    evtloopfcn.SetParList(&parlistfcn);
    *fLog << "Event loop for fcnDisp has been setup" << endl;

    // address of evtloopfcn is used in CallMinuit()


    //-----------------------------------------------------------------------
    //
    //----------   Start of minimization part   --------------------
    //
    // Do the minimization with MINUIT
    //
    // Be careful: This is not thread safe
    //
    *fLog << "========================================================" << endl;
    *fLog << "Start minimization for Disp" << endl;


    // -------------------------------------------
    // prepare call to MINUIT
    //

    // get initial values of parameters 
    fVinit = dispparams->GetParameters();
    fStep  = dispparams->GetStepsizes();

    TString name[fVinit.GetSize()];
    fStep.Set(fVinit.GetSize());
    fLimlo.Set(fVinit.GetSize());
    fLimup.Set(fVinit.GetSize());
    fFix.Set(fVinit.GetSize());

    fNpar = fVinit.GetSize();

    // define names, step sizes, lower and upper limits of the parameters
    // fFix[] = 0;  means the parameter is not fixed
    for (UInt_t i=0; i<fNpar; i++)
    {
        name[i]   = "p";
        name[i]  += i+1;
        //fStep[i]  = TMath::Abs(fVinit[i]/10.0);
        fLimlo[i] = -100.0;
        fLimup[i] =  100.0;
        fFix[i]   =     0;
    }

    // fix some parameters
    //for (UInt_t i=0; i<fNpar; i++)
    //{
    //  if (i == 1  ||  i==3)
    //	{
    //      fStep[i] = 0.0;
    //      fFix[i]  =   1;
    //	}
    //}
 

    // -------------------------------------------
    // call MINUIT

    TStopwatch clock;
    clock.Start();

    *fLog << "before calling CallMinuit" << endl;

    MMinuitInterface inter;               
    Bool_t rc = inter.CallMinuit(fcnDisp, name,
                                 fVinit, fStep, fLimlo, fLimup, fFix,
                                 &evtloopfcn, "MIGRAD", kFALSE);
 
    *fLog << "after calling CallMinuit" << endl;

    *fLog << "Time spent for the minimization in MINUIT :   " << endl;;
    clock.Stop();
    clock.Print();


    if (!rc)
        return kFALSE;

    *fLog << "Minimization for Disp finished" << endl;
    *fLog << "========================================================" << endl;


    // -----------------------------------------------------------------
    // in 'fcnDisp' the optimum parameter values were put into MDisp

    // write optimum parameter values onto file fFilenameParam
    
    TFile outparam(fFilenameParam, "RECREATE"); 
    dispparams->Write();
    outparam.Close();

    *fLog << "Optimum parameter values for Disp were written onto file '"
              << fFilenameParam << "' :" << endl;

    const TArrayD &check = dispparams->GetParameters();
    for (Int_t i=0; i<check.GetSize(); i++)
        *fLog << check[i] << ",  ";
    *fLog << endl;



    //-------------------------------------------
    // draw the plots
    fHDispTrain->Draw();

    *fLog << "End of Disp optimization part" << endl;
    *fLog << "======================================================" << endl;

    return kTRUE;
}


// -----------------------------------------------------------------------
//
// Test the result of the Disp optimization on the test sample
//

Bool_t MFindDisp::TestParams()
{
    if (fMatrixTestCalc == NULL || fMatrixTestHists == NULL)
    {
      *fLog << "MFindDisp::TestParams; test matrices are not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTestCalc->GetM().GetNrows() <= 0  ||  fMatrixTestHists->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindDisp::TestParams; test matrices have no entries" 
              << endl;
        return kFALSE;
    }

    // -------------   TEST the Disp optimization    -------------------
    //
    *fLog << "Test the Disp optimization on the test sample" << endl;

    // -----------------------------------------------------------------
    // read optimum parameter values from file fFilenameParam
    // into array 'dispPar'

    TFile inparam(fFilenameParam);
    MDispParameters dispin; 
    dispin.Read("MDispParameters");
    inparam.Close();

    *fLog << "Optimum parameter values for Disp were read from file '";
    *fLog << fFilenameParam << "' :" << endl;

    const TArrayD &dispPar = dispin.GetParameters();
    for (Int_t i=0; i<dispPar.GetSize(); i++)
        *fLog << dispPar[i] << ",  ";
    *fLog << endl;
    //---------------------------


    // -----------------------------------------------------------------
    MParList  parlist2;
    MTaskList tasklist2;

    MDispParameters *dispparams = fDispCalcTest->GetDispParameters();
    dispparams->SetParameters(dispPar);

    MMatrixLoop loopcalc(fMatrixTestCalc);
    MMatrixLoop loophists(fMatrixTestHists);

    MHMatrix *imatrix = NULL;
    TArrayI  imap;
    fHDispTest->GetMatrixMap(imatrix,imap);

    // attention: argument of MHDisp is the name of MImageParDisp container, that should
    // be the same than the name given to it when creating MDispCalc object at the MFindDisp
    // constructor:  fDispCalcTrain = new MDispCalc("DispTest");
    // fill the plots for Disp and sum up the Minimization parameter contributions
    MHDisp hdisp1("DispTest");
    hdisp1.SetName("MHDispCorr");
    hdisp1.SetSelectedPos(1);
    hdisp1.SetMatrixMap(imatrix,imap);
    MFillH filldispplots1("MHDispCorr[MHDisp]", "");
    
    MHDisp hdisp2("DispTest");
    hdisp2.SetName("MHDispWrong");
    hdisp2.SetSelectedPos(2);
    hdisp2.SetMatrixMap(imatrix,imap);
    MFillH filldispplots2("MHDispWrong[MHDisp]", "");
    
    MHDisp hdisp3("DispTest");
    hdisp3.SetName("MHDispM3Long");
    hdisp3.SetSelectedPos(3);
    hdisp3.SetMatrixMap(imatrix,imap);
    MFillH filldispplots3("MHDispM3Long[MHDisp]", "");
    
    MHDisp hdisp4("DispTest");
    hdisp4.SetName("MHDispAsym");
    hdisp4.SetSelectedPos(4);
    hdisp4.SetMatrixMap(imatrix,imap);
    MFillH filldispplots4("MHDispAsym[MHDisp]", "");
    

    //******************************
    // entries in MParList

    parlist2.AddToList(&tasklist2);
    parlist2.AddToList(dispparams);
    parlist2.AddToList(&hdisp1);
    parlist2.AddToList(&hdisp2);
    parlist2.AddToList(&hdisp3);
    parlist2.AddToList(&hdisp4);
    parlist2.AddToList(fCam);
    parlist2.AddToList(fMatrixTestCalc);
    parlist2.AddToList(fMatrixTestHists);

    //******************************
    // entries in MTaskList

    tasklist2.AddToList(&loopcalc);
    tasklist2.AddToList(&loophists);
    tasklist2.AddToList(fDispCalcTest);
    tasklist2.AddToList(&filldispplots1);
    tasklist2.AddToList(&filldispplots2);
    tasklist2.AddToList(&filldispplots3);
    tasklist2.AddToList(&filldispplots4);

    //******************************

    MProgressBar bar2;
    MEvtLoop evtloop2;
    evtloop2.SetParList(&parlist2);
    evtloop2.SetName("EvtLoopTestParams");
    evtloop2.SetProgressBar(&bar2);
    Int_t maxevents2 = -1;
    if (!evtloop2.Eventloop(maxevents2))
        return kFALSE;

    tasklist2.PrintStatistics(0, kTRUE);


    //-------------------------------------------
    // draw the plots

    hdisp1.Draw();
    hdisp2.Draw();
    hdisp3.Draw();
    hdisp4.Draw();

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


    *fLog << "Test of Disp otimization finished" << endl;
    *fLog << "======================================================" << endl;

    return kTRUE;
}
