/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Wolfgang Wittek, 7/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Marcos Lopez, 05/2004 <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

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

#include <math.h>            // fabs 

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

#include "MBinning.h"
#include "MContinue.h"
#include "MMySuperCuts.h"
#include "MMySuperCutsCalc.h"
#include "MDataElement.h"
#include "MDataMember.h"

#include "MEvtLoop.h"
#include "MFSelFinal.h"
#include "MF.h"
#include "MFEventSelector.h"
#include "MFEventSelector2.h"
#include "MFillH.h"
//#include "MGeomCamCT1Daniel.h"
#include "MFEventSelector.h"
#include "MGeomCamMagic.h"
#include "MH3.h"
#include "MHSupercuts.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"

#include "TVector2.h"
#include "MSrcPosCam.h"
#include "MHillasSrcCalc.h"
#include "MSrcPosCalc.h"
#include "MObservatory.h"

ClassImp(MMyFindSuperCuts);

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)
//
// -------------------------
//   - Read the Matrix
//   - Fill a histogram with the alpha (6th column)
//   - Fit the alpha plot
//   - Get result of fit (significance)
//
static void fcnSuperCuts(Int_t &npar, Double_t *gin, Double_t &f, 
                         Double_t *par, Int_t iflag)
{

    //
    // Set by hand the paramters range
    //
    if ( par[0]<0 || par[0]>1.8 || par[4]>par[0] || par[4]<0 || par[0]<par[4] ||
	 par[2]>1 || par[3]>1 ||  par[4]>1 || par[1]>1   || par[6]>1 || par[7]>1   ||
	 par[8]<0 || par[8]>1.8 || par[12]>par[0] || par[12]<0 || par[8]<par[12] ||
	 par[13]<0 || par[13]>1. || par[14]<0 || par[14]>1 || par[15]<0 || par[15]>1 || 
	 par[16]<0 || par[16]>1.8 || par[20]<0 || par[20]>1.8 || par[16]<par[20] ||
	 // par[4]>0.18 ||
	 par[24]>200 || par[28]<-200 || 
	 par[32]>1 || par[36]<0 || 
  	 par[40]>1
	 )
    {
	gLog << "Parameter out of limit"  << endl;
	f=1e10;
	cout << f << endl;
	return;
    }

    //-------------------------------------------------------------
    // 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();
  
    MMySuperCuts *super = (MMySuperCuts*)plistfcn->FindObject("MMySuperCuts");
    if (!super)
    {
        gLog << "fcnSupercuts : MMySuperCuts object '" << "MMySuperCuts"
            << "' not found... aborting" << endl;
        return;
    }

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

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

    
    //$$$$$$$$$$$$$$$$$$$$$
    // for testing
    gLog.Separator("Current params");
      TArrayD checkparameters = super->GetParameters();
      gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ",  " 
           << fNparx  << endl;
    //  gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
  //   for (Int_t i=0; i<fNparx; i++)
//     {
//       gLog << i << ",  " << par[i] << ",  " << checkparameters[i] << endl;
//     }
  super->Print();
    
    //$$$$$$$$$$$$$$$$$$$$$


    //
    // EventLoop: Fill the alpha plot with these cuts
    //
    gLog.SetNullOutput(kTRUE);
    evtloopfcn->Eventloop();

    //tasklistfcn->PrintStatistics(0, kTRUE);



    //
    // Calculate significance
    //
    MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3");
    if (!alpha)
        return;

    TH1 &alphaHist = alpha->GetHist();
    alphaHist.SetName("alpha-fcnSupercuts");

    //-------------------------------------------
    // 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
    //
    const Double_t alphasig = 10.0;
    const Double_t alphamin = 30.0;
    const Double_t alphamax = 90.0;
    const Int_t    degree   =    0;

    Bool_t drawpoly;
    Bool_t fitgauss;
    TCanvas* c; 
    if (iflag == 3)
    {
	c= new TCanvas();
        drawpoly  = kTRUE;
        fitgauss  = kTRUE;
    }
    else
    {
        drawpoly  = kFALSE;
        fitgauss  = kFALSE;
    }
    fitgauss  = kTRUE;
    const Bool_t print = kFALSE;

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


    const Bool_t rc = findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
                                        alphasig, drawpoly, fitgauss, print);
    gLog.SetNullOutput(kFALSE);

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

    // 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
    MHSupercuts *plotsuper = (MHSupercuts*)plistfcn->FindObject("MHSupercuts");
    if (plotsuper)
        plotsuper->Fill(&findsig);


    //
    // Print info
    //
    gLog << "(Non, Nbg, Nex, S/N, sigma, Significane) = " 
	 << findsig.GetNon() << " " 
	 << findsig.GetNbg()<< " "
	 << findsig.GetNex() << " "
	 << findsig.GetNex()/  findsig.GetNbg()<< " "
	 << findsig.GetSigmaGauss() << " " 
	 << findsig.GetSignificance() << endl;

    const Double_t sigmagauss = findsig.GetSigmaGauss();
   //   if (sigmagauss>10)
