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

/* ======================================================================== *\
!
! *
! * 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@uni-sw.gwdg.de>
!   Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Abelardo Moralejo 4/2003 <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MMcEnergyEst                                                            //
//                                                                         //
// Class for otimizing the parameters of the energy estimator              //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MMcEnergyEst.h"

#include <TStopwatch.h>
#include <TVirtualFitter.h>
#include <TMinuit.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"

extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
TMinuit *minuit;

ClassImp(MMcEnergyEst);

// --------------------------------------------------------------------------
//
// 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");
    SetHadronnessName("MHadronness");
}

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

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

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

  MFCT1SelFinal filterhadrons;
  filterhadrons.SetHadronnessName(fHadronnessName);
  filterhadrons.SetCuts(fMaxHadronness, fMaxAlpha, fMaxDist);
  filterhadrons.SetInverted();

  cout << "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)
  {
    cout << "colenergy, colimpact = " << colenergy << ",  " 
         << colimpact << endl;
    return;
  }

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

  //-----------------------------------------------------------------------
  //
  // Setup the eventloop which will be executed in the fcn of MINUIT 
  //
  cout << "--------------------------------------" << endl;
  cout << "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);


  cout << "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
  //
  cout << "--------------------------------------" << endl;
  cout << "Start minimization" << endl;

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

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

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

  if (minuit->SetErrorDef(1))
    {
      cout << "SetErrorDef failed." << endl;
      return;
    }

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

  //
  // 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 = minuit->DefineParameter(i, name, vinit, step, limlo, limup);
      if (!rc)
        continue;

      cout << "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 = minuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
      if (!rc)
        continue;

      cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
      return;
    }

  TStopwatch clock;
  clock.Start();

  // Now ready for minimization step:

  gLog.SetNullOutput(kTRUE);
  Bool_t rc = minuit->Migrad();
  gLog.SetNullOutput(kFALSE);
  
  if (rc)
    {
      cout << "Migrad failed." << endl;
      return;
    }

  cout << endl;
  clock.Stop();
  clock.Print();
  cout << endl;

  cout << "Resulting Chisq: " << minuit->fAmin << endl;

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

  //    eest.Print();
  eest.StopMapping();
  cout << "Minimization finished" << endl;

}

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

  for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++)
    cout << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl; 

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

}

//------------------------------------------------------------------------
//
// fcn calculates the function to be minimized (using TMinuit::Migrad)
//
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  MEvtLoop *evtloop = (MEvtLoop*)minuit->GetObjectFit();
 
  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();

  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
  MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");

  eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));

  evtloop->Eventloop();

  f = eval->GetChisq();
}