//
//
/*
  class SobIter:
  root implementation of the S/B Iterative method
  for selection tuning

  Version 2.0 Released 30 April 09

  Further Information/Inquiries:
  Nikolaos.Rompotis@Cern.ch

  Nikolaos Rompotis - 23 March 2009 - initial release
  26 March 09: Detector specific tuning implemented
  03 April 09: bug corrected in TChain implementation that was stopping the
               code when more than one files for bkg were specified
  30 April 09: version 2.0 This implementation gives identical results with
               the previous version, but has some more flexibility in the
               definition of signal and bkg. Instead of giving a signal and
               bkg files a signal filter and a bkg filter is defined in the
               input configuration and the electrons are read from a single
               dataset
  11 May 09:   version 2.1 enables to use an artificially lower cut in any
               variable that we want
  22 June 09:  version 3.1  bug  is corrected in sig/bkg mc values
               code extended to receive many different types as input
  22 Sept 09:  modification for the october exersize
  15 Marc 10:  met type chosen in cfg, preselection option in cfg
               some refinements in order to be able to calculate the S/B init
               old cfg still work, but the numbers may not have the same
               interpretation
               TO BE DONE: proper integration of Zee tuning, inverted cuts
*/

#include <iostream>
#include <vector>
#include "TString.h"
#include "TTree.h"
#include "TBranch.h"
#include "TFile.h"
#include "TGraph.h"
#include "TMath.h"
#include "TH1F.h"
#include "TChain.h"
#include "TDatime.h"
#include "TCanvas.h"
#include "TROOT.h"


using namespace std;

struct Filter {
  // this structure allows the application of a filter to an event
  // the filter contains ET cut, MET cut and particular type of event cut
  // see function PassFilter in SobIter class
  double ET;
  bool ET_LESS;

  double MET;
  bool MET_LESS;

  vector<Short_t> TYPE;
  bool TYPE_EQUAL;

  bool EXCLUDE;
  TString NAME;
};

class SobIter
{
public:
  SobIter( double step, bool relative, double lowest_eff,
	   double min_step, 
	   vector<TString> barrel_vars, vector<TString> endcaps_vars,
	   vector<TString> barrel_names, vector<TString> endcaps_names,
	   vector<double> barrel_initcuts,  vector<double> endcaps_initcuts, 
	   vector<double> barrel_mins,  vector<double> endcaps_mins,
	   vector<int> barrel_points,   vector<int> endcaps_points,
	   double *weights, int nw,  vector<TString> files, 
	   Filter signalfilter, Filter bkgfilter,
	   vector<double> lowerBarrel, vector<double> lowerEndcaps,
	   TString metType, int preSelType,
	   vector< vector<double> > selections, vector<double> labels);
  void Analysis();
private:
  bool RELATIVE_STEP_;
  double LOWEST_EFF_;
  vector<double> InitCuts_barrel_;
  vector<double> InitCuts_endcaps_;
  vector<double> Mins_barrel_;
  vector<double> Mins_endcaps_;
  vector<int> Points_barrel_; 
  vector<int> Points_endcaps_; 
  vector<TString> Branches_barrel_;
  vector<TString> Branches_endcaps_;
  vector<TString> Branches_;
  vector<int> BranchToVariable_;
  vector<TString> Names_barrel_;
  vector<TString> Names_endcaps_;
  vector<double> LowerAllowedBarrel_;
  vector<double> LowerAllowedEndcaps_;


  vector<TString> DATAFILES_;
  double minStep_;
  double *w_;
  int Nw_;
  double Step_;
  bool isAlive_;
  Filter SignalFilter_;
  Filter BkgFilter_;
  //
  //vector<double> CurrentCutsBarrel_;
  //vector<double> CurrentCutsEndcaps_;
  int n_barrel_;
  int n_endcaps_;
  int n_variables_;
  TString metType_;
  bool UseExtraPreselection_;
  int preSelType_;
  //
  //  bool CheckCuts( double *values, double eta);
  bool PassFilter(Filter myfilter, double et, double met, 
		  short int type, short int exclusion);
  void  PrintFilter(Filter filter);
  int FindName(TString name, vector<TString> vec);
  vector<double>  ReturnSeffForSoB
  (double target, TH1F* h_signal, TH1F* h_bkg, double rest_signal, 
   double rest_bkg, double SIGNAL_INIT, double tolerance);
  double SobIter::Bisection
  (TGraph *g, double a_, double b_, double target, double tolerance);
  // extra
  vector<double> makeRoundCuts(vector<double> vec, vector<TString> names);
  double RoundNumber(double n);
  vector< vector<double> > CurrentCutsBarrel_;
  vector< vector<double> > CurrentCutsEndcaps_;
  vector< vector<double> > CurrentCutsBarrelRound_;
  vector< vector<double> > CurrentCutsEndcapsRound_;
  vector < vector<double> > selections_;
  vector<double> labels_;
  bool CheckCutsRound( double *values, double eta, int j);
  bool CheckCuts( double *values, double eta, int j);
};


