PageRenderTime 154ms CodeModel.GetById 16ms app.highlight 127ms RepoModel.GetById 1ms app.codeStats 0ms

/ElectroWeakAnalysis/WENu/macros/PlotCombiner.cc

https://github.com/aivanov-cern/cmssw
C++ | 1161 lines | 855 code | 52 blank | 254 comment | 122 complexity | 60f55b9a6a6ce00a85788ee2418c21e2 MD5 | raw file
   1/*
   2     Macro to make the plots .......................................
   3
   4     Instructions:
   5     a. set up an input file that looks like the following:
   6     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   7     # zee or wenu
   8     wenu
   9     # file name, type (sig, qcd, bce, gje, ewk), weight
  10     histos_wenu.root     sig     1.46
  11     histos_q20_30.root   qcd     0
  12     histos_q30_80.root   qcd     100.
  13     histos_q80_170.root  qcd     0
  14     histos_b20_30.root   bce     0
  15     histos_b30_80.root   bce     0
  16     histos_b80_170.root  bce     0
  17     histos_zee.root      ewk     0
  18     histos_wtaunu.root   ewk     0
  19     histos_ztautau.root  ewk     0
  20     histos_gj15.root     gje     0
  21     histos_gj20.root     gje     0
  22     histos_gj25.root     gje     10.12
  23     histos_gj30.root     gje     0
  24     histos_gj35.root     gje     0
  25     histos_wmunu.root    ewk     0
  26     histos_ttbar.root    ewk     0
  27     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  28     lines that start with # are considered to be comments
  29     line 2 has wenu or zee. From line 4 the list of the histo files are listed
  30     (first word) then a type that could be sig,qcd,bce, gj or ewk in order to
  31     discriminate among different sources of bkgs and finally the weight that we
  32     want to weight the histogram entries. This particular example is for Wenu. For
  33     Zee one has to put type sig in the zee file and ewk in the Wenu file. The order
  34     of the files is arbitrary. Files with weight 0 will be ignored.
  35     After you have set up this code you run a root macro to combine the plots.
  36     You can do (not recommended - it actually crushes - to be debugged)
  37     root -b PlotCombiner.cc 
  38     or to compile it within root (recommended)
  39     root -b
  40     root [1] .L PlotCombiner.cc++
  41     root [2] PlotCombiner()
  42     
  43     and you finally get the plots.
  44
  45     For the ABCD method:
  46     ^^^^^^^^^^^^^^^^^^^^
  47     you have to insert in the 2nd line instead of wenu or zee the keyword abcd(...)
  48     The files should contain ewk samples, sig samples and qcd samples (but also read
  49     later). The only absolutely necessary files are the sig ones.
  50     Example:
  51     abcd(I=0.95,dI=0.01,Fz=0.6,dFz=0.01,FzP=0.56, dFzP=0.2,ewkerror=0.1,METCut=30.,mc)
  52     These parameters keep the same notation as in the note. The last parameter (data)
  53     can take 3 values:
  54     data: calculate in ABCD as in data. This means that the histograms denoted with
  55           sig,qcd,bce,gje  are used as of the same kind and ewk as the MC ewk. 
  56           The background is substructed as in data
  57     mcOnly: here we ignore all the input parameters I, dI etc. All parameters are taken
  58           from MC by forcing Fqcd=1
  59     mc:   <THIS IS WHAT ONE NORMALLY USES> input mc samples, calculation of statistical
  60           and systematics as in CMS AN 2009/004, systematic and statistic error 
  61	   calculation. This option also creates the plots of the variation of the
  62           signal prediction vs the parameter variation. In order to set the limits of
  63           the desired variation you have to edit the values in line 113 of this code
  64           (they are hardwired in the code)
  65     TO DO:
  66     functionalities to plot more kind of plots, e.g. efficiencies
  67     
  68     
  69     Further Questions/Contact:
  70     
  71         nikolaos.rompotis @ cern.ch
  72
  73
  74
  75	 Nikolaos Rompotis - 29 June 09
  76	 18 Sept 09:  1st updgrade: input files in a text file
  77	 28 May  10:  bug in IMET corrected, thanks to Sadia Khalil
  78	 Imperial College London
  79	 
  80	 
  81*/
  82
  83
  84#include <iostream>
  85#include <fstream>
  86#include <vector>
  87#include <string>
  88#include "TString.h"
  89#include "TROOT.h"
  90#include "TStyle.h"
  91#include "TH1F.h"
  92#include "TFile.h"
  93#include "TCanvas.h"
  94#include "TGraph.h"
  95#include "TLegend.h"
  96
  97void plotMaker(TString histoName, TString typeOfplot,
  98	       vector<TString> file, vector<TString> type, 
  99	       vector<double> weight,  TString xtitle);
 100
 101void abcd(vector<TString> file, vector<TString> type, vector<double> weight,
 102	  double METCut, double I, double dI, double Fz, double dFz, 
 103	  double FzP, double dFzP, double ewkerror,
 104	  double data, double mc, double mcOnly);
 105double  searchABCDstring(TString abcdString, TString keyword);
 106double Trionym(double a, double b, double c, double sum);
 107double CalcABCD
 108(double I, double Fz, double FzP, double K, double ewk,
 109 double Na_, double Nb_, double Nc_, double Nd_,
 110 double Ea_, double Eb_, double Ec_, double Ed_);
 111
 112// values for systematics plots: it is fraction of the MC value
 113const double EWK_SYST_MIN = 0.3;
 114const double EWK_SYST_MAX = 0.3;
 115//
 116const double I_SYST_MIN = 0.05;
 117const double I_SYST_MAX = 0.05;
 118//
 119const double FZ_SYST_MIN = 0.1;
 120const double FZ_SYST_MAX = 0.1;
 121//
 122const double FZP_SYST_MIN = 0.1;
 123const double FZP_SYST_MAX = 0.1;
 124//
 125const double K_SYST_MIN = 0.8;
 126const double K_SYST_MAX = 0.8;
 127
 128
 129using namespace std;
 130
 131void PlotCombiner()
 132{
 133  // read the file
 134  ifstream input("inputFiles");
 135  int i = 0;
 136  TString typeOfplot = "";
 137  vector<TString> types;
 138  vector<TString> files;
 139  vector<double> weights;
 140
 141  if (input.is_open()) {
 142    std::string myline;
 143    while (! input.eof()) {
 144      getline(input, myline);
 145      TString line(myline);
 146      TString c('#');
 147      TString empty(' ');
 148      if (line[0] != c) {
 149	++i;
 150	if (i==1) typeOfplot=line;
 151	else {
 152	  // read until you find 3 words
 153	  TString fname("");
 154	  TString ftype("");
 155	  TString fw("");
 156	  int lineSize = (int) line.Length();
 157	  int j=0;
 158	  while (j<lineSize) {
 159	    if(line[j] != empty) fname += line[j];
 160	    else break;
 161	    ++j;
 162	  }
 163	  while (j<lineSize) {
 164	    if(line[j] != empty) ftype += line[j];
 165	    else if(ftype.Length()==3) break;
 166	    ++j;
 167	  }
 168	  while (j<lineSize) {
 169	    if(line[j] != empty) fw += line[j];
 170	    else{ if(fw.Length()>0) break;}
 171	    ++j;
 172	  }
 173	  if (fname.Length() == 0) break;
 174	  files.push_back(fname);
 175	  types.push_back(ftype);
 176	  double w = fw.Atof();
 177	  weights.push_back(w);
 178	  if (w>0)
 179	    std::cout << fname << ", " << ftype << ", "<< w << std::endl;
 180	}
 181      }
 182    }
 183    input.close();
 184  }
 185  else {
 186    std::cout << "File with name inputFile was not found" << std::endl;
 187    return;
 188  }
 189
 190  // now you can launch the jobs
 191  if (typeOfplot == "wenu") {
 192    cout << "wenu plot maker" << endl;
 193    //        ====================
 194    // =====> WHICH HISTOS TO PLOT
 195    //        ====================
 196    plotMaker("h_met", typeOfplot, files, types, weights, "MET (GeV)");
 197  }
 198  else if (typeOfplot == "zee"){
 199    cout << "zee plot maker" << endl;
 200    //        ====================
 201    // =====> WHICH HISTOS TO PLOT
 202    //        ====================
 203    plotMaker("h_mee", typeOfplot, files, types, weights, "M_{ee} (GeV)");
 204  }
 205  else if (typeOfplot(0,4) == "abcd") {
 206    // now read the parameters of the ABCD method
 207    // look for parameter I and dI
 208    double I = searchABCDstring(typeOfplot, "I");
 209    double dI= searchABCDstring(typeOfplot, "dI");
 210    // look for parameter Fz
 211    double Fz = searchABCDstring(typeOfplot, "Fz");
 212    double dFz= searchABCDstring(typeOfplot, "dFz");
 213    // look for parameter FzP
 214    double FzP = searchABCDstring(typeOfplot, "FzP");
 215    double dFzP= searchABCDstring(typeOfplot, "dFzP");
 216    // look for the MET cut
 217    double METCut =searchABCDstring(typeOfplot, "METCut");
 218    // do you want data driven only?
 219    double data = searchABCDstring(typeOfplot, "data");
 220    double mc = searchABCDstring(typeOfplot, "mc");
 221    double mcOnly = searchABCDstring(typeOfplot, "mcOnly");
 222    // what is the ewk error?
 223    double ewkerror = searchABCDstring(typeOfplot, "ewkerror");
 224    // sanity check:
 225    if (METCut<0 || (data<-0.7 && mc<-0.7 && mcOnly<-0.7)) {
 226      cout << "Error in your configurtion!" << endl;
 227      if (METCut <0) cout << "Error in MET Cut" << endl;
 228      else cout << "You need to specify one mc or data or mcOnly"
 229		<< endl;
 230      abort();
 231    }
 232    if (mc>-0.7 && mc <0 && ewkerror<0) {
 233      cout << "You have specified mc option, but you have forgotten"
 234	   << " to set the ewkerror!" << endl;
 235      abort();
 236    }
 237    //        ===============================
 238    // =====> ABCD METHOD FOR BKG SUBTRACTION
 239    //        ===============================
 240    cout << "doing ABCD with input: " << typeOfplot << endl;
 241    abcd(files, types, weights, METCut, I, dI, Fz, dFz, FzP, dFzP,
 242	 ewkerror, data, mc, mcOnly);
 243
 244  }
 245  // force the program to abort in order to clear the memory
 246  // and avoid further use of the interpreter after
 247  abort();
 248
 249}
 250
 251void abcd( vector<TString> file, vector<TString> type, vector<double> weight, 
 252	   double METCut, double I, double dI, double Fz, double dFz, 
 253	   double FzP, double dFzP, double ewkerror,
 254	   double data, double mc, double mcOnly)
 255{
 256  gROOT->Reset();
 257  gROOT->ProcessLine(".L tdrstyle.C"); 
 258  gROOT->ProcessLine("setTDRStyle()");
 259  //
 260  std::cout << "Trying ABCD method for Background subtration" << std::endl;
 261  //
 262  // histogram names to use:
 263  TString histoName_Ba("h_met_EB");
 264  TString histoName_Bb("h_met_inverse_EB");
 265  TString histoName_Ea("h_met_EE");
 266  TString histoName_Eb("h_met_inverse_EE");
 267  //
 268  // find one file and get the dimensions of your histogram
 269  int fmax = (int) file.size();
 270  int NBins = 0; double min = 0; double max = -1;
 271  for (int i=0; i<fmax; ++i) {
 272    if (weight[i]>0) {
 273      //      cout << "Loading file " << file[i] << endl;
 274      TFile f(file[i]);
 275      TH1F *h = (TH1F*) f.Get(histoName_Ba);
 276      NBins = h->GetNbinsX();
 277      min = h->GetBinLowEdge(1);
 278      max = h->GetBinLowEdge(NBins+1);
 279      break;
 280    }
 281  }
 282  if (NBins ==0 || (max<min)) {
 283    std::cout << "PlotCombiner::abcd error: Could not find valid histograms in file"
 284	      << std::endl;
 285    abort();
 286  }
 287  cout << "Histograms with "<< NBins <<" bins  and range " << min << "-" << max  << endl;
 288  //
 289  // Wenu Signal .......................................................
 290  TH1F h_wenu("h_wenu", "h_wenu", NBins, min, max);
 291  TH1F h_wenu_inv("h_wenu_inv", "h_wenu_inv", NBins, min, max);
 292  for (int i=0; i<fmax; ++i) {
 293    if (type[i] == "sig" && weight[i]>0 ) {
 294      TFile f(file[i]);
 295      //
 296      TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
 297      h_wenu.Add(h_ba, weight[i]);
 298      TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
 299      h_wenu.Add(h_ea, weight[i]);
 300      //
 301      TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
 302      h_wenu_inv.Add(h_bb, weight[i]);
 303      TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
 304      h_wenu_inv.Add(h_eb, weight[i]);
 305    }
 306  }
 307  // QCD Bkgs
 308  TH1F h_qcd("h_qcd", "h_qcd", NBins, min, max);
 309  TH1F h_qcd_inv("h_qcd_inv", "h_qcd_inv", NBins, min, max);
 310  for (int i=0; i<fmax; ++i) {
 311    if ((type[i] == "qcd" || type[i] == "bce" 
 312	 || type[i] == "gje") && weight[i]>0) {
 313      TFile f(file[i]);
 314      //
 315      TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
 316      h_qcd.Add(h_ba, weight[i]);
 317      TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
 318      h_qcd.Add(h_ea, weight[i]);
 319      //
 320      TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
 321      h_qcd_inv.Add(h_bb, weight[i]);
 322      TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
 323      h_qcd_inv.Add(h_eb, weight[i]);
 324    }
 325  }
 326  //
 327  TH1F h_ewk("h_ewk", "h_ewk", NBins, min, max);
 328  TH1F h_ewk_inv("h_ewk_inv", "h_ewk_inv", NBins, min, max);
 329  for (int i=0; i<fmax; ++i) {
 330    if ( type[i] == "ewk" && weight[i]>0) {
 331      TFile f(file[i]);
 332      //
 333      TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
 334      h_ewk.Add(h_ba, weight[i]);
 335      TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
 336      h_ewk.Add(h_ea, weight[i]);
 337      //
 338      TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
 339      h_ewk_inv.Add(h_bb, weight[i]);
 340      TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
 341      h_ewk_inv.Add(h_eb, weight[i]);
 342    }
 343  }
 344  //
 345  // calculate the METCut position
 346  //
 347  // this is calculated as a low edge bin of your input histogram
 348  // METCut = min + (max-min)*IMET/NBins
 349  int IMET = int ((METCut - min)/(max-min) * double(NBins)); 
 350  // check whether it is indeed a low egde position
 351  double metCalc = min + (max-min)*double(IMET)/double(NBins);
 352  if (metCalc < METCut || metCalc > METCut) {
 353    std::cout << "PlotCombiner:abcd: your MET Cut is not in low egde bin position"
 354	      << std::endl;
 355  }
 356  cout << "MET Cut in " << METCut << "GeV corresponds to bin #" << IMET << endl;
 357  // Calculate the population in the ABCD Regions now
 358  // signal
 359  double a_sig = h_wenu.Integral(IMET,NBins+1);
 360  double b_sig = h_wenu.Integral(0,IMET-1);
 361  double c_sig = h_wenu_inv.Integral(0,IMET-1);
 362  double d_sig = h_wenu_inv.Integral(IMET,NBins+1);
 363  // qcd
 364  double a_qcd = h_qcd.Integral(IMET,NBins+1);
 365  double b_qcd = h_qcd.Integral(0,IMET-1);
 366  double c_qcd = h_qcd_inv.Integral(0,IMET-1);
 367  double d_qcd = h_qcd_inv.Integral(IMET,NBins+1);
 368  // ewk
 369  double a_ewk = h_ewk.Integral(IMET,NBins+1);
 370  double b_ewk = h_ewk.Integral(0,IMET-1);
 371  double c_ewk = h_ewk_inv.Integral(0,IMET-1);
 372  double d_ewk = h_ewk_inv.Integral(IMET,NBins+1);
 373  ////////////////////////////////////////////////
 374
 375  //
 376  // now the parameters of the method
 377  if (data < 0 && data >-0.75) {  // select value -0.5 that gives the
 378                                  // string parser
 379    // now everything is done from data + input
 380    std::cout << "Calculating ABCD Result and Stat Error Assuming DATA"
 381	      << std::endl << "Summary: in this implementation we have assumed"
 382	      << " that what real 'data' appear with type sig in the input"
 383	      << std::endl << "No systematics available with this type of"
 384	      << " calculation. If you need systematics try one of the other"
 385	      << " options" << std::endl;
 386    double A = (1.0-I)*(FzP-Fz);
 387    double B = I*(FzP+1.0)*(FzP*(c_sig-c_ewk)-(d_sig-d_ewk)) + 
 388      (1+Fz)*(1-I)*((a_sig-a_ewk)-dFzP*(b_sig-b_ewk));
 389    double C = I*(1.+Fz)*(1.+FzP)*((d_sig-d_ewk)*(b_sig-b_ewk) - (a_sig-a_ewk)*(c_sig-c_ewk));
 390    //
 391    // signal calculation:
 392    double S = Trionym(A,B,C, a_sig+b_sig);
 393
 394    // the errors now:
 395    // calculate the statistical error now:
 396    double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
 397    double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
 398    double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
 399    double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
 400    //
 401    double Na = a_sig, Nb = b_sig, Nc=c_sig, Nd = d_sig;
 402    double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
 403    if (A != 0) {
 404      
 405      ApI   = -(FzP-Fz);
 406      ApFz  = -(1.0-I);
 407      ApFzP = (1.0-I);
 408      ApNa  = 0.0;
 409      ApNb  = 0.0;
 410      ApNc  = 0.0;
 411      ApNd  = 0.0;
 412      
 413      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 414      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 415      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 416      BpNa  =  (1.0-I)*(1.0+Fz);
 417      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 418      BpNc  = I*(FzP+1.0)*Fz;
 419      BpNd  = -I*(FzP+1.0);
 420      
 421      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 422      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 423      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 424      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 425      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 426      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 427      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 428      
 429      SpI   = (-BpI   + (B*BpI   -2.0*ApI*C   -2.0*A*CpI)  /fabs(2.0*A*S+B)- 2.0*ApI*S)  /(2.0*A);
 430      SpFz  = (-BpFz  + (B*BpFz  -2.0*ApFz*C  -2.0*A*CpFz) /fabs(2.0*A*S+B)- 2.0*ApFz*S) /(2.0*A);
 431      SpFzP = (-BpFzP + (B*BpFzP -2.0*ApFzP*C -2.0*A*CpFzP)/fabs(2.0*A*S+B)- 2.0*ApFzP*S)/(2.0*A);
 432      SpNa  = (-BpNa  + (B*BpNa  -2.0*ApNa*C  -2.0*A*CpNa) /fabs(2.0*A*S+B)- 2.0*ApNa*S) /(2.0*A);
 433      SpNb  = (-BpNb  + (B*BpNb  -2.0*ApNb*C  -2.0*A*CpNb) /fabs(2.0*A*S+B)- 2.0*ApNb*S) /(2.0*A);
 434      SpNc  = (-BpNc  + (B*BpNc  -2.0*ApNc*C  -2.0*A*CpNc) /fabs(2.0*A*S+B)- 2.0*ApNc*S) /(2.0*A);
 435      SpNd  = (-BpNd  + (B*BpNd  -2.0*ApNd*C  -2.0*A*CpNd) /fabs(2.0*A*S+B)- 2.0*ApNd*S) /(2.0*A);
 436    }
 437    else {
 438      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 439      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 440      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 441      BpNa  =  (1.0-I)*(1.0+Fz);
 442      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 443      BpNc  = I*(FzP+1.0)*Fz;
 444      BpNd  = -I*(FzP+1.0);
 445      
 446      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 447      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 448      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 449      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 450      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 451      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 452      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 453      
 454      SpI   = -CpI/B+C*BpI/(B*B);
 455      SpFz  = -CpFz/B+C*BpFz/(B*B);
 456      SpFzP = -CpFzP/B+C*BpFzP/(B*B);
 457      SpNa  = -CpNa/B+C*BpNa/(B*B);
 458      SpNb  = -CpNb/B+C*BpNb/(B*B);
 459      SpNc  = -CpNc/B+C*BpNc/(B*B);
 460      SpNd  = -CpNd/B+C*BpNd/(B*B);
 461    }
 462    double  DS;
 463    DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
 464	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
 465    // warning: S here denotes the method prediction ..........
 466    cout << "********************************************************" << endl;
 467    cout << "Signal Prediction: " << S << "+-" << DS << "(stat)" << endl;
 468    cout << "********************************************************" << endl;
 469    cout << "Parameters used in calculation: " << endl;
 470    cout << "I=  " << I << "+-" << dI << endl;
 471    cout << "Fz= " << Fz << "+-" << dFz << endl;
 472    cout << "FzP=" << FzP << "+-" << dFzP << endl;
 473    cout << endl;
 474    cout << "ABCD Regions population:" << endl;
 475    cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
 476         << ", ewk=" << a_ewk << endl;
 477    cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
 478         << ", ewk=" << b_ewk << endl;
 479    cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
 480         << ", ewk=" << c_ewk << endl;
 481    cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
 482         << ", ewk=" << d_ewk << endl;
 483    cout << endl;
 484    //
 485    cout << "Statistical Error Summary: " << endl;
 486    cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
 487    cout << "due to FzP= "<< SpFzP*dFzP
 488	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
 489    cout << "due to  I = "<< SpI*dI
 490	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
 491    cout << "due to Na = "<< SpNa*sqrt(Na)
 492	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
 493    cout << "due to Nb = "<< SpNb*sqrt(Nb)
 494	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
 495    cout << "due to Nc = "<< SpNc*sqrt(Nc)
 496	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
 497    cout << "due to Nd = "<< SpNd*sqrt(Nd)
 498	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
 499    cout << "Total Statistical Error: " 
 500	 << DS << ", (" << DS*100./S << "%)"<< endl;
 501    cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
 502  }
 503  //
 504  //
 505  //  this is the main option of the algorithm: the one implemented in the 
 506  //  Analysis Note
 507  //
 508  if (mc < 0 && mc >-0.75) {  // select value -0.5 that gives the 
 509                              // string parser
 510    
 511    //////// STATISTICAL ERROR CALCULATION /////////////////////////////
 512    double A = (1.0-I)*(FzP-Fz);
 513    double B = I*(FzP+1.0)*(FzP*(c_sig+c_qcd)-(d_sig+d_qcd)) + 
 514      (1+Fz)*(1-I)*((a_sig+a_qcd)-dFzP*(b_sig+b_qcd));
 515    double C = I*(1.+Fz)*(1.+FzP)*((d_sig+d_qcd)*(b_sig+b_qcd) - 
 516				   (a_sig+a_qcd)*(c_sig+c_qcd));
 517    //
 518    // signal calculation:
 519    double S = Trionym(A,B,C, a_sig+b_sig);
 520    //
 521    double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
 522    double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
 523    double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
 524    double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
 525    //
 526    double Na = a_sig+a_qcd+a_ewk, Nb = b_sig+b_qcd+b_ewk;
 527    double Nc=c_sig+c_qcd+c_ewk,   Nd = d_sig+d_qcd+d_ewk;
 528    double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
 529    if (A != 0) {
 530      
 531      ApI   = -(FzP-Fz);
 532      ApFz  = -(1.0-I);
 533      ApFzP = (1.0-I);
 534      ApNa  = 0.0;
 535      ApNb  = 0.0;
 536      ApNc  = 0.0;
 537      ApNd  = 0.0;
 538      
 539      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 540      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 541      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 542      BpNa  =  (1.0-I)*(1.0+Fz);
 543      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 544      BpNc  = I*(FzP+1.0)*Fz;
 545      BpNd  = -I*(FzP+1.0);
 546      
 547      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 548      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 549      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 550      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 551      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 552      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 553      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 554      
 555      SpI   = (-BpI   + (B*BpI   -2.0*ApI*C   -2.0*A*CpI)  /fabs(2.0*A*S+B)- 2.0*ApI*S)  /(2.0*A);
 556      SpFz  = (-BpFz  + (B*BpFz  -2.0*ApFz*C  -2.0*A*CpFz) /fabs(2.0*A*S+B)- 2.0*ApFz*S) /(2.0*A);
 557      SpFzP = (-BpFzP + (B*BpFzP -2.0*ApFzP*C -2.0*A*CpFzP)/fabs(2.0*A*S+B)- 2.0*ApFzP*S)/(2.0*A);
 558      SpNa  = (-BpNa  + (B*BpNa  -2.0*ApNa*C  -2.0*A*CpNa) /fabs(2.0*A*S+B)- 2.0*ApNa*S) /(2.0*A);
 559      SpNb  = (-BpNb  + (B*BpNb  -2.0*ApNb*C  -2.0*A*CpNb) /fabs(2.0*A*S+B)- 2.0*ApNb*S) /(2.0*A);
 560      SpNc  = (-BpNc  + (B*BpNc  -2.0*ApNc*C  -2.0*A*CpNc) /fabs(2.0*A*S+B)- 2.0*ApNc*S) /(2.0*A);
 561      SpNd  = (-BpNd  + (B*BpNd  -2.0*ApNd*C  -2.0*A*CpNd) /fabs(2.0*A*S+B)- 2.0*ApNd*S) /(2.0*A);
 562    }
 563    else {
 564      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 565      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 566      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 567      BpNa  =  (1.0-I)*(1.0+Fz);
 568      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 569      BpNc  = I*(FzP+1.0)*Fz;
 570      BpNd  = -I*(FzP+1.0);
 571      
 572      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 573      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 574      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 575      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 576      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 577      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 578      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 579      
 580      SpI   = -CpI/B+C*BpI/(B*B);
 581      SpFz  = -CpFz/B+C*BpFz/(B*B);
 582      SpFzP = -CpFzP/B+C*BpFzP/(B*B);
 583      SpNa  = -CpNa/B+C*BpNa/(B*B);
 584      SpNb  = -CpNb/B+C*BpNb/(B*B);
 585      SpNc  = -CpNc/B+C*BpNc/(B*B);
 586      SpNd  = -CpNd/B+C*BpNd/(B*B);
 587    }
 588    double  DS;
 589    DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
 590	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
 591
 592    ////////////////////////////////////////////////////////////////////////
 593    // SYSTEMATICS CALCULATION /////////////////////////////////////////////
 594    ////////////////////////////////////////////////////////////////////////
 595    // recalculate the basic quantities
 596    double Imc = (a_sig + b_sig) / (a_sig + b_sig + c_sig + d_sig);
 597    double dImc = sqrt(Imc*(1-Imc)/(a_sig + b_sig + c_sig + d_sig));
 598    double Fzmc = a_sig/b_sig;
 599    double e =a_sig/(a_sig + b_sig);
 600    double de = sqrt(e*(1-e)/(a_sig + b_sig));
 601    double alpha = de/(2.*Fzmc-e);
 602    double dFzmc = alpha/(1-alpha);
 603    double FzPmc = d_sig/c_sig;
 604    double ep =d_sig/(c_sig + d_sig);
 605    double dep = sqrt(ep*(1-ep)/(c_sig + d_sig));
 606    double alphap = dep/(2.*FzPmc-ep);
 607    double dFzPmc = alphap/(1-alphap);
 608    //
 609    // calculate the K parameter as it is in MC:
 610    double KMC = (d_qcd/c_qcd)/(a_qcd/b_qcd);
 611    double SMC = a_sig + b_sig;
 612    //
 613    double dfz = Fz -Fzmc;
 614    double di = I - Imc;
 615    double dfzp = FzP - FzPmc;
 616    double fk = fabs(1-KMC);
 617    ////////////////////////////////////////////////////////////////////////
 618    // ewk error: this error has to be inserted by hand
 619    double fm = 1.-ewkerror;
 620    double fp = 1.+ewkerror;
 621    double S_EWK_PLUS = CalcABCD(Imc, Fzmc, FzPmc, KMC, fp, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 622    double S_EWK_MINUS = CalcABCD(Imc, Fzmc, FzPmc, KMC, fm, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 623    // error in K
 624    double S_K= CalcABCD(Imc, Fzmc, FzPmc, 1., 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 625    // error in Fz
 626    double S_FZ= CalcABCD(Imc, Fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 627    // error in FzP
 628    double S_FZP= CalcABCD(Imc, Fzmc, FzP, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 629    // error in I
 630    double S_I = CalcABCD(I, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 631    //
 632    // sanity tets
 633    //cout << "Smc=" << SMC<< ", " << CalcABCD(Imc, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed)
 634    // << endl;
 635    //abort();
 636    //
 637    // ************ plots for the systematics calculation ****************
 638    // ewk plot
 639    int const POINTS = 10;
 640    int const allPOINTS = 2*POINTS;
 641    TGraph g_ewk(allPOINTS);
 642    TGraph g_fz(allPOINTS);
 643    TGraph g_fzp(allPOINTS);
 644    TGraph g_k(allPOINTS);
 645    TGraph g_i(allPOINTS);
 646    //
 647    double ewk_syst_min = EWK_SYST_MIN; // because this is just fraction
 648    double i_syst_min = Imc*(1.-I_SYST_MIN);
 649    double fz_syst_min = Fzmc*(1.-FZ_SYST_MIN);
 650    double fzp_syst_min = FzPmc*(1.-FZP_SYST_MIN);
 651    double k_syst_min = KMC*(1.-K_SYST_MIN);
 652    //
 653    double ewk_syst_max = EWK_SYST_MAX; // because this is just fraction
 654    double i_syst_max = Imc*I_SYST_MAX;
 655    double fz_syst_max = Fzmc*FZ_SYST_MAX;
 656    double fzp_syst_max = FzPmc*FZP_SYST_MAX;
 657    double k_syst_max = KMC*K_SYST_MAX;
 658    //
 659    // negative points
 660    for (int i=0; i<POINTS; ++i) {
 661      double x_ewk = 1.-ewk_syst_min + (ewk_syst_min)*i/POINTS;
 662      double x_fz  = fz_syst_min + (Fzmc-fz_syst_min)*i/POINTS;
 663      double x_fzp = fzp_syst_min + (FzPmc-fzp_syst_min)*i/POINTS;
 664      double x_k   = k_syst_min + (KMC-k_syst_min)*i/POINTS;
 665      double x_i   = i_syst_min + (Imc-i_syst_min)*i/POINTS;
 666      //
 667      double y_ewk= CalcABCD(Imc, Fzmc, FzPmc, KMC, x_ewk, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 668      double y_fz = CalcABCD(Imc, x_fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 669      double y_fzp= CalcABCD(Imc, Fzmc, x_fzp, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 670      double y_k  = CalcABCD(Imc, Fzmc, FzPmc, x_k, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 671      double y_i  = CalcABCD(x_i, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 672      //
 673      g_ewk.SetPoint(i,(x_ewk-1.)*100., 100.*fabs(y_ewk-SMC)/SMC);
 674      g_fz.SetPoint(i,(x_fz-Fzmc)*100./Fzmc, 100.*fabs(y_fz-SMC)/SMC);
 675      g_fzp.SetPoint(i,(x_fzp-FzPmc)*100./FzPmc, 100.*fabs(y_fzp-SMC)/SMC);
 676      g_i.SetPoint(i,(x_i-Imc)*100./Imc, 100.*fabs(y_i-SMC)/SMC);
 677      g_k.SetPoint(i,(x_k-KMC)*100./KMC, 100.*fabs(y_k-SMC)/SMC);
 678    }
 679    //
 680    // positive points
 681    for (int i=0; i<=POINTS; ++i) {
 682      double x_ewk = 1.+ewk_syst_max*i/POINTS; 
 683      double x_fz  = Fzmc+fz_syst_max*i/POINTS;
 684      double x_fzp = FzPmc+fzp_syst_max*i/POINTS;
 685      double x_k   = KMC+k_syst_max*i/POINTS;
 686      double x_i   = Imc+i_syst_max*i/POINTS;
 687      //
 688      double y_ewk= CalcABCD(Imc, Fzmc, FzPmc, KMC, x_ewk, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 689      double y_fz = CalcABCD(Imc, x_fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 690      double y_fzp= CalcABCD(Imc, Fzmc, x_fzp, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 691      double y_k  = CalcABCD(Imc, Fzmc, FzPmc, x_k, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 692      double y_i  = CalcABCD(x_i, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
 693      //
 694      g_ewk.SetPoint(i+POINTS,(x_ewk-1.)*100., 100.*fabs(y_ewk-SMC)/SMC);
 695      g_fz.SetPoint(i+POINTS,(x_fz-Fzmc)*100./Fzmc, 100.*fabs(y_fz-SMC)/SMC);
 696      g_fzp.SetPoint(i+POINTS,(x_fzp-FzPmc)*100./FzPmc, 100.*fabs(y_fzp-SMC)/SMC);
 697      g_i.SetPoint(i+POINTS,(x_i-Imc)*100./Imc, 100.*fabs(y_i-SMC)/SMC);
 698      g_k.SetPoint(i+POINTS,(x_k-KMC)*100./KMC, 100.*fabs(y_k-SMC)/SMC);
 699    }
 700    TString yaxis("(S-S_{mc})/S_{mc} (%)");
 701    TCanvas c;
 702    g_ewk.SetLineWidth(2);
 703    g_ewk.GetXaxis()->SetTitle("EWK Variation (%)");
 704    g_ewk.GetYaxis()->SetTitle(yaxis);
 705    g_ewk.Draw("AL");
 706    c.Print("ewk_syst_variation.C");
 707    //
 708    g_fz.SetLineWidth(2);
 709    g_fz.GetXaxis()->SetTitle("F_{z} Variation (%)");
 710    g_fz.GetYaxis()->SetTitle(yaxis);
 711    g_fz.Draw("AL");
 712    c.Print("fz_syst_variation.C");
 713    //
 714    g_fzp.SetLineWidth(2);
 715    g_fzp.GetXaxis()->SetTitle("F_{z}' Variation (%)");
 716    g_fzp.GetYaxis()->SetTitle(yaxis);
 717    g_fzp.Draw("AL");
 718    c.Print("fzp_syst_variation.C");
 719    //
 720    g_i.SetLineWidth(2);
 721    g_i.GetXaxis()->SetTitle("I Variation (%)");
 722    g_i.GetYaxis()->SetTitle(yaxis);
 723    g_i.Draw("AL");
 724    c.Print("i_syst_variation.C");
 725    //
 726    g_k.SetLineWidth(2);
 727    g_k.GetXaxis()->SetTitle("K Variation (%)");
 728    g_k.GetYaxis()->SetTitle(yaxis);
 729    g_k.Draw("AL");
 730    c.Print("k_syst_variation.C");
 731    //
 732    // ******************************************************************
 733    //
 734    //
 735    // 
 736    double err_ewk = std::max(fabs(SMC-S_EWK_PLUS),fabs(SMC-S_EWK_MINUS));
 737    double err_fz = fabs(SMC-S_FZ);
 738    double err_fzp = fabs(SMC-S_FZP);
 739    double err_i  = fabs(SMC-S_I);
 740    double err_k = fabs(SMC-S_K);
 741    //
 742    double DS_syst = sqrt(err_ewk*err_ewk + err_fz*err_fz + err_fzp*err_fzp+
 743			  err_i*err_i + err_k*err_k);
 744    //
 745    cout << "********************************************************" << endl;
 746    cout << "Signal Prediction: " << S << "+-" << DS << "(stat) +-"
 747	 << DS_syst << "(syst)"  << endl;
 748    cout << "stat error: " << 100.*DS/S <<"%" << endl;
 749    cout << "syt  error: " << 100.*DS_syst/S<< "%"  << endl;
 750    cout << "********************************************************" << endl;
 751    cout << "Parameters used in calculation: " << endl;
 752    cout << "I=  " << I << "+-" << dI << endl;
 753    cout << "Fz= " << Fz << "+-" << dFz << endl;
 754    cout << "FzP=" << FzP << "+-" << dFzP << endl;
 755    cout << "EWK error assumed to be: " << ewkerror << endl;
 756    cout << endl;
 757    cout << "ABCD Regions population:" << endl;
 758    cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
 759         << ", ewk=" << a_ewk << endl;
 760    cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
 761         << ", ewk=" << b_ewk << endl;
 762    cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
 763         << ", ewk=" << c_ewk << endl;
 764    cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
 765         << ", ewk=" << d_ewk << endl;
 766    cout << endl;
 767    cout << "Parameters from MC: " << endl;
 768    cout << "I=  " << Imc << "+-" << dImc << endl;
 769    cout << "Fz= " << Fzmc << "+-" << dFzmc << endl;
 770    cout << "FzP=" << FzPmc << "+-" << dFzPmc << endl;
 771    cout << endl;
 772    cout << "Real value of K=" << KMC << endl;
 773    cout << "Real value of Signal=" << SMC << endl;
 774    cout << endl;
 775    cout << "Difference Measured - MC value (% wrt MC value except K=1): " 
 776	 << endl;
 777    cout << "Fz : " << dfz  << ", (" << dfz*100./Fzmc << "%)" << endl;
 778    cout << "FzP: " << dfzp << ", (" << dfzp*100./FzPmc << "%)"  << endl;
 779    cout << "I  : " << di   << ", (" << di*100./Imc << "%)"  << endl;
 780    cout << "K  : " << fk   << ", (" << fk*100./1. << "%)"  << endl;
 781    cout << endl;
 782    //
 783    cout << "DETAILS OF THE CALCULATION" << endl;
 784    cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl;
 785    cout << "Statistical Error Summary: " << endl;
 786    cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
 787    cout << "due to FzP= "<< SpFzP*dFzP
 788	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
 789    cout << "due to  I = "<< SpI*dI
 790	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
 791    cout << "due to Na = "<< SpNa*sqrt(Na)
 792	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
 793    cout << "due to Nb = "<< SpNb*sqrt(Nb)
 794	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
 795    cout << "due to Nc = "<< SpNc*sqrt(Nc)
 796	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
 797    cout << "due to Nd = "<< SpNd*sqrt(Nd)
 798	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
 799    cout << "Total Statistical Error: " 
 800	 << DS << ", (" << DS*100./S << "%)"<< endl;
 801    cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
 802    cout << endl;
 803    cout << "Systematic Error Summary:" << endl;
 804    cout << "due to k   = " << err_k << " ( " << err_k*100./S  << "%)" << endl;
 805    cout << "due to Fz  = " << err_fz << " ( " << err_fz*100./S  << "%)" << endl;
 806    cout << "due to FzP = " << err_fzp << " ( " << err_fzp*100./S  << "%)" << endl;
 807    cout << "due to I   = " << err_i << " ( " << err_i*100./S  << "%)" << endl;
 808    cout << "due to EWK = " << err_ewk << " ( " << err_ewk*100./S  << "%)" << endl;
 809
 810    cout << "Syst Error percentages are wrt S prediction, not S mc" << endl;
 811  }
 812  //
 813  //
 814  if (mcOnly < 0 && mcOnly >-0.75) {  // select value -0.5 that gives the
 815                                      // string parser
 816    cout << "=======================================================" << endl;
 817    cout << "Calculating ABCD Result and Stat Error Assuming MC ONLY"  << endl;
 818    cout << "=======================================================" << endl;
 819    cout << "All input parameters that the user have inserted will be "
 820	 << "ignored and recalculated from MC" << endl;
 821    cout << "This option will not give you systematics estimation" << endl;
 822    // recalculate the basic quantities
 823    I = (a_sig + b_sig) / (a_sig + b_sig + c_sig + d_sig);
 824    dI = sqrt(I*(1-I)/(a_sig + b_sig + c_sig + d_sig));
 825    Fz = a_sig/b_sig;
 826    double e =a_sig/(a_sig + b_sig);
 827    double de = sqrt(e*(1-e)/(a_sig + b_sig));
 828    double alpha = de/(2.*Fz-e);
 829    dFz = alpha/(1-alpha);
 830    FzP = d_sig/c_sig;
 831    double ep =d_sig/(c_sig + d_sig);
 832    double dep = sqrt(ep*(1-ep)/(c_sig + d_sig));
 833    double alphap = dep/(2.*FzP-ep);
 834    dFzP = alphap/(1-alphap);
 835    //
 836    double KMC = (d_qcd/c_qcd)/(a_qcd/b_qcd);
 837    //
 838    // now everything is done from data + input
 839    double A = (1.0-I)*(FzP-Fz);
 840    double B = I*(FzP+1.0)*(FzP*(c_sig+c_qcd)-(d_sig+d_qcd)) + 
 841      (1+Fz)*(1-I)*((a_sig+a_qcd)-dFzP*(b_sig+b_qcd));
 842    double C = I*(1.+Fz)*(1.+FzP)*((d_sig+d_qcd)*(b_sig+b_qcd) - 
 843				   (a_sig+a_qcd)*(c_sig+c_qcd));
 844    //
 845    // signal calculation:
 846    double S = Trionym(A,B,C, a_sig+b_sig);
 847
 848    // the errors now:
 849    // calculate the statistical error now:
 850    double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
 851    double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
 852    double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
 853    double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
 854    //
 855    double Na = a_sig+a_qcd+a_ewk, Nb = b_sig+b_qcd+b_ewk;
 856    double Nc=c_sig+c_qcd+c_ewk,   Nd = d_sig+d_qcd+d_ewk;
 857    double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
 858    if (A != 0) {
 859      
 860      ApI   = -(FzP-Fz);
 861      ApFz  = -(1.0-I);
 862      ApFzP = (1.0-I);
 863      ApNa  = 0.0;
 864      ApNb  = 0.0;
 865      ApNc  = 0.0;
 866      ApNd  = 0.0;
 867      
 868      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 869      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 870      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 871      BpNa  =  (1.0-I)*(1.0+Fz);
 872      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 873      BpNc  = I*(FzP+1.0)*Fz;
 874      BpNd  = -I*(FzP+1.0);
 875      
 876      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 877      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 878      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 879      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 880      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 881      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 882      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 883      
 884      SpI   = (-BpI   + (B*BpI   -2.0*ApI*C   -2.0*A*CpI)  /fabs(2.0*A*S+B)- 2.0*ApI*S)  /(2.0*A);
 885      SpFz  = (-BpFz  + (B*BpFz  -2.0*ApFz*C  -2.0*A*CpFz) /fabs(2.0*A*S+B)- 2.0*ApFz*S) /(2.0*A);
 886      SpFzP = (-BpFzP + (B*BpFzP -2.0*ApFzP*C -2.0*A*CpFzP)/fabs(2.0*A*S+B)- 2.0*ApFzP*S)/(2.0*A);
 887      SpNa  = (-BpNa  + (B*BpNa  -2.0*ApNa*C  -2.0*A*CpNa) /fabs(2.0*A*S+B)- 2.0*ApNa*S) /(2.0*A);
 888      SpNb  = (-BpNb  + (B*BpNb  -2.0*ApNb*C  -2.0*A*CpNb) /fabs(2.0*A*S+B)- 2.0*ApNb*S) /(2.0*A);
 889      SpNc  = (-BpNc  + (B*BpNc  -2.0*ApNc*C  -2.0*A*CpNc) /fabs(2.0*A*S+B)- 2.0*ApNc*S) /(2.0*A);
 890      SpNd  = (-BpNd  + (B*BpNd  -2.0*ApNd*C  -2.0*A*CpNd) /fabs(2.0*A*S+B)- 2.0*ApNd*S) /(2.0*A);
 891    }
 892    else {
 893      BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
 894      BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
 895      BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
 896      BpNa  =  (1.0-I)*(1.0+Fz);
 897      BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
 898      BpNc  = I*(FzP+1.0)*Fz;
 899      BpNd  = -I*(FzP+1.0);
 900      
 901      CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
 902      CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 903      CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
 904      CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
 905      CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
 906      CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
 907      CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
 908      
 909      SpI   = -CpI/B+C*BpI/(B*B);
 910      SpFz  = -CpFz/B+C*BpFz/(B*B);
 911      SpFzP = -CpFzP/B+C*BpFzP/(B*B);
 912      SpNa  = -CpNa/B+C*BpNa/(B*B);
 913      SpNb  = -CpNb/B+C*BpNb/(B*B);
 914      SpNc  = -CpNc/B+C*BpNc/(B*B);
 915      SpNd  = -CpNd/B+C*BpNd/(B*B);
 916    }
 917    double  DS;
 918    DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
 919	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
 920    // warning: S here denotes the method prediction ..........
 921    cout << "********************************************************" << endl;
 922    cout << "Signal Prediction: " << S << "+-" << DS << "(stat)" << endl;
 923    cout << "********************************************************" << endl;
 924    cout << "Parameters used in calculation: " << endl;
 925    cout << "I=  " << I << "+-" << dI << endl;
 926    cout << "Fz= " << Fz << "+-" << dFz << endl;
 927    cout << "FzP=" << FzP << "+-" << dFzP << endl;
 928    cout << endl;
 929    cout << "ABCD Regions population:" << endl;
 930    cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
 931	 << ", ewk=" << a_ewk << endl;
 932    cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
 933	 << ", ewk=" << b_ewk << endl;
 934    cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
 935	 << ", ewk=" << c_ewk << endl;
 936    cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
 937	 << ", ewk=" << d_ewk << endl;
 938    cout << "K value from MC: " << KMC << endl;
 939    cout << endl;
 940    cout << "Statistical Error Summary: " << endl;
 941    cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
 942    cout << "due to FzP= "<< SpFzP*dFzP
 943	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
 944    cout << "due to  I = "<< SpI*dI
 945	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
 946    cout << "due to Na = "<< SpNa*sqrt(Na)
 947	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
 948    cout << "due to Nb = "<< SpNb*sqrt(Nb)
 949	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
 950    cout << "due to Nc = "<< SpNc*sqrt(Nc)
 951	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
 952    cout << "due to Nd = "<< SpNd*sqrt(Nd)
 953	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
 954    cout << "Total Statistical Error: " 
 955	 << DS << ", (" << DS*100./S << "%)"<< endl;
 956    cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
 957  }
 958
 959
 960}
 961
 962void plotMaker(TString histoName, TString wzsignal,
 963	       vector<TString> file, vector<TString> type, 
 964	       vector<double> weight, TString xtitle)
 965{
 966  gROOT->Reset();
 967  gROOT->ProcessLine(".L tdrstyle.C"); 
 968  gROOT->ProcessLine("setTDRStyle()");
 969  // automatic recognition of histogram dimension
 970  int NBins = 0; double min = 0; double max = -1;
 971  for (int i=0; i< (int) file.size(); ++i) {
 972    if (weight[i]>0) {
 973      TFile f(file[i]);
 974      TH1F *h = (TH1F*) f.Get(histoName);
 975      NBins = h->GetNbinsX();
 976      min = h->GetBinLowEdge(1);
 977      max = h->GetBinLowEdge(NBins+1);
 978      break;
 979    }
 980  }
 981  if (NBins ==0 || (max<min)) {
 982    std::cout << "PlotCombiner::abcd error: Could not find valid histograms in file"
 983              << std::endl;
 984    abort();
 985  }
 986  cout << "Histograms with "<< NBins <<" bins  and range " << min << "-" << max  << endl;
 987  // Wenu Signal .......................................................
 988  TH1F h_wenu("h_wenu", "h_wenu", NBins, min, max);
 989  int fmax = (int) file.size();
 990  for (int i=0; i<fmax; ++i) {
 991    if (type[i] == "sig" && weight[i]>0) {
 992      TFile f(file[i]);
 993      TH1F *h = (TH1F*) f.Get(histoName);
 994      h_wenu.Add(h, weight[i]);
 995    }
 996  }
 997  // Bkgs ..............................................................
 998  //
 999  // QCD light flavor
1000  TH1F h_qcd("h_qcd", "h_qcd", NBins, min, max);
1001  for (int i=0; i<fmax; ++i) {
1002    if (type[i] == "qcd" && weight[i]>0) {
1003      TFile f(file[i]);
1004      TH1F *h = (TH1F*) f.Get(histoName);
1005      h_qcd.Add(h, weight[i]);
1006    }
1007  }
1008  // QCD heavy flavor
1009  TH1F h_bce("h_bce", "h_bce", NBins, min, max);
1010  for (int i=0; i<fmax; ++i) {
1011    if (type[i] == "bce" && weight[i]>0) {
1012      TFile f(file[i]);
1013      TH1F *h = (TH1F*) f.Get(histoName);
1014      h_bce.Add(h, weight[i]);
1015    }
1016  }
1017  // QCD Gjets
1018  TH1F h_gj("h_gj", "h_gj", NBins, min, max);
1019  for (int i=0; i<fmax; ++i) {
1020    if (type[i] == "gje" && weight[i]>0) {
1021      TFile f(file[i]);
1022      TH1F *h = (TH1F*) f.Get(histoName);
1023      h_gj.Add(h, weight[i]);
1024    }
1025  }
1026  // Other EWK bkgs
1027  TH1F h_ewk("h_ewk", "h_ewk", NBins, min, max);
1028  for (int i=0; i<fmax; ++i) {
1029    if (type[i] == "ewk" && weight[i]>0) {
1030      TFile f(file[i]);
1031      TH1F *h = (TH1F*) f.Get(histoName);
1032      h_ewk.Add(h, weight[i]);
1033    }
1034  }
1035  //
1036  // ok now decide how to plot them:
1037  // first the EWK bkgs
1038  h_ewk.SetFillColor(3);
1039  //
1040  // then the gjets
1041  h_gj.Add(&h_ewk);
1042  h_gj.SetFillColor(1);
1043  // thent the QCD dijets
1044  h_bce.Add(&h_qcd);
1045  h_bce.Add(&h_gj);
1046  h_bce.SetFillColor(2);
1047  // and the signal at last
1048  TH1F h_tot("h_tot", "h_tot", NBins, min, max);
1049  h_tot.Add(&h_bce);
1050  h_tot.Add(&h_wenu);
1051  h_wenu.SetLineColor(4);  h_wenu.SetLineWidth(2);
1052  //
1053  TCanvas c;
1054  h_tot.GetXaxis()->SetTitle(xtitle);
1055  h_tot.Draw("PE");
1056  h_bce.Draw("same");
1057  h_gj.Draw("same");
1058  h_ewk.Draw("same");
1059  h_wenu.Draw("same");
1060
1061  // the Legend
1062  TLegend  leg(0.6,0.65,0.95,0.92);
1063  if (wzsignal == "wenu")
1064    leg.AddEntry(&h_wenu, "W#rightarrow e#nu","l");
1065  else
1066    leg.AddEntry(&h_wenu, "Z#rightarrow ee","l");
1067  leg.AddEntry(&h_tot, "Signal + Bkg","p");
1068  leg.AddEntry(&h_bce, "dijets","f");
1069  leg.AddEntry(&h_gj, "#gamma + jets","f");
1070  leg.AddEntry(&h_ewk,  "EWK+t#bar t", "f");
1071  leg.Draw("same");
1072
1073  c.Print("test.png");
1074
1075
1076
1077}
1078
1079
1080//
1081// reads the ABCD string and returns its value
1082// value is whatever it exists after the = sign and
1083// before the comma or the parenthesis
1084//
1085// if the string is not contained returns -1
1086// if there is no value, but the string is contained -0.5
1087// if there is an error in the algorithm return -99 and print error
1088// else returns its value
1089double searchABCDstring(TString abcdString, TString keyword)
1090{
1091  int size = keyword.Sizeof();
1092  int existsEntry = abcdString.Index(keyword);
1093  //
1094  if (existsEntry==-1) return -1.;
1095  //
1096  TString previousVal = abcdString(existsEntry-1);
1097  if (!(previousVal == "," || previousVal == " " || 
1098	previousVal == "(" )) return -1.;
1099  //
1100  TString afterVal = abcdString(existsEntry+size-1);
1101  //std::cout << "afterVal=" << afterVal << std::endl;
1102  if (afterVal =="," || afterVal==")") return -0.5;
1103  else if (afterVal != "=") return -1.;
1104  //
1105  // now find the comma or the parenthesis after the = sign
1106  int comma = abcdString.Index(",",existsEntry);
1107  //std::cout << "first comma=" << comma << std::endl;
1108  if (comma<0) comma = abcdString.Index(")",existsEntry);
1109  if (comma<0) {
1110    std::cout << "Error in parcing abcd line, chech syntax " 
1111	      << std::endl;
1112    return -99.;
1113  }
1114  TString svalue=abcdString(existsEntry+size,comma-existsEntry-size);
1115  std::cout << "Found ABCD parameter "<< keyword
1116	    << " with value "  << svalue << endl;
1117  // convert this to a float
1118  double value = svalue.Atof();
1119  return value;
1120
1121}
1122
1123
1124double Trionym(double a, double b, double c, double sum)
1125{
1126  if (a==0) {
1127    return -c/b;
1128  }
1129  double D2 = b*b - 4.*a*c;
1130  //return (-b + sqrt(D2)) / (2.*a);
1131  if (D2 > 0) {
1132    double s1 = (-b + sqrt(D2)) / (2.*a);
1133    double s2 = (-b - sqrt(D2)) / (2.*a);
1134    double solution =   fabs(s1-sum)<fabs(s2-sum)?s1:s2;
1135    return solution;
1136  }
1137  else  {
1138    return -1.;  
1139  }
1140}
1141
1142//
1143// the naming of the variables and the order is in this way for historical
1144// reasons
1145// your complains for different Nd_ and nd to G.Daskalakis :P
1146//
1147double CalcABCD
1148(double I, double Fz, double FzP, double K, double ewk,
1149 double Na_, double Nb_, double Nc_, double Nd_ ,
1150 double Ea_, double Eb_, double Ec_, double Ed_)
1151{
1152  double A, B, C;
1153  A = (1.0-I)*(FzP-K*Fz);
1154  B = I*(FzP+1.0)*(K*Fz*(Nc_-ewk*Ec_)-(Nd_-ewk*Ed_))+
1155    (1.0-I)*(1.0+Fz)*(K*(Na_-ewk*Ea_)-FzP*(Nb_-ewk*Eb_));
1156  C = I*(1.0+Fz)*(1.0+FzP)*((Nd_-ewk*Ed_)*(Nb_-ewk*Eb_)-
1157			    K*(Na_-ewk*Ea_)*(Nc_-ewk*Ec_));
1158  //
1159  return Trionym(A, B, C, Na_+Nb_-Ea_-Eb_);
1160
1161}