// ----------------------------------------------------------------------------------------------------- //
//                                                                                                       //
// Compute expected numbers of WW->dilepton events in different measurement channels, based on           //
// Monte Carlo numbers, cross sections, and luminosities.                                                //
//                                                                                                       //
// Instructions :                                                                                        //
// ==============                                                                                        //
// (1) Set the flag "apply_syst", depending on whether systematic error estimates are to be included     //
//     in the grand summary tables.                                                                      //
//     WARNING : in general, systematic errors are assumed uncorrelated between different dilepton       //
//     channels and background sources (the WW signal systematic is treated differently, because we      //
//     know that the dominant systematic error sources (e.g. pT(WW) distribution) are common between     //
//     e-e, e-mu and mu-mu. Correlations should be studied in detail and the appropriate modifications   //
//     made to the error combination code.                                                               //
//                                                                                                       //
// (2) For each background type, enter :                                                                 //
//     (i)   the number of events analysed                                                               //
//     (ii)  the cross section in pb                                                                     //
//     (iii) the filter efficiency (e.g. reported by HEPG filter in simulation log file).                //
//           Set to 1.0 if no filter was applied. For the WW signal, the branching ratio is applied      //
//           here.                                                                                       //
//     (iv)  any k-factor that should be applied to the cross section in (ii). Note that if (ii) is      //
//           already the best estimate (NLO or NNLO) of the cross section, no k-factor should be         //
//           specified here.                                                                             //
//     (v)   systematic error estimate.                                                                  //
//     (vi)  the input ROOT file containing the cut counter 2-d histogram.                               //
//                                                                                                       //
// (3) Specificy the variable "luminosity_systematic", currently 6%.                                     //
//                                                                                                       //
// (4) Specify the luminosities (in pb-1) for the different run ranges used in the analysis :            //
//     (i)   lum_nosi : good run list without silicon requirement (usually the largest number).          //
//     (ii)  lum_nosicmx : CMX good run list without silicon requirement.                                //
//     (iii) lum_si  : good run list with silicon requirement.                                           //
//     (iv)  lum_sicmx : CMX good run list with silicon requirement (usually the smallest number).       //
//                                                                                                       //
// (5) Specify the trigger scale factors. Relevant CDF notes listed below.                               //
//                                                                                                       //
// (6) Specify the (data/MC) efficiency scale factors. Relevant CDF notes listed below.                  //
//                                                                                                       //
// (7) Run :  setup CDF software, then > gmake clean; gmake; ./ww_dilepton.exe                           //
//                                                                                                       //
//                                                                                                       //
//                                                                                                       //
// Date        Person               Notes                                                                //
// ============================================================================================          //
// 29.07.03    Dave Waters          First version.                                                       //
//                                                                                                       //
// ----------------------------------------------------------------------------------------------------- //

#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <string>
#include <math.h>
#include <vector>
#include <algorithm>

#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom3.h"
#include "TVector3.h"
#include "TMath.h"

using std::cout;
using std::hex;
using std::dec;
using std::endl;

typedef std::vector<double> rowd;
typedef std::vector<rowd> tabled;

typedef std::vector<int> rowi;
typedef std::vector<rowi> tablei;

typedef std::vector<std::string> label;
typedef std::vector<label> table_label;

enum {
  null = -1,
  ee   = 0,
  emu  = 1,
  mumu = 2
};


// Functions and classes required for Feldman & Cousins interval calculation :
// ---------------------------------------------------------------------------
static void feldman_cousins(int n0, double bg, double cl);
static void feldman_cousins(int n0, double bg, double cl, double& signal, double& up_error, double& down_error);
static double poisson(double mean, double bg, int n);
void cbelt(double mu, double bg, double cl, double& low, double& high);
class ratio{
public:
  int n;
  double v;
};
static bool largerRatio(ratio e1, ratio e2) {
  return (e1.v>e2.v);
}
double rangauss(double mean, double sigma, TRandom3* ranGen);