SobIter::SobIter( double step, bool relative, double lowest_eff,
	   double min_step, 
	   vector<TString> barrel_vars, vector<TString> endcaps_vars,
 	   vector<TString> barrel_names, vector<TString> endcaps_names,
	   vector<double> barrel_initcuts,  vector<double> endcaps_initcuts, 
	   vector<double> barrel_mins,  vector<double> endcaps_mins,
	   vector<int> barrel_points,   vector<int> endcaps_points,
	   double *weights, int nw,  vector<TString> files, 
	   Filter signalfilter, Filter bkgfilter, 
	   vector<double> lowerBarrel, vector<double> lowerEndcaps,
	   TString metType, int preSelType,
	   vector< vector<double> > selections, vector<double> labels)
{
  isAlive_ =true;
  selections_ = selections; labels_ = labels;
  // check for consistency
  int n = (int) barrel_vars.size();
  if (n != (int) barrel_initcuts.size() || n != (int) barrel_mins.size()
      || n != (int)  barrel_points.size() 
      || n != (int) barrel_names.size() ) {
    cout << "Inconsistent Number: Barrel Variables " << endl;
    cout << "branches: " << n << endl;
    cout << "init cuts: " << (int) barrel_initcuts.size() << endl;
    cout << "min cuts: " << (int) barrel_mins.size() << endl;
    cout << "points cuts: " << (int) barrel_points.size() << endl;
    isAlive_ = false;
    return;
  }
  //
  int n1 = (int) endcaps_vars.size();
  if (n1 != (int) endcaps_initcuts.size() || n1 != (int) endcaps_mins.size()
      || n1 != (int)  endcaps_points.size() 
      || n1 != (int) endcaps_names.size()) {
    cout << "Inconsistent Number: Endcap Variables " << endl;
    cout << "branches: " << n1 << endl;
    cout << "init cuts: " << (int) endcaps_initcuts.size() << endl;
    cout << "min cuts: " << (int) endcaps_mins.size() << endl;
    cout << "points cuts: " << (int) endcaps_points.size() << endl;

    isAlive_ = false;
    return;
  }
  //
  if ( (int) files.size() == 0){
    cout << "You must have at least one sample" << endl;
    isAlive_ = false;
    return;
  }
  if (nw == 0) {
    cout << "You should insert at least one weight" << endl;
    isAlive_ = false;
    return;
  }
  if (n1 == 0 && n ==0) {
    cout << "You should have at least one variable" << endl;
    isAlive_ = false;
    return;
  }
  InitCuts_barrel_ = barrel_initcuts;
  InitCuts_endcaps_ = endcaps_initcuts;
  Mins_barrel_ = barrel_mins;
  Mins_endcaps_ = endcaps_mins;
  Points_barrel_ = barrel_points;
  Points_endcaps_ = endcaps_points;
  Branches_barrel_ = barrel_vars;
  Branches_endcaps_ = endcaps_vars;
  Names_barrel_ = barrel_names;
  Names_endcaps_ = endcaps_names;
  LowerAllowedBarrel_ = lowerBarrel;
  LowerAllowedEndcaps_ = lowerEndcaps;
  metType_ = metType;
  preSelType_ = preSelType;
  if (preSelType_ > 0) UseExtraPreselection_ = true;
  else  UseExtraPreselection_ = false;
  w_ = weights;
  Nw_ = nw;
  DATAFILES_ = files;
  Step_ = step;
  RELATIVE_STEP_ = relative;
  LOWEST_EFF_ = lowest_eff;
  minStep_ = min_step;
  n_barrel_ = (int) Names_barrel_.size();
  n_endcaps_ = (int) Names_endcaps_.size();
  /// now find the common set of branches
  for (int i=0; i< n_barrel_; ++i) {
    Branches_.push_back(Branches_barrel_[i]);
    BranchToVariable_.push_back(i);
  }
  int counter = 0;
  for (int i=0; i<n_endcaps_; ++i) {
    // is branch endcaps i the same as one of the barrel branches?
    int ind = FindName(Branches_endcaps_[i],Branches_barrel_);
    if (ind == -1) {
      Branches_.push_back(Branches_endcaps_[i]);
      BranchToVariable_.push_back(n_barrel_ + counter);
      ++counter;
      cout << "didn't find " << Branches_endcaps_[i] << " asign ind: "
	   << n_barrel_ + counter << endl;
    }
    else
      BranchToVariable_.push_back(ind);
  }
  n_variables_ = (int)  Branches_.size();
  SignalFilter_ = signalfilter;
  BkgFilter_ = bkgfilter;
  PrintFilter(SignalFilter_);
  PrintFilter(BkgFilter_);
  cout << "SobIter object was created..." << endl;

}

