//
//
/*
  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 = "event_MET", int preSelType = 0);
  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);
};


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)
{
  isAlive_ =true;
  // 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]);


  //
  //
  // 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;
  //return;
  // / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 
  /////////////////////////////////////////////////////////////////
  //                                         //////////////////////
  //   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.;
    //
    // LOOP OVER THE TREE ENTRIES ........................................
    //
    // Loop over the signal 
    for (int iT = 0; iT < N_signal; ++iT) { // loop over the signal
      tree.GetEntry(iT);
      bool pass = CheckCuts(vars, probe_sc_eta);
      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)) {
	/////////////////////////////////////////////////////////////
	//cout << "met=" << event_MET << ", et=" << probe_sc_et 
	if (Iter == 0) { 
	  if (in_barrel || in_endcaps) {
	    signal_init2 += (w_[type]*w_[type]);
	    if (in_barrel)
	      signal_init_barrel += w_[type];
	    if (in_endcaps)
	      signal_init_endcaps += w_[type];
	  }
	}
	//
	if (pass && in_barrel) {
	  current_signal_barrel += w_[type];
	  dS2 +=  w_[type]*w_[type];
	  if (type == 0) current_signal_in_sigSample_eb += w_[type];
	  for (int i=0; i<n_barrel_; ++i) {
	    h_barrel[i].Fill(vars[i], w_[type]);
	  }
	}
	else if (pass && in_endcaps) {
	  current_signal_endcaps += w_[type];
	  dS2 +=  w_[type]*w_[type];
	  if (type == 0) current_signal_in_sigSample_ee += w_[type];
	  for (int i=0; i<n_endcaps_; ++i)
	    h_endcaps[i].Fill(vars[ BranchToVariable_[i+n_barrel_] ],w_[type]);
	}
	/////////////////////////////////////////////////////////////////////
      }
      //// bkg related ......................................................
      if (PassFilter(BkgFilter_, probe_sc_et, event_MET, type, exclusion)) {
	if (Iter == 0 && somewhere) {
	  bkg_init += w_[type]; bkg_init2 += (w_[type]*w_[type]);
	}
	//
	if (pass && in_barrel) {
	  current_bkg_barrel += w_[type];
	  dB2 +=  w_[type]*w_[type];
	  if (type == 0) current_signal_in_bkgSample_eb += w_[type];
	  for (int i=0; i<n_barrel_; ++i)
	    b_barrel[i].Fill(vars[i], w_[type]);
	}
	else if (pass && in_endcaps) {
	  current_bkg_endcaps += w_[type];
	  dB2 +=  w_[type]*w_[type];
	  if (type == 0) current_signal_in_bkgSample_ee += w_[type];
	  for (int i=0; i<n_endcaps_; ++i)
	    b_endcaps[i].Fill(vars[ BranchToVariable_[i+n_barrel_] ],w_[type]);
	}
      }
      // MC values for signal and bkg -- to be changed -- this is for wenu
      // only
      bool Ex =  true;
      if (SignalFilter_.EXCLUDE) Ex = exclusion != 1;
      if (probe_sc_et> SignalFilter_.ET && Ex && somewhere) {
	if (type == 0) {
	  if (Iter==0) {
	    signal_init_mc += w_[type]; signal_init_mc2+=(w_[type]*w_[type]);
	  }
	  if (pass) {
	    current_actual_signal += w_[type]; dS2_mc += w_[type]*w_[type];
	  }
	}
	//	if (type >=1  && type <=13) {
	else { // all the rest are bkg
	  if (Iter==0) {
	    bkg_init_mc += w_[type];  bkg_init_mc2 += (w_[type]*w_[type]);
	  }
	  if (pass) {
	    current_actual_bkg += w_[type]; dB2_mc += w_[type]*w_[type];
	  }
	}

      }
      //////////////////////////////////////////////////////////////////////
    } // end loop over the samples
    // ...................................................................
    // -------------------------------------------------------------------
    // List of the counters that you have available now ------------------
    // ...................................................................
    // signal_init_barrel           signal sample contents at Iter == 0 
    // signal_init_endcaps, signal_init          
    // ...................................................................
    // signal_init_mc             number of events with type == 0 that
    //                            pass the kinematic restrictions of the
    //                            signal sample (ET, fiducial, exclusion)
    //                            at Iter == 0
    // bkg_init_mc                number of events with type != 0 + 
    //                            signal sample kinematic restrictions
    //                         >> may need to be modified for Zee-driven
    // ...................................................................
    // current_signal_barrel      number of events in signal sample that
    // current_signal_endcaps     pass kinematics+selection
    // current_bkg_barrel         number of events in bkg sample that pass
    // current_bkg_endcaps        kinematics+selection
    // ...................................................................
    // current_actual_signal      number of events that pass the kinematics
    // current_actual_bkg         + selection of the signal sample
    // ...................................................................
    // current_signal_in_sigSample  number of events with type == 0 in the
    // current_signal_in_bkgSample  signal and bkg sample
    // ...................................................................
    // to be done - slight redefinitions needed for the Zee-driven tuning
    // -------------------------------------------------------------------
    if (Iter==0) {
      signal_init = signal_init_barrel + signal_init_endcaps;
      sobInit = signal_init/bkg_init;
      sobInitError = sobInit * sqrt(signal_init2/(signal_init*signal_init)+
				    bkg_init2/(bkg_init*bkg_init));
      sobInitMC = signal_init_mc/bkg_init_mc;
      sobInitMCError = sobInitMC * 
	sqrt(signal_init_mc2/(signal_init_mc*signal_init_mc)+
	     bkg_init_mc2/(bkg_init_mc*bkg_init_mc));
    }
    current_signal = current_signal_barrel + current_signal_endcaps;
    current_bkg = current_bkg_barrel + current_bkg_endcaps;
    if (current_bkg == 0 || current_signal == 0) break;
    double sob = current_signal / current_bkg;
    double sob_error = sob * sqrt(dS2/(current_signal*current_signal)+
				  dB2/(current_bkg*current_bkg));
    double sob_barrel = 
      (current_bkg_barrel>0)?(current_signal_barrel/current_bkg_barrel):0;
    double sob_endcaps = 
      (current_bkg_endcaps>0)?(current_signal_endcaps/current_bkg_endcaps):0;
    if (Iter==0) {
      // this is in order to be able to have available the S/B after presel
      cout << "initial S/B= " << sobInit << " +- " << sobInitError
	   << " fromMC: " << sobInitMC << " +- " << sobInitMCError
	   << endl;
      cout << "Current signal = " << current_signal
	   << " Current bkg = " << current_bkg
	   << ", init=" << signal_init << " sigMC=" << current_actual_signal
	   << ", bkgMC=" << bkg_init_mc << ", sigInitMC="<<signal_init_mc
	   << endl;
    }
     // debugging...
    double current_signal_in_sigSample =current_signal_in_sigSample_eb +
      current_signal_in_sigSample_ee;
    double current_signal_in_bkgSample =current_signal_in_bkgSample_eb +
      current_signal_in_bkgSample_ee;
    cout << "Current signal = " << current_signal 
	 << " Current bkg = " << current_bkg
	 << ", init=" << signal_init << " sigMC=" << current_actual_signal
         << ", bkgMC=" << current_actual_bkg << ", sigInitMC="<<signal_init_mc
	 << ", CurrentSigInSig=" << current_signal_in_sigSample
	 << ", CurrentSigInBkg=" << current_signal_in_bkgSample
	 << endl;
    
    double seff = current_signal/signal_init; 
    double seff_barrel = 
      (signal_init_barrel>0)?(current_signal_barrel/signal_init_barrel):0;
    double seff_endcaps = 
      (signal_init_endcaps>0)?(current_signal_endcaps/signal_init_endcaps):0;
    double sob_mc = 
      (current_actual_bkg>0)?(current_actual_signal/current_actual_bkg):0;
    double sob_mc_error = 
      sob_mc * sqrt(dS2_mc/(current_actual_signal*current_actual_signal)+
		    dB2_mc/(current_actual_bkg*current_actual_bkg));

    double seff_mc = 
      (signal_init_mc>0)?(current_actual_signal/signal_init_mc):0;
    //
    TString selection = "(";
    for (int i=0; i<n_barrel_-1; ++i) {
      selection +=  CurrentCutsBarrel_[i];
      selection += ", ";
    }
    if (n_barrel_>0)
      selection +=  CurrentCutsBarrel_[n_barrel_-1];
    else
      selection +=  "(no barrel variables)";
    selection += " // ";
    for (int i=0; i<n_endcaps_-1; ++i) {
      selection += CurrentCutsEndcaps_[i];
      selection += ", ";
    }
    if (n_endcaps_>0)
      selection += CurrentCutsEndcaps_[n_endcaps_-1];
    else
      selection += "(no endcap variables)";
    selection += ")";
    //
    if (Iter == 0) {
      cout << "Initial Iteration: S/B=" << sob << "  Initial Signal="
	   << signal_init << " in EB: " << signal_init_barrel << " in EE: "
	   << signal_init_endcaps << ", Initial Conditions: "
	   << selection << endl;
    }
    else {
      cout << "Iteration #" << Iter << " S/B=" << sob << "+/-" << sob_error
	   << "; seff=" << seff  << "; S="<< current_signal << "; B="
	   << current_bkg << ";  BARREL: S/B=" << sob_barrel << "; seff="
	   << seff_barrel << "; ENDCAPS: S/B=" << sob_endcaps << "; seff="
	   << seff_endcaps << ";      Cuts: " << selection << endl;
    }
    SOB.push_back(sob);    SEFF.push_back(seff);
    SOB_ERROR.push_back(sob_error);
    SELECTION.push_back(selection);
    SOB_MC.push_back(sob_mc);    SEFF_MC.push_back(seff_mc);
    SOB_MC_ERROR.push_back(sob_mc_error);
    // -----------------------------------------------------------------------
    // variable  tuning  to be performed here --------------------------------
    //
    // set the S/B target ....................................................
    double abs_step;
    if (RELATIVE_STEP_) {
      abs_step = Step_ * sob;
      ///if (Iter == 0) abs_step = 0.01;
      if (abs_step > minStep_) abs_step = minStep_;
    }
    else abs_step = Step_;
    double target = sob + abs_step;
    double tolerance = abs_step/100.;
    cout << "Targeting S/B " << target << ", S/B step is " << abs_step << endl;
    //
    // for each of the variables calculate the seff for the target S/B
    double *tmp_cuts  = new double[n_barrel_+n_endcaps_];
    double *tmp_seffs = new double[n_barrel_+n_endcaps_];
    double *tmp_sobs  = new double[n_barrel_+n_endcaps_];
    double error_flag = 0;
    // run over the barrel variables
    for (int i=0; i< n_barrel_; ++i) {
      cout << "Working on Barrel " << Names_barrel_[i] << endl;
      vector<double> a = ReturnSeffForSoB
	(target, &h_barrel[i], &b_barrel[i], 
	 current_signal_endcaps, current_bkg_endcaps, signal_init, tolerance);
      tmp_cuts[i] = a[0];    tmp_seffs[i] = a[1];    tmp_sobs[i] = a[2];
      error_flag += a[3];
      if (a[0]< LowerAllowedBarrel_[i]) {
	cout << "Cut Value Lower than Limit: cut will not be moved..."<< endl;
	tmp_seffs[i] = 0.;
      }
    }
    // run over the endcap variables
    for (int i=0; i< n_endcaps_; ++i) {
      cout << "Working on Endcaps " << Names_endcaps_[i] << endl;
      vector<double> a = ReturnSeffForSoB
	(target, &h_endcaps[i], &b_endcaps[i], 
	 current_signal_barrel, current_bkg_barrel, signal_init, tolerance);
      tmp_cuts[i+n_barrel_] = a[0];    
      tmp_seffs[i+n_barrel_] = a[1];    tmp_sobs[i+n_barrel_] = a[2];
      error_flag += a[3];     
      // check whether you are below a certain limit
      if (a[0]< LowerAllowedEndcaps_[i]) {
	cout << "Cut Value Lower than Limit: cut will not be moved..."<< endl;
	tmp_seffs[i+n_barrel_] = 0.;
      }
    }
    // find the maximim efficiency
    int imax = TMath::LocMax(n_barrel_ + n_endcaps_, tmp_seffs);

    double movement;
    if (imax<n_barrel_) 
      movement = CurrentCutsBarrel_[imax] - tmp_cuts[imax];      
    else 
      movement = CurrentCutsEndcaps_[imax-n_barrel_] - tmp_cuts[imax];      
    if (tmp_seffs[imax] < -1) {
      cout << "S/B cannot increased more ..." << endl;
      break;
    }
    // update the variable cuts now
    if (imax < n_barrel_) {
      vector <double> tmp;
      tmp = CurrentCutsBarrel_;
      CurrentCutsBarrel_.clear();
      for (int i=0; i< n_barrel_; ++i) {
	double cut =  tmp[i];
	if (i == imax) cut = tmp_cuts[imax];
	CurrentCutsBarrel_.push_back(cut);
      }
      cout << "Iter #" << Iter << ": Variable Chosen: EB " 
	   << Names_barrel_[imax] << " to " << tmp_cuts[imax] 
	   << " moved by " << movement << endl;
    }
    else {
      vector<double> tmp;
      tmp = CurrentCutsEndcaps_;
      CurrentCutsEndcaps_.clear();
      for (int i=0; i< n_endcaps_; ++i) {
	double cut =  tmp[i];
	if (i == (imax-n_barrel_)) cut = tmp_cuts[imax];
	CurrentCutsEndcaps_.push_back(cut);
      }
      cout << "Iter #" << Iter << ": Variable Chosen: EE " 
	   << Names_endcaps_[imax-n_barrel_] << " to " << tmp_cuts[imax] 
	   << " moved by " << movement << endl;
    }
    //for (int i=0; i< n_endcaps_; ++i)
    //  cout << CurrentCutsEndcaps_[i] << endl;
    // release memory:
    delete [] tmp_cuts;    delete [] tmp_seffs;    delete [] tmp_sobs;
    // -----------------------------------------------------------------------
    //cout << "Current Signal is: " << current_signal << endl
    //	 << "Current Bkg is: " << current_bkg << endl;
    if (seff < LOWEST_EFF_)  break;
    if (movement == 0) break;
    if (sob_previous == sob || sob_previous > sob || error_flag != 0) break;
    sob_previous = sob;
  } // end S/B Iteration loop
  //
  // print out the summary of the method:
  const int Niter = (int) SOB.size();
  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)
    cout << "sob[" << i << "]=" << SOB[i] << "; seff[" << i << "]=" << SEFF[i]
	 << "; sob_mc[" << i << "]=" << SOB_MC[i] 
	 << "; seff_mc[" << i << "]=" << SEFF_MC[i] 
	 << "; sob_error[" << i << "]=" << SOB_ERROR[i] 
	 << ";  sob_mc_error[" << i << "]=" << SOB_MC_ERROR[i] 
	 << ";     //  "  << SELECTION[i] << endl;


  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;

}

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

}