//        {
//          gLog << "fcnSupercuts : Alpha plot too width = " << sigmagauss << endl;
//           f = 1.e10;
//         return;
//      }


    //
    // Get significance
    //
    const Double_t significance = findsig.GetSignificance();
    //const Double_t ratio = findsig.GetNex()/findsig.GetNbg();

  //      f = significance>0 ? -significance : 0;
//     f = -TMath::Abs(significance)/sqrt(sigmagauss);
//     f = - significance*sqrt(findsig.GetNex())/sqrt(sigmagauss); 
    f = - significance*significance/sqrt(sigmagauss);

 // f = - significance*findsig.GetProb();
   cout << f << endl; 

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

//      //f = -ratio; 
//       if(significance < 5) 
//    	f=1e10;

//      f=f*ratio*ratio;
//      cout << f << " " << significance << " " << ratio  << endl; 

}



// --------------------------------------------------------------------------
//
// Default constructor.
//
MMyFindSuperCuts::MMyFindSuperCuts(const char *name, const char *title)
    : fHowManyTrain(10000), fHowManyTest(10000), fMatrixFilter(NULL),fSizeCutLow(2000.0),fSizeCutUp(10000000), fOptimizationMode(0)
{
    fName  = name  ? name  : "MMyFindSuperCuts";
    fTitle = title ? title : "Optimizer of the supercuts";

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

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

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


    // Define columns of matrices
    *fLog << inf << "Init Mapping of Train Matrix" << endl;
    fCalcHadTrain->InitMapping(fMatrixTrain); 
    *fLog << inf << "Init Mapping of Test Matrix" << endl;
    fCalcHadTest->InitMapping(fMatrixTest);
}


