source: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc@ 5618

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  1/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Abelardo Moralejo 4/2003 <mailto:moralejo@pd.infn.it>
!   Author(s): Marcos Lopez 5/2004 <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MMcEnergyEst
//
// Class for otimizing the parameters of the energy estimator
//
// FIXME: the class must be made more flexible, allowing for different
// parametrizations to be used. Also a new class is needed, a container
// for the parameters of the energy estimator.
// FIXME: the starting values of the parameters are now fixed.
//
/////////////////////////////////////////////////////////////////////////////
#include "MMcEnergyEst.h"

#include <math.h>            // fabs on Alpha

#include <TMinuit.h>
#include <TStopwatch.h>
#include <TVirtualFitter.h>

#include "MParList.h"
#include "MTaskList.h"
#include "MGeomCamCT1.h"
#include "MFEventSelector.h"
#include "MReadTree.h"
#include "MFCT1SelFinal.h"
#include "MHMatrix.h"
#include "MEnergyEstParam.h"
#include "MMatrixLoop.h"
#include "MChisqEval.h"
#include "MEvtLoop.h"
#include "MDataElement.h"
#include "MDataMember.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MParameters.h"

ClassImp(MMcEnergyEst);

using namespace std;

//------------------------------------------------------------------------
//
// fcn calculates the function to be minimized (using TMinuit::Migrad)
//
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
    MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();

    MParList  *plist  = evtloop->GetParList();
    MTaskList *tlist = (MTaskList*)plist->FindObject("MTaskList");
 

    // Pass current minuit parameters to the energy estimation class 
    MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
    eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));

    
    evtloop->Eventloop();


    // Get result of the ChiSquare 
    MParameterD *eval = (MParameterD*)plist->FindObject("MFitResult", "MParameterD");
    
    f = eval->GetVal();
}


// --------------------------------------------------------------------------
//
// Default constructor.
//
MMcEnergyEst::MMcEnergyEst(const char *name, const char *title)
{
    fName  = name  ? name  : "MMcEnergyEst";
    fTitle = title ? title : "Optimizer of the energy estimator";

    SetHillasName("MHillas");
    SetHillasSrcName("MHillasSrc");

    //
    // Set initial values of the parameters (close enough to the final 
    // ones, taken from previous runs of the procedure). Parameter 
    // fA[4] is not used in the default energy estimation model (from 
    // D. Kranich).
    //
    fA.Set(5);
    fB.Set(7);

    fA[0] =  21006.2;
    fA[1] = -43.2648;
    fA[2] = -690.403;
    fA[3] = -0.0428544;
    fA[4] =  1.;
    fB[0] = -3391.05;
    fB[1] =  136.58;
    fB[2] =  0.253807;
    fB[3] =  254.363;
    fB[4] =  61.0386;
    fB[5] = -0.0190984;
    fB[6] = -421695;
}

Bool_t MMcEnergyEst::SetCoeff(TArrayD &coeff)
{
  if (coeff.GetSize() != fA.GetSize()+fB.GetSize())
    {
      *fLog << err << dbginf << "Wrong number of parameters!" << endl;
      return(kFALSE);
    }

  for (Int_t i = 0; i < fA.GetSize(); i++)
    fA[i] = coeff[i];
  for (Int_t i = 0; i < fB.GetSize(); i++)
    fB[i] = coeff[i+fA.GetSize()];

  return(kTRUE);

}

