/* ======================================================================== *\
!
! *
! * 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    6/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Wolfgang Wittek 7/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MCT1FindSupercuts                                                       //
//                                                                         //
// Class for otimizing the parameters of the supercuts                     //
//                                                                         //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCT1FindSupercuts.h"

#include <math.h>            // fabs 

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

#include "MBinning.h"
#include "MContinue.h"
#include "MCT1Supercuts.h"
#include "MCT1SupercutsCalc.h"
#include "MDataElement.h"
#include "MDataMember.h"

#include "MEvtLoop.h"
#include "MFCT1SelFinal.h"
#include "MF.h"
#include "MFEventSelector.h"
#include "MFEventSelector2.h"
#include "MFillH.h"
#include "MGeomCamCT1Daniel.h"
#include "MH3.h"
#include "MHCT1Supercuts.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 "MTaskList.h"


ClassImp(MCT1FindSupercuts);

using namespace std;


//------------------------------------------------------------------------
//
// fcnSupercuts 
//
// - calculates the quantity to be minimized (using TMinuit)
//
// - the quantity to be minimized is (-1)*significance of the gamma signal
//   in the alpha distribution (after cuts)
//
// - the parameters to be varied in the minimization are the cut parameters
//   (par)
//
static void fcnSupercuts(Int_t &npar, Double_t *gin, Double_t &f, 
                         Double_t *par, Int_t iflag)
{
    //-------------------------------------------------------------
    // save pointer to the MINUIT object for optimizing the supercuts
    // because it will be overwritten 
    // when fitting the alpha distribution in MHFindSignificance
    TMinuit *savePointer = gMinuit;
    //-------------------------------------------------------------


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

    MParList *plistfcn  = 
              (MParList*)evtloopfcn->GetParList(); 

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

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


    //-------------------------------------------------------
    // transfer current parameter values to MCT1Supercuts
    //
    Double_t fMin,   fEdm,  fErrdef;
    Int_t    fNpari, fNparx, fIstat;
    gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);

    //gLog << "fcnSupercuts : npar, fNpari, fNparx = " << npar << ",  " 
    //     << fNpari  << ",  " << fNparx << endl;

    TArrayD d(fNparx, par);
    super->SetParams(fNparx, d);


    //-------------------------------------------
    // plot alpha with the current cuts
    evtloopfcn->Eventloop();

    TString  mh3Name = "AlphaFcn";
    MH3* alpha = (MH3*)plistfcn->FindObject(mh3Name, "MH3");

    TH1 &alphaHist = alpha->GetHist();
    alphaHist.SetXTitle("|alpha|  [\\circ]");  

    //-------------------------------------------
    // set Minuit pointer to zero in order not to destroy the TMinuit
    // object for optimizing the supercuts
    gMinuit = NULL;

    //=================================================================
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    Double_t alphasig = 20.0;
    Double_t alphamin = 30.0;
    Double_t alphamax = 90.0;
    Int_t    degree   =    4;

    Bool_t   drawpoly;
    Bool_t   fitgauss;
    if ( iflag == 3)
    {
      drawpoly  = kTRUE;
      fitgauss  = kTRUE;
    }
    else
    {
      drawpoly  = kFALSE;
      fitgauss  = kFALSE;
    }
    drawpoly  = kFALSE;
    fitgauss  = kFALSE;

    Bool_t   print     = kTRUE;

    MHFindSignificance findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);
    
    Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
				  alphasig, drawpoly, fitgauss, print);

    //=================================================================

    // reset gMinuit to the MINUIT object for optimizing the supercuts 
    gMinuit = savePointer;
    //-------------------------------------------

    if (!rc)
    {
      gLog << "fcnSupercuts : FindSigma() failed" << endl;
      f = 1.e10;
      return;
    }

    // plot some quantities during the optimization
    plotsuper->Fill(&findsig);


    //------------------------
    // get significance
    Double_t significance = findsig.GetSignificance();
    f = significance>0.0 ? -significance : 0.0;


    //------------------------
    // optimize signal/background ratio
    //Double_t ratio = findsig.GetNbg()>0.0 ? 
    //                 findsig.GetNex()/findsig.GetNbg() : 0.0; 
    //f = -ratio;

    //-------------------------------------------
    // final calculations
    //if (iflag == 3)
    //{
    //
    //}    

    //-------------------------------------------------------------
}


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

    //---------------------------
    // default values 
    fFilenameTrain = "";
    fFilenameTest  = "";
    fFilenameParam = "";

    fHowManyTrain = 10000;
    fHowManyTest  = 10000;

    fHadronnessName = "";

    fMatrixFilter = NULL;

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

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

    // objects of MCT1SupercutsCalc to which these matrices are attached
    fCalcHadTrain = new MCT1SupercutsCalc("SupercutsCalcTrain");
    fCalcHadTest  = new MCT1SupercutsCalc("SupercutsCalcTest");

    // Define columns of matrices
    fCalcHadTrain->InitMapping(fMatrixTrain);
    fCalcHadTest->InitMapping(fMatrixTest);
}

// --------------------------------------------------------------------------
//
// Default destructor.
//
MCT1FindSupercuts::~MCT1FindSupercuts()
{
    delete fCam;
    delete fMatrixTrain;
    delete fMatrixTest;
    delete fCalcHadTrain;
    delete fCalcHadTest;
}

// --------------------------------------------------------------------------
//
// 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 MCT1FindSupercuts::DefineTrainMatrix(const TString &nametrain,
	                  const Int_t howmanytrain, MH3 &hreftrain)
{
    if ( nametrain == ""  ||  howmanytrain <= 0 )
      return kFALSE;

    *fLog << "=============================================" << endl;
    *fLog << "fill training matrix from file '" << nametrain 
          << "',   select " << howmanytrain 
          << " events according to a distribution given by the MH3 object '"
          << hreftrain.GetName() << "'" << endl;


    MParList  plist;
    MTaskList tlist;

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

    //MFEventSelector2 seltrain(hreftrain);
    //seltrain.SetNumMax(howmanytrain);
    MFEventSelector seltrain;
    seltrain.SetNumSelectEvts(howmanytrain);

    MFillH filltrain(fMatrixTrain);
    filltrain.SetFilter(&seltrain);

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

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

    tlist.AddToList(&read);
    tlist.AddToList(&seltrain);
    tlist.AddToList(&filltrain);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrain");
    evtloop.SetProgressBar(&bar);

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");
    Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
    if ( fabs(howmanygenerated-howmanytrain) > 3.0*sqrt(howmanytrain) )
    {
      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
	    << howmanygenerated 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    *fLog << "training matrix was filled" << endl;
    *fLog << "=============================================" << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Define the matrix for the test sample
//
// alltogether 'howmanytest' events are read from file 'nametest'
//
// the events are selected according to a target distribution 'hreftest'
//
//
Bool_t MCT1FindSupercuts::DefineTestMatrix(const TString &nametest,
	                  const Int_t howmanytest, MH3 &hreftest)
{
    if ( nametest == ""  ||  howmanytest <= 0 )
      return kFALSE;

    *fLog << "=============================================" << endl;
    *fLog << "fill test matrix from file '" << nametest 
          << "',   select " << howmanytest 
          << " events according to a distribution given by the MH3 object '"
          << hreftest.GetName() << "'" << endl;


    MParList  plist;
    MTaskList tlist;

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

    //MFEventSelector2 seltest(hreftest);
    //seltest.SetNumMax(howmanytest);
    MFEventSelector seltest;
    seltest.SetNumSelectEvts(howmanytest);

    MFillH filltest(fMatrixTest);
    filltest.SetFilter(&seltest);

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

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

    tlist.AddToList(&read);
    tlist.AddToList(&seltest);
    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTest");
    evtloop.SetProgressBar(&bar);

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTest->Print("SizeCols");
    Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
    if ( fabs(howmanygenerated-howmanytest) > 3.0*sqrt(howmanytest) )
    {
      *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events ("
	    << howmanygenerated 
            << ") is incompatible with the no.of requested events ("
            << howmanytest << ")" << endl;
    }

    *fLog << "test matrix was filled" << endl;
    *fLog << "=============================================" << endl;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Define the matrices for the training and test sample respectively
//
//
//
Bool_t MCT1FindSupercuts::DefineTrainTestMatrix(
                          const TString &name,
	                  const Int_t howmanytrain, MH3 &hreftrain,
	                  const Int_t howmanytest,  MH3 &hreftest)
{
    *fLog << "=============================================" << endl;
    *fLog << "fill training matrix from file '" << name 
          << "',   select " << howmanytrain 
          << " events for the training according to a distribution given by the MH3 object '"
          << hreftrain.GetName() << "'" << endl;

    *fLog << "fill test matrix from the same file '" << name 
          << "',   select " << howmanytest 
          << " events for the training according to a distribution given by the MH3 object '"
          << hreftest.GetName() << "'" << endl;


    MParList  plist;
    MTaskList tlist;

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

    //MFEventSelector2 seltrain(hreftrain);
    //seltrain.SetNumMax(howmanytrain);
    MFEventSelector seltrain;
    seltrain.SetName("SelTrain");
    seltrain.SetNumSelectEvts(howmanytrain);

    MFillH filltrain(fMatrixTrain);
    filltrain.SetName("fillMatrixTrain");
    filltrain.SetFilter(&seltrain);

    // consider this event as candidate for a test event 
    // only if event was not accepted as a training event
    MContinue cont(&seltrain);

    Float_t fcorr =  (Float_t)read.GetEntries() / 
                    ((Float_t)read.GetEntries()-(Float_t)howmanytrain);

    *fLog << "entries, howmanytrain, fcorr = " << read.GetEntries()
          << ",  " << howmanytrain << ",  " << fcorr << endl;

    //MFEventSelector2 seltest(hreftest);
    //seltrain.SetNumMax(howmanytest*fcorr);
    MFEventSelector seltest;
    seltest.SetName("SelTest");
    seltest.SetNumSelectEvts(howmanytest*fcorr);

    MFillH filltest(fMatrixTest);
    filltest.SetName("fillMatrixTest");
    filltest.SetFilter(&seltest);

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

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

    tlist.AddToList(&read);
    tlist.AddToList(&seltrain);
    tlist.AddToList(&filltrain);

    tlist.AddToList(&cont);

    tlist.AddToList(&seltest);
    tlist.AddToList(&filltest);

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

    MProgressBar bar;
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetName("EvtLoopMatrixTrainTest");
    evtloop.SetProgressBar(&bar);

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");
    Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
    if ( fabs(generatedtrain-howmanytrain) > 3.0*sqrt(howmanytrain) )
    {
      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
	    << generatedtrain 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    fMatrixTest->Print("SizeCols");
    Int_t generatedtest = fMatrixTest->GetM().GetNrows();
    if ( fabs(generatedtest-howmanytest) > 3.0*sqrt(howmanytest) )
    {
      *fLog << "MCT1FindSupercuts::DefineTestMatrix; 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;


  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read training and test matrices from files
//
//
// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
//

Bool_t MCT1FindSupercuts::ReadMatrix(TString &filetrain, TString &filetest)
{
  //--------------------------
  // read in training matrix

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

  *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
        << filetrain << "'" << endl;


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

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

  *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
        << filetest << "'" << endl;


  return kTRUE;  
}

// --------------------------------------------------------------------------
//
// Write training and test matrices onto files
//
//
// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
//
//
Bool_t MCT1FindSupercuts::WriteMatrix(TString &filetrain, TString &filetest)
{
  //--------------------------
  // write out training matrix

  TFile filetr(filetrain, "RECREATE", "");

  *fLog << "nach TFile" << endl;

  fMatrixTrain->Write();

  *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '"
        << filetrain << "'" << endl;


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

  TFile filete(filetest, "RECREATE", "");
  fMatrixTest->Print("SizeCols");
  fMatrixTest->Write();

  *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '"
        << filetest << "'" << endl;


  return kTRUE;  

}

//------------------------------------------------------------------------
//
// Steering program for optimizing the supercuts
// ---------------------------------------------
//
//      the criterion for the 'optimum' is 
//
//          - a maximum significance of the gamma signal in the alpha plot, 
//            in which the supercuts have been applied
//
// The various steps are :
//
// - setup the event loop to be executed for each call to fcnSupercuts 
// - call TMinuit to do the minimization of (-significance) :
//        the fcnSupercuts function calculates the significance 
//                                             for the current cut values
//        for this - the alpha plot is produced in the event loop, 
//                   in which the cuts have been applied
//                 - the significance of the gamma signal in the alpha plot 
//                   is calculated
//
// Needed as input : (to be set by the Set functions)
//
// - fFilenameParam      name of file to which optimum values of the 
//                       parameters are written
//
// - fHadronnessName     name of container where MCT1SupercutsCalc will
//                       put the supercuts hadronness
//
// - for the minimization, the starting values of the parameters are taken as 
//   MCT1Supercuts::GetParams(fVinit)
//
//----------------------------------------------------------------------
Bool_t MCT1FindSupercuts::FindParams()
{
    // Setup the event loop which will be executed in the 
    //                 fcnSupercuts function  of MINUIT
    //

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

    if ( fHadronnessName == "")
    {
      *fLog << "MCT1FindSupercuts::FindParams; hadronness name is not defined... aborting"
            << endl;
      return kFALSE;
    }

    if ( fMatrixTrain == NULL)
    {
      *fLog << "MCT1FindSupercuts::FindParams; training matrix is not defined... aborting"
            << endl;
      return kFALSE;
    }

    MParList  parlistfcn;
    MTaskList tasklistfcn;

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

    // create container for the supercut parameters
    // and set them to their default values
    MCT1Supercuts super;
    super.InitParams();

    // calculate supercuts hadronness
    fCalcHadTrain->SetHadronnessName(fHadronnessName);

    // apply the supercuts
    MF scfilter(fHadronnessName+".fHadronness>0.5");
    MContinue supercuts(&scfilter);

    // plot |alpha|
    TString  mh3Name = "AlphaFcn";
    MBinning binsalpha("Binning"+mh3Name);
    binsalpha.SetEdges(54, -12.0, 96.0);

    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alpha("MatrixTrain[7]");
    alpha.SetName(mh3Name);

    MFillH fillalpha(&alpha);

    // book histograms to be filled during the optimization
    //                              ! not in the event loop !
    MHCT1Supercuts plotsuper;
    plotsuper.SetupFill(&parlistfcn);

    //******************************
    // entries in MParList (extension of old MParList)
    
    parlistfcn.AddToList(&tasklistfcn);
    parlistfcn.AddToList(&super);
    parlistfcn.AddToList(fCam);
    parlistfcn.AddToList(fMatrixTrain);

    parlistfcn.AddToList(&binsalpha);
    parlistfcn.AddToList(&alpha);

    parlistfcn.AddToList(&plotsuper);

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

    tasklistfcn.AddToList(&loop);
    tasklistfcn.AddToList(fCalcHadTrain);
    tasklistfcn.AddToList(&supercuts);
    tasklistfcn.AddToList(&fillalpha);


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

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

    // address of evtloopfcn is used in CallMinuit()
    fObjectFit = &evtloopfcn;


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


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

    fNpar = 48;

    // get initial values of parameters from MCT1SupercutsCalc
    TArrayD vinit(fNpar);
    super.GetParams(fNpar, vinit);    

    for (UInt_t i=0; i<fNpar; i++)
    {
      fVinit[i] = vinit[i];
    }

    for (UInt_t i=0; i<fNpar; i++)
    {
      sprintf(&fParName[i][0], "p%d", i+1);
      fStep[i]  = fabs(fVinit[i]/10.0);
      fLimlo[i] = -10.0;
      fLimup[i] =  10.0;
      fFix[i]   =     0;
    }

    fMethod = "SIMPLEX";        
    Bool_t nulloutput = kFALSE;   

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

    TStopwatch clock;
    clock.Start();

    MMinuitInterface inter;               
    Bool_t rc = inter.CallMinuit( fcnSupercuts, fNpar,  fParName,   
                                  fVinit, fStep, fLimlo, fLimup, fFix,   
                                  fObjectFit, fMethod, nulloutput);
 
    *fLog << "Time spent for the minimization in MINUIT :   " << endl;;
    clock.Stop();
    clock.Print();

    plotsuper.DrawClone();

    if (!rc) return kFALSE;

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


    // -----------------------------------------------------------------
    // in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put 
    //                    into MCT1Supercuts

    // write optimum parameter values onto file filenameParam
    
    TFile outparam(fFilenameParam, "RECREATE"); 
    super.Write();
    outparam.Close();

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

    TArrayD check(72);
    super.GetParams(72, check);
    for (Int_t i=0; i<72; i++)
    {
      *fLog << check[i] << ",  ";
    }
    *fLog << endl;



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

    return kTRUE;
}


// -----------------------------------------------------------------------
//
// Test the supercuts on the test sample
//

Bool_t MCT1FindSupercuts::TestParams()
{
    Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();

    if ( howmanygenerated <= 0 )
    {
      *fLog << "MCT1FindSupercuts::TestParams; test matrix is empty... aborting"
            << endl;
      return kFALSE;
    }

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


    // -----------------------------------------------------------------
    // read optimum parameter values from file filenameParam
    // into array 'supercutsPar'

    TFile inparam(fFilenameParam);
    MCT1Supercuts scin; 
    scin.Read("MCT1Supercuts");
    inparam.Close();

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

    TArrayD supercutsPar(72);;
    scin.GetParams(72, supercutsPar);
    for (Int_t i=0; i<72; i++)
    {
      *fLog << supercutsPar[i] << ",  ";
    }
    *fLog << endl;
    //---------------------------


    // -----------------------------------------------------------------
    if ( fHadronnessName == "")
    {
      *fLog << "MCT1FindSupercuts::TestParams; hadronness name is not defined... aborting"
            << endl;
      return kFALSE;
    }


    MParList  parlist2;
    MTaskList tasklist2;

    MCT1Supercuts supercuts;
    supercuts.SetParams(72, supercutsPar);

    fCalcHadTest->SetHadronnessName(fHadronnessName);


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

    MMatrixLoop loop(fMatrixTest);

    // plot alpha before applying the supercuts
    TString  mh3NameB = "AlphaBefore";
    MBinning binsalphabef("Binning"+mh3NameB);
    binsalphabef.SetEdges(54, -12.0, 96.0);

    // |alpha| is assumed to be in column 7 of the matrix
    MH3 alphabefore("MatrixTest[7]");
    alphabefore.SetName(mh3NameB);
    MFillH fillalphabefore(&alphabefore);
    fillalphabefore.SetName("FillAlphaBefore");

    // apply the supercuts
    MF scfilter(fHadronnessName+".fHadronness>0.5");
    MContinue applysupercuts(&scfilter);

    // plot alpha after applying the supercuts
    TString  mh3NameA = "AlphaAfter";
    MBinning binsalphaaft("Binning"+mh3NameA);
    binsalphaaft.SetEdges(54, -12.0, 96.0);

    MH3 alphaafter("MatrixTest[7]");
    alphaafter.SetName(mh3NameA);
    MFillH fillalphaafter(&alphaafter);
    fillalphaafter.SetName("FillAlphaAfter");

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

    parlist2.AddToList(&tasklist2);

    parlist2.AddToList(&supercuts);

    parlist2.AddToList(fCam);
    parlist2.AddToList(fMatrixTest);

    parlist2.AddToList(&binsalphabef);
    parlist2.AddToList(&alphabefore);

    parlist2.AddToList(&binsalphaaft);
    parlist2.AddToList(&alphaafter);

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

    tasklist2.AddToList(&loop);
    tasklist2.AddToList(&fillalphabefore);

    tasklist2.AddToList(fCalcHadTest);
    tasklist2.AddToList(&applysupercuts);

    tasklist2.AddToList(&fillalphaafter);

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

    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 alpha plots

    MH3* alphaBefore = (MH3*)(parlist2.FindObject(mh3NameB, "MH3"));
    TH1  &alphaHist1 = alphaBefore->GetHist();
    alphaHist1.SetXTitle("|alpha|  [\\circ]");

    MH3* alphaAfter = (MH3*)(parlist2.FindObject(mh3NameA, "MH3"));
    TH1  &alphaHist2 = alphaAfter->GetHist();
    alphaHist2.SetXTitle("|alpha|  [\\circ]");

    TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
    c->Divide(2,1);

    gROOT->SetSelectedPad(NULL);

    c->cd(1);
    alphaHist1.DrawCopy();

    c->cd(2);
    alphaHist2.DrawCopy();



    //-------------------------------------------
    // fit alpha distribution to get the number of excess events and
    // calculate significance of gamma signal in the alpha plot
  
    Double_t alphasig = 20.0;
    Double_t alphamin = 30.0;
    Double_t alphamax = 90.0;
    Int_t    degree   =    4;
    Double_t significance = -99.0;
    Bool_t   drawpoly  = kTRUE;
    Bool_t   fitgauss  = kFALSE;
    Bool_t   print     = kTRUE;

    MHFindSignificance findsig;
    findsig.SetRebin(kTRUE);
    findsig.SetReduceDegree(kFALSE);

    findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, 
                      alphasig, drawpoly, fitgauss, print);
    significance = findsig.GetSignificance();
    Double_t alphasi = findsig.GetAlphasi();

    *fLog << "Significance of gamma signal after supercuts in test sample : "
         << significance << " (for |alpha| < " << alphasi << " degrees)" 
         << endl;
    //-------------------------------------------


    *fLog << "Test of supercuts part finished" << endl;
    *fLog << "======================================================" << endl;

    return kTRUE;
}


//===========================================================================