// --------------------------------------------------------------------------
//
// Default destructor.
//
MMyFindSuperCuts::~MMyFindSuperCuts()
{
    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 MMyFindSuperCuts::DefineTrainMatrix(
			 const TString &filenametrain, MH3 &hreftrain,
	                 const Int_t howmanytrain, const TString &outfiletrain)
{
    if (filenametrain.IsNull() || howmanytrain <= 0)
        return kFALSE;

    *fLog << endl
	  << "***************************************************************"
	  << endl
	  << " Filling training matrix from file: '" << filenametrain << "'" 
	  << endl 
          << " Selecting: " << howmanytrain << " events";
   
    if (!hreftrain.GetHist().GetEntries()==0)
    {
	*fLog << ", according to a distribution given by the MH3 object '"
	      << hreftrain.GetName() << "'" << endl;
    }
    else
    {
	*fLog << ", Randomly" << endl;
    }

    *fLog << "***************************************************************"
	  << endl;


    //
    // ParList
    //
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);

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


    //
    // TaskList
    //
    MReadTree read("Events", filenametrain); 
    //MReadTree read("Parameters", filenametrain);
    read.DisableAutoScheme();

    //  MFEventSelector2 seltrain(hreftrain);
//      seltrain.SetNumMax(howmanytrain);
//      seltrain.SetName("selectTrain");

    // Random selection of events
    MFEventSelector seltrain; 
    seltrain.SetNumSelectEvts(howmanytrain);
    seltrain.SetName("selectTrain");

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

//      // ReCalculate hillass
//      MObservatory obs;
//      plist.AddToList(&obs);
    
//      MSrcPosCalc srccalc;
//      srccalc.SetPositionXY(3.46945e-17, 0.0487805);
//      MHillasSrcCalc hsrc;     
  

   
    // entries in MTaskList 
    tlist.AddToList(&read);
    tlist.AddToList(&seltrain);
   //   tlist.AddToList(&srccalc);
//      tlist.AddToList(&hsrc);
    tlist.AddToList(&filltrain);

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

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

  //    if (!evtloop.PreProcess())
//        return kFALSE;

//      Int_t i=0;
//      while(tlist.Process())
//      {
//  	i++;
//  	if (i>howmanytrain)
//  	    break;
//      }

//      evtloop.PostProcess();

    tlist.PrintStatistics(0, kTRUE);


    //
    // Print
    //
    fMatrixTrain->Print("SizeCols");
    Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
    if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MMyFindSuperCuts::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 (outfiletrain != "")
    {
      TFile filetr(outfiletrain, "RECREATE");
      fMatrixTrain->Write();
      filetr.Close();

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

    //   fMatrixTrain->Print("data");

    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 MMyFindSuperCuts::DefineTestMatrix(
			  const TString &nametest, MH3 &hreftest,
	                  const Int_t howmanytest, const TString &filetest)
{
    if (nametest.IsNull() || howmanytest<=0)
        return kFALSE;

    *fLog << "=============================================" << endl;
    *fLog << "fill test matrix from file '" << nametest 
          << "',   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;
    }


    //
    // ParList
    //
    MParList  plist;
    MTaskList tlist;

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


    //
    // TaskList
    //
    MReadTree read("Parameters", nametest);
    read.DisableAutoScheme();

   //   MFEventSelector2 seltest(hreftest);
//      seltest.SetNumMax(howmanytest);
//      seltest.SetName("selectTest");
 
    MFillH filltest(fMatrixTest);
    //filltest.SetFilter(&seltest);
    
    // entries in MTaskList 
    tlist.AddToList(&read);
    //tlist.AddToList(&seltest);
    tlist.AddToList(&filltest);


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

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

    Int_t i=0;
    while(tlist.Process())
    {
	i++;
	if (i>100000)
	    break;
    }

    evtloop.PostProcess();


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

    tlist.PrintStatistics(0, kTRUE);


    //
    // Print
    //
    fMatrixTest->Print("SizeCols");
    const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
    if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MMyFindSuperCuts::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 << "MMyFindSuperCuts::DefineTestMatrix; Test matrix was written onto file '"
            << filetest << "'" << endl;
    }


  return kTRUE;
}



// --------------------------------------------------------------------------
//
// Define the matrices for the training and test sample respectively
//
//
//
Bool_t MMyFindSuperCuts::DefineTrainTestMatrix(
			  const TString &name, MH3 &href,
	                  const Int_t howmanytrain, const Int_t howmanytest,
                          const TString &filetrain, const TString &filetest)
{
    *fLog << "=============================================" << endl;
    *fLog << "fill training and test matrix from file '" << name 
          << "',   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;
    }


    //
    // ParList
    //
    MParList  plist;
    MTaskList tlist;

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



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

    MFEventSelector2 selector(href);
    selector.SetNumMax(howmanytrain+howmanytest);
    selector.SetName("selectTrainTest");
    selector.SetInverted();

    MContinue cont(&selector);
    cont.SetName("ContTrainTest");

    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 MTaskList 
    tlist.AddToList(&read);
    tlist.AddToList(&cont);
    tlist.AddToList(&split);
    tlist.AddToList(&filltrain);
    tlist.AddToList(&conttrain);
    tlist.AddToList(&filltest);

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

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

    tlist.PrintStatistics(0, kTRUE);

    fMatrixTrain->Print("SizeCols");
    const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
    if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
    {
      *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; no.of generated events ("
	    << generatedtrain 
            << ") is incompatible with the no.of requested events ("
            << howmanytrain << ")" << endl;
    }

    fMatrixTest->Print("SizeCols");
    const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
    if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
    {
      *fLog << "MMyFindSuperCuts::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 << "MMyFindSuperCuts::DefineTrainTestMatrix; Training matrix was written onto file '"
            << filetrain << "'" << endl;
    }

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

      *fLog << "MMyFindSuperCuts::DefineTrainTestMatrix; Test matrix was written onto file '"
            << filetest << "'" << endl;
    }

  return kTRUE;
}