//------------------------------------------------------------------------
//
// Optimization (via migrad minimization) of parameters of energy estimation.
//
void MMcEnergyEst::FindParams()
{
  MParList parlist;

  MFEventSelector selector;
  selector.SetNumSelectEvts(fNevents);
  *fLog << inf << "Read events from file '" << fInFile << "'" << endl;    

  MReadTree read("Events");
  read.AddFile(fInFile);
  read.DisableAutoScheme();
  read.SetSelector(&selector);

  *fLog << inf << "Define columns of matrix" << endl;
  MHMatrix matrix;
  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
  Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact");

  if (colenergy < 0  ||  colimpact < 0)
  {
    *fLog << err << dbginf << "colenergy, colimpact = " << colenergy << ",  " 
          << colimpact << endl;
    return;
  }

  MEnergyEstParam eest(fHillasName);
  eest.Add(fHillasSrcName);
  eest.InitMapping(&matrix);
  
  *fLog << inf << "--------------------------------------" << endl;
  *fLog << inf << "Fill events into the matrix" << endl;
  if ( !matrix.Fill(&parlist, &read, fEventFilter) )
    return;
  *fLog << inf << "Matrix was filled with " << matrix.GetNumRows() 
        << inf << " events" << endl;  

  //-----------------------------------------------------------------------
  //
  // Setup the eventloop which will be executed in the fcn of MINUIT 
  //
  *fLog << inf << "--------------------------------------" << endl;
  *fLog << inf << "Setup event loop to be used in 'fcn'" << endl;

  MTaskList tasklist;

  MMatrixLoop loop(&matrix);

  MChisqEval eval;
  eval.SetY1(new MDataElement(&matrix, colenergy));
  eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
  eval.SetOwner();

  //
  // Entries in MParList

  parlist.AddToList(&tasklist);

  //
  // Entries in MTaskList

  tasklist.AddToList(&loop);
  tasklist.AddToList(&eest);
  tasklist.AddToList(&eval);


  *fLog << inf << "Event loop was setup" << endl;
  MEvtLoop evtloop;
  evtloop.SetParList(&parlist);

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

  gMinuit = new TMinuit(12);
  gMinuit->SetPrintLevel(-1);

  gMinuit->SetFCN(fcn);
  gMinuit->SetObjectFit(&evtloop);

  // Ready for: gMinuit->mnexcm("SET ERR", arglist, 1, ierflg)

  if (gMinuit->SetErrorDef(1))
    {
      *fLog << err << dbginf << "SetErrorDef failed." << endl;
      return;
    }

  //
  // Set starting values and step sizes for parameters
  //
  for (Int_t i=0; i<fA.GetSize(); i++)
    {
      TString name = "fA[";
      name += i;
      name += "]";
      Double_t vinit = fA[i];
      Double_t step  = fabs(fA[i]/3);

      Double_t limlo = 0; // limlo=limup=0: no limits
      Double_t limup = 0; 

      Bool_t rc = gMinuit->DefineParameter(i, name, vinit, step, limlo, limup);
      if (!rc)
        continue;

      *fLog << err << dbginf << "Error in defining parameter #" << i << endl;
      return;
    }

  for (Int_t i=0; i<fB.GetSize(); i++)
    {
      TString name = "fB[";
      name += i;
      name += "]";
      Double_t vinit = fB[i];
      Double_t step  = fabs(fB[i]/3);

      Double_t limlo = 0; // limlo=limup=0: no limits
      Double_t limup = 0;

      Bool_t rc = gMinuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
      if (!rc)
        continue;

      *fLog << err << dbginf << "Error in defining parameter #" << i+fA.GetSize() << endl;
      return;
    }

  TStopwatch clock;
  clock.Start();

  // Now ready for minimization step:

  gLog.SetNullOutput(kTRUE);
  Bool_t rc = gMinuit->Migrad();
  gLog.SetNullOutput(kFALSE);
  
  if (rc)
    {
      *fLog << err << dbginf << "Migrad failed." << endl;
      return;
    }

  *fLog << inf << endl;
  clock.Stop();
  clock.Print();
  *fLog << inf << endl;

  *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl;

  //
  // Update values of fA, fB:
  //
  for (Int_t i = 0; i < fA.GetSize(); i++)
    {
      Double_t x1, x2;
      gMinuit->GetParameter(i,x1,x2);
      fA[i] = x1;
    }
  for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++)
    {
      Double_t x1, x2;
      gMinuit->GetParameter(i,x1,x2);
      fB[i-fA.GetSize()] = x1;
    }

  //    eest.Print();
  eest.StopMapping();
  *fLog << inf << "Minimization for the energy estimator finished" << endl;

}

//------------------------------------------------------------------------
//
// Print current values of parameters
//
void MMcEnergyEst::Print(Option_t *o) const
{
  for (Int_t i = 0; i < fA.GetSize(); i++)
    *fLog << inf << "Parameter " << i << ": " << const_cast<TArrayD&>(fA)[i] << endl;

  for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++)
    *fLog << inf << "Parameter " << i << ": " << const_cast<TArrayD&>(fB)[i-fA.GetSize()] << endl;

  /*
    // Print results
    Double_t amin, edm, errdef;
    Int_t nvpar, nparx, icstat;
    gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
    gMinuit->mnprin(3, amin);
  */

}