source: trunk/MagicSoft/Mars/macros/CT1EgyEst.C@ 2153

/* ======================================================================== *\
!
! *
! * 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,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Abelardo Moralejo, 03/2003 <mailto:moralejo@pd.infn.it>
!              Wolfgang Wittek,   04/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

#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 "MLog.h"
#include "MLogManip.h"

void CT1EgyEst()
{
  //  TString inPath        = "~magican/ct1test/wittek/marsoutput/optionA/";
  TString inPath        = "./";
  TString fileNameIn    = "MC1.root";

  //  TString outPath       = "~magican/ct1test/wittek/marsoutput/optionA/";
  TString outPath       = "./";
  TString energyParName = "energyest.root";
 
  TString hilName    = "MHillas";
  TString hilSrcName = "MHillasSrc";

  //  TString hadronnessName = "MHadronness";
  TString hadronnessName = "HadNN";

  //  Int_t howMany = 10000;
  Int_t howMany = 2000;

  Float_t maxhadronness = 0.3;
  Float_t maxalpha      = 20.0;
  Float_t maxdist       = 1.;

  CT1EEst(inPath,  fileNameIn, 
            outPath, energyParName, 
            hilName, hilSrcName,    hadronnessName, 
            howMany, maxhadronness, maxalpha, maxdist);
}

//------------------------------------------------------------------------
//
// Arguments of CT1EEst :
// ------------------------
//
// inPath,  fileNameIn       path and name of input file (MC gamma events)
// outPath, energyParName    path and name of file containing the parameters
//                           of the energy estimator
// hilName, hilSrcName       names of "MHillas" and "MHillasSrc" containers
// hadronnessName            name of container holding the hadronness
// howMany                   no.of events to be read from input file
// maxhadronness             maximum value of hadronness to be accepted
// maxalpha                  maximum value of |ALPHA| to be accepted
// maxdist                   maximum value of DIST (degrees) to be accepted
//

void CT1EEst(TString inPath,  TString fileNameIn, 
               TString outPath, TString energyParName,
               TString hilName, TString hilSrcName, TString hadronnessName,
               Int_t howMany, 
               Float_t maxhadronness, Float_t maxalpha, Float_t maxdist)
{
  cout << "================================================================"
       << endl;
  cout << "Start of energy estimation part" << endl;
  cout << "" << endl;
  cout << "CT1EEst input values : " << endl;
  cout << "---------------------- " << endl;
  cout << "inPath, fileNameIn = " 
       <<  inPath << ",  " << fileNameIn << endl;
  cout << "outPath, energyParName = " 
       <<  outPath << ",  " << energyParName << endl;
  cout << "hilName, hilSrcName, hadronnessName = " << hilName << ",  "
       << hilSrcName << ",  " << hadronnessName << endl;
  cout << "howMany, maxhadronness, maxalpha, maxdist = " << howMany << ",  "
       << maxhadronness << ",  " << maxalpha << ",  " << maxdist << endl;
  cout << "" << endl;


  //==========================================================================
  // Start of Part 1 (determination of the parameters of the energy estimator)
  //

  //-----------------------------------------------------------------------
  // Fill events into a MHMatrix, 
  // to be used later in the minimization by MINUIT
  //

  MMcEnergyEst Optimize;

  TString inputfile(outPath);
  inputfile += fileNameIn;
  Optimize.SetInFile(inputfile);

  TString paramout(outPath);
  paramout += energyParName;
  Optimize.SetOutFile(paramout);

  MFCT1SelFinal filterhadrons;
  filterhadrons.SetHadronnessName(hadronnessName);
  filterhadrons.SetCuts(maxhadronness, maxalpha, maxdist);
  filterhadrons.SetInverted();
  Optimize.SetEventFilter(&filterhadrons);

  Optimize.SetNevents(howMany);

  //
  // Find the optimal parameters for the energy estimator:
  //
  Optimize.FindParams();

  cout << "--------------------------------------" << endl;
  cout << "Write parameters of energy estimator onto file '" 
       << paramout << endl;
  cout << "--------------------------------------" << endl;

  Optimize.Print();

  TFile out(paramout, "recreate");
  Optimize.SetReadyToSave();
  Optimize.Write();
  out.Close();

  //
  // END of Part 1
  //==========================================================================



  //==========================================================================
  // Start of Part 2 (test quality of energy estimation)
  //
  //

  //--------------------------------------------
  // Read the parameters of the energy estimator 
  //

  TString paramin(outPath);
  paramin += energyParName;

  TFile enparam(paramin);
  MMcEnergyEst mcest;
  mcest.Read("MMcEnergyEst");
  enparam.Close();

  cout << "--------------------------------------" << endl;
  cout << "Read parameters of energy estimator from file '" 
       << paramin << endl;

  TArrayD parA(5);
  TArrayD parB(7);

  for (Int_t i=0; i<parA.GetSize(); i++)
    parA[i] = mcest.GetCoeff(i);
  for (Int_t i=0; i<parB.GetSize(); i++)
    parB[i] = mcest.GetCoeff(i+parA.GetSize());

  //-----------------------------------------------
  // Create energy estimator task, add necessary containers and 
  // initialize parameters read from file:
  //

  MEnergyEstParam eest2(hilName);
  eest2.Add(hilSrcName);

  eest2.SetCoeffA(parA);
  eest2.SetCoeffB(parB);


  //=======================================================================
  // Setup the event loop
  //
  cout << "--------------------------------------" << endl;
  cout << "Setup event loop for checking quality of energy estimation" << endl;


  MTaskList tlist2;
  MParList  parlist2;  

  //-----------------------------------------------
  // Read events
  //

  TString inputfile(inPath);
  inputfile += fileNameIn;

  cout << "Read events from file '" << inputfile << "'" << endl;    
  MReadTree read2("Events");
  read2.AddFile(inputfile);
  read2.DisableAutoScheme();


  //-----------------------------------------------
  // Select events
  //
  cout << "Select events with hadronness < " << maxhadronness 
        << " and |alpha| < " << maxalpha << endl; 
  MFCT1SelFinal hcut2;
  hcut2.SetHadronnessName(hadronnessName); hcut2;
  hcut2.SetCuts(maxhadronness, maxalpha, maxdist);

  MContinue cont(&hcut2);

  parlist2.AddToList(&tlist2);

  //********************************
  // Entries in MTaskList

  tlist2.AddToList(&read2);
  tlist2.AddToList(&cont);
  tlist2.AddToList(&eest2);

  //
  // Create Object MHMcEnergyMigration containing useful histograms,
  // and task MHMcEnergyMigration to fill them:
  //

  MHMcEnergyMigration mighist;

  parlist2.AddToList(&mighist);

  MFillH migfill(&mighist, "MMcEvt");

  tlist2.AddToList(&migfill);

  MBinning BinningE("BinningE");
  MBinning BinningTheta("BinningTheta");
  BinningE.SetEdgesLog(50, 500, 50000.);
  BinningTheta.SetEdges(5, 0., 50.);
  
  parlist2.AddToList(&BinningE);
  parlist2.AddToList(&BinningTheta);

  MBinning BinningDE("BinningDE");
  MBinning BinningImpact("BinningImpact");
 
  BinningDE.SetEdges(60, -1.2, 1.2);
  BinningImpact.SetEdges(50, 0., 400.);
  parlist2.AddToList(&BinningDE);
  parlist2.AddToList(&BinningImpact);

  cout << "Event loop was setup" << endl; 
  MProgressBar bar;

  MEvtLoop evtloop2;
  evtloop2.SetProgressBar(&bar);
  evtloop2.SetParList(&parlist2);

  if (!evtloop2.Eventloop())
    return;

  TString paramout(outPath);
  paramout += energyParName;

  TFile out2(paramout, "UPDATE");
  mighist.Write();
  out2.Close();

  mighist.Draw();

  cout << "Quality histograms were added onto the file '" << paramout << endl;
  cout << endl;
  cout << "End of energy estimation part" << endl;
  cout << "================================================================" << endl;
  //
  // End of Part 2 (test quality of energy estimation)
  //==========================================================================
  //

  return;

}