// --------------------------------------------------------------------------
//
//  Read training and test matrices from files
//
Bool_t MMyFindSuperCuts::ReadMatrix(const TString &filetrain, const TString &filetest)
{
    //
    // Read in training matrix
    //
    TFile filetr(filetrain);
    fMatrixTrain->Read("MatrixTrain");
    fMatrixTrain->Print("SizeCols");
    
    *fLog << "MMyFindSuperCuts::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 <<"MMyFindSuperCuts::ReadMatrix; Test matrix was read in from file '"
	  << filetest << "'" << endl;
    filete.Close();
    

    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 MMySuperCutsCalc will
//                       put the supercuts hadronness
//
// - for the minimization, the starting values of the parameters are taken  
//     - from file parSCinit (if it is != "")
//     - or from the arrays params and/or steps 
//
//----------------------------------------------------------------------
Bool_t MMyFindSuperCuts::FindParams(TString parSCinit,
                                     TArrayD &params, TArrayD &steps )
{

    fCalcHadTrain->SetSizeCuts(fSizeCutLow,fSizeCutUp);
    fCalcHadTest->SetSizeCuts(fSizeCutLow,fSizeCutUp);


    // Setup the event loop which will be executed in the 
    //                 fcnSupercuts function  of MINUIT
    //
    // parSCinit is the name of the file containing the initial values 
    // of the parameters; 
    // if parSCinit = ""   'params' and 'steps' are taken as initial values

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


    //
    // Check that all the necessary containers already exists
    //
    if (fHadronnessName.IsNull())
    {
      *fLog << "MyFindSuperCuts::FindParams; hadronness name is not defined... aborting"
            << endl;
      return kFALSE;
    }

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

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



    //
    // ParList
    //
    MParList  parlistfcn;
    MTaskList tasklistfcn;
    parlistfcn.AddToList(&tasklistfcn);
    
    
   
    //
    // create container for the supercut parameters
    // and set them to their initial values
    //
    MMySuperCuts super;

    // take initial values from file parSCinit
    if (parSCinit != "")
    {
	TFile inparam(parSCinit);
	super.Read("MMySuperCuts");
	inparam.Close();
	*fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from file " << parSCinit << endl;
	super.Print();
    }

    // take initial values from 'params' and/or 'steps'
    //else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
    if (params.GetSize() != 0  || steps.GetSize()  != 0 )
    {
      if (params.GetSize()  != 0)
      {
        *fLog << "MMyFindSuperCuts::FindParams; initial values of parameters are taken from 'params'" << endl;
        super.SetParameters(params);
      }
      if (steps.GetSize()  != 0)
      {
        *fLog << "MMyFindSuperCuts::FindParams; initial step sizes are taken from 'steps'"
              << endl;
        super.SetStepsizes(steps);
      }
    }
    else
    {
        *fLog << "MMyFindSuperCuts::FindParams; initial values and step sizes are taken from the MMySuperCuts constructor"
              << endl;
    }
 

    // Alpha binning
    const TString  mh3Name = "AlphaFcn";
    MBinning binsalpha("Binning"+mh3Name);
    //    binsalpha.SetEdges(54, -12.0, 96.0);
    binsalpha.SetEdges(45, 0.0, 90.0);

    // Alpha plot |alpha|
    MH3 alpha("MatrixTrain[6]");  // |alpha| assumed to be in column 6
        alpha.SetName(mh3Name);

    *fLog << warn << "WARNING----->ALPHA IS ASSUMED TO BE IN COLUMN 6!!!!!!!!"
	  << endl;


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

   
    // Add to ParList
    parlistfcn.AddToList(&super);
    parlistfcn.AddToList(fCam);
    parlistfcn.AddToList(fMatrixTrain);
    parlistfcn.AddToList(&binsalpha);
    parlistfcn.AddToList(&alpha);
    parlistfcn.AddToList(&plotsuper);



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

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

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

    // Fill Alpha Plot
    MFillH fillalpha(&alpha);

   
    // Add to TaskList
    tasklistfcn.AddToList(&loop);
    tasklistfcn.AddToList(fCalcHadTrain);
    tasklistfcn.AddToList(&supercuts);
    tasklistfcn.AddToList(&fillalpha);


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

  


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

    // get initial values of parameters 
    fVinit = super.GetParameters();
    fStep  = super.GetStepsizes();

    TString name[fVinit.GetSize()];
    fStep.Set(fVinit.GetSize());
    fLimlo.Set(fVinit.GetSize());
    fLimup.Set(fVinit.GetSize());
    fFix.Set(fVinit.GetSize());

    fNpar = fVinit.GetSize();



    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; 
	fLimlo[i] = 0.0;
        fLimup[i] = 0.0;
        fFix[i]   =   0;
    }

    // these parameters make no sense, fix them at 0.0
    //  for (UInt_t i=32; i<fNpar; i++)
//      {
//  	//if( i!=3 && i!=7 && i!=11 && i!=15 && i!=19 && i!=23 )
//  	{

//  	    //fVinit[i] = 0.0;
//  	    fStep[i]  = 0.0;
//  	    fFix[i]   = 1;
//  	}
//  	//   fVinit[3+4*i] = 0.0;
//  	//  	    fStep[3+4*i]  = 0.0;
//  	//  	    fFix[3+4*i]   = 1;
	
//      }
    // these parameters make no sense, fix them at 0.0
   //   for (UInt_t i=0; i<24; i++)
//      {
//  	//fVinit[i] = 0.0;
//  	fStep[i]  = 0.0;
//  	fFix[i]   = 1;
//      }


    //-----------------------------------------------------------------
    //
    // 0=supercuts, 1=dyn from supercuts, 2=supercuts, 3=all
    //
    gLog << "*********************************" <<   fOptimizationMode << endl;
   
    switch (fOptimizationMode)
    {
	case 0:   
	    for (UInt_t i=0; i<fNpar; i++)
	    {
		if(i%4==0)      // Only optimaze par0, par4, par8,...
		    continue;
		
		fStep[i]  = 0.0;
		fFix[i]   = 1;
	    }	
	    break;

	case 1:   
	    for (UInt_t i=0; i<fNpar; i++) // Optimize all but par0, par4, par8
	    {
		if(i%4!=0 && i!=19 && i!=23) 
		    continue;

		fStep[i]  = 0.0;
		fFix[i]   = 1;

	    }	

	    break;
    }

    //
    // DistLo/Up doens't dipend on Dist2
    // 
    fStep[19]  = 0.0;
    fFix[19]   = 1;
    fStep[23]  = 0.0;
    fFix[23]   = 1;



    //
    // Fix all parameters from 24 to the end
    //
    for (UInt_t i=24; i<fNpar; i++)
    {
 	//fVinit[i] = 0.0;
 	fStep[i]  = 0.0;
 	fFix[i]   = 1;
    }

    
  //    // Static Cuts
//      //
//      for (UInt_t i=0; i<fNpar; i++)
//      { 
// 	 if( i!=0 &&  i!=4 && i!=8 &&  i!=12  &&  i!=16 &&  i!=20 && 
//   	     i!=24 &&  i!=28 && i!=32 &&  i!=36  &&  i!=40  ) 
//  	{
//  	    fStep[i]  = 0.0;
//  	    fFix[i]   = 1;
//  	}
//      }

  //
  //    // 
    //
   //   for (UInt_t i=0; i<fNpar; i++)
//      {
//  	if( i==0 || i==4 || i==8 ||  i==12  ||  i==16 ||  i==20) 
//  	{
//  	    //fVinit[i] = 0.0;
//  	    fStep[i]  = 0.0;
//  	    fFix[i]   = 1;
//  	}
//      }

 
    // Fix
    //  for (UInt_t i=0; i<36; i++)
//      {
//  	    fStep[i]  = 0.0;
//  	    fFix[i]   = 1;	
//      } 
  //   for (UInt_t i=44; i<fNpar; i++)
//     {
// 	    fStep[i]  = 0.0;
// 	    fFix[i]   = 1;	
//     }

//   for (UInt_t i=37; i<=39; i++)
//      {
//  	    //fVinit[i] = 0.0;
//  	    fStep[i]  = 0.0;
//  	    fFix[i]   = 1;	
//      }


    //      fStep[40]  = .01;  
    //          fStep[28]  = 10; 
    //  fStep[24]  = 10;  
    //     fVinit[24] = 100; 
    //    fVinit[28] = -74;
     //    fVinit[32] = 0.16;
//         fVinit[36] = 0.03;
    //         fVinit[40] = 0.21;
//         fVinit[44] = 0.0;



  //fVinit[28] = 0;
   //   fVinit[40] = 0.5;

//        fVinit[44] = -1e6;
      //       fStep[44]  = 1;
//      fVinit[36] = 0.;
//      //fStep[36]  = 1.0;


    // vary only first 48 parameters
    //for (UInt_t i=0; i<fNpar; i++)
    //{
    //    if (i >= 48)
    //	{
    //      fStep[i] = 0.0;
    //      fFix[i]  =   1;
    //	}
    //

   //  //}
//     // vary only first 23 parameters
//     //
//     for (UInt_t i=0; i<fNpar; i++)
//     {
//        if (i >= 24)
//     	{
// 	    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(fcnSuperCuts, name,
                                 fVinit, fStep, fLimlo, fLimup, fFix,
                                 &evtloopfcn, "SIMPLEX", kFALSE);
 
    *fLog << "after calling CallMinuit" << endl;

    *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 MMySuperCuts

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

    const TArrayD &check = super.GetParameters();
    for (Int_t i=0; i<check.GetSize(); 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 MMyFindSuperCuts::TestParams()
{
    if (fMatrixTest->GetM().GetNrows() <= 0)
    {
        *fLog << "MMyFindSuperCuts::TestParams; test matrix has no entries" 
              << 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);
    MMySuperCuts scin; 
    scin.Read("MMySuperCuts");
    inparam.Close();

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

    const TArrayD &supercutsPar = scin.GetParameters();
    for (Int_t i=0; i<supercutsPar.GetSize(); i++)
        *fLog << supercutsPar[i] << ",  ";
    *fLog << endl;
    //---------------------------


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

    MParList  parlist2;
    MTaskList tasklist2;

    MMySuperCuts supercuts;
    supercuts.SetParameters(supercutsPar);

    fCalcHadTest->SetHadronnessName(fHadronnessName);


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

    MMatrixLoop loop(fMatrixTest);

    // plot alpha before applying the supercuts
    const 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[6]");
    alphabefore.SetName(mh3NameB);

    TH1 &alphahistb = alphabefore.GetHist();
    alphahistb.SetName("alphaBefore-TestParams");

    MFillH fillalphabefore(&alphabefore);
    fillalphabefore.SetName("FillAlphaBefore");

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

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

    MH3 alphaafter("MatrixTest[6]");
    alphaafter.SetName(mh3NameA);

    TH1 &alphahista = alphaafter.GetHist();
    alphahista.SetName("alphaAfter-TestParams");

    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]");
    alphaHist2.SetName("alpha-TestParams");

    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
  
    const Double_t alphasig = 20.0;
    const Double_t alphamin = 30.0;
    const Double_t alphamax = 90.0;
    const Int_t    degree   =    2;
    const Bool_t   drawpoly  = kTRUE;
    const Bool_t   fitgauss  = kTRUE;
    const Bool_t   print     = kTRUE;

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

    findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree, 
                      alphasig, drawpoly, fitgauss, print);

    const Double_t significance = findsig.GetSignificance();
    const 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;
}