void SobIter::Analysis()
{
  if (not isAlive_) return;
  TString relative = "yes";
  if (not RELATIVE_STEP_) relative = "no";
  cout << "==============================" << endl;
  cout << "Iterative Method starts ...." << endl;
  cout << "Step = " << Step_ << ", relative? "<< relative
       << " Min step = " << minStep_       << endl;
  cout << "Variables and initial values: " << endl;
  cout << "BARREL: " << endl;
  for (int i=0; i<n_barrel_; ++i) {
    cout << i << ". "<<Names_barrel_[i] << ": "<<InitCuts_barrel_[i] << endl;
  }
  cout << "ENDCAPS: " << endl;
  for (int i=0; i<n_endcaps_; ++i) {
    cout << i << ". "<<Names_endcaps_[i] << ": " <<InitCuts_endcaps_[i]<<endl;
  }
  cout << "WEIGHTS IN USE: " << endl;
  for (int i=0; i<Nw_; ++i) cout << "w_[" << i << "]=" << w_[i] << endl;
  cout << "==============================" << endl;
  //
  //
  // Define your Current Cut Values:
  //for (int i=0; i<n_barrel_; ++i)
  //  CurrentCutsBarrel_.push_back(InitCuts_barrel_[i]);
  //for (int i=0; i<n_endcaps_; ++i)
  //  CurrentCutsEndcaps_.push_back(InitCuts_endcaps_[i]);
  // make the rounded cuts
  for (int j=0; j< (int) selections_.size(); ++j) {
    vector<double> CurrentCutsBarrel;
    vector<double> CurrentCutsEndcaps;
    for (int i=0; i<n_barrel_; ++i) {
      CurrentCutsBarrel.push_back(selections_[j][i]);
    }
    for (int i=0; i<n_endcaps_; ++i)
      CurrentCutsEndcaps.push_back(selections_[j][i+n_barrel_]);
    // make the rounded cuts
    vector<double> CurrentCutsBarrelRound  = makeRoundCuts(CurrentCutsBarrel,
							   Names_barrel_);
    vector<double> CurrentCutsEndcapsRound = makeRoundCuts(CurrentCutsEndcaps,
							   Names_endcaps_);
    CurrentCutsBarrel_.push_back(CurrentCutsBarrel);
    CurrentCutsEndcaps_.push_back(CurrentCutsEndcaps);

    CurrentCutsBarrelRound_.push_back(CurrentCutsBarrelRound);
    CurrentCutsEndcapsRound_.push_back(CurrentCutsEndcapsRound);
  }


  //
  //
  // for book keeping
  vector<double> SOB;   vector<double> SOB_MC;
  vector<double> SEFF;  vector<double> SEFF_MC;
  vector<double> SOB_ERROR;  vector<double> SOB_MC_ERROR;
  vector<TString> SELECTION;
  //
  // prepare the input files
  TChain  tree("probe_tree");
  vector<TString>::iterator it = DATAFILES_.begin();
  for ( ; it !=  DATAFILES_.end(); ++it) {
    tree.Add(*it); 
  }
  
  //
  // 
  const int N_signal = tree.GetEntries();
  //cout << "Signal sample contains " << N_signal << " electrons" << endl;
  //
  // / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 
  /////////////////////////////////////////////////////////////////
  //                                         //////////////////////
  //   Sig/Bkg Branches   ...........................//////////////
  // ..................................................////////////
  // Common branch:
  double probe_sc_eta;  // eta of the supercluster
  double probe_sc_et;
  double event_MET;
  short int type;  // type of the event
  short int exclusion; // ==1 when electron belongs to event with 2 elecs
  //
  tree.SetBranchAddress("probe_sc_eta", &probe_sc_eta);  
  tree.SetBranchAddress("probe_sc_et", &probe_sc_et);  
  // met default case:
  // tree.SetBranchAddress("event_MET", &event_MET);  
  tree.SetBranchAddress(metType_, &event_MET);  
  // extra preselection
  // preSelType_ = 0   nothing extra
  //             = 1   Valid first PXB
  //             = 2   Dcot Dist - prechosen values
  //             = 3   expected missing hits (Not yet implemented)
  //             = 12  both 1 and 2
  // to be replaced with biwise operators
  // 
  int probe_ele_validHitInFirstPixelBarrel;
  double probe_ele_dist, probe_ele_dcot;
  if (UseExtraPreselection_) {
    if (preSelType_ == 1 || preSelType_ == 12)
      tree.SetBranchAddress("probe_ele_validHitInFirstPixelBarrel",
			    &probe_ele_validHitInFirstPixelBarrel);
    if (preSelType_ == 2 || preSelType_ == 12) {
      tree.SetBranchAddress("probe_ele_dist", &probe_ele_dist);
      tree.SetBranchAddress("probe_ele_dcot", &probe_ele_dcot);
    }
  }
  tree.SetBranchAddress("type", &type);  
  tree.SetBranchAddress("exclude", &exclusion);  
  //
  // Signal and Bkg branches
  double *vars = new double[n_variables_];
  for (int i =0; i< n_variables_; ++i) { 
    tree.SetBranchAddress(Branches_[i], &vars[i]);
  }
  //////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////
  ///////////////////////////////////////////
  ////////////////////////////////////
  ////////////////////////
  //
  //   ITERATION LOOP
  //   ^^^^^^^^^^^^^^
  //
  //
  double signal_init = 0;    // needed for the efficiency definition
  //double signal_init_barrel = 0;
  //double signal_init_endcaps = 0;
  //double signal_init_mc = 0;
  //double bkg_init_mc = 0;
  //double bkg_init    = 0;
  //double signal_init2 = 0, bkg_init2 = 0;
  //double signal_init_mc2 = 0, bkg_init_mc2 = 0;
  //double sob_previous=-1;
  //double sobInit = 0, sobInitMC = 0;
  //double sobInitError = 0, sobInitMCError = 0;
  // .....................
  //for (int Iter = 0; Iter<100000; ++Iter) {  // S/B Iteration loop
  // define your histograms .........................................
  // barrel
  //vector <TH1F> h_barrel; // signal
  //vector <TH1F> b_barrel; // bkg
  //for (int i=0; i<n_barrel_; ++i) {
  //  TString name = "h_";
  //  name += i;
  //  TH1F h(name, name, Points_barrel_[i], Mins_barrel_[i], 
  //	   CurrentCutsBarrel_[i]);
  //  h_barrel.push_back(h);
  //  TString b_name = "b_";
  //  b_name += i;
  //  TH1F b(b_name, b_name, Points_barrel_[i], Mins_barrel_[i], 
  //	   CurrentCutsBarrel_[i]);
  // b_barrel.push_back(b);
  //}
  // endcaps
  //vector <TH1F> h_endcaps; // signal
  //vector <TH1F> b_endcaps; // bkg
  //for (int i=0; i<n_endcaps_; ++i) {
  //  TString name = "h_ee_";
  //  name += i;
  //  TH1F h(name, name, Points_endcaps_[i], Mins_endcaps_[i], 
  //	   CurrentCutsEndcaps_[i]);
  //  h_endcaps.push_back(h);
  //  TString b_name = "b_ee_";
  //  b_name += i;
  //  TH1F b(b_name, b_name, Points_endcaps_[i], Mins_endcaps_[i], 
  //	   CurrentCutsEndcaps_[i]);
  //  b_endcaps.push_back(b);
  //}
  //double current_signal = 0;
  //double current_bkg = 0;
  //double current_signal_barrel = 0;    double current_bkg_barrel = 0;
  //double current_signal_endcaps = 0;   double current_bkg_endcaps = 0;
  //double current_actual_signal = 0; // assuming type 0 is the signal
  //double current_actual_bkg = 0; // assuming type!= 0 is the bkg
  //double current_signal_in_sigSample_eb = 0;
  //double current_signal_in_sigSample_ee = 0;
  //double current_signal_in_bkgSample_eb = 0;
  //double current_signal_in_bkgSample_ee = 0;
  // statistical errors on s/b
  //double dS2 = 0. , dB2 = 0.;
  //double dS2_mc = 0. , dB2_mc = 0.;
  // ROUNDING specific ................................................
  // initializations: for each selection you keep the S/B and the seff
  int number_of_selections = (int) selections_.size();
  //  double signal_init = 0;
  double *old_S  = new double[number_of_selections];
  double *old_B = new double[number_of_selections];
  //
  double *round_S  = new double[number_of_selections];
  double *round_B = new double[number_of_selections];
  //
  // initial value to 0
  for (int i=0; i< number_of_selections; ++i) {
    old_S[i] = 0;    old_B[i] = 0;    round_S[i] = 0;    round_B[i] = 0;
  }
  // ...................................................................
  //
  // LOOP OVER THE TREE ENTRIES ........................................
  //
  // Loop over the signal 
  for (int iT = 0; iT < N_signal; ++iT) { // loop over the signal
    //if (iT%100 == 0) 
    //cout << "Event # " << iT  << endl;
    tree.GetEntry(iT);
    //    bool pass = CheckCuts(vars, probe_sc_eta);
    bool pass = true;
    if (UseExtraPreselection_) {
      if (preSelType_ == 1)
	pass = pass && probe_ele_validHitInFirstPixelBarrel == 1;
      if (preSelType_ >= 2) {
	bool pass_ConvPartner = true;
	if (probe_ele_dcot < 990) {
	  pass_ConvPartner = fabs(probe_ele_dcot) > 0.02 &&
	    fabs(probe_ele_dist) > 0.02;
	}
	if (preSelType_ == 2)
	  pass = pass && pass_ConvPartner;
	else if (preSelType_ == 12) {
	  pass = pass && pass_ConvPartner && 
	    probe_ele_validHitInFirstPixelBarrel == 1;
	}
      }
    }
    bool in_barrel = fabs(probe_sc_eta) < 1.479&& n_barrel_>0;
    bool in_endcaps = fabs(probe_sc_eta) > 1.479&& n_endcaps_>0;
    bool somewhere = in_barrel || in_endcaps;
    //
    ///////// this is signal related
    if (PassFilter(SignalFilter_, probe_sc_et, event_MET, type, exclusion)) {
      if (somewhere) {
	signal_init += w_[type];
	// now let us check the selections
	for (int jSel =0 ; jSel < number_of_selections; ++jSel) {
	  pass = pass && CheckCuts(vars, probe_sc_eta,jSel);
	  bool passRound = pass && CheckCutsRound(vars, probe_sc_eta,jSel);
	  if (pass) {
	    old_S[jSel] += w_[type];
	  }
	  if (passRound) {
	    round_S[jSel] += w_[type]; 
	  }
	}
      }
    }
    //// bkg related ......................................................
    if (PassFilter(BkgFilter_, probe_sc_et, event_MET, type, exclusion)) {
      if (somewhere) {
	for (int jSel =0 ; jSel < number_of_selections; ++jSel) {
	  pass = pass && CheckCuts(vars, probe_sc_eta,jSel);
	  bool passRound = pass && CheckCutsRound(vars, probe_sc_eta,jSel);
	  if (pass) {
	    old_B[jSel] += w_[type];
	  }
	  if (passRound) {
	    round_B[jSel] += w_[type]; 
	  }
	}
      }
    }
    //////////////////////////////////////////////////////////////////////
  } // end loop over the samples
  
  //
  // now we have to print out the output
  //
  cout << "----------------- DEBUGGING --------------------------" << endl;
  cout << "For Debugging purposes I write down the original seffs" << endl;
  cout << "Sould be same for original and the old selection" << endl;
  for (int i=0; i<number_of_selections; ++i) 
    cout << labels_[i] << " - " << old_S[i]/signal_init << endl;
  cout << "----------------- DEBUGGING --------------------------" << endl;
  //
  cout << "Comparison Rounded and UnRounded Values" << endl;
  cout << "---------------------------------------" << endl;
  cout << " UnRounded (seff-s/b)       Rounded (seff-s/b)\n ";
  cout << "BARREL   ";
  for (int j=0; j<number_of_selections; ++j) {
    printf(" (%4.3f-%5.2f)  (%4.3f-%5.2f) ", 
	   old_S[j]/signal_init, old_S[j]/old_B[j],
	   round_S[j]/signal_init, round_S[j]/round_B[j]);
  }
  printf("\n");
  for (int i=0; i<n_barrel_; ++i) {
    printf("%8s ", Names_barrel_[i].Data());
    for (int j=0; j<number_of_selections; ++j ) {
      if (Names_barrel_[i].Contains("Tdist") || 
	  Names_barrel_[i].Contains("Tdcot")) {
	printf("   %10.6f     %3.1e    |",  
	   (1.-CurrentCutsBarrel_[j][i])/CurrentCutsBarrel_[j][i], 
	   (1.-CurrentCutsBarrelRound_[j][i])/CurrentCutsBarrelRound_[j][i]);
      } else {
	printf("   %10.6f     %3.1e    |",  
	       CurrentCutsBarrel_[j][i], CurrentCutsBarrelRound_[j][i]);
      }
    }
    printf("\n");
  }
  cout << "ENDCAPS" << endl;
  for (int i=0; i<n_endcaps_; ++i) {
    printf("%8s ", Names_endcaps_[i].Data());
    for (int j=0; j<number_of_selections; ++j ) {
      if (Names_endcaps_[i].Contains("Tdist") || 
	  Names_endcaps_[i].Contains("Tdcot")) {
	printf("   %10.6f     %3.1e    |",  
	   (1.-CurrentCutsEndcaps_[j][i])/CurrentCutsEndcaps_[j][i], 
	   (1.-CurrentCutsEndcapsRound_[j][i])/CurrentCutsEndcapsRound_[j][i]);
      } else {
	printf("   %10.6f     %3.1e    |",  
	       CurrentCutsEndcaps_[j][i], CurrentCutsEndcapsRound_[j][i]);
      }      
    }
    printf("\n");
  }
  // print out in table format the values only
  cout << "Rounded Values only " << endl;
  cout << "-------------------"  << endl;
  cout << "BARREL   ";
  for (int j=0; j<number_of_selections; ++j) {
    printf("    %4.3f   ", round_S[j]/signal_init);
  }  
  printf("\n");
  for (int i=0; i<n_barrel_; ++i) {
    printf("%8s ", Names_barrel_[i].Data());
    for (int j=0; j<number_of_selections; ++j ) {
      if (Names_barrel_[i].Contains("Tdist") || 
	  Names_barrel_[i].Contains("Tdcot")) {
	printf("  %3.1e  |", 
	  (1-CurrentCutsBarrelRound_[j][i])/CurrentCutsBarrelRound_[j][i]);
      } else {
	printf("  %3.1e  |", CurrentCutsBarrelRound_[j][i]);
      }
    }
    printf("\n");
  }
  cout << "ENDCAPS" << endl;
  for (int i=0; i<n_endcaps_; ++i) {
    printf("%8s ", Names_endcaps_[i].Data());
    for (int j=0; j<number_of_selections; ++j ) {
      if (Names_endcaps_[i].Contains("Tdist") || 
	  Names_endcaps_[i].Contains("Tdcot")) {
	printf("  %3.1e  |", 
	  (1-CurrentCutsEndcapsRound_[j][i])/CurrentCutsEndcapsRound_[j][i]);
      } else {
	printf("  %3.1e  |",  CurrentCutsEndcapsRound_[j][i]);
      }
    }
    printf("\n");
  }


  // standard iterative output format
  // printing only the rounded values
  cout << "******************************************************" << endl;
  cout << "RESULTS WITH THE STANDARD ITERATIVE ALGO OUTPUT FORMAT" << endl;
  cout << "******************************************************" << endl;
  const int Niter = number_of_selections;
  cout << "Tuning Summary: " << endl;
  cout << "Cuts appear in the following order (BARREL/ENDCAPS): " << endl;
  for (int i=0; i<n_barrel_-1; ++i)
    cout << Names_barrel_[i] << ",";
  if (n_barrel_>0)
    cout << Names_barrel_[n_barrel_-1] << " / ";
  else
    cout << "(no barrel variables)" << " / ";
  for (int i=0; i<n_endcaps_-1; ++i)
    cout << Names_endcaps_[i] << ",";
  if (n_endcaps_>0)
    cout << Names_endcaps_[n_endcaps_-1] << endl;
  else
    cout << "(no endcap variables)" << endl;
  cout << "double sob[" << Niter << "], seff[" << Niter << "], sob_mc["
       << Niter << "], seff_mc[" << Niter << "]," << "sob_error["
       << Niter << "], sob_mc_error[" << Niter << "];"   << endl;
  for (int i=0;i<Niter; ++i) {
    double SOB = round_S[i]/round_B[i];
    double SEFF = round_S[i]/signal_init;
    // calculate the selection
    TString selection("");
    for (int ii=0; ii<n_barrel_-1; ++ii) {
     selection +=  CurrentCutsBarrelRound_[i][ii];
     selection += ", ";
    }
    if (n_barrel_>0)
     selection +=  CurrentCutsBarrelRound_[i][n_barrel_-1];
    else
     selection +=  "(no barrel variables)";
    selection += " // ";
    for (int ii=0; ii<n_endcaps_-1; ++ii) {
     selection += CurrentCutsEndcapsRound_[i][ii];
     selection += ", ";
    }
    if (n_endcaps_>0)
     selection += CurrentCutsEndcapsRound_[i][n_endcaps_-1];
    else
     selection += "(no endcap variables)";
    selection += ")";
    // end selection calculation
    cout << "sob[" << i << "]=" << SOB << "; seff[" << i << "]=" << SEFF
	 << "; sob_mc[" << i << "]=" << SOB
	 << "; seff_mc[" << i << "]=" << SEFF
	 << "; sob_error[" << i << "]=" << 0
	 << ";  sob_mc_error[" << i << "]=" << 0
	 << ";     //  "  << selection << endl;
  }



  delete [] old_S;    delete [] old_B;
  delete [] round_S;  delete [] round_B;


  delete [] vars;

}

