CMS 3D CMS Logo
PlotCombiner.cc
Go to the documentation of this 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 
97 void plotMaker(TString histoName, TString typeOfplot,
98  vector<TString> file, vector<TString> type,
99  vector<double> weight, TString xtitle);
100 
101 void 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);
105 double searchABCDstring(TString abcdString, TString keyword);
106 double Trionym(double a, double b, double c, double sum);
107 double 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
113 const double EWK_SYST_MIN = 0.3;
114 const double EWK_SYST_MAX = 0.3;
115 //
116 const double I_SYST_MIN = 0.05;
117 const double I_SYST_MAX = 0.05;
118 //
119 const double FZ_SYST_MIN = 0.1;
120 const double FZ_SYST_MAX = 0.1;
121 //
122 const double FZP_SYST_MIN = 0.1;
123 const double FZP_SYST_MAX = 0.1;
124 //
125 const double K_SYST_MIN = 0.8;
126 const double K_SYST_MAX = 0.8;
127 
128 
129 using namespace std;
130 
131 void 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 
251 void 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);
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 
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 
593  // SYSTEMATICS CALCULATION /////////////////////////////////////////////
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);
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 
962 void 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
1089 double 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 
1124 double 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 //
1147 double 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 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
HGCalRecHit_cfi.weights
weights
Definition: HGCalRecHit_cfi.py:6
CaloTowersParam_cfi.mc
mc
Definition: CaloTowersParam_cfi.py:8
type
type
Definition: HCALResponse.h:21
funct::D2
Divides< B, C > D2
Definition: Factorize.h:147
alpha
float alpha
Definition: AMPTWrapper.h:95
Trionym
double Trionym(double a, double b, double c, double sum)
Definition: PlotCombiner.cc:1124
mps_fire.i
i
Definition: mps_fire.py:338
plotMaker
void plotMaker(TString histoName, TString typeOfplot, vector< TString > file, vector< TString > type, vector< double > weight, TString xtitle)
Definition: PlotCombiner.cc:962
PlotCombiner
void PlotCombiner()
Definition: PlotCombiner.cc:131
w
const double w
Definition: UKUtility.cc:23
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
I_SYST_MAX
const double I_SYST_MAX
Definition: PlotCombiner.cc:117
EWK_SYST_MIN
const double EWK_SYST_MIN
Definition: PlotCombiner.cc:113
EWK_SYST_MAX
const double EWK_SYST_MAX
Definition: PlotCombiner.cc:114
weight
Definition: weight.py:1
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
types
types
Definition: AlignPCLThresholds_PayloadInspector.cc:30
input
static std::string const input
Definition: EdmProvDump.cc:48
I_SYST_MIN
const double I_SYST_MIN
Definition: PlotCombiner.cc:116
patCaloMETCorrections_cff.C
C
Definition: patCaloMETCorrections_cff.py:45
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
FZ_SYST_MAX
const double FZ_SYST_MAX
Definition: PlotCombiner.cc:120
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
mps_splice.line
line
Definition: mps_splice.py:76
abcd
void abcd(vector< TString > file, vector< TString > type, vector< double > weight, double METCut, double I, double dI, double Fz, double dFz, double FzP, double dFzP, double ewkerror, double data, double mc, double mcOnly)
Definition: PlotCombiner.cc:251
Exhume::I
const std::complex< double > I
Definition: I.h:8
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
create_public_lumi_plots.leg
leg
Definition: create_public_lumi_plots.py:1166
value
Definition: value.py:1
min
T min(T a, T b)
Definition: MathUtil.h:58
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
CalcABCD
double CalcABCD(double I, double Fz, double FzP, double K, double ewk, double Na_, double Nb_, double Nc_, double Nd_, double Ea_, double Eb_, double Ec_, double Ed_)
Definition: PlotCombiner.cc:1148
FZP_SYST_MIN
const double FZP_SYST_MIN
Definition: PlotCombiner.cc:122
K_SYST_MIN
const double K_SYST_MIN
Definition: PlotCombiner.cc:125
b
double b
Definition: hdecay.h:120
FZP_SYST_MAX
const double FZP_SYST_MAX
Definition: PlotCombiner.cc:123
S
double S(const TLorentzVector &, const TLorentzVector &)
Definition: Particle.cc:99
alignmentValidation.fname
string fname
main script
Definition: alignmentValidation.py:959
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
K_SYST_MAX
const double K_SYST_MAX
Definition: PlotCombiner.cc:126
a
double a
Definition: hdecay.h:121
FZ_SYST_MIN
const double FZ_SYST_MIN
Definition: PlotCombiner.cc:119
patCaloMETCorrections_cff.A
A
Definition: patCaloMETCorrections_cff.py:43
gather_cfg.cout
cout
Definition: gather_cfg.py:144
searchABCDstring
double searchABCDstring(TString abcdString, TString keyword)
Definition: PlotCombiner.cc:1089
Generated for CMSSW Reference Manual by doxygen 1.8.11