int main()
{
  
  // counter[i][j] : 0<=i<=7  : cut number
  //                 0<=j<=26 : dilepton category 
  //                        0 : TCE-TCE
  //                        1 : TCE-PHX
  //                        2 : TCE-PEM 
  //                        3 : PHX-PHX
  //                        4 : PHX-PEM
  //                        5 : TCE-CMUP
  //                        6 : TCE-CMU
  //                        7 : TCE-CMP
  //                        8 : TCE-CMX
  //                        9 : TCE-CMIO
  //                       10 : PHX-CMUP
  //                       11 : PHX-CMU
  //                       12 : PHX-CMP
  //                       13 : PHX-CMX
  //                       14 : PHX-CMIO
  //                       15 : PEM-CMUP
  //                       16 : PEM-CMX
  //                       17 : CMUP-CMUP
  //                       18 : CMUP-CMU
  //                       19 : CMUP-CMP
  //                       20 : CMUP-CMX
  //                       21 : CMUP-CMIO
  //                       22 : CMX-CMX
  //                       23 : CMX-CMU
  //                       24 : CMX-CMP
  //                       25 : CMX-CMIO
  //                       26 : TRI-LEPTONS

  // Systematics on or off ?
  // -----------------------
  bool apply_syst = true;

  // Drell-Yan e+e-
  // ==============
  double dy_ee_nevt_total  = 456986;     // Total number of MC events processed at ntuple level.
  double dy_ee_cs          = 236.0;      // Pythia cross section, pb.
  double dy_ee_filter_eff  = 0.5991667;  // HEPG filter efficiency.
  double dy_ee_k_factor    = 1.40;       // NNLO/LO from CDF-6202
  double dy_ee_systematic     = 0.00;
  if (apply_syst) {
    dy_ee_systematic     = 0.30; // From DY e+e- studies
  }
  tabled dy_ee_results(8,rowd(27));
  // TFile dy_ee("root_files_26.07.2003/wwana_4111_072503_ZEE.root"); // pre-blessing
  // TFile dy_ee("root_files_26.07.2003/wwana_4111_010803_ZEE_test.root");
  TFile dy_ee("root_files_26.07.2003/wwana_4111_040803_ZEE.root");
  TH2F* dy_ee_counter = (TH2F*)(dy_ee.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      dy_ee_results[i][j] = dy_ee_counter->GetBinContent(i,j);
    }
  }

  // Drell-Yan mu+mu-
  // ================
  double dy_mumu_nevt_total  = 539443;     // Total number of MC events processed at ntuple level.
  double dy_mumu_cs          = 236.0;      // Pythia cross section, pb.
  double dy_mumu_filter_eff  = 0.6109434;  // HEPG filter efficiency.
  double dy_mumu_k_factor    = 1.40;       // NNLO/LO from CDF-6202
  double dy_mumu_systematic   = 0.00;
  if (apply_syst) {
    dy_mumu_systematic   = 0.35; // From DY mu+mu- studies
  }
  tabled dy_mumu_results(8,rowd(27)); // No TCE-TCE
  // TFile dy_mumu("root_files_26.07.2003/wwana_4111_072503_ZMUMU.root"); // pre-blessing
  TFile dy_mumu("root_files_26.07.2003/wwana_4111_040803_ZMUMU.root");
  TH2F* dy_mumu_counter = (TH2F*)(dy_mumu.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      dy_mumu_results[i][j] = dy_mumu_counter->GetBinContent(i,j);
    }
  }

  // Drell-Yan tau+tau-
  // ==================
  double dy_tautau_nevt_total  = 418269;     // Total number of MC events processed at ntuple level.
  double dy_tautau_cs          = 236.0;      // Pythia cross section, pb.
  double dy_tautau_filter_eff  = 0.0088750;  // HEPG filter efficiency.
  double dy_tautau_k_factor    = 1.40;       // NNLO/LO from CDF-6202
  double dy_tautau_systematic = 0.00;
  if (apply_syst) {
    dy_tautau_systematic = 0.35; // Maximum of the e+e- and mu+mu- errors.
  }
  tabled dy_tautau_results(8,rowd(27)); 
  // TFile dy_tautau("root_files_26.07.2003/wwana_4111_072503_ZTAUTAU.root"); //pre-blessing
  TFile dy_tautau("root_files_26.07.2003/wwana_4111_040803_ZTAUTAU.root");
  TH2F* dy_tautau_counter = (TH2F*)(dy_tautau.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      dy_tautau_results[i][j] = dy_tautau_counter->GetBinContent(i,j);
    }
  }
  
  // ttbar :
  // =======
  double ttbar_nevt_total  = 380653;     // Total number of MC events processed at ntuple level.
  double ttbar_cs          = 7.0;        // Besth theoretical cross section value, pb.  // CHANGE TO 6.7 ???
  double ttbar_filter_eff  = 1.00;       // HEPG filter efficiency.
  double ttbar_k_factor    = 1.00;       
  double ttbar_systematic  = 0.0;
  if (apply_syst) ttbar_systematic  = 0.30;       // Systematic uncertainty on ttbar Monte Carlo expectation.
  tabled ttbar_results(8,rowd(27)); 
  // TFile ttbar("root_files_26.07.2003/wwana_4111_072503_TTBAR.root"); //pre-blessing
  TFile ttbar("root_files_26.07.2003/wwana_4111_040803_TTBAR.root");
  TH2F* ttbarcounter = (TH2F*)(ttbar.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      ttbar_results[i][j] = ttbarcounter->GetBinContent(i,j);
    }
  }

  // WZ :
  // ====
  double wz_nevt_total  = 80744;     // Total number of MC events processed at ntuple level.
  double wz_cs          = 1.98*2.;      // Best theoretical cross section value, pb.
  double wz_filter_eff  = 0.0580;    // HEPG filter efficiency.
  double wz_k_factor    = 1.00;       
  double wz_systematic  = 0.0;
  if (apply_syst) wz_systematic  = 0.07;      // Systematic uncertainty on wz Monte Carlo expectation.
  tabled wz_results(8,rowd(27)); 
  // TFile wz("root_files_26.07.2003/wwana_4111_072503_WZ.root"); // pre-blessing
  TFile wz("root_files_26.07.2003/wwana_4111_040803_WZ.root"); // blessing
  TH2F* wzcounter = (TH2F*)(wz.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      wz_results[i][j] = wzcounter->GetBinContent(i,j);
    }
  }

  // WW :
  // ====
  double ww_nevt_total  = 825081;     // Total number of MC events processed at ntuple level.
  double ww_cs          = 13.25;      // Best theoretical cross section value, pb.
  // Note about W->lv branching ratios :
  // -----------------------------------
  // 0.1068  = PDG measured value. We used this for Winter 2003 analysis and Summer 2003 pre-blessing.
  // 0.11043 = value used by Campbell & Ellis in their calculations (hep-ph/9905386). We should use the same value.
  // double ww_filter_eff  = pow(0.1068*3,2);
  double ww_filter_eff  = pow(0.11043*3,2);

  double ww_k_factor    = 1.00;       
  double ww_systematic  = 0.0;
  double ww_correlated  = 0.0;
  if (apply_syst) {
    ww_systematic  = 0.22;       // Systematic uncertainty on ww Monte Carlo expectation.
    ww_correlated  = 0.21;       // The fraction of the systematic error correlated amongst different channels.
  }
  tabled ww_results(8,rowd(27));   
  // TFile ww("root_files_26.07.2003/wwana_4111_072503_WW.root"); // pre-blessing
  TFile ww("root_files_26.07.2003/wwana_4111_040803_WW.root");
  TH2F* wwcounter = (TH2F*)(ww.Get("global_counter"));
  for (int i=0; i<8; ++i) {
    for (int j=0; j<27; ++j) {
      ww_results[i][j] = wwcounter->GetBinContent(i,j);
    }
  }

  // *****************************
  // *** FINAL RESULTS SECTION ***
  // *****************************
  cout << "\n" << " WW->Dilepton Summary of Results" << endl;
  cout << " ===============================\n" << endl ;

  // Sort out the labels :
  // =====================
  label tags(27);
  label region(27);

  tags[0] ="TCE-TCE   ";  region[0]="CC";
  tags[1] ="TCE-PHX   ";  region[1]="CP";
  tags[2] ="TCE-PEM   ";  region[2]="CP";
  tags[3] ="PHX-PHX   ";  region[3]="PP";
  tags[4] ="PHX-PEM   ";  region[4]="PP";
  tags[5] ="TCE-CMUP  ";  region[5]="CC";
  tags[6] ="TCE-CMU   ";  region[6]="CC";
  tags[7] ="TCE-CMP   ";  region[7]="CC";
  tags[8] ="TCE-CMX   ";  region[8]="CC";
  tags[9] ="TCE-CMIO  ";  region[9]="CC";
  tags[10]="PHX-CMUP  ";  region[10]="CP";
  tags[11]="PHX-CMU   ";  region[11]="CP";
  tags[12]="PHX-CMP   ";  region[12]="CP";
  tags[13]="PHX-CMX   ";  region[13]="CP";
  tags[14]="PHX-CMIO  ";  region[14]="CP";
  tags[15]="PEM-CMUP  ";  region[15]="CP";
  tags[16]="PEM-CMX   ";  region[16]="CP";
  tags[17]="CMUP-CMUP ";  region[17]="CC";
  tags[18]="CMUP-CMU  ";  region[18]="CC";
  tags[19]="CMUP-CMP  ";  region[19]="CC";
  tags[20]="CMUP-CMX  ";  region[20]="CC";
  tags[21]="CMUP-CMIO ";  region[21]="CC";
  tags[22]="CMX-CMX   ";  region[22]="CC";
  tags[23]="CMX-CMU   ";  region[23]="CC";
  tags[24]="CMX-CMP   ";  region[24]="CC";
  tags[25]="CMX-CMIO  ";  region[25]="CC";
  tags[26]="TRILEPTON ";  region[26]="NONE";
  
  label cuts(8);
  cuts[0]=" ID           : ";
  cuts[1]=" ISO          : ";
  cuts[2]=" Conv+Cosm    : ";
  cuts[3]=" Z-veto       : ";
  cuts[4]=" Missing-Et   : ";
  cuts[5]=" Delta-Phi    : ";
  cuts[6]=" 0-jet        : ";
  cuts[7]=" OS           : ";
  
  // Luminosities :
  // ==============
  rowd luminosity(27);
  double luminosity_systematic = 0.00;
  if (apply_syst) {
    luminosity_systematic = 0.06; // 6%
  }

  rowi allowed_combination(27);
  rowd trigger_acceptance(27);
  rowd eff_fac(27);

  // Winter acceptance = Winter categories (no PHX or CMIO)
  // Summer acceptance = Summer categories (PHX and CMIO)
  // Winter data = Winter run list and corresponding luminosities
  // Summer data = Summer run list and corresponding luminosities
  bool winter_acceptance = false;
  bool summer_acceptance = true;
  bool winter_data = false;
  bool summer_data = true;

  double lum_nosi=0;
  double lum_nosicmx=0;
  double lum_si=0;
  double lum_sicmx=0;
  
  int data_ee = 0;
  int data_mm = 0;
  int data_em = 0;

  // SCALE FACTOR : 
  // --------------
  // An overall scale factor applied to each dilepton category equally, via the luminosity.
  double overall_factor = 0.951/0.973; // data/MC z0 efficiency

  if(winter_data){
    lum_nosi=(71995.2/1000.)*overall_factor;
    lum_nosicmx=(55919.5/1000.)*overall_factor;
    lum_si=(43617.4/1000.)*overall_factor;
    lum_sicmx=(33862.3/1000.)*overall_factor;

    data_ee = 0;  // After losing TCE-TCE event (CES Chi2 bug)
    data_mm = 0;
    data_em = 1;
  }
  if(summer_data){
    lum_nosi=125.823*overall_factor;
    lum_nosicmx=109.442*overall_factor;
    lum_si=94.4042*overall_factor;
    // lum_si=125.823;
    lum_sicmx=84.4638*overall_factor;

    data_ee = 2;
    data_mm = 1;
    data_em = 2;
  }

  // TRIGGER SCALE FACTORS :
  // -----------------------
  // Numbers used for Summer 2003 pre-blessing :
  // double trig_tce  = 0.96;   // CDF-6234
  // double trig_cmup = 0.932;  // CDF-6293
  // double trig_cmx  = 0.914;  // CDF-6293
  // double trig_phx = 0.91  // AND NO SCALE FACTOR

  // Numbers used for Summer 2003 blessing :
  double trig_tce  = 0.96;   // CDF-6234
  double trig_cmup = 0.890;  // CDF-6610
  double trig_cmx  = 0.966;  // CDF-6610

  // MUON SCALE FACTORS (difference in DATA/MC efficiencies) :
  // ---------------------------------------------------------
  // Numbers used for Summer 2003 pre-blessing :
  //double cmu_fac  = 0.986;  // CDF-6347
  //double cmp_fac  = 0.974;  // CDF-6347
  //double cmup_fac = 0.962;  // CDF-6347
  //double cmx_fac  = 1.000;  // CDF-6347
  //double cmio_fac = 1.000;  // CDF-6517
  // TCE : no scale factor (CDF-6262)
  // PHX : no scale factor (Chris Hill)

  // Numbers used for Summer 2003 blessing :
  double cmu_fac  = 0.898; // CDF-6612 Table 12.
  double cmp_fac  = 0.920; // CDF-6612 Table 13.
  double cmup_fac = 0.927; // CDF-6612 Table 4, option 1.
  double cmx_fac  = 0.992; // CDF-6612 Table 10, 
  double cmio_fac = 1.000;

  if(summer_acceptance){
    // E-E
    // ===
    luminosity[0]  = lum_nosi;      allowed_combination[0] = ee;    trigger_acceptance[0] = 1.0-(1.0-trig_tce)*(1.0-trig_tce)     ; eff_fac[0]=1.0;  // TCE-TCE    pb^-1	       
    luminosity[1]  = lum_si;        allowed_combination[1] = ee;    trigger_acceptance[1] = trig_tce                              ; eff_fac[1]=1.0;  // TCE-PHX    pb^-1	       
    luminosity[2]  = 0.0;           allowed_combination[2] = null;  trigger_acceptance[2] = 0.0                                   ; eff_fac[2]=1.0;  // TCE-PEM    pb^-1	       
    luminosity[3]  = 0.0;           allowed_combination[3] = null;  trigger_acceptance[3] = 0.0                                   ; eff_fac[3]=1.0;  // PHX-PHX    pb^-1	       
    luminosity[4]  = 0.0;           allowed_combination[4] = null;  trigger_acceptance[4] = 0.0                                   ; eff_fac[4]=1.0;  // PHX-PEM    pb^-1	       
    								                                                                                         			       
    // E-MU							                                                                                        			       
    // ====							                                                                                        			       
    luminosity[5]   = lum_nosi;     allowed_combination[5]  = emu;  trigger_acceptance[5] = 1.0-(1.0-trig_tce)*(1.0-trig_cmup)    ; eff_fac[5]=cmup_fac;   // TCE-CMUP   pb^-1	       
    luminosity[6]   = lum_nosi;     allowed_combination[6]  = emu;  trigger_acceptance[6] = trig_tce                              ; eff_fac[6]=cmu_fac;    // TCE-CMU    pb^-1	       
    luminosity[7]   = lum_nosi;     allowed_combination[7]  = emu;  trigger_acceptance[7] = trig_tce                              ; eff_fac[7]=cmp_fac;    // TCE-CMP    pb^-1	       
    luminosity[8]   = lum_nosicmx;  allowed_combination[8]  = emu;  trigger_acceptance[8] = 1.0-(1.0-trig_tce)*(1.0-trig_cmx)     ; eff_fac[8]=cmx_fac;    // TCE-CMX    pb^-1	       
    luminosity[9]   = lum_nosi;     allowed_combination[9]  = emu;  trigger_acceptance[9] = trig_tce                              ; eff_fac[9]=cmio_fac;   // TCE-CMIO   pb^-1	       
    luminosity[10]  = lum_si;       allowed_combination[10] = emu;  trigger_acceptance[10] = trig_cmup                            ; eff_fac[10]=cmup_fac;  // PHX-CMUP   pb^-1	       
    luminosity[11]  = 0.0;          allowed_combination[11] = null; trigger_acceptance[11] = 0.0                                  ; eff_fac[11]=cmu_fac;   // PHX-CMU    pb^-1	       
    luminosity[12]  = 0.0;          allowed_combination[12] = null; trigger_acceptance[12] = 0.0                                  ; eff_fac[12]=cmp_fac;   // PHX-CMP    pb^-1	       
    luminosity[13]  = lum_sicmx;    allowed_combination[13] = emu;  trigger_acceptance[13] = trig_cmx                             ; eff_fac[13]=cmx_fac;   // PHX-CMX    pb^-1	       
    luminosity[14]  = 0.0;          allowed_combination[14] = null; trigger_acceptance[14] = 0.0                                  ; eff_fac[14]=cmio_fac;  // PHX-CMIO   pb^-1	       
    luminosity[15]  = 0.0;          allowed_combination[15] = null; trigger_acceptance[15] = 0.0                                  ; eff_fac[15]=cmup_fac;  // PEM-CMUP   pb^-1	       
    luminosity[16]  = 0.0;          allowed_combination[16] = null; trigger_acceptance[16] = 0.0                                  ; eff_fac[16]=cmx_fac;   // PEM-CMX    pb^-1	       
    								                                                                                        			       
    // MU-MU 							                                                                                        			       
    // =====							                                                                                        			       
    luminosity[17]  = lum_nosi;     allowed_combination[17] = mumu;  trigger_acceptance[17] = 1.0-(1.0-trig_cmup)*(1.0-trig_cmup) ; eff_fac[17]=cmup_fac*cmup_fac; // CMUP-CMUP   pb^-1
    luminosity[18]  = lum_nosi;     allowed_combination[18] = mumu;  trigger_acceptance[18] = trig_cmup                           ; eff_fac[18]=cmup_fac*cmu_fac;  // CMUP-CMU    pb^-1
    luminosity[19]  = lum_nosi;     allowed_combination[19] = mumu;  trigger_acceptance[19] = trig_cmup                           ; eff_fac[19]=cmup_fac*cmp_fac;  // CMUP-CMP    pb^-1
    luminosity[20]  = lum_nosicmx;  allowed_combination[20] = mumu;  trigger_acceptance[20] = 1.0-(1.0-trig_cmup)*(1.0-trig_cmx)  ; eff_fac[20]=cmup_fac*cmx_fac;  // CMUP-CMX    pb^-1
    luminosity[21]  = lum_nosi;     allowed_combination[21] = mumu;  trigger_acceptance[21] = trig_cmup                           ; eff_fac[21]=cmup_fac*cmio_fac; // CMUP-CMIO   pb^-1
    luminosity[22]  = lum_nosicmx;  allowed_combination[22] = mumu;  trigger_acceptance[22] = 1.0-(1.0-trig_cmx)*(1.0-trig_cmx)   ; eff_fac[22]=cmx_fac*cmx_fac;   // CMX-CMX     pb^-1
    luminosity[23]  = lum_nosicmx;  allowed_combination[23] = mumu;  trigger_acceptance[23] = trig_cmx                            ; eff_fac[23]=cmx_fac*cmu_fac;   // CMX-CMU     pb^-1
    luminosity[24]  = lum_nosicmx;  allowed_combination[24] = mumu;  trigger_acceptance[24] = trig_cmx                            ; eff_fac[24]=cmx_fac*cmp_fac;   // CMX-CMP     pb^-1
    luminosity[25]  = lum_nosicmx;  allowed_combination[25] = mumu;  trigger_acceptance[25] = trig_cmx                            ; eff_fac[25]=cmx_fac*cmio_fac;  // CMX-CMIO    pb^-1
    
  }
  
  if(winter_acceptance){
    // E-E
    // ===
    luminosity[0]  = lum_nosi;      allowed_combination[0] = ee;   trigger_acceptance[0] = 1.0-(1.0-trig_tce)*(1.0-trig_tce)      ; eff_fac[0]=1.0;  // TCE-TCE    pb^-1	       
    luminosity[1]  = 0.0;           allowed_combination[1] = null; trigger_acceptance[1] = 0.0					  ; eff_fac[1]=1.0;  // TCE-PHX    pb^-1	       
    luminosity[2]  = 0.0;           allowed_combination[2] = null; trigger_acceptance[2] = 0.0					  ; eff_fac[2]=1.0;  // TCE-PEM    pb^-1	       
    luminosity[3]  = 0.0;           allowed_combination[3] = null; trigger_acceptance[3] = 0.0					  ; eff_fac[3]=1.0;  // PHX-PHX    pb^-1	       
    luminosity[4]  = 0.0;           allowed_combination[4] = null; trigger_acceptance[4] = 0.0					  ; eff_fac[4]=1.0;  // PHX-PEM    pb^-1	       
    								                                                                                       			                            
    // E-MU							                                                                                      			                            
    // ====							                                                                                      			                            
    luminosity[5]   = lum_nosi;     allowed_combination[5]  = emu; trigger_acceptance[5] = 1.0-(1.0-trig_tce)*(1.0-trig_cmup)	  ; eff_fac[5]=cmup_fac;   // TCE-CMUP   pb^-1	       
    luminosity[6]   = lum_nosi;     allowed_combination[6]  = emu; trigger_acceptance[6] = trig_tce				  ; eff_fac[6]=cmu_fac;    // TCE-CMU    pb^-1	       
    luminosity[7]   = lum_nosi;     allowed_combination[7]  = emu; trigger_acceptance[7] = trig_tce				  ; eff_fac[7]=cmp_fac;    // TCE-CMP    pb^-1	       
    luminosity[8]   = lum_nosicmx;  allowed_combination[8]  = emu; trigger_acceptance[8] = 1.0-(1.0-trig_tce)*(1.0-trig_cmx)	  ; eff_fac[8]=cmx_fac;    // TCE-CMX    pb^-1	       
    luminosity[9]   = 0.0;          allowed_combination[9]  = null;trigger_acceptance[9] = 0.0					  ; eff_fac[9]=cmio_fac;   // TCE-CMIO   pb^-1	       
    luminosity[10]  = 0.0;          allowed_combination[10] = null;trigger_acceptance[10] = 0.0					  ; eff_fac[10]=cmup_fac;  // PHX-CMUP   pb^-1	       
    luminosity[11]  = 0.0;          allowed_combination[11] = null;trigger_acceptance[11] = 0.0					  ; eff_fac[11]=cmu_fac;   // PHX-CMU    pb^-1	       
    luminosity[12]  = 0.0;          allowed_combination[12] = null;trigger_acceptance[12] = 0.0					  ; eff_fac[12]=cmp_fac;   // PHX-CMP    pb^-1	       
    luminosity[13]  = 0.0;          allowed_combination[13] = null;trigger_acceptance[13] = 0.0					  ; eff_fac[13]=cmx_fac;   // PHX-CMX    pb^-1	       
    luminosity[14]  = 0.0;          allowed_combination[14] = null;trigger_acceptance[14] = 0.0					  ; eff_fac[14]=cmio_fac;  // PHX-CMIO   pb^-1	       
    luminosity[15]  = 0.0;          allowed_combination[15] = null;trigger_acceptance[15] = 0.0					  ; eff_fac[15]=cmup_fac;  // PEM-CMUP   pb^-1	       
    luminosity[16]  = 0.0;          allowed_combination[16] = null;trigger_acceptance[16] = 0.0					  ; eff_fac[16]=cmx_fac;   // PEM-CMX    pb^-1	       
    								                                                                                      			                            
    // MU-MU 							                                                                                      			                            
    // =====							                                                                                      			                            
    luminosity[17]  = lum_nosi;     allowed_combination[17] = mumu; trigger_acceptance[17] = 1.0-(1.0-trig_cmup)*(1.0-trig_cmup)  ; eff_fac[17]=cmup_fac*cmup_fac; // CMUP-CMUP   pb^-1
    luminosity[18]  = lum_nosi;     allowed_combination[18] = mumu; trigger_acceptance[18] = trig_cmup				  ; eff_fac[18]=cmup_fac*cmu_fac;  // CMUP-CMU    pb^-1
    luminosity[19]  = lum_nosi;     allowed_combination[19] = mumu; trigger_acceptance[19] = trig_cmup				  ; eff_fac[19]=cmup_fac*cmp_fac;  // CMUP-CMP    pb^-1
    luminosity[20]  = lum_nosicmx;  allowed_combination[20] = mumu; trigger_acceptance[20] = 1.0-(1.0-trig_cmup)*(1.0-trig_cmx)	  ; eff_fac[20]=cmup_fac*cmx_fac;  // CMUP-CMX    pb^-1
    luminosity[21]  = 0.0;          allowed_combination[21] = null; trigger_acceptance[21] = 0.0				  ; eff_fac[21]=cmup_fac*cmio_fac; // CMUP-CMIO   pb^-1
    luminosity[22]  = lum_nosicmx;  allowed_combination[22] = mumu; trigger_acceptance[22] = 1.0-(1.0-trig_cmx)*(1.0-trig_cmx)	  ; eff_fac[22]=cmx_fac*cmx_fac;   // CMX-CMX     pb^-1
    luminosity[23]  = lum_nosicmx;  allowed_combination[23] = mumu; trigger_acceptance[23] = trig_cmx				  ; eff_fac[23]=cmx_fac*cmu_fac;   // CMX-CMU     pb^-1
    luminosity[24]  = lum_nosicmx;  allowed_combination[24] = mumu; trigger_acceptance[24] = trig_cmx				  ; eff_fac[24]=cmx_fac*cmp_fac;   // CMX-CMP     pb^-1
    luminosity[25]  = 0.0;          allowed_combination[25] = null; trigger_acceptance[25] = 0.0				  ; eff_fac[25]=cmx_fac*cmio_fac;  // CMX-CMIO    pb^-1
    
    
  }




  // For 100pb^-1 tables :
  // ---------------------
  // for (int i=0; i<27; ++i) {
  //   luminosity[i] = 100.0;
  // }  
  
  // Print ee, mumu and tautau summary tables :
  // ==========================================
  cout << " *** Drell-Yan e+e- *** " << endl; 
  tabled dy_ee_cut_totals(8,rowd(3));
  tabled dy_ee_cut_fractions(8,rowd(3));
  tabled dy_ee_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    dy_ee_cut_totals[i][0] = 0;
    dy_ee_cut_totals[i][1] = 0;
    dy_ee_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	dy_ee_cut_totals[i][0] += dy_ee_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	dy_ee_cut_totals[i][1] += dy_ee_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	dy_ee_cut_totals[i][2] += dy_ee_results[i][j];  // MU-MU
      }
    }

    if (i>0) {
      if (dy_ee_cut_totals[i-1][0] > 0) {
	dy_ee_cut_fractions[i][0]=dy_ee_cut_totals[i][0]/dy_ee_cut_totals[i-1][0];
	dy_ee_cut_errors[i][0]=sqrt((dy_ee_cut_fractions[i][0]*(1-dy_ee_cut_fractions[i][0]))/dy_ee_cut_totals[i-1][0]);
      }
      if (dy_ee_cut_totals[i-1][1] > 0) {
	dy_ee_cut_fractions[i][1]=dy_ee_cut_totals[i][1]/dy_ee_cut_totals[i-1][1];
	dy_ee_cut_errors[i][1]=sqrt((dy_ee_cut_fractions[i][1]*(1-dy_ee_cut_fractions[i][1]))/dy_ee_cut_totals[i-1][1]);
      }
      if (dy_ee_cut_totals[i-1][2] > 0) {
	dy_ee_cut_fractions[i][2]=dy_ee_cut_totals[i][2]/dy_ee_cut_totals[i-1][2];
	dy_ee_cut_errors[i][2]=sqrt((dy_ee_cut_fractions[i][2]*(1-dy_ee_cut_fractions[i][2]))/dy_ee_cut_totals[i-1][2]);
      }
      
    }

    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",dy_ee_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",dy_ee_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",dy_ee_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",dy_ee_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",dy_ee_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",dy_ee_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",dy_ee_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",dy_ee_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",dy_ee_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }



  cout << endl << " *** Drell-Yan mu+mu- *** " << endl; 
  tabled dy_mumu_cut_totals(8,rowd(3));
  tabled dy_mumu_cut_fractions(8,rowd(3));
  tabled dy_mumu_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    dy_mumu_cut_totals[i][0] = 0;
    dy_mumu_cut_totals[i][1] = 0;
    dy_mumu_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	dy_mumu_cut_totals[i][0] += dy_mumu_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	dy_mumu_cut_totals[i][1] += dy_mumu_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	dy_mumu_cut_totals[i][2] += dy_mumu_results[i][j];  // MU-MU
      }
    }

    if (i>0) {
      if (dy_mumu_cut_totals[i-1][0] > 0) {
	dy_mumu_cut_fractions[i][0]=dy_mumu_cut_totals[i][0]/dy_mumu_cut_totals[i-1][0];
	dy_mumu_cut_errors[i][0]=sqrt((dy_mumu_cut_fractions[i][0]*(1-dy_mumu_cut_fractions[i][0]))/dy_mumu_cut_totals[i-1][0]);
      }
      if (dy_mumu_cut_totals[i-1][1] > 0) {
	dy_mumu_cut_fractions[i][1]=dy_mumu_cut_totals[i][1]/dy_mumu_cut_totals[i-1][1];
	dy_mumu_cut_errors[i][1]=sqrt((dy_mumu_cut_fractions[i][1]*(1-dy_mumu_cut_fractions[i][1]))/dy_mumu_cut_totals[i-1][1]);
      }
      if (dy_mumu_cut_totals[i-1][2] > 0) {
	dy_mumu_cut_fractions[i][2]=dy_mumu_cut_totals[i][2]/dy_mumu_cut_totals[i-1][2];
	dy_mumu_cut_errors[i][2]=sqrt((dy_mumu_cut_fractions[i][2]*(1-dy_mumu_cut_fractions[i][2]))/dy_mumu_cut_totals[i-1][2]);
      }

    }

    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",dy_mumu_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",dy_mumu_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",dy_mumu_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",dy_mumu_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",dy_mumu_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",dy_mumu_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",dy_mumu_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",dy_mumu_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",dy_mumu_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }

  cout << endl << " *** Drell-Yan tau+tau- *** " << endl; 
  tabled dy_tautau_cut_totals(8,rowd(3));
  tabled dy_tautau_cut_fractions(8,rowd(3));
  tabled dy_tautau_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    dy_tautau_cut_totals[i][0] = 0;
    dy_tautau_cut_totals[i][1] = 0;
    dy_tautau_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	dy_tautau_cut_totals[i][0] += dy_tautau_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	dy_tautau_cut_totals[i][1] += dy_tautau_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	dy_tautau_cut_totals[i][2] += dy_tautau_results[i][j];  // MU-MU
      }
    }

    if (i>0) {
      if (dy_tautau_cut_totals[i-1][0] > 0) {
	dy_tautau_cut_fractions[i][0]=dy_tautau_cut_totals[i][0]/dy_tautau_cut_totals[i-1][0];
	dy_tautau_cut_errors[i][0]=sqrt((dy_tautau_cut_fractions[i][0]*(1-dy_tautau_cut_fractions[i][0]))/dy_tautau_cut_totals[i-1][0]);
      }
      if (dy_tautau_cut_totals[i-1][1] > 0) {
	dy_tautau_cut_fractions[i][1]=dy_tautau_cut_totals[i][1]/dy_tautau_cut_totals[i-1][1];
	dy_tautau_cut_errors[i][1]=sqrt((dy_tautau_cut_fractions[i][1]*(1-dy_tautau_cut_fractions[i][1]))/dy_tautau_cut_totals[i-1][1]);
      }
      if (dy_tautau_cut_totals[i-1][2] > 0) {
	dy_tautau_cut_fractions[i][2]=dy_tautau_cut_totals[i][2]/dy_tautau_cut_totals[i-1][2];
	dy_tautau_cut_errors[i][2]=sqrt((dy_tautau_cut_fractions[i][2]*(1-dy_tautau_cut_fractions[i][2]))/dy_tautau_cut_totals[i-1][2]);
      }

    }

    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",dy_tautau_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",dy_tautau_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",dy_tautau_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",dy_tautau_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",dy_tautau_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",dy_tautau_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",dy_tautau_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",dy_tautau_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",dy_tautau_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }
  

  cout << endl << " *** ttbar *** " << endl; 
  tabled ttbar_cut_totals(8,rowd(3));
  tabled ttbar_cut_fractions(8,rowd(3));
  tabled ttbar_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    ttbar_cut_totals[i][0] = 0;
    ttbar_cut_totals[i][1] = 0;
    ttbar_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	ttbar_cut_totals[i][0] += ttbar_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	ttbar_cut_totals[i][1] += ttbar_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	ttbar_cut_totals[i][2] += ttbar_results[i][j];  // MU-MU
      }
    }

    if (i>0) {
      if (ttbar_cut_totals[i-1][0] > 0) {
	ttbar_cut_fractions[i][0]=ttbar_cut_totals[i][0]/ttbar_cut_totals[i-1][0];
	ttbar_cut_errors[i][0]=sqrt((ttbar_cut_fractions[i][0]*(1-ttbar_cut_fractions[i][0]))/ttbar_cut_totals[i-1][0]);
      }
      if (ttbar_cut_totals[i-1][1] > 0) {
	ttbar_cut_fractions[i][1]=ttbar_cut_totals[i][1]/ttbar_cut_totals[i-1][1];
	ttbar_cut_errors[i][1]=sqrt((ttbar_cut_fractions[i][1]*(1-ttbar_cut_fractions[i][1]))/ttbar_cut_totals[i-1][1]);
      }
      if (ttbar_cut_totals[i-1][2] > 0) {
	ttbar_cut_fractions[i][2]=ttbar_cut_totals[i][2]/ttbar_cut_totals[i-1][2];
	ttbar_cut_errors[i][2]=sqrt((ttbar_cut_fractions[i][2]*(1-ttbar_cut_fractions[i][2]))/ttbar_cut_totals[i-1][2]);
      }

    }

    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",ttbar_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",ttbar_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",ttbar_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",ttbar_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",ttbar_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",ttbar_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",ttbar_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",ttbar_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",ttbar_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }
  
  
  cout << endl << " *** WZ *** " << endl; 
  tabled wz_cut_totals(8,rowd(3));
  tabled wz_cut_fractions(8,rowd(3));
  tabled wz_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    wz_cut_totals[i][0] = 0;
    wz_cut_totals[i][1] = 0;
    wz_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	wz_cut_totals[i][0] += wz_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	wz_cut_totals[i][1] += wz_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	wz_cut_totals[i][2] += wz_results[i][j];  // MU-MU
      }
    }
    
    if (i>0) {
      if (wz_cut_totals[i-1][0] > 0) {
	wz_cut_fractions[i][0]=wz_cut_totals[i][0]/wz_cut_totals[i-1][0];
	wz_cut_errors[i][0]=sqrt((wz_cut_fractions[i][0]*(1-wz_cut_fractions[i][0]))/wz_cut_totals[i-1][0]);
      }
      if (wz_cut_totals[i-1][1] > 0) {
	wz_cut_fractions[i][1]=wz_cut_totals[i][1]/wz_cut_totals[i-1][1];
	wz_cut_errors[i][1]=sqrt((wz_cut_fractions[i][1]*(1-wz_cut_fractions[i][1]))/wz_cut_totals[i-1][1]);
      }
      if (wz_cut_totals[i-1][2] > 0) {
	wz_cut_fractions[i][2]=wz_cut_totals[i][2]/wz_cut_totals[i-1][2];
	wz_cut_errors[i][2]=sqrt((wz_cut_fractions[i][2]*(1-wz_cut_fractions[i][2]))/wz_cut_totals[i-1][2]);
      }
      
    }
    
    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",wz_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",wz_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",wz_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",wz_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",wz_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",wz_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",wz_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",wz_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",wz_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }
  
  
  cout << endl << " *** WW *** " << endl; 
  tabled ww_cut_totals(8,rowd(3));
  tabled ww_cut_fractions(8,rowd(3));
  tabled ww_cut_errors(8,rowd(3));
  for (int i=0; i<8; ++i) {

    ww_cut_totals[i][0] = 0;
    ww_cut_totals[i][1] = 0;
    ww_cut_totals[i][2] = 0;

    for (int j=0; j<26; ++j) {
      if (allowed_combination[j]==ee) {
	ww_cut_totals[i][0] += ww_results[i][j];  // E-E
      } else if (allowed_combination[j]==emu) {
	ww_cut_totals[i][1] += ww_results[i][j];  // E-MU
      } else if (allowed_combination[j]==mumu) {
	ww_cut_totals[i][2] += ww_results[i][j];  // MU-MU
      }
    }

    if (i>0) {
      if (ww_cut_totals[i-1][0] > 0) {
	ww_cut_fractions[i][0]=ww_cut_totals[i][0]/ww_cut_totals[i-1][0];
	ww_cut_errors[i][0]=sqrt((ww_cut_fractions[i][0]*(1-ww_cut_fractions[i][0]))/ww_cut_totals[i-1][0]);
      }
      if (ww_cut_totals[i-1][1] > 0) {
	ww_cut_fractions[i][1]=ww_cut_totals[i][1]/ww_cut_totals[i-1][1];
	ww_cut_errors[i][1]=sqrt((ww_cut_fractions[i][1]*(1-ww_cut_fractions[i][1]))/ww_cut_totals[i-1][1]);
      }
      if (ww_cut_totals[i-1][2] > 0) {
	ww_cut_fractions[i][2]=ww_cut_totals[i][2]/ww_cut_totals[i-1][2];
	ww_cut_errors[i][2]=sqrt((ww_cut_fractions[i][2]*(1-ww_cut_fractions[i][2]))/ww_cut_totals[i-1][2]);
      }
    }

    char formatted_1[25];
    char formatted_2[25];
    char formatted_3[25];
    sprintf(formatted_1,"%12.5f",ww_cut_totals[i][0]);
    sprintf(formatted_2,"%7.4f",ww_cut_fractions[i][0]*100.0);
    sprintf(formatted_3,"%7.4f",ww_cut_errors[i][0]*100.0);
    cout << cuts[i] 
	 << " : ee = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) ";
    sprintf(formatted_1,"%12.5f",ww_cut_totals[i][2]);
    sprintf(formatted_2,"%7.4f",ww_cut_fractions[i][2]*100.0);
    sprintf(formatted_3,"%7.4f",ww_cut_errors[i][2]*100.0);
    cout << "   mm = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "; 
    sprintf(formatted_1,"%12.5f",ww_cut_totals[i][1]);
    sprintf(formatted_2,"%7.4f",ww_cut_fractions[i][1]*100.0);
    sprintf(formatted_3,"%7.4f",ww_cut_errors[i][1]*100.0);
    cout << "   em = " << formatted_1 << " ( " << formatted_2 << " +- " << formatted_3 << " %) "
	 << endl;
  }


  // Calculate expectations :
  // ========================
  int final_cut = 7; // position in cut array

  double dy_ee_final[4];
  double dy_ee_stat_error[4];
  double dy_ee_final_error_uncorrelated[4];
  double dy_ee_final_error[4];
  double dy_ee_mc_factor = 1./(dy_ee_nevt_total/(dy_ee_filter_eff*dy_ee_cs*dy_ee_k_factor));
  dy_ee_final[0] = 0; dy_ee_stat_error[0] = 0;
  dy_ee_final[1] = 0; dy_ee_stat_error[1] = 0;
  dy_ee_final[2] = 0; dy_ee_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      dy_ee_final[0]      += dy_ee_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor;
      dy_ee_stat_error[0] += pow(sqrt(dy_ee_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      dy_ee_final[1]      += dy_ee_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor;
      dy_ee_stat_error[1] += pow(sqrt(dy_ee_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      dy_ee_final[2]      += dy_ee_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor;
      dy_ee_stat_error[2] += pow(sqrt(dy_ee_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_ee_mc_factor,2);
    }
  }
  dy_ee_stat_error[0]  = sqrt(dy_ee_stat_error[0]);
  dy_ee_stat_error[1]  = sqrt(dy_ee_stat_error[1]);
  dy_ee_stat_error[2]  = sqrt(dy_ee_stat_error[2]);
  dy_ee_final_error_uncorrelated[0] = sqrt( pow((dy_ee_stat_error[0]),2) +
					    pow((dy_ee_final[0]*dy_ee_systematic),2) ); // assumes uncorrelated systematic errors. Probably wrong.
  dy_ee_final_error[0]              = sqrt( pow((dy_ee_stat_error[0]),2) +
					    pow((dy_ee_final[0]*dy_ee_systematic),2) +
					    pow((dy_ee_final[0]*luminosity_systematic),2) ); // luminosity piece
  dy_ee_final_error_uncorrelated[1] = sqrt( pow((dy_ee_stat_error[1]),2) +
					    pow((dy_ee_final[1]*dy_ee_systematic),2) );
  dy_ee_final_error[1]              = sqrt( pow((dy_ee_stat_error[1]),2) +
					    pow((dy_ee_final[1]*dy_ee_systematic),2) +
					    pow((dy_ee_final[1]*luminosity_systematic),2) ); // luminosity piece            
  dy_ee_final_error_uncorrelated[2] = sqrt( pow((dy_ee_stat_error[2]),2) +
					    pow((dy_ee_final[2]*dy_ee_systematic),2) );
  dy_ee_final_error[2]              = sqrt( pow((dy_ee_stat_error[2]),2) +
					    pow((dy_ee_final[2]*dy_ee_systematic),2) +
					    pow((dy_ee_final[2]*luminosity_systematic),2) ); // luminosity piece            
  
  dy_ee_final[3]=dy_ee_final[0]+dy_ee_final[1]+dy_ee_final[2];
  dy_ee_final_error_uncorrelated[3]=sqrt( pow(dy_ee_final_error_uncorrelated[0],2) + 
					  pow(dy_ee_final_error_uncorrelated[1],2) +
					  pow(dy_ee_final_error_uncorrelated[2],2) );
  dy_ee_final_error[3]=sqrt( pow(dy_ee_final_error_uncorrelated[3],2) +
			     pow(dy_ee_final[3]*luminosity_systematic,2) );
  

  double dy_mumu_final[4];
  double dy_mumu_stat_error[4];
  double dy_mumu_final_error_uncorrelated[4];
  double dy_mumu_final_error[4];
  double dy_mumu_mc_factor = 1./(dy_mumu_nevt_total/(dy_mumu_filter_eff*dy_mumu_cs*dy_mumu_k_factor));
  dy_mumu_final[0] = 0; dy_mumu_stat_error[0] = 0;
  dy_mumu_final[1] = 0; dy_mumu_stat_error[1] = 0;
  dy_mumu_final[2] = 0; dy_mumu_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      dy_mumu_final[0]      += dy_mumu_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor;
      dy_mumu_stat_error[0] += pow(sqrt(dy_mumu_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      dy_mumu_final[1]      += dy_mumu_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor;
      dy_mumu_stat_error[1] += pow(sqrt(dy_mumu_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      dy_mumu_final[2]      += dy_mumu_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor;
      dy_mumu_stat_error[2] += pow(sqrt(dy_mumu_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_mumu_mc_factor,2);
    }
  }
  dy_mumu_stat_error[0]  = sqrt(dy_mumu_stat_error[0]);
  dy_mumu_stat_error[1]  = sqrt(dy_mumu_stat_error[1]);
  dy_mumu_stat_error[2]  = sqrt(dy_mumu_stat_error[2]);
  dy_mumu_final_error_uncorrelated[0] = sqrt( pow((dy_mumu_stat_error[0]),2) +
					      pow((dy_mumu_final[0]*dy_mumu_systematic),2) );
  dy_mumu_final_error[0]              = sqrt( pow((dy_mumu_stat_error[0]),2) +
					      pow((dy_mumu_final[0]*dy_mumu_systematic),2) +
					      pow((dy_mumu_final[0]*luminosity_systematic),2) ); // luminosity piece
  dy_mumu_final_error_uncorrelated[1] = sqrt( pow((dy_mumu_stat_error[1]),2) +
					      pow((dy_mumu_final[1]*dy_mumu_systematic),2) );
  dy_mumu_final_error[1]              = sqrt( pow((dy_mumu_stat_error[1]),2) +
					      pow((dy_mumu_final[1]*dy_mumu_systematic),2) +
					      pow((dy_mumu_final[1]*luminosity_systematic),2) ); // luminosity piece            
  dy_mumu_final_error_uncorrelated[2] = sqrt( pow((dy_mumu_stat_error[2]),2) +
					      pow((dy_mumu_final[2]*dy_mumu_systematic),2) );
  dy_mumu_final_error[2]              = sqrt( pow((dy_mumu_stat_error[2]),2) +
					      pow((dy_mumu_final[2]*dy_mumu_systematic),2) +
					      pow((dy_mumu_final[2]*luminosity_systematic),2) ); // luminosity piece            
  
  dy_mumu_final[3]=dy_mumu_final[0]+dy_mumu_final[1]+dy_mumu_final[2];
  dy_mumu_final_error_uncorrelated[3]=sqrt( pow(dy_mumu_final_error_uncorrelated[0],2) + 
					    pow(dy_mumu_final_error_uncorrelated[1],2) +
					    pow(dy_mumu_final_error_uncorrelated[2],2) );
  dy_mumu_final_error[3]=sqrt( pow(dy_mumu_final_error_uncorrelated[3],2) +
			       pow(dy_mumu_final[3]*luminosity_systematic,2) );
  
  

  double dy_tautau_final[4];
  double dy_tautau_stat_error[4];
  double dy_tautau_final_error_uncorrelated[4];
  double dy_tautau_final_error[4];
  double dy_tautau_mc_factor = 1./(dy_tautau_nevt_total/(dy_tautau_filter_eff*dy_tautau_cs*dy_tautau_k_factor));
  dy_tautau_final[0] = 0; dy_tautau_stat_error[0] = 0;
  dy_tautau_final[1] = 0; dy_tautau_stat_error[1] = 0;
  dy_tautau_final[2] = 0; dy_tautau_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      dy_tautau_final[0]      += dy_tautau_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor;
      dy_tautau_stat_error[0] += pow(sqrt(dy_tautau_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      dy_tautau_final[1]      += dy_tautau_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor;
      dy_tautau_stat_error[1] += pow(sqrt(dy_tautau_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      dy_tautau_final[2]      += dy_tautau_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor;
      dy_tautau_stat_error[2] += pow(sqrt(dy_tautau_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*dy_tautau_mc_factor,2);
    }
  }
  dy_tautau_stat_error[0]  = sqrt(dy_tautau_stat_error[0]);
  dy_tautau_stat_error[1]  = sqrt(dy_tautau_stat_error[1]);
  dy_tautau_stat_error[2]  = sqrt(dy_tautau_stat_error[2]);
  dy_tautau_final_error_uncorrelated[0] = sqrt( pow((dy_tautau_stat_error[0]),2) +
						pow((dy_tautau_final[0]*dy_tautau_systematic),2) );
  dy_tautau_final_error[0]              = sqrt( pow((dy_tautau_stat_error[0]),2) +
						pow((dy_tautau_final[0]*dy_tautau_systematic),2) +
						pow((dy_tautau_final[0]*luminosity_systematic),2) ); // luminosity piece
  dy_tautau_final_error_uncorrelated[1] = sqrt( pow((dy_tautau_stat_error[1]),2) +
						pow((dy_tautau_final[1]*dy_tautau_systematic),2) );
  dy_tautau_final_error[1]              = sqrt( pow((dy_tautau_stat_error[1]),2) +
						pow((dy_tautau_final[1]*dy_tautau_systematic),2) +
						pow((dy_tautau_final[1]*luminosity_systematic),2) ); // luminosity piece            
  dy_tautau_final_error_uncorrelated[2] = sqrt( pow((dy_tautau_stat_error[2]),2) +
						pow((dy_tautau_final[2]*dy_tautau_systematic),2) );
  dy_tautau_final_error[2]              = sqrt( pow((dy_tautau_stat_error[2]),2) +
						pow((dy_tautau_final[2]*dy_tautau_systematic),2) +
						pow((dy_tautau_final[2]*luminosity_systematic),2) ); // luminosity piece            
  
  dy_tautau_final[3]=dy_tautau_final[0]+dy_tautau_final[1]+dy_tautau_final[2];
  dy_tautau_final_error_uncorrelated[3]=sqrt( pow(dy_tautau_final_error_uncorrelated[0],2) + 
					      pow(dy_tautau_final_error_uncorrelated[1],2) +
					      pow(dy_tautau_final_error_uncorrelated[2],2) );
  dy_tautau_final_error[3]=sqrt( pow(dy_tautau_final_error_uncorrelated[3],2) +
				 pow(dy_tautau_final[3]*luminosity_systematic,2) );
  
  

  double ttbar_final[4];
  double ttbar_stat_error[4];
  double ttbar_final_error_uncorrelated[4];
  double ttbar_final_error[4];
  double ttbar_mc_factor = 1./(ttbar_nevt_total/(ttbar_filter_eff*ttbar_cs*ttbar_k_factor));
  ttbar_final[0] = 0; ttbar_stat_error[0] = 0;
  ttbar_final[1] = 0; ttbar_stat_error[1] = 0;
  ttbar_final[2] = 0; ttbar_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      ttbar_final[0]      += ttbar_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor;
      ttbar_stat_error[0] += pow(sqrt(ttbar_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      ttbar_final[1]      += ttbar_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor;
      ttbar_stat_error[1] += pow(sqrt(ttbar_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      ttbar_final[2]      += ttbar_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor;
      ttbar_stat_error[2] += pow(sqrt(ttbar_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ttbar_mc_factor,2);
    }
  }
  ttbar_stat_error[0]  = sqrt(ttbar_stat_error[0]);
  ttbar_stat_error[1]  = sqrt(ttbar_stat_error[1]);
  ttbar_stat_error[2]  = sqrt(ttbar_stat_error[2]);
  ttbar_final_error_uncorrelated[0] = sqrt( pow((ttbar_stat_error[0]),2) +
					    pow((ttbar_final[0]*ttbar_systematic),2) );
  ttbar_final_error[0]              = sqrt( pow((ttbar_stat_error[0]),2) +
					    pow((ttbar_final[0]*ttbar_systematic),2) +
					    pow((ttbar_final[0]*luminosity_systematic),2) ); // luminosity piece
  ttbar_final_error_uncorrelated[1] = sqrt( pow((ttbar_stat_error[1]),2) +
					    pow((ttbar_final[1]*ttbar_systematic),2) );
  ttbar_final_error[1]              = sqrt( pow((ttbar_stat_error[1]),2) +
					    pow((ttbar_final[1]*ttbar_systematic),2) +
					    pow((ttbar_final[1]*luminosity_systematic),2) ); // luminosity piece            
  ttbar_final_error_uncorrelated[2] = sqrt( pow((ttbar_stat_error[2]),2) +
					    pow((ttbar_final[2]*ttbar_systematic),2) );
  ttbar_final_error[2]              = sqrt( pow((ttbar_stat_error[2]),2) +
					    pow((ttbar_final[2]*ttbar_systematic),2) +
					    pow((ttbar_final[2]*luminosity_systematic),2) ); // luminosity piece            
  
  ttbar_final[3]=ttbar_final[0]+ttbar_final[1]+ttbar_final[2];
  ttbar_final_error_uncorrelated[3]=sqrt( pow(ttbar_final_error_uncorrelated[0],2) + 
					  pow(ttbar_final_error_uncorrelated[1],2) +
					  pow(ttbar_final_error_uncorrelated[2],2) );
  ttbar_final_error[3]=sqrt( pow(ttbar_final_error_uncorrelated[3],2) +
			     pow(ttbar_final[3]*luminosity_systematic,2) );
  

  double wz_final[4];
  double wz_stat_error[4]; 
  double wz_final_error_uncorrelated[4];
  double wz_final_error[4];
  double wz_mc_factor = 1./(wz_nevt_total/(wz_filter_eff*wz_cs*wz_k_factor));
  wz_final[0] = 0; wz_stat_error[0] = 0;
  wz_final[1] = 0; wz_stat_error[1] = 0;
  wz_final[2] = 0; wz_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      wz_final[0]      += wz_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor;
      wz_stat_error[0] += pow(sqrt(wz_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      wz_final[1]      += wz_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor;
      wz_stat_error[1] += pow(sqrt(wz_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      wz_final[2]      += wz_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor;
      wz_stat_error[2] += pow(sqrt(wz_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*wz_mc_factor,2);
    }
  }
  wz_stat_error[0]  = sqrt(wz_stat_error[0]);
  wz_stat_error[1]  = sqrt(wz_stat_error[1]);
  wz_stat_error[2]  = sqrt(wz_stat_error[2]);
  wz_final_error_uncorrelated[0] = sqrt( pow((wz_stat_error[0]),2) +
					 pow((wz_final[0]*wz_systematic),2) );
  wz_final_error[0]              = sqrt( pow((wz_stat_error[0]),2) +
					 pow((wz_final[0]*wz_systematic),2) +
					 pow((wz_final[0]*luminosity_systematic),2) ); // luminosity piece
  wz_final_error_uncorrelated[1] = sqrt( pow((wz_stat_error[1]),2) +
					 pow((wz_final[1]*wz_systematic),2) );
  wz_final_error[1]              = sqrt( pow((wz_stat_error[1]),2) +
					 pow((wz_final[1]*wz_systematic),2) +
					 pow((wz_final[1]*luminosity_systematic),2) ); // luminosity piece            
  wz_final_error_uncorrelated[2] = sqrt( pow((wz_stat_error[2]),2) +
					 pow((wz_final[2]*wz_systematic),2) );
  wz_final_error[2]              = sqrt( pow((wz_stat_error[2]),2) +
					 pow((wz_final[2]*wz_systematic),2) +
					 pow((wz_final[2]*luminosity_systematic),2) ); // luminosity piece            
  
  wz_final[3]=wz_final[0]+wz_final[1]+wz_final[2];
  wz_final_error_uncorrelated[3]=sqrt( pow(wz_final_error_uncorrelated[0],2) + 
				       pow(wz_final_error_uncorrelated[1],2) +
				       pow(wz_final_error_uncorrelated[2],2) );
  wz_final_error[3]=sqrt( pow(wz_final_error_uncorrelated[3],2) +
			  pow(wz_final[3]*luminosity_systematic,2) );
  
  double ww_final[4];
  double ww_stat_error[4];
  double ww_final_error[4];
  double ww_mc_factor = 1./(ww_nevt_total/(ww_filter_eff*ww_cs*ww_k_factor));
  ww_final[0] = 0; ww_stat_error[0] = 0;
  ww_final[1] = 0; ww_stat_error[1] = 0;
  ww_final[2] = 0; ww_stat_error[2] = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i]==ee) {
      ww_final[0]      += ww_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor;
      ww_stat_error[0] += pow(sqrt(ww_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor,2);
    } else if (allowed_combination[i]==emu) {
      ww_final[1]      += ww_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor;
      ww_stat_error[1] += pow(sqrt(ww_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor,2);
    } else if (allowed_combination[i]==mumu) {
      ww_final[2]      += ww_results[final_cut][i]*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor;
      ww_stat_error[2] += pow(sqrt(ww_results[final_cut][i])*luminosity[i]*trigger_acceptance[i]*eff_fac[i]*ww_mc_factor,2);
    }
  }
  ww_stat_error[0]  = sqrt(ww_stat_error[0]);
  ww_stat_error[1]  = sqrt(ww_stat_error[1]);
  ww_stat_error[2]  = sqrt(ww_stat_error[2]);
  ww_final_error[0] = sqrt( pow((ww_stat_error[0]),2) +
			    pow((ww_final[0]*ww_systematic),2) +
			    pow((ww_final[0]*luminosity_systematic),2) ); // luminosity piece
  ww_final_error[1] = sqrt( pow((ww_stat_error[1]),2) +
			    pow((ww_final[1]*ww_systematic),2) +
			    pow((ww_final[1]*luminosity_systematic),2) ); // luminosity piece            
  ww_final_error[2] = sqrt( pow((ww_stat_error[2]),2) +
			    pow((ww_final[2]*ww_systematic),2) +
			    pow((ww_final[2]*luminosity_systematic),2) ); // luminosity piece            

  ww_final[3]=ww_final[0]+ww_final[1]+ww_final[2];
  // Take into account correlated signal acceptance systematic~:
  double ww_uncorrelated = sqrt(pow(ww_systematic,2)-pow(ww_correlated,2));
  ww_final_error[3]=sqrt( pow((ww_correlated*(ww_final[0]+ww_final[1]+ww_final[2])),2) +     // correlated piece
			  pow((ww_uncorrelated*(ww_final[0])),2) +                           // uncorrelated pieces
			  pow((ww_uncorrelated*(ww_final[1])),2) +
			  pow((ww_uncorrelated*(ww_final[2])),2) +                           
			  pow((ww_stat_error[0]),2) +                                        // statistical pieces
			  pow((ww_stat_error[1]),2) +
			  pow((ww_stat_error[2]),2) +                                        
			  pow((ww_final[0]+ww_final[1]+ww_final[2])*luminosity_systematic,2) ); // luminosity piece
  
  // Uncorrelated error calculation :
  // ww_final_error[3]=sqrt( pow(ww_final_error[0],2) + 
  //			  pow(ww_final_error[1],2) +
  //			  pow(ww_final_error[2],2) );



  // QCD Fake Background - treated separately from Monte Carlo based background estimates :
  // ======================================================================================
  double fake_final[4];
  double fake_final_stat[4];
  double fake_final_syst[4];
  double fake_final_error[4];

  // double fake_systematic = 0.85;
  // fake_final_error[0] = sqrt ( pow(0.033,2) + pow((fake_systematic*fake_final[0]),2) );
  // fake_final_error[1] = sqrt ( pow(0.374,2) + pow((fake_systematic*fake_final[1]),2) );;
  // fake_final_error[2] = sqrt ( pow(0.076,2) + pow((fake_systematic*fake_final[2]),2) );
  // fake_final[3]=fake_final[0]+fake_final[1]+fake_final[2];
  // fake_final_error[3]=sqrt( pow(fake_final_error[0],2) + 
  //	   		       pow(fake_final_error[1],2) +
  //			       pow(fake_final_error[2],2) );  

  if (winter_data && winter_acceptance) {
    fake_final[0]      = 0.041; // E-E
    fake_final_stat[0] = 0.004;
    fake_final_syst[0] = 0.001;
    fake_final[1]      = 0.113; // E-M
    fake_final_stat[1] = 0.065;
    fake_final_syst[1] = 0.020;
    fake_final[2]      = 0.049; // M-M
    fake_final_stat[2] = 0.039;
    fake_final_syst[2] = 0.012;
    fake_final[3]      = 0.204; // TOTAL
    fake_final_stat[3] = 0.076;
    fake_final_syst[3] = 0.023;
    //     fake_final[0] = 0.018; // E-E
    //     fake_final_error[0] = 0.005;
    //     fake_final[1] = 0.086; // E-MU
    //     fake_final_error[1] = 0.065;
    //     fake_final[2] = 0.049; // MU-MU
    //     fake_final_error[2] = 0.039;
    //     fake_final[3]=0.154;
    //     fake_final_error[3] = sqrt( pow(0.077,2) + pow(0.032,2) );
  } else if (summer_data && winter_acceptance) {
    fake_final[0]      = 0.069; // E-E
    fake_final_stat[0] = 0.007;
    fake_final_syst[0] = 0.002;
    fake_final[1]      = 0.251; // E-M
    fake_final_stat[1] = 0.091;
    fake_final_syst[1] = 0.034;
    fake_final[2]      = 0.162; // M-M
    fake_final_stat[2] = 0.087;
    fake_final_syst[2] = 0.034;
    fake_final[3]      = 0.481; // TOTAL
    fake_final_stat[3] = 0.126;
    fake_final_syst[3] = 0.049;
    //     fake_final[0] = 0.031; // E-E
    //     fake_final_error[0] = 0.009;
    //     fake_final[1] = 0.178; // E-MU
    //     fake_final_error[1] = 0.084;
    //     fake_final[2] = 0.162; // MU-MU
    //     fake_final_error[2] = 0.087;
    //     fake_final[3]=0.370;
    //     fake_final_error[3] = sqrt( pow(0.130,2) + pow(0.072,2) );
  } else if (winter_data && summer_acceptance) {
    fake_final[0]      = 0.079; // E-E
    fake_final_stat[0] = 0.006;
    fake_final_syst[0] = 0.019;
    fake_final[1]      = 0.175; // E-M
    fake_final_stat[1] = 0.069;
    fake_final_syst[1] = 0.025;
    fake_final[2]      = 0.049; // M-M
    fake_final_stat[2] = 0.039;
    fake_final_syst[2] = 0.012;
    fake_final[3]      = 0.302; // TOTAL
    fake_final_stat[3] = 0.080;
    fake_final_syst[3] = 0.033;
    //     fake_final[0] = 0.072; // E-E
    //     fake_final_error[0] = 0.008;
    //     fake_final[1] = 0.151; // E-MU
    //     fake_final_error[1] = 0.069;
    //     fake_final[2] = 0.049; // MU-MU
    //     fake_final_error[2] = 0.039;
    //     fake_final[3]=0.272;
    //     fake_final_error[3] = sqrt( pow(0.081,2) + pow(0.121,2) );
  } else if (summer_data && summer_acceptance) {
    fake_final[0]      = 0.155; // E-E
    fake_final_stat[0] = 0.009;
    fake_final_syst[0] = 0.043;
    fake_final[1]      = 0.399; // E-M
    fake_final_stat[1] = 0.102;
    fake_final_syst[1] = 0.047;
    fake_final[2]      = 0.162; // M-M
    fake_final_stat[2] = 0.087;
    fake_final_syst[2] = 0.034;
    fake_final[3]      = 0.716; // TOTAL
    fake_final_stat[3] = 0.134;
    fake_final_syst[3] = 0.072;
    //     fake_final[0] = 0.129; // E-E
    //     fake_final_error[0] = 0.014;
    //     fake_final[1] = 0.316; // E-MU
    //     fake_final_error[1] = 0.096;
    //     fake_final[2] = 0.162; // MU-MU
    //     fake_final_error[2] = 0.087;
    //     fake_final[3]=0.607;
    //     fake_final_error[3] = sqrt( pow(0.139,2) + pow(0.221,2) );
  }

  if (apply_syst) {
    for (int i=0; i<4; ++i) {
      fake_final_error[i] = sqrt( pow(fake_final_stat[i],2) +
				  pow(fake_final_syst[i],2) );
    }
  } else {
    for (int i=0; i<4; ++i) {
      fake_final_error[i] = fake_final_stat[i];
    }
  }
  

  // Background totals, assuming uncorrelated errors apart from luminosity :
  // =======================================================================
  double background_final[4];
  double background_final_error_uncorrelated[4];
  double background_final_error[4];
  for(int i=0; i<4; ++i) {
    background_final[i] = 
      dy_ee_final[i] +
      dy_mumu_final[i] +
      dy_tautau_final[i] + 
      wz_final[i] +
      fake_final[i] +
      ttbar_final[i];
    background_final_error_uncorrelated[i] = sqrt ( pow(dy_ee_final_error_uncorrelated[i],2) +
						    pow(dy_mumu_final_error_uncorrelated[i],2) +
						    pow(dy_tautau_final_error_uncorrelated[i],2) + 
						    pow(wz_final_error_uncorrelated[i],2) +
						    pow(fake_final_error[i],2) +
						    pow(ttbar_final_error_uncorrelated[i],2) );	
    background_final_error[i] = sqrt ( pow(background_final_error_uncorrelated[i],2) +
				       pow( (dy_ee_final[i] +
					     dy_mumu_final[i] +
					     dy_tautau_final[i] +
					     wz_final[i] +
					     ttbar_final[i])*luminosity_systematic,2) ); // luminosity systematic
  }


  // Print-out expectations :
  // ========================
  cout << "\n\n";
  cout << " E-E RESULTS : " << endl;
  cout << " ============= " << endl;
  cout << "   DY-ee            : " << dy_ee_final[0] << " +- " << dy_ee_final_error[0] << endl;
  cout << "   DY-mumu          : " << dy_mumu_final[0] << " +- " << dy_mumu_final_error[0]<< endl;
  cout << "   DY-tautau        : " << dy_tautau_final[0] << " +- " << dy_tautau_final_error[0] << endl;
  cout << "   WZ               : " << wz_final[0] << " +- " << wz_final_error[0] << endl;
  cout << "   Fake             : " << fake_final[0] << " +- " << fake_final_error[0] << endl;
  cout << "   ttbar            : " << ttbar_final[0] << " +- " << ttbar_final_error[0] << endl;
  cout << "   TOTAL BKGND.     : " << background_final[0] << " +- " << background_final_error[0] << endl;
  cout << "\n";
  cout << "   WW               : " << ww_final[0] << " +- " << ww_final_error[0] << endl;

  cout << "\n";

  cout << " MU-MU RESULTS : " << endl;
  cout << " =============== " << endl;
  cout << "   DY-ee            : " << dy_ee_final[2] << " +- " << dy_ee_final_error[2] << endl;
  cout << "   DY-mumu          : " << dy_mumu_final[2] << " +- " << dy_mumu_final_error[2] << endl;
  cout << "   DY-tautau        : " << dy_tautau_final[2] << " +- " << dy_tautau_final_error[2] << endl;
  cout << "   WZ               : " << wz_final[2] << " +- " << wz_final_error[2] << endl;
  cout << "   Fake             : " << fake_final[2] << " +- " << fake_final_error[2] << endl;
  cout << "   ttbar            : " << ttbar_final[2] << " +- " << ttbar_final_error[2] << endl;
  cout << "   TOTAL BKGND.     : " << background_final[2] << " +- " << background_final_error[2] << endl;
  cout << "\n";
  cout << "   WW               : " << ww_final[2] << " +- " << ww_final_error[2] << endl;

  cout << "\n";

  cout << " E-MU RESULTS : " << endl;
  cout << " ============== " << endl;
  cout << "   DY-ee            : " << dy_ee_final[1] << " +- " << dy_ee_final_error[1] << endl;
  cout << "   DY-mumu          : " << dy_mumu_final[1] << " +- " << dy_mumu_final_error[1] << endl;
  cout << "   DY-tautau        : " << dy_tautau_final[1] << " +- " << dy_tautau_final_error[1] << endl;
  cout << "   WZ               : " << wz_final[1] << " +- " << wz_final_error[1] << endl;
  cout << "   Fake             : " << fake_final[1] << " +- " << fake_final_error[1] << endl;
  cout << "   ttbar            : " << ttbar_final[1] << " +- " << ttbar_final_error[1] << endl;
  cout << "   TOTAL BKGND.     : " << background_final[1] << " +- " << background_final_error[1] << endl;
  cout << "\n";
  cout << "   WW               : " << ww_final[1] << " +- " << ww_final_error[1] << endl;

  cout << "\n";

  cout << " TOTAL RESULTS : " << endl;
  cout << " =============== " << endl;
  cout << "   DY-ee            : " << dy_ee_final[3] << " +- " << dy_ee_final_error[3] << endl;
  cout << "   DY-mumu          : " << dy_mumu_final[3] << " +- " << dy_mumu_final_error[3] << endl;
  cout << "   DY-tautau        : " << dy_tautau_final[3] << " +- " << dy_tautau_final_error[3] << endl;
  cout << "   WZ               : " << wz_final[3] << " +- " << wz_final_error[3] << endl;
  cout << "   Fake             : " << fake_final[3] << " +- " << fake_final_error[3] << endl;
  cout << "   ttbar            : " << ttbar_final[3] << " +- " << ttbar_final_error[3] << endl;
  cout << "   TOTAL BKGND.     : " << background_final[3] << " +- " << background_final_error[3] << endl;
  cout << "\n";
  cout << "   WW               : " << ww_final[3] << " +- " << ww_final_error[3] << endl;

  
  label cuts_latex(8);
  cuts_latex[0]=" Lepton ID";
  cuts_latex[1]=" Isolation";
  cuts_latex[2]=" Conv+Cosmic";
  cuts_latex[3]=" $Z$-veto"; 
  cuts_latex[4]=" $\\not\\!\\!{E}_{T}>$25 GeV   : ";
  cuts_latex[5]=" $\\Delta\\Phi>20^{\\circ}$ if $\\not\\!\\!{E}_{T}<$50";
  cuts_latex[6]=" 0 jets";
  cuts_latex[7]=" Opposite Sign";
  char formatted_1_ee[50];
  char formatted_2_ee[50];  
  char formatted_3_ee[50];  
  char formatted_1_mm[50];
  char formatted_2_mm[50];  
  char formatted_3_mm[50];  
  char formatted_1_em[50];
  char formatted_2_em[50];  
  char formatted_3_em[50];  

  cout << "\n\n WW SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{$WW$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(ww_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",ww_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",ww_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(ww_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",ww_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",ww_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(ww_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",ww_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",ww_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  cout << "\n\n WW ACCEPTANCE SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c|c||c|}" << endl;
  cout << "\\hline" << endl;
  cout << "$A$ & $A_{ee}$ & $A_{{\\mu}{\\mu}}$ & $A_{e{\\mu}}$ & $A_{ee+{\\mu}{\\mu}+e{\\mu}}$ \\\\ \\hline" << endl;
  double acc_ee = ww_cut_totals[final_cut][0]/(ww_nevt_total);
  double acc_mm = ww_cut_totals[final_cut][2]/(ww_nevt_total);
  double acc_em = ww_cut_totals[final_cut][1]/(ww_nevt_total);
  double acc_tot = (ww_cut_totals[final_cut][0] + ww_cut_totals[final_cut][2] + ww_cut_totals[final_cut][1])/(ww_nevt_total);
  // Binomial errors :
  double acc_ee_error = sqrt((acc_ee*(1.0-acc_ee))/ww_nevt_total);
  double acc_mm_error = sqrt((acc_mm*(1.0-acc_mm))/ww_nevt_total);
  double acc_em_error = sqrt((acc_em*(1.0-acc_em))/ww_nevt_total);
  double acc_tot_error = sqrt((acc_tot*(1.0-acc_tot))/ww_nevt_total);
  // Scale by (ww_filter_eff) to get total WW acceptance estimates, and multiply by 100 to get percentages :
  acc_ee *= 100.0*ww_filter_eff;
  acc_mm *= 100.0*ww_filter_eff;
  acc_em *= 100.0*ww_filter_eff;
  acc_tot *= 100.0*ww_filter_eff;
  acc_ee_error *= 100.0*ww_filter_eff;
  acc_mm_error *= 100.0*ww_filter_eff;
  acc_em_error *= 100.0*ww_filter_eff;
  acc_tot_error *= 100.0*ww_filter_eff;
  // Relative acceptances :  
  double rel_ee = 100.0*(acc_ee/acc_tot);
  double rel_mm = 100.0*(acc_mm/acc_tot);
  double rel_em = 100.0*(acc_em/acc_tot);
  double rel_tot = 100.0*(acc_tot/acc_tot);
  
  char formatted_acc_ee[50];
  char formatted_acc_mm[50];  
  char formatted_acc_em[50];  
  char formatted_acc_tot[50];  
  char formatted_acc_ee_error[50];
  char formatted_acc_mm_error[50];  
  char formatted_acc_em_error[50];  
  char formatted_acc_tot_error[50];  
  char formatted_rel_ee[50];
  char formatted_rel_mm[50];  
  char formatted_rel_em[50];  
  char formatted_rel_tot[50];  
  sprintf(formatted_acc_ee,"%5.3f",acc_ee);
  sprintf(formatted_acc_mm,"%5.3f",acc_mm);
  sprintf(formatted_acc_em,"%5.3f",acc_em);
  sprintf(formatted_acc_tot,"%5.3f",acc_tot);
  sprintf(formatted_acc_ee_error,"%5.3f",acc_ee_error);
  sprintf(formatted_acc_mm_error,"%5.3f",acc_mm_error);
  sprintf(formatted_acc_em_error,"%5.3f",acc_em_error);
  sprintf(formatted_acc_tot_error,"%5.3f",acc_tot_error);
  sprintf(formatted_rel_ee,"%4.1f",rel_ee);
  sprintf(formatted_rel_mm,"%4.1f",rel_mm);
  sprintf(formatted_rel_em,"%4.1f",rel_em);
  sprintf(formatted_rel_tot,"%4.1f",rel_tot);

  cout << " & (" << formatted_acc_ee << " $\\pm$ " << formatted_acc_ee_error << ") \\%" 
       << " & (" << formatted_acc_mm << " $\\pm$ " << formatted_acc_mm_error << ") \\%" 
       << " & (" << formatted_acc_em << " $\\pm$ " << formatted_acc_em_error << ") \\%" 
       << " & (" << formatted_acc_tot << " $\\pm$ " << formatted_acc_tot_error << ") \\%"  
       << " \\\\ \\hline" << endl;
  cout << " $A_{r}$ & " << formatted_rel_ee << " \\%" 
       << " & " << formatted_rel_mm << " \\%" 
       << " & " << formatted_rel_em << " \\%" 
       << " & " << formatted_rel_tot << " \\%" 
       << " \\\\ \\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;
 
  // Calculate relative acceptance by region :
  // -----------------------------------------
  int cctot = 0;
  int cptot = 0;
  int pptot = 0;
  for (int i=0; i<26; ++i) {
    if (allowed_combination[i] != null) {
      if (region[i]=="CC") cctot += int(ww_results[final_cut][i]);
      if (region[i]=="CP") cptot += int(ww_results[final_cut][i]);
      if (region[i]=="PP") pptot += int(ww_results[final_cut][i]);
    }
  }
  cout << "WW ACCEPTANCE : CCTOT = " << cctot << endl;
  cout << "WW ACCEPTANCE : CPTOT = " << cptot << endl;
  cout << "WW ACCEPTANCE : PPTOT = " << pptot << endl;

  cout << "\n\n EE SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{Drell-Yan $e^{+}e^{-}$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(dy_ee_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",dy_ee_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",dy_ee_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(dy_ee_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",dy_ee_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",dy_ee_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(dy_ee_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",dy_ee_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",dy_ee_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;

  cout << "\n\n MUMU SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{Drell-Yan $\\mu^{+}\\mu^{-}$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(dy_mumu_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",dy_mumu_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",dy_mumu_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(dy_mumu_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",dy_mumu_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",dy_mumu_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(dy_mumu_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",dy_mumu_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",dy_mumu_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  cout << "\n\n TAUTAU SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{Drell-Yan $\\tau^{+}\\tau^{-}$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(dy_tautau_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",dy_tautau_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",dy_tautau_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(dy_tautau_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",dy_tautau_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",dy_tautau_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(dy_tautau_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",dy_tautau_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",dy_tautau_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  cout << "\n\n WZ SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{$WZ$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(wz_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",wz_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",wz_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(wz_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",wz_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",wz_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(wz_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",wz_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",wz_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  cout << "\n\n TTBAR SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.30}" << endl;
  cout << "\\large" << endl;
  cout << "\\begin{tabular}{|l||c|c||c|c||c|c|}" << endl;
  cout << "\\hline" << endl;
  cout << "&  \\multicolumn{6}{|c|}{$t\\bar{t}$} \\\\" << endl;
  cout << "\\cline{2-7} & \\multicolumn{2}{|c||}{ee} & \\multicolumn{2}{|c||}{${\\mu}{\\mu}$} & \\multicolumn{2}{|c|}{$e{\\mu}$} \\\\" << endl;
  cout << "\\cline{2-7} Cut & & $\\%$  &  & $\\%$  &  & $\\%$ \\\\ \\hline \\hline" << endl;
  for (int i=0; i<8; ++i) {  
    sprintf(formatted_1_ee,"%7i",int(ttbar_cut_totals[i][0]));
    sprintf(formatted_2_ee,"%6.2f",ttbar_cut_fractions[i][0]*100.0);
    sprintf(formatted_3_ee,"%6.2f",ttbar_cut_errors[i][0]*100.0);
    
    sprintf(formatted_1_mm,"%7i",int(ttbar_cut_totals[i][2]));
    sprintf(formatted_2_mm,"%6.2f",ttbar_cut_fractions[i][2]*100.0);
    sprintf(formatted_3_mm,"%6.2f",ttbar_cut_errors[i][2]*100.0);

    sprintf(formatted_1_em,"%7i",int(ttbar_cut_totals[i][1]));
    sprintf(formatted_2_em,"%6.2f",ttbar_cut_fractions[i][1]*100.0);
    sprintf(formatted_3_em,"%6.2f",ttbar_cut_errors[i][1]*100.0);

    cout << cuts_latex[i] << " & "
	 << formatted_1_ee << " & " << formatted_2_ee << " $\\pm$ " << formatted_3_ee << " & " 
	 << formatted_1_mm << " & " << formatted_2_mm << " $\\pm$ " << formatted_3_mm << " & "
	 << formatted_1_em << " & " << formatted_2_em << " $\\pm$ " << formatted_3_em << " \\\\ " << endl; 
  }
  cout << "\\hline" << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  // Grand Summary latex table :
  // ===========================
  cout << "\n\n GRAND SUMMARY TABLE \n\n" << endl;
  cout << "\\begin{table}[h]" << endl;
  cout << "\\renewcommand{\\arraystretch}{1.40}" << endl;
  cout << "\\large" << endl;
  cout << "\\centering" << endl;
  cout << "\\begin{tabular}{|l||c|c|c||c|}" << endl;
  cout << "\\hline" << endl;
  cout << " & \\multicolumn{4}{|c|}{CDF Run~${\\rm II}$ Summer 2003 Preliminary} \\\\ \\hline \\hline" << endl;
  cout << "Source & $ee$ & ${\\mu}{\\mu}$ & $e{\\mu}$ & ${\\ell}{\\ell}$ \\\\ \\hline" << endl;

  char formatted_central[4][50];
  char formatted_error[4][50];
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",dy_ee_final[i]);
    sprintf(formatted_error[i],"%7.4f",dy_ee_final_error[i]);
  }
  cout << "Drell-Yan $e^{+}e^{-}$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",dy_mumu_final[i]);
    sprintf(formatted_error[i],"%7.4f",dy_mumu_final_error[i]);
  }
  cout << "Drell-Yan $\\mu^{+}\\mu^{-}$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",dy_tautau_final[i]);
    sprintf(formatted_error[i],"%7.4f",dy_tautau_final_error[i]);
  }
  cout << "Drell-Yan $\\tau^{+}\\tau^{-}$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",wz_final[i]);
    sprintf(formatted_error[i],"%7.4f",wz_final_error[i]);
  }
  cout << "$WZ$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",fake_final[i]);
    sprintf(formatted_error[i],"%7.4f",fake_final_error[i]);
  }
  cout << "Fake " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",ttbar_final[i]);
    sprintf(formatted_error[i],"%7.4f",ttbar_final_error[i]);
  }
  cout << "$t\\bar{t}$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ \\hline " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",background_final[i]);
    sprintf(formatted_error[i],"%7.4f",background_final_error[i]);
  }
  cout << "Total Background " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ \\hline \\hline " << endl;
  for (int i=0; i<4; ++i) {
    sprintf(formatted_central[i],"%7.4f",ww_final[i]);
    sprintf(formatted_error[i],"%7.4f",ww_final_error[i]);
  }
  cout << "$WW\\rightarrow {\\rm dileptons}$ " 
       << " & " << formatted_central[0] << " $\\pm$ " << formatted_error[0] 
       << " & " << formatted_central[2] << " $\\pm$ " << formatted_error[2] 
       << " & " << formatted_central[1] << " $\\pm$ " << formatted_error[1] 
       << " & " << formatted_central[3] << " $\\pm$ " << formatted_error[3] 
       << " \\\\ \\hline \\hline " << endl;
  cout << "{\\bf Run 2 Data} "
       << " & " << data_ee 
       << " & " << data_mm 
       << " & " << data_em
       << " & " << (data_ee+data_mm+data_em)
       << " \\\\ \\hline  " << endl;
  cout << "\\end{tabular}" << endl;
  cout << "\\end{table}" << endl;


  // Calculate Feldman-Cousins confidence interval :
  // -----------------------------------------------
  int nobs          = data_ee+data_mm+data_em;
  double background = background_final[3];
  double confidence = 0.95; // "1-sigma"
  double signal     = 0.0;
  double up_error   = 0.0;
  double down_error = 0.0;
  feldman_cousins(nobs,background,confidence,signal,up_error,down_error);
  // Turn this into a cross section estimate :
  double ww_cs_sf   = ww_cs/ww_final[3];
  double ww_cs_meas = signal*ww_cs_sf;
  double ww_cs_up   = up_error*ww_cs_sf;
  double ww_cs_down = down_error*ww_cs_sf;

  cout << "\n  Feldman-Cousins cross section estimate : " << endl;
  cout << "  ---------------------------------------- " << endl;
  cout << "  sigma(WW) = " << ww_cs_meas << " + " << ww_cs_up << " - " << ww_cs_down << "\n" << endl;

  // Now propagate the errors in a simple minded way :
  // -------------------------------------------------
  // Efficiency : propagate error through equation : sigma = (N_obs-N_bk)/eff*L
  //              this gives a systematic error = sigma * (d-eff/eff)
  double ww_cs_eff_syst = ww_cs_meas*ww_systematic;
  // Luminosity : propagates the same way as the efficiency :
  double ww_cs_lum_syst = ww_cs_meas*luminosity_systematic;
  // Background uncertainty : propagation gives error = (d-N_bk)/eff*L
  // "background_final_error_uncorrelated" doesn't contain luminosity uncertainty
  double ww_cs_bk_syst = background_final_error_uncorrelated[3]*ww_cs_sf;  
  
  cout << "  Efficiency systematic = " << ww_cs_eff_syst << endl;
  cout << "  Background systematic = " << ww_cs_bk_syst << endl;
  cout << "  Luminosity systematic = " << ww_cs_lum_syst << endl;
  
  double ww_cs_syst_combined = sqrt(pow(ww_cs_eff_syst,2)+pow(ww_cs_bk_syst,2));

  cout << "  Final result : " << endl;
  cout << "  -------------- " << endl;
  cout << "  sigma(WW) = " << ww_cs_meas << " + " << ww_cs_up << " - " << ww_cs_down
       << " (stat) +- " << ww_cs_syst_combined << " (syst) +- " << ww_cs_lum_syst << " (lum) pb " << endl;

  // feldman_cousins(0,3.0,0.90,signal,up_error,down_error);


  cout << "  \n\n  Observation Significance : " << endl;
  cout << "  -------------------------- " << endl;
  

  // Find CL at which CI become one sided :
  // --------------------------------------
  if (down_error < signal) {
    double new_cl = confidence;
    double old_cl = confidence;
    double cl_hi = 1.0;
    double cl_lo = new_cl;
    double signal_new = signal;
    double up_error_new = up_error;
    double down_error_new = down_error; 
    while (true) {
      if (down_error_new < signal_new) {
	cl_lo  = new_cl;
	new_cl = (new_cl+cl_hi)/2.0;
      } else {
	cl_hi  = new_cl;
	new_cl = (new_cl+cl_lo)/2.0; // down_error_new == signal_new
      }
      feldman_cousins(nobs,background,new_cl,signal_new,up_error_new,down_error_new);      
      if (fabs(new_cl-old_cl) < 0.0001) break;
      old_cl = new_cl;
    }

    cout << "  CI at which CL becomes one sided : " << new_cl << endl;
  }
    
  // Probability that background gives >= nobs events :
  // --------------------------------------------------
  double bkprob = 1.0;
  for (int i=0; i<nobs; ++i) {
    bkprob -= poisson(0.0,background,i);
  }
  cout << "  Poisson probability of fixed background (mu_b = " << background << ") giving >= nobs (" << nobs << ") : " << bkprob << "\n" << endl;
  // Allow the background to vary :
  TRandom3 ranGen;
  ranGen.SetSeed(123456);
  double rmssum = 0.0;
  double bkprob_exp_sum = 0.0;
  Int_t nexp = 100000;
  for (int iexp = 0; iexp < nexp; ++iexp) {
    double bk_ran = rangauss(background, background_final_error[3], &ranGen);
    rmssum += pow((bk_ran-background),2);
    // For this experiment, calculate the probability of the background giving >= nobs events :
    double bkprob_exp = 1.0;
    for (int i=0; i<nobs; ++i) {
	bkprob_exp -= poisson(0.0,bk_ran,i);
      }
    bkprob_exp_sum += bkprob_exp;
  }
  cout << "  Background uncertainty                                                          = " << background_final_error[3] << endl;
  cout << "  Background RMS from pseudo-experiments                                          = " << sqrt(rmssum/nexp) << endl;
  cout << "  Average probability of background giving >= nobs events from pseudo-experiments = " << (bkprob_exp_sum/nexp) << "\n" << endl;

  return 1;

}

static void feldman_cousins(int n0, double bg, double cl, double& signal, double& up_error, double& down_error) {
  
  // Following the prescription in "A Unified Approach to the Classical Statistical Analysis of Small Signals",
  // Feldman & Cousins, physics/9711021.
  // WARNINGS :
  // ----------
  // (1) Code not tested under all circumstances
  // (2) There are "mild pathologies" referred to by Feldman & Cousins. ALWAYS check the results
  //     against tables 2-9 in the above publication.


  cout << "\n\n  ********************************************************" << endl;
  cout << "  * Feldman-Cousins interval calculation                 *" << endl;
  cout << "  ********************************************************" << endl;
  cout << "  n-observed          : " << n0 << endl;
  cout << "  background          : " << bg << endl;
  cout << "  confidence level    : " << cl << endl;
  
  //   Brute force scanning method :
  //   -----------------------------
  //   double mu_max = n0+2.0*sqrt(n0);
  //   // double mu_min = max(0.0,n0-2.0*sqrt(n0));
  //   double mu_min = 0.0;
  //   double mu_lo  = 5000.0;
  //   double mu_hi  = -5000.0;
  //   int nmu = 2000;
  //   for (int imu=0; imu<nmu; imu++) {
  //     double mu = mu_min+((mu_max-mu_min)/nmu)*imu;
  //     double lo = 0.0;
  //     double hi = 0.0;
  //     cbelt(mu, bg, cl, lo, hi);
  //     if (n0>=lo && n0<=hi) {
  //       // This mu must be included in the range :
  //       if (mu<mu_lo) mu_lo=mu;
  //       if (mu>mu_hi) mu_hi=mu;
  //     }
  //   }

  // A coarse check to find which values of mu have a belt covering n0 :
  double mu_lo    = 0.0;
  double mu_hi    = 0.0;
  double low_out  = 0.0;   // low-mu outside the CI
  double low_in   = -1.0;  // low-mu inside the CI
  double high_in  = 0.0;   // high-mu inside the CI
  double high_out = -1.0;  // high-mu outside the CI
  for (int imu=0; imu<100; ++imu) {
    double mu = 1.0*imu;
    double lo = 0.0;
    double hi = 0.0;
    cbelt(mu, bg, cl, lo, hi);
    if (n0>hi) low_out=mu;
    if (n0>=lo && n0<=hi && low_in==-1.0) low_in=mu;
    if (n0>=lo && n0<=hi) high_in=mu;
    if (n0<lo && high_out==-1.0) high_out=mu;
  }

  // Detailed search for the mu-boundaries :
  bool found_lo = false;
  bool found_hi = false;
  double thresh = 0.0001;
  double mu_guess_old = 10000.0;
  while(!found_lo) {
    // new best estimate for low mu :
    double mu_guess = (low_out+low_in)/2.0;
    double lo = 0.0;
    double hi = 0.0;
    cbelt(mu_guess, bg, cl, lo, hi);
    if (n0>hi) {
      low_out = mu_guess;
    } else if (n0<=hi) {
      low_in = mu_guess;
    }
    if (fabs(mu_guess-mu_guess_old)<thresh) {
      mu_lo = mu_guess;
      found_lo = true;
    }    
    mu_guess_old = mu_guess;
  }

  mu_guess_old = 10000.0;
  while(!found_hi) {
    // new best estimate for high mu :
    double mu_guess = (high_out+high_in)/2.0;
    double lo = 0.0;
    double hi = 0.0;
    cbelt(mu_guess, bg, cl, lo, hi);
    if (n0>=lo) {
      high_in = mu_guess;
    } else if (n0<lo) {
      high_out = mu_guess;
    }
    if (fabs(mu_guess-mu_guess_old)<thresh) {
      mu_hi = mu_guess;
      found_hi = true;
    }
    mu_guess_old = mu_guess;
  }
  
  signal     = TMath::Max(0.0,n0-bg);
  up_error   = mu_hi-signal;
  down_error = signal-mu_lo;
  cout << "  confidence interval : [ " << mu_lo << " , " << mu_hi << " ]" << endl;
  cout << "  signal              : " << signal << " + " << up_error << " - " << down_error << endl;
  cout << "  WARNING : ALWAYS CHECK AGAINST TABLES IN physics/9711021 " << endl;
  cout << "  ********************************************************\n\n" << endl;
}

static double poisson(double mean, double bg, int n) {
  return pow((mean+bg),n)*exp(-1.0*(mean+bg))/TMath::Factorial(n);
}

void cbelt(double mu, double bg, double cl, double& low, double& high) {
  // Calculate the likelihood ratio R :
  std::vector<ratio> ratios;
  for (int i=0; i<1000; ++i) {
    double mu_best   = TMath::Max(0.0, i-bg);
    double pnmu      = poisson(mu, bg, i);
    double pnmu_best = poisson(mu_best, bg, i);
    ratio e;
    e.n = i;
    if (pnmu_best > 0.0) {
      e.v = pnmu/pnmu_best;
    }
    else {
      e.v = 0.0;
    }
    ratios.push_back(e);
  }
  // Rank the ratio elements :
  std::sort(ratios.begin(), ratios.end(), largerRatio);

  double ptot = 0.0;
  low  = 5000.0;
  high = -5000.0;
  // Now find the smallest possible conservative range :
  for (std::vector<ratio>::iterator ir = ratios.begin(); ir != ratios.end(); ++ir) {
    ptot += poisson(mu, bg, (*ir).n);
    if ((*ir).n < low)  low  = (*ir).n;
    if ((*ir).n > high) high = (*ir).n;
    if (ptot >= cl) break;
  }

  return;
}

double rangauss(double mean, double sigma, TRandom3* ranGen) {

  // From Numerical Recipes.
  // Two Gaussian random numbers are generated at once,
  // hence the following static state information :
  static int iset=0;
  static double gset;
  double fac,rsq,v1,v2;
  
  if (iset == 0) {
    // Must regenerate :
    do {
      v1 = 2.0*(ranGen->Rndm())-1.0;
      v2 = 2.0*(ranGen->Rndm())-1.0;
      rsq = v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    
    fac=sqrt(-2.0*log(rsq)/rsq);
    gset=(v1*fac)*sigma+mean;
    iset=1;
    return (v2*fac)*sigma+mean;
  } else {
    iset=0;
    return gset;
  }
}