/*
bool SobIter::CheckCuts( double *values, double eta)
{
  bool res = true;
  if (fabs(eta) < 1.479) {
    for (int i=0; i< n_barrel_; ++i) {
      bool cut = (fabs(values[i])< CurrentCutsBarrel_[i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_barrel_[i] << " val=" << values[i]
	//     << " current cut =" << CurrentCutsBarrel_[i]
	//     << endl;
	break;
      }
    }
  }
  else {
    for (int i=0; i< n_endcaps_; ++i) {
      bool cut = (fabs(values[ BranchToVariable_[i+n_barrel_] ])< 
		  CurrentCutsEndcaps_[i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_endcaps_[i] << " val=" 
	//     << fabs(values[ BranchToVariable_[i+n_barrel_] ])
	//     << " current cut =" << CurrentCutsEndcaps_[i]
	//    << endl;
	break;
      }
    }
  }
  return res;

}
*/

bool SobIter::CheckCuts( double *values, double eta, int j)
{
  bool res = true;
  if (fabs(eta) < 1.479) {
    for (int i=0; i< n_barrel_; ++i) {
      bool cut = (fabs(values[i])< CurrentCutsBarrel_[j][i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_barrel_[i] << " val=" << values[i]
	//     << " current cut =" << CurrentCutsBarrel_[i]
	//     << endl;
	break;
      }
    }
  }
  else {
    for (int i=0; i< n_endcaps_; ++i) {
      bool cut = (fabs(values[ BranchToVariable_[i+n_barrel_] ])< 
		  CurrentCutsEndcaps_[j][i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_endcaps_[i] << " val=" 
	//     << fabs(values[ BranchToVariable_[i+n_barrel_] ])
	//     << " current cut =" << CurrentCutsEndcaps_[i]
	//    << endl;
	break;
      }
    }
  }
  return res;

}



bool SobIter::CheckCutsRound( double *values, double eta, int j)
{
  bool res = true;
  if (fabs(eta) < 1.479) {
    for (int i=0; i< n_barrel_; ++i) {
      bool cut = (fabs(values[i])< CurrentCutsBarrelRound_[j][i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_barrel_[i] << " val=" << values[i]
	//     << " current cut =" << CurrentCutsBarrel_[i]
	//     << endl;
	break;
      }
    }
  }
  else {
    for (int i=0; i< n_endcaps_; ++i) {
      bool cut = (fabs(values[ BranchToVariable_[i+n_barrel_] ])< 
		  CurrentCutsEndcapsRound_[j][i]);
      res = res && cut;
      if ( not cut) { 
	//cout << "failed #" << i << Names_endcaps_[i] << " val=" 
	//     << fabs(values[ BranchToVariable_[i+n_barrel_] ])
	//     << " current cut =" << CurrentCutsEndcaps_[i]
	//    << endl;
	break;
      }
    }
  }
  return res;

}

//
// try to find the string name in a vector of strings vec
// return the index of vec where you have found it, or -1
// if it is not there
int SobIter::FindName(TString name, vector<TString> vec)
{
  int res = -1;
  for (int i=0; i< (int) vec.size(); ++i)
    if (name == vec[i]) { res = i; break;}

  return res;
}


vector<double>  SobIter::ReturnSeffForSoB
(double target, TH1F* h_signal, TH1F* h_bkg,  
 double rest_signal, double rest_bkg, double SIGNAL_INIT, double tolerance)
{
   int points = h_signal->GetNbinsX();
  //cout << "DEBUG: S/B =" 
  //     << h_signal->Integral(0,nbins)/h_bkg->Integral(0,nbins) << endl;
  //
  double *sob = new double[points];
  double *seff= new double[points];
  double *cut = new double[points];
  double max_sob=0; double max_cut=0; double min_sob = 10000;
  int RealPoints = 0;
  for (int i=1; i < points+1; ++i) {
    double signal = h_signal->Integral(0,i) + rest_signal;
    double bkg = h_bkg->Integral(0,i)+rest_bkg;
    //cout << "signal " << signal << " signal_init " << signal_init << endl;
    seff[i-1] = signal/SIGNAL_INIT;
    if (bkg == 0) {
      cout << "ERROR!!!! Bkg is zero!!!" << endl;
      continue;
    }
    RealPoints++;
    sob[RealPoints-1] = signal/bkg;  //cout << signal/bkg << endl;
    cut[RealPoints-1] = h_signal->GetBinLowEdge(i) + h_signal->GetBinWidth(i);
    if (sob[RealPoints-1] > max_sob && sob[RealPoints-1] < 1000. ) { 
      // for more stable results
      max_sob = sob[RealPoints-1]; max_cut = cut[RealPoints-1];}
    if (sob[RealPoints-1] < min_sob) { min_sob = sob[RealPoints-1];}
  }
  TGraph g_sob(RealPoints, cut, sob); // this is to store the sob
  TGraph g_seff(RealPoints, cut, seff); // this is to store the sob

  //
  // bisection now to find the needed cut value
  int error_flag = 0;
  if (max_sob > 1000){
    cout << "ERROR: Your lower limit in this variable"
			  << " is too low: max sob > 1000!!!!!" << endl;
    error_flag = 1;
  }
  //
  if (sob[0]<sob[RealPoints-1]) {
    min_sob = sob[RealPoints-1];
  }
  //
  //
  // target too small
  if (min_sob > target) {
    cout << "ERROR:  min S/B is "<< min_sob <<" > target "<< target  << endl;
    cout << "Last Bin: Cut["<< RealPoints-1 <<"] = " << cut[RealPoints-1] 
	 << " Seff["<< RealPoints-1 <<"] = " << seff[RealPoints-1]
	 << " S/B["<< RealPoints-1  <<"] = " << sob[RealPoints-1]  << endl;
    error_flag = 2;
  }
  //
  //
  //
  // what happens if the target is too big
  if (max_sob < target) {
    if ( sob[0] < max_sob ) {
      cout << "NOTICE: S/B maximum "<< max_sob <<" is located at cut " 
	   << max_cut
	   << " and it is smaller than the target " << target << endl;
    }
    else {
      cout << "ERROR: S/B max is "<< max_sob  <<" < target " << target << endl;
      cout << "First Bin: Cut[0] = " << cut[0] << " Seff[0] = " 
	   << seff[0] << " S/B[0] = "  << sob[0]  << endl;
    }
    error_flag = 3;
  }

  //  double tolerance = target / 1000.; ==> inserted as an argument
  double working_point = 
    (error_flag == 2)?cut[RealPoints-1]:Bisection(&g_sob, max_cut, 
						  cut[RealPoints-1], 
						  target, tolerance);
  double working_seff = g_seff.Eval(working_point,0,"S");
  if (working_seff>1. && error_flag == 0) {
    cout << "SERIOUS ERROR: Fitting/Bisection gives efficiency " 
	 << working_seff	 << " manual replacement to 1" << endl;
    working_seff = 1.;
  }
  double working_sob = g_sob.Eval(working_point,0,"S");
  //
  // release memory:
  delete [] sob;   delete [] seff;   delete [] cut; 
  vector<double> res;
  res.push_back(working_point); // first entry
  if (error_flag == 2) {
    cout  << "Serious Error detected: program will exit..." << endl;
    res.push_back(-1.);
  }
  else if (error_flag == 3) {
    res.push_back(-3.);
    cout << "No active cut here " << endl;
  }
  else if (working_point > 999.) {
    cout << "Your actual values is seff = " << working_seff 
	 << " but this working point will be excluded because" 
	 <<  " the output is 1000 - Bisection Failure!!!!"
	 << endl;
    res.push_back(-1.);
  }
  else {
    res.push_back(working_seff); // sendond entry
    cout << "Working cut "<< working_point
	 << " Working seff is " << working_seff 
	 << " Working Sob is " << working_sob 
	 << endl;
  }
  res.push_back(working_sob); // 3rd entry
  if (error_flag == 2) {
    res.push_back(1);  // 4th entry
  } 
  else
    res.push_back(0);
  return res;


}
//
// this function is Bisection_Method2 from SelectionOptimizer class
// written by NR - December 08
//
// the meaning of tolerance  is wrt the target value
// i.e. if |f(c) - Target| < tolerance => stop
double SobIter::Bisection
(TGraph *g, double a_, double b_, double target, double tolerance)
{
  double a = a_;
  double b = b_;
  double c = (a+b)/2.;
  double fa = g->Eval(a,0,"S") - target;
  double fb = g->Eval(b,0,"S") - target;
  if (fa*fb < 0) {
    for (int i =0; i < 100000; ++i) {
      c = (a + b)/2.;
      double fc = g->Eval(c,0,"S") - target;
      if (fabs(fc) < tolerance) break;
      if (fc*fa < 0) {	fb = fc; b =c;      }
      else if (fc*fb < 0) {fa = fc; a=c;}
      else break;
    }

  }
  else if (fa*fb == 0) return (fa == 0.0 ? a:b);
  else return 10000;

  return c;
}

//
//
//
bool SobIter::PassFilter
(Filter myfilter, double et, double met, short int t, short int exclusion)
{
  // check the ET cut
  if (myfilter.ET>=0) {
    if  (myfilter.ET_LESS) { // passes if et < ET
      if ((et < myfilter.ET) == false) return false;
    }
    else {  // passes if  et > ET
      if ((et > myfilter.ET) == false) return false;
    }
  }
  // check the MET cut
  if (myfilter.MET>=0) {
    if  (myfilter.MET_LESS) { // passes if met < MET
      if ((met < myfilter.MET) == false) return false;
    }
    else {  // passes if  met > MET
      if ((met > myfilter.MET) == false) return false;
    }
  }
  // check the type cut
  vector<Short_t> mytypes = myfilter.TYPE;
  int ntypes = mytypes.size();
  bool skip = false;
  if (ntypes == 1) { if (mytypes[0] <0) skip = true; }
  if (ntypes>=0 && (not skip)) {
    if (myfilter.TYPE_EQUAL) { // passes if t == TYPE
      //      for at least one of the specified types
      bool res = false;
      // use logical OR: if one equality is true, then
      // the variable res will become true
      for (int i=0; i<ntypes; ++i) {
	res = res || (t == mytypes[i]);
      }
      if (res == false) return false;
    }
    else {  // passes if t != TYPE
      // for each of the specified types
      for (int i=0; i<ntypes; ++i) {
	if (t == mytypes[i]) return false;
      }
    }
  }
  if (myfilter.EXCLUDE) { // fails if exclusion == 1
    if (exclusion == 1) return false;
  }
  return true;
}


void SobIter::PrintFilter(Filter filter)
{
  cout << "Filter Summary of filter: " << filter.NAME << endl;
  if (filter.ET >=0) {
    if (filter.ET_LESS)
      cout << "Cuts: et<" << filter.ET << endl;
    else
      cout << "Cuts: et>" << filter.ET << endl;
  }
  if (filter.MET >=0) {
    if (filter.MET_LESS)
      cout << "Cuts: met<" << filter.MET << endl;
    else
      cout << "Cuts: met>" << filter.MET << endl;
  }
  if (filter.TYPE_EQUAL) 
    cout << "Event type =";// << filter.TYPE << endl;
  else
    cout << "Event type!=";// << filter.TYPE << endl;
  vector<Short_t> mytypes= filter.TYPE;
  for (int i=0; i< (int) mytypes.size(); ++i)
    {
      cout << mytypes[i] << "  ";
    }
  cout << endl;
  if (filter.EXCLUDE)
    cout << "Events with a second electron are excluded" << endl;
  else
    cout << "All electrons included "  << endl;
  cout << "End Summary for filter " << filter.NAME << endl;

}



vector<double> SobIter::makeRoundCuts( vector<double> vec, 
				       vector<TString> names) {
  vector<double> res;
  vector<TString>::iterator n= names.begin();
  for (vector<double>::iterator it = vec.begin(); it != vec.end(); ++it) {
    //cout << "Variable " <<  *n << endl;
    double theCut = *it;
    double theRoundedCut = theCut;
    // check whether the name is dist or dcot
    TString theName = *n;
    if (theName.Contains("Tdist") || theName.Contains("Tdcot")) {
      //cout << "found name: " << theName << endl;
      //////theRoundedCut = theCut;
      // then transform the cut, then round it
      theCut = (1-theCut) / theCut;
      theRoundedCut = RoundNumber(theCut);
      // and finally return the transform of the rounded:
      theRoundedCut = 1/(1.+ theRoundedCut);
    }
    else {
      theRoundedCut = RoundNumber(theCut);
    }
    //
    res.push_back(theRoundedCut);
    ++n;
  }
  return res;

}

//
// Chris' algorithm
// 
// round to 3 significant digits
//     put it in the from  0.XXX x 10^exp
//     hence the smaller is 0.100 and the biggest 0.999
// choose the closest value of the following:
// 0.10 0.12 0.15 0.20 0.25 0.30 0.40 0.50 0.60  0.70  0.80  0.90  1.0
// if there are 2 that are close together, choose the highest value
// i.e. for 0.55 select 0.60 instead of 0.50
double SobIter::RoundNumber(double n) {
  // bring number in 0.XXX format
  //cout << " -----> entering RoundNumber with: "<< n  << endl;
  if (n == 0) return n;
  double factor = 1;
  if (n>1.)
    factor = pow(10, 1+(unsigned) log10(n));
  else {
    factor = 1.;
    double nbk = factor*n;
    while ( 1>0 ) {
      if (nbk < 1 && 10*nbk >1) break;
      factor = factor*10;
      nbk = factor*n;
    }
    factor = 1./factor;
  }
  double num = n/factor;
  unsigned p = 3-(unsigned)log10(num);
  double theNum = 0;
  if ((num * pow(10.,p) - int(num * pow(10.,p))) >= 0.5) 
    theNum = int(num * pow(10.,p)) +1;
  else
    theNum = int(num * pow(10.,p));
  double format3 = theNum / pow(10.,p);
  //cout << "\n >>>>> format3: " << format3 << " -- " << factor  << endl;
  // find now which is closer
  vector<double> points;
  points.push_back(0.10);  points.push_back(0.12);  points.push_back(0.15);
  points.push_back(0.20);  points.push_back(0.25);  
  points.push_back(0.30);  points.push_back(0.40);
  points.push_back(0.50);  points.push_back(0.60);  points.push_back(0.70);
  points.push_back(0.80);  points.push_back(0.90);  points.push_back(1.00);
  //
  int nPoints = (int) points.size();
  //
  for (int i=0; i<nPoints-1; ++i) {
    //cout << " >>>>> point #" << i << endl;
    double diff0 = format3 - points[i];
    double diff1 = points[i+1] -format3;
    //cout << " >>>>>  diif01 " << diff0 << " - " << diff1 << endl;
    if (diff0 == 0) return points[i]*factor;
    if (diff1 == 0) return points[i+1]*factor;
    if (diff0 > 0 && diff1 >0) {
      if (diff0 > diff1) return points[i+1]*factor;
      else  if (diff0 < diff1) return points[i]*factor;
      else  return points[i+1]*factor;
    }
  }
  // if you are here there is a mistake in the algo
  cout << "ERROR!!!!! in RoundNumber" << endl;
  return -999999.;

}