/* ======================================================================== *\
!
! *
! * 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 "MDisp.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 "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 (MDisp::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 oject (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
    MDisp *disp = (MDisp*)plistfcn->FindObject("MDisp");
    if (!disp)
    {
        gLog << "fcnDisp : MDisp object '" << "MDisp"
            << "' not found... aborting" << endl;
        return;
    }

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

    //
    // transfer current parameter values to MDisp
    //
    // 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);

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

    // checking that parameters have been properly set
    TArrayD checkparameters = disp->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 matrix to compute the Chi^2
    // for the current parameter values
    evtloopfcn->Eventloop();

    tasklistfcn->PrintStatistics(0, kTRUE);

    //-------------------------------------------
    // get the Chi^2 value for the current values of the parameters
    f = hdisp->GetChi2();
}


// --------------------------------------------------------------------------
//
// 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;
    *fLog << inf << "MFindDisp::MFindDisp; address of MFDisp object = "
          << fDispFilter << endl;    

    // matrices to contain the training/test samples
    fMatrixTrain = new MHMatrix("MatrixTrain");
    fMatrixTest  = new MHMatrix("MatrixTest");

    // objects of MDispCalc to which these matrices are attached
    fDispCalcTrain = new MDispCalc("DispTrain");
    fDispCalcTest  = new MDispCalc("DispTest");

    // define columns of matrices
    fDispCalcTrain->InitMapping(fMatrixTrain);
    fDispCalcTest->InitMapping(fMatrixTest);
}

// --------------------------------------------------------------------------
//
// Default destructor.
//
MFindDisp::~MFindDisp()
{
    delete fCam;
    delete fMatrixTrain;
    delete fMatrixTest;
    delete fDispCalcTrain;
    delete fDispCalcTest;
}


// --------------------------------------------------------------------------
//
// Define the matrix 'fMatrixTrain' 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 Matrix 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 filltrain(fMatrixTrain);
    filltrain.SetFilter(&seltrain);
    filltrain.SetName("fillMatrixTrain");

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

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

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

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

    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 Matrix
    fMatrixTrain->Print("SizeCols");

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

    *fLog << "TRAINING Matrix was filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out training matrix

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

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

    if (display != NULL)
      delete display;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Define the matrix 'fMatrixTest' 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 Matrix 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 filltest(fMatrixTest);
    filltest.SetFilter(&seltest);
    filltest.SetName("fillMatrixTest");

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

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

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

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

    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 Matrix
    fMatrixTest->Print("SizeCols");

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

    *fLog << "TEST Matrix was filled" << endl;
    *fLog << "=============================================" << endl;


    //--------------------------
    // write out test matrix

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

      *fLog << "MFindDisp::DefineTestMatrix; Test matrix was 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 Matrix 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 filltrain(fMatrixTrain);
    filltrain.SetFilter(&split);
    filltrain.SetName("fillMatrixTrain");

    // 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 filltest(fMatrixTest);
    filltest.SetName("fillMatrixTest");


    //******************************
    // entries in MParList 
    
    plist.AddToList(&tlist);
    plist.AddToList(fCam);
    plist.AddToList(fMatrixTrain);
    plist.AddToList(fMatrixTest);

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

    tlist.AddToList(&read);
    if (fDispFilter != NULL)
      tlist.AddToList(&contdisp);
    tlist.AddToList(&cont);
    tlist.AddToList(&filltrain);
    tlist.AddToList(&conttrain);
    tlist.AddToList(&filltest);

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

    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 Matrix
    fMatrixTrain->Print("SizeCols");

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


    // print the filled Test Matrix
    fMatrixTest->Print("SizeCols");

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


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


    //--------------------------
    // write out training matrix

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

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

    //--------------------------
    // write out test matrix

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

      *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrix was 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 matrix

  TFile filetr(filetrain);
  fMatrixTrain->Read("MatrixTrain");
  fMatrixTrain->Print("SizeCols");

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


  //--------------------------
  // read in test matrix

  TFile filete(filetest);
  fMatrixTest->Read("MatrixTest");
  fMatrixTest->Print("SizeCols");

  *fLog << "MFindDisp::ReadMatrix; Test matrix was 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 Chi^2 for the current 
//                                                parameter values;
//        for this - Disp is calculated in the event loop by calling 
//                   MDispCalc::Process() ==> MDisp::Calc() 
//                 - the Chi^2 contributions are summed up in the event loop 
//                   by calling MHDisp::Fill()
//                 - after the event loop the Chi^2 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 Disp 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 (fMatrixTrain == NULL)
    {
      *fLog << "MFindDisp::FindParams; training matrix is not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTrain->GetM().GetNrows() <= 0)
    {
      *fLog << "MFindDisp::FindParams; training matrix has no entries"
            << endl;
      return kFALSE;
    }

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

    // loop over rows of matrix
    MMatrixLoop loop(fMatrixTrain);

    //--------------------------------
    // create container for the Disp parameters
    // and set them to their initial values
    MDisp disp;

    // take initial values from file parDispInit
    if (parDispInit != "")
    {
      TFile inparam(parDispInit);
      disp.Read("MDisp");
      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;
        disp.SetParameters(params);
      }
      if (steps.GetSize()  != 0)
      {
        *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'"
              << endl;
        disp.SetStepsizes(steps);
      }
    }
    else
    {
        *fLog << "MFindDisp::FindParams; initial values and step sizes are taken "
	      << "from the MDisp constructor" << endl;
    }

    // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
    TArrayI map;
    MHMatrix *tmpmatrix = NULL;
    tmpmatrix = fDispCalcTrain->GetMatrixMap(map);

    // 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("DispTrain");
    MHDisp hdisp("DispTrain");
    hdisp.SetMatrixMap(tmpmatrix,map);
    // fill the plots for Disp and sum up the Chi^2 contributions
    MFillH filldispplots("MHDisp", "");

    //******************************
    // entries in MParList
    
    parlistfcn.AddToList(&tasklistfcn);
    parlistfcn.AddToList(&disp);
    parlistfcn.AddToList(&hdisp);
    parlistfcn.AddToList(fCam);
    parlistfcn.AddToList(fMatrixTrain);

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

    tasklistfcn.AddToList(&loop);
    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 = disp.GetParameters();
    fStep  = disp.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"); 
    disp.Write();
    outparam.Close();

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

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



    //-------------------------------------------
    // draw the plots
    hdisp.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 (fMatrixTest == NULL)
    {
      *fLog << "MFindDisp::TestParams; test matrix is not defined... aborting"
            << endl;
      return kFALSE;
    }

    if (fMatrixTest->GetM().GetNrows() <= 0)
    {
        *fLog << "MFindDisp::TestParams; test matrix has 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);
    MDisp dispin; 
    dispin.Read("MDisp");
    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;

    MDisp disp;
    disp.SetParameters(dispPar);

    MMatrixLoop loop(fMatrixTest);

    // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
    TArrayI map;
    MHMatrix *tmpmatrix = NULL;
    tmpmatrix = fDispCalcTest->GetMatrixMap(map);

    // 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 Chi^2 contributions
    MHDisp hdisp1("DispTest");
    hdisp1.SetName("MHDispCorr");
    hdisp1.SetSelectedPos(1);
    hdisp1.SetMatrixMap(tmpmatrix,map);
    MFillH filldispplots1("MHDispCorr[MHDisp]", "");
    
    MHDisp hdisp2("DispTest");
    hdisp2.SetName("MHDispWrong");
    hdisp2.SetSelectedPos(2);
    hdisp2.SetMatrixMap(tmpmatrix,map);
    MFillH filldispplots2("MHDispWrong[MHDisp]", "");
    
    MHDisp hdisp3("DispTest");
    hdisp3.SetName("MHDispM3Long");
    hdisp3.SetSelectedPos(3);
    hdisp3.SetMatrixMap(tmpmatrix,map);
    MFillH filldispplots3("MHDispM3Long[MHDisp]", "");
    
    MHDisp hdisp4("DispTest");
    hdisp4.SetName("MHDispAsym");
    hdisp4.SetSelectedPos(4);
    hdisp4.SetMatrixMap(tmpmatrix,map);
    MFillH filldispplots4("MHDispAsym[MHDisp]", "");
    

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

    parlist2.AddToList(&tasklist2);
    parlist2.AddToList(&disp);
    parlist2.AddToList(&hdisp1);
    parlist2.AddToList(&hdisp2);
    parlist2.AddToList(&hdisp3);
    parlist2.AddToList(&hdisp4);
    parlist2.AddToList(fCam);
    parlist2.AddToList(fMatrixTest);

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

    tasklist2.AddToList(&loop);
    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;
}








