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

  2.      Macro to make the plots .......................................
    3. 

    4. 

  4.      Instructions:
    5. 

  5.      a. set up an input file that looks like the following:
    6. 

  6.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    7. 

  7.      # zee or wenu
    8. 

  8.      wenu
    9. 

  9.      # file name, type (sig, qcd, bce, gje, ewk), weight
    10. 

  10.      histos_wenu.root     sig     1.46
    11. 

  11.      histos_q20_30.root   qcd     0
    12. 

  12.      histos_q30_80.root   qcd     100.
    13. 

  13.      histos_q80_170.root  qcd     0
    14. 

  14.      histos_b20_30.root   bce     0
    15. 

  15.      histos_b30_80.root   bce     0
    16. 

  16.      histos_b80_170.root  bce     0
    17. 

  17.      histos_zee.root      ewk     0
    18. 

  18.      histos_wtaunu.root   ewk     0
    19. 

  19.      histos_ztautau.root  ewk     0
    20. 

  20.      histos_gj15.root     gje     0
    21. 

  21.      histos_gj20.root     gje     0
    22. 

  22.      histos_gj25.root     gje     10.12
    23. 

  23.      histos_gj30.root     gje     0
    24. 

  24.      histos_gj35.root     gje     0
    25. 

  25.      histos_wmunu.root    ewk     0
    26. 

  26.      histos_ttbar.root    ewk     0
    27. 

  27.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    28. 

  28.      lines that start with # are considered to be comments
    29. 

  29.      line 2 has wenu or zee. From line 4 the list of the histo files are listed
    30. 

  30.      (first word) then a type that could be sig,qcd,bce, gj or ewk in order to
    31. 

  31.      discriminate among different sources of bkgs and finally the weight that we
    32. 

  32.      want to weight the histogram entries. This particular example is for Wenu. For
    33. 

  33.      Zee one has to put type sig in the zee file and ewk in the Wenu file. The order
    34. 

  34.      of the files is arbitrary. Files with weight 0 will be ignored.
    35. 

  35.      After you have set up this code you run a root macro to combine the plots.
    36. 

  36.      You can do (not recommended - it actually crushes - to be debugged)
    37. 

  37.      root -b PlotCombiner.cc 
    38. 

  38.      or to compile it within root (recommended)
    39. 

  39.      root -b
    40. 

  40.      root [1] .L PlotCombiner.cc++
    41. 

  41.      root [2] PlotCombiner()
    42. 

  42.      
    43. 

  43.      and you finally get the plots.
    44. 

    45. 

  45.      For the ABCD method:
    46. 

  46.      ^^^^^^^^^^^^^^^^^^^^
    47. 

  47.      you have to insert in the 2nd line instead of wenu or zee the keyword abcd(...)
    48. 

  48.      The files should contain ewk samples, sig samples and qcd samples (but also read
    49. 

  49.      later). The only absolutely necessary files are the sig ones.
    50. 

  50.      Example:
    51. 

  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. 

  52.      These parameters keep the same notation as in the note. The last parameter (data)
    53. 

  53.      can take 3 values:
    54. 

  54.      data: calculate in ABCD as in data. This means that the histograms denoted with
    55. 

  55.            sig,qcd,bce,gje  are used as of the same kind and ewk as the MC ewk. 
    56. 

  56.            The background is substructed as in data
    57. 

  57.      mcOnly: here we ignore all the input parameters I, dI etc. All parameters are taken
    58. 

  58.            from MC by forcing Fqcd=1
    59. 

  59.      mc:   <THIS IS WHAT ONE NORMALLY USES> input mc samples, calculation of statistical
    60. 

  60.            and systematics as in CMS AN 2009/004, systematic and statistic error 
    61. 

  61. 	   calculation. This option also creates the plots of the variation of the
    62. 

  62.            signal prediction vs the parameter variation. In order to set the limits of
    63. 

  63.            the desired variation you have to edit the values in line 113 of this code
    64. 

  64.            (they are hardwired in the code)
    65. 

  65.      TO DO:
    66. 

  66.      functionalities to plot more kind of plots, e.g. efficiencies
    67. 

  67.      
    68. 

  68.      
    69. 

  69.      Further Questions/Contact:
    70. 

  70.      
    71. 

  71.          nikolaos.rompotis @ cern.ch
    72. 

    73. 

    74. 

    75. 

  75. 	 Nikolaos Rompotis - 29 June 09
    76. 

  76. 	 18 Sept 09:  1st updgrade: input files in a text file
    77. 

  77. 	 28 May  10:  bug in IMET corrected, thanks to Sadia Khalil
    78. 

  78. 	 Imperial College London
    79. 

  79. 	 
    80. 

  80. 	 
    81. 

  81. */
    82. 

  82. 

  83. 

  84. #include <iostream>
    85. 

  85. #include <fstream>
    86. 

  86. #include <vector>
    87. 

  87. #include <string>
    88. 

  88. #include "TString.h"
    89. 

  89. #include "TROOT.h"
    90. 

  90. #include "TStyle.h"
    91. 

  91. #include "TH1F.h"
    92. 

  92. #include "TFile.h"
    93. 

  93. #include "TCanvas.h"
    94. 

  94. #include "TGraph.h"
    95. 

  95. #include "TLegend.h"
    96. 

  96. 

  97. void plotMaker(TString histoName, TString typeOfplot,
    98. 

  98. 	       vector<TString> file, vector<TString> type, 
    99. 

  99. 	       vector<double> weight,  TString xtitle);
    100. 

  100. 

  101. void abcd(vector<TString> file, vector<TString> type, vector<double> weight,
    102. 

  102. 	  double METCut, double I, double dI, double Fz, double dFz, 
    103. 

  103. 	  double FzP, double dFzP, double ewkerror,
    104. 

  104. 	  double data, double mc, double mcOnly);
    105. 

  105. double  searchABCDstring(TString abcdString, TString keyword);
    106. 

  106. double Trionym(double a, double b, double c, double sum);
    107. 

  107. double CalcABCD
    108. 

  108. (double I, double Fz, double FzP, double K, double ewk,
    109. 

  109.  double Na_, double Nb_, double Nc_, double Nd_,
    110. 

  110.  double Ea_, double Eb_, double Ec_, double Ed_);
    111. 

  111. 

  112. // values for systematics plots: it is fraction of the MC value
    113. 

  113. const double EWK_SYST_MIN = 0.3;
    114. 

  114. const double EWK_SYST_MAX = 0.3;
    115. 

  115. //
    116. 

  116. const double I_SYST_MIN = 0.05;
    117. 

  117. const double I_SYST_MAX = 0.05;
    118. 

  118. //
    119. 

  119. const double FZ_SYST_MIN = 0.1;
    120. 

  120. const double FZ_SYST_MAX = 0.1;
    121. 

  121. //
    122. 

  122. const double FZP_SYST_MIN = 0.1;
    123. 

  123. const double FZP_SYST_MAX = 0.1;
    124. 

  124. //
    125. 

  125. const double K_SYST_MIN = 0.8;
    126. 

  126. const double K_SYST_MAX = 0.8;
    127. 

  127. 

  128. 

  129. using namespace std;
    130. 

  130. 

  131. void PlotCombiner()
    132. 

  132. {
    133. 

  133.   // read the file
    134. 

  134.   ifstream input("inputFiles");
    135. 

  135.   int i = 0;
    136. 

  136.   TString typeOfplot = "";
    137. 

  137.   vector<TString> types;
    138. 

  138.   vector<TString> files;
    139. 

  139.   vector<double> weights;
    140. 

  140. 

  141.   if (input.is_open()) {
    142. 

  142.     std::string myline;
    143. 

  143.     while (! input.eof()) {
    144. 

  144.       getline(input, myline);
    145. 

  145.       TString line(myline);
    146. 

  146.       TString c('#');
    147. 

  147.       TString empty(' ');
    148. 

  148.       if (line[0] != c) {
    149. 

  149. 	++i;
    150. 

  150. 	if (i==1) typeOfplot=line;
    151. 

  151. 	else {
    152. 

  152. 	  // read until you find 3 words
    153. 

  153. 	  TString fname("");
    154. 

  154. 	  TString ftype("");
    155. 

  155. 	  TString fw("");
    156. 

  156. 	  int lineSize = (int) line.Length();
    157. 

  157. 	  int j=0;
    158. 

  158. 	  while (j<lineSize) {
    159. 

  159. 	    if(line[j] != empty) fname += line[j];
    160. 

  160. 	    else break;
    161. 

  161. 	    ++j;
    162. 

  162. 	  }
    163. 

  163. 	  while (j<lineSize) {
    164. 

  164. 	    if(line[j] != empty) ftype += line[j];
    165. 

  165. 	    else if(ftype.Length()==3) break;
    166. 

  166. 	    ++j;
    167. 

  167. 	  }
    168. 

  168. 	  while (j<lineSize) {
    169. 

  169. 	    if(line[j] != empty) fw += line[j];
    170. 

  170. 	    else{ if(fw.Length()>0) break;}
    171. 

  171. 	    ++j;
    172. 

  172. 	  }
    173. 

  173. 	  if (fname.Length() == 0) break;
    174. 

  174. 	  files.push_back(fname);
    175. 

  175. 	  types.push_back(ftype);
    176. 

  176. 	  double w = fw.Atof();
    177. 

  177. 	  weights.push_back(w);
    178. 

  178. 	  if (w>0)
    179. 

  179. 	    std::cout << fname << ", " << ftype << ", "<< w << std::endl;
    180. 

  180. 	}
    181. 

  181.       }
    182. 

  182.     }
    183. 

  183.     input.close();
    184. 

  184.   }
    185. 

  185.   else {
    186. 

  186.     std::cout << "File with name inputFile was not found" << std::endl;
    187. 

  187.     return;
    188. 

  188.   }
    189. 

  189. 

  190.   // now you can launch the jobs
    191. 

  191.   if (typeOfplot == "wenu") {
    192. 

  192.     cout << "wenu plot maker" << endl;
    193. 

  193.     //        ====================
    194. 

  194.     // =====> WHICH HISTOS TO PLOT
    195. 

  195.     //        ====================
    196. 

  196.     plotMaker("h_met", typeOfplot, files, types, weights, "MET (GeV)");
    197. 

  197.   }
    198. 

  198.   else if (typeOfplot == "zee"){
    199. 

  199.     cout << "zee plot maker" << endl;
    200. 

  200.     //        ====================
    201. 

  201.     // =====> WHICH HISTOS TO PLOT
    202. 

  202.     //        ====================
    203. 

  203.     plotMaker("h_mee", typeOfplot, files, types, weights, "M_{ee} (GeV)");
    204. 

  204.   }
    205. 

  205.   else if (typeOfplot(0,4) == "abcd") {
    206. 

  206.     // now read the parameters of the ABCD method
    207. 

  207.     // look for parameter I and dI
    208. 

  208.     double I = searchABCDstring(typeOfplot, "I");
    209. 

  209.     double dI= searchABCDstring(typeOfplot, "dI");
    210. 

  210.     // look for parameter Fz
    211. 

  211.     double Fz = searchABCDstring(typeOfplot, "Fz");
    212. 

  212.     double dFz= searchABCDstring(typeOfplot, "dFz");
    213. 

  213.     // look for parameter FzP
    214. 

  214.     double FzP = searchABCDstring(typeOfplot, "FzP");
    215. 

  215.     double dFzP= searchABCDstring(typeOfplot, "dFzP");
    216. 

  216.     // look for the MET cut
    217. 

  217.     double METCut =searchABCDstring(typeOfplot, "METCut");
    218. 

  218.     // do you want data driven only?
    219. 

  219.     double data = searchABCDstring(typeOfplot, "data");
    220. 

  220.     double mc = searchABCDstring(typeOfplot, "mc");
    221. 

  221.     double mcOnly = searchABCDstring(typeOfplot, "mcOnly");
    222. 

  222.     // what is the ewk error?
    223. 

  223.     double ewkerror = searchABCDstring(typeOfplot, "ewkerror");
    224. 

  224.     // sanity check:
    225. 

  225.     if (METCut<0 || (data<-0.7 && mc<-0.7 && mcOnly<-0.7)) {
    226. 

  226.       cout << "Error in your configurtion!" << endl;
    227. 

  227.       if (METCut <0) cout << "Error in MET Cut" << endl;
    228. 

  228.       else cout << "You need to specify one mc or data or mcOnly"
    229. 

  229. 		<< endl;
    230. 

  230.       abort();
    231. 

  231.     }
    232. 

  232.     if (mc>-0.7 && mc <0 && ewkerror<0) {
    233. 

  233.       cout << "You have specified mc option, but you have forgotten"
    234. 

  234. 	   << " to set the ewkerror!" << endl;
    235. 

  235.       abort();
    236. 

  236.     }
    237. 

  237.     //        ===============================
    238. 

  238.     // =====> ABCD METHOD FOR BKG SUBTRACTION
    239. 

  239.     //        ===============================
    240. 

  240.     cout << "doing ABCD with input: " << typeOfplot << endl;
    241. 

  241.     abcd(files, types, weights, METCut, I, dI, Fz, dFz, FzP, dFzP,
    242. 

  242. 	 ewkerror, data, mc, mcOnly);
    243. 

  243. 

  244.   }
    245. 

  245.   // force the program to abort in order to clear the memory
    246. 

  246.   // and avoid further use of the interpreter after
    247. 

  247.   abort();
    248. 

  248. 

  249. }
    250. 

  250. 

  251. void abcd( vector<TString> file, vector<TString> type, vector<double> weight, 
    252. 

  252. 	   double METCut, double I, double dI, double Fz, double dFz, 
    253. 

  253. 	   double FzP, double dFzP, double ewkerror,
    254. 

  254. 	   double data, double mc, double mcOnly)
    255. 

  255. {
    256. 

  256.   gROOT->Reset();
    257. 

  257.   gROOT->ProcessLine(".L tdrstyle.C"); 
    258. 

  258.   gROOT->ProcessLine("setTDRStyle()");
    259. 

  259.   //
    260. 

  260.   std::cout << "Trying ABCD method for Background subtration" << std::endl;
    261. 

  261.   //
    262. 

  262.   // histogram names to use:
    263. 

  263.   TString histoName_Ba("h_met_EB");
    264. 

  264.   TString histoName_Bb("h_met_inverse_EB");
    265. 

  265.   TString histoName_Ea("h_met_EE");
    266. 

  266.   TString histoName_Eb("h_met_inverse_EE");
    267. 

  267.   //
    268. 

  268.   // find one file and get the dimensions of your histogram
    269. 

  269.   int fmax = (int) file.size();
    270. 

  270.   int NBins = 0; double min = 0; double max = -1;
    271. 

  271.   for (int i=0; i<fmax; ++i) {
    272. 

  272.     if (weight[i]>0) {
    273. 

  273.       //      cout << "Loading file " << file[i] << endl;
    274. 

  274.       TFile f(file[i]);
    275. 

  275.       TH1F *h = (TH1F*) f.Get(histoName_Ba);
    276. 

  276.       NBins = h->GetNbinsX();
    277. 

  277.       min = h->GetBinLowEdge(1);
    278. 

  278.       max = h->GetBinLowEdge(NBins+1);
    279. 

  279.       break;
    280. 

  280.     }
    281. 

  281.   }
    282. 

  282.   if (NBins ==0 || (max<min)) {
    283. 

  283.     std::cout << "PlotCombiner::abcd error: Could not find valid histograms in file"
    284. 

  284. 	      << std::endl;
    285. 

  285.     abort();
    286. 

  286.   }
    287. 

  287.   cout << "Histograms with "<< NBins <<" bins  and range " << min << "-" << max  << endl;
    288. 

  288.   //
    289. 

  289.   // Wenu Signal .......................................................
    290. 

  290.   TH1F h_wenu("h_wenu", "h_wenu", NBins, min, max);
    291. 

  291.   TH1F h_wenu_inv("h_wenu_inv", "h_wenu_inv", NBins, min, max);
    292. 

  292.   for (int i=0; i<fmax; ++i) {
    293. 

  293.     if (type[i] == "sig" && weight[i]>0 ) {
    294. 

  294.       TFile f(file[i]);
    295. 

  295.       //
    296. 

  296.       TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
    297. 

  297.       h_wenu.Add(h_ba, weight[i]);
    298. 

  298.       TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
    299. 

  299.       h_wenu.Add(h_ea, weight[i]);
    300. 

  300.       //
    301. 

  301.       TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
    302. 

  302.       h_wenu_inv.Add(h_bb, weight[i]);
    303. 

  303.       TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
    304. 

  304.       h_wenu_inv.Add(h_eb, weight[i]);
    305. 

  305.     }
    306. 

  306.   }
    307. 

  307.   // QCD Bkgs
    308. 

  308.   TH1F h_qcd("h_qcd", "h_qcd", NBins, min, max);
    309. 

  309.   TH1F h_qcd_inv("h_qcd_inv", "h_qcd_inv", NBins, min, max);
    310. 

  310.   for (int i=0; i<fmax; ++i) {
    311. 

  311.     if ((type[i] == "qcd" || type[i] == "bce" 
    312. 

  312. 	 || type[i] == "gje") && weight[i]>0) {
    313. 

  313.       TFile f(file[i]);
    314. 

  314.       //
    315. 

  315.       TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
    316. 

  316.       h_qcd.Add(h_ba, weight[i]);
    317. 

  317.       TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
    318. 

  318.       h_qcd.Add(h_ea, weight[i]);
    319. 

  319.       //
    320. 

  320.       TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
    321. 

  321.       h_qcd_inv.Add(h_bb, weight[i]);
    322. 

  322.       TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
    323. 

  323.       h_qcd_inv.Add(h_eb, weight[i]);
    324. 

  324.     }
    325. 

  325.   }
    326. 

  326.   //
    327. 

  327.   TH1F h_ewk("h_ewk", "h_ewk", NBins, min, max);
    328. 

  328.   TH1F h_ewk_inv("h_ewk_inv", "h_ewk_inv", NBins, min, max);
    329. 

  329.   for (int i=0; i<fmax; ++i) {
    330. 

  330.     if ( type[i] == "ewk" && weight[i]>0) {
    331. 

  331.       TFile f(file[i]);
    332. 

  332.       //
    333. 

  333.       TH1F *h_ba = (TH1F*) f.Get(histoName_Ba);
    334. 

  334.       h_ewk.Add(h_ba, weight[i]);
    335. 

  335.       TH1F *h_ea = (TH1F*) f.Get(histoName_Ea);
    336. 

  336.       h_ewk.Add(h_ea, weight[i]);
    337. 

  337.       //
    338. 

  338.       TH1F *h_bb = (TH1F*) f.Get(histoName_Bb);
    339. 

  339.       h_ewk_inv.Add(h_bb, weight[i]);
    340. 

  340.       TH1F *h_eb = (TH1F*) f.Get(histoName_Eb);
    341. 

  341.       h_ewk_inv.Add(h_eb, weight[i]);
    342. 

  342.     }
    343. 

  343.   }
    344. 

  344.   //
    345. 

  345.   // calculate the METCut position
    346. 

  346.   //
    347. 

  347.   // this is calculated as a low edge bin of your input histogram
    348. 

  348.   // METCut = min + (max-min)*IMET/NBins
    349. 

  349.   int IMET = int ((METCut - min)/(max-min) * double(NBins)); 
    350. 

  350.   // check whether it is indeed a low egde position
    351. 

  351.   double metCalc = min + (max-min)*double(IMET)/double(NBins);
    352. 

  352.   if (metCalc < METCut || metCalc > METCut) {
    353. 

  353.     std::cout << "PlotCombiner:abcd: your MET Cut is not in low egde bin position"
    354. 

  354. 	      << std::endl;
    355. 

  355.   }
    356. 

  356.   cout << "MET Cut in " << METCut << "GeV corresponds to bin #" << IMET << endl;
    357. 

  357.   // Calculate the population in the ABCD Regions now
    358. 

  358.   // signal
    359. 

  359.   double a_sig = h_wenu.Integral(IMET,NBins+1);
    360. 

  360.   double b_sig = h_wenu.Integral(0,IMET-1);
    361. 

  361.   double c_sig = h_wenu_inv.Integral(0,IMET-1);
    362. 

  362.   double d_sig = h_wenu_inv.Integral(IMET,NBins+1);
    363. 

  363.   // qcd
    364. 

  364.   double a_qcd = h_qcd.Integral(IMET,NBins+1);
    365. 

  365.   double b_qcd = h_qcd.Integral(0,IMET-1);
    366. 

  366.   double c_qcd = h_qcd_inv.Integral(0,IMET-1);
    367. 

  367.   double d_qcd = h_qcd_inv.Integral(IMET,NBins+1);
    368. 

  368.   // ewk
    369. 

  369.   double a_ewk = h_ewk.Integral(IMET,NBins+1);
    370. 

  370.   double b_ewk = h_ewk.Integral(0,IMET-1);
    371. 

  371.   double c_ewk = h_ewk_inv.Integral(0,IMET-1);
    372. 

  372.   double d_ewk = h_ewk_inv.Integral(IMET,NBins+1);
    373. 

  373.   ////////////////////////////////////////////////
    374. 

  374. 

  375.   //
    376. 

  376.   // now the parameters of the method
    377. 

  377.   if (data < 0 && data >-0.75) {  // select value -0.5 that gives the
    378. 

  378.                                   // string parser
    379. 

  379.     // now everything is done from data + input
    380. 

  380.     std::cout << "Calculating ABCD Result and Stat Error Assuming DATA"
    381. 

  381. 	      << std::endl << "Summary: in this implementation we have assumed"
    382. 

  382. 	      << " that what real 'data' appear with type sig in the input"
    383. 

  383. 	      << std::endl << "No systematics available with this type of"
    384. 

  384. 	      << " calculation. If you need systematics try one of the other"
    385. 

  385. 	      << " options" << std::endl;
    386. 

  386.     double A = (1.0-I)*(FzP-Fz);
    387. 

  387.     double B = I*(FzP+1.0)*(FzP*(c_sig-c_ewk)-(d_sig-d_ewk)) + 
    388. 

  388.       (1+Fz)*(1-I)*((a_sig-a_ewk)-dFzP*(b_sig-b_ewk));
    389. 

  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. 

  390.     //
    391. 

  391.     // signal calculation:
    392. 

  392.     double S = Trionym(A,B,C, a_sig+b_sig);
    393. 

  393. 

  394.     // the errors now:
    395. 

  395.     // calculate the statistical error now:
    396. 

  396.     double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
    397. 

  397.     double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
    398. 

  398.     double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
    399. 

  399.     double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
    400. 

  400.     //
    401. 

  401.     double Na = a_sig, Nb = b_sig, Nc=c_sig, Nd = d_sig;
    402. 

  402.     double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
    403. 

  403.     if (A != 0) {
    404. 

  404.       
    405. 

  405.       ApI   = -(FzP-Fz);
    406. 

  406.       ApFz  = -(1.0-I);
    407. 

  407.       ApFzP = (1.0-I);
    408. 

  408.       ApNa  = 0.0;
    409. 

  409.       ApNb  = 0.0;
    410. 

  410.       ApNc  = 0.0;
    411. 

  411.       ApNd  = 0.0;
    412. 

  412.       
    413. 

  413.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    414. 

  414.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    415. 

  415.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    416. 

  416.       BpNa  =  (1.0-I)*(1.0+Fz);
    417. 

  417.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    418. 

  418.       BpNc  = I*(FzP+1.0)*Fz;
    419. 

  419.       BpNd  = -I*(FzP+1.0);
    420. 

  420.       
    421. 

  421.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    422. 

  422.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    423. 

  423.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    424. 

  424.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    425. 

  425.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    426. 

  426.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    427. 

  427.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    428. 

  428.       
    429. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  436.     }
    437. 

  437.     else {
    438. 

  438.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    439. 

  439.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    440. 

  440.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    441. 

  441.       BpNa  =  (1.0-I)*(1.0+Fz);
    442. 

  442.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    443. 

  443.       BpNc  = I*(FzP+1.0)*Fz;
    444. 

  444.       BpNd  = -I*(FzP+1.0);
    445. 

  445.       
    446. 

  446.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    447. 

  447.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    448. 

  448.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    449. 

  449.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    450. 

  450.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    451. 

  451.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    452. 

  452.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    453. 

  453.       
    454. 

  454.       SpI   = -CpI/B+C*BpI/(B*B);
    455. 

  455.       SpFz  = -CpFz/B+C*BpFz/(B*B);
    456. 

  456.       SpFzP = -CpFzP/B+C*BpFzP/(B*B);
    457. 

  457.       SpNa  = -CpNa/B+C*BpNa/(B*B);
    458. 

  458.       SpNb  = -CpNb/B+C*BpNb/(B*B);
    459. 

  459.       SpNc  = -CpNc/B+C*BpNc/(B*B);
    460. 

  460.       SpNd  = -CpNd/B+C*BpNd/(B*B);
    461. 

  461.     }
    462. 

  462.     double  DS;
    463. 

  463.     DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
    464. 

  464. 	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
    465. 

  465.     // warning: S here denotes the method prediction ..........
    466. 

  466.     cout << "********************************************************" << endl;
    467. 

  467.     cout << "Signal Prediction: " << S << "+-" << DS << "(stat)" << endl;
    468. 

  468.     cout << "********************************************************" << endl;
    469. 

  469.     cout << "Parameters used in calculation: " << endl;
    470. 

  470.     cout << "I=  " << I << "+-" << dI << endl;
    471. 

  471.     cout << "Fz= " << Fz << "+-" << dFz << endl;
    472. 

  472.     cout << "FzP=" << FzP << "+-" << dFzP << endl;
    473. 

  473.     cout << endl;
    474. 

  474.     cout << "ABCD Regions population:" << endl;
    475. 

  475.     cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
    476. 

  476.          << ", ewk=" << a_ewk << endl;
    477. 

  477.     cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
    478. 

  478.          << ", ewk=" << b_ewk << endl;
    479. 

  479.     cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
    480. 

  480.          << ", ewk=" << c_ewk << endl;
    481. 

  481.     cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
    482. 

  482.          << ", ewk=" << d_ewk << endl;
    483. 

  483.     cout << endl;
    484. 

  484.     //
    485. 

  485.     cout << "Statistical Error Summary: " << endl;
    486. 

  486.     cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
    487. 

  487.     cout << "due to FzP= "<< SpFzP*dFzP
    488. 

  488. 	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
    489. 

  489.     cout << "due to  I = "<< SpI*dI
    490. 

  490. 	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
    491. 

  491.     cout << "due to Na = "<< SpNa*sqrt(Na)
    492. 

  492. 	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
    493. 

  493.     cout << "due to Nb = "<< SpNb*sqrt(Nb)
    494. 

  494. 	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
    495. 

  495.     cout << "due to Nc = "<< SpNc*sqrt(Nc)
    496. 

  496. 	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
    497. 

  497.     cout << "due to Nd = "<< SpNd*sqrt(Nd)
    498. 

  498. 	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
    499. 

  499.     cout << "Total Statistical Error: " 
    500. 

  500. 	 << DS << ", (" << DS*100./S << "%)"<< endl;
    501. 

  501.     cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
    502. 

  502.   }
    503. 

  503.   //
    504. 

  504.   //
    505. 

  505.   //  this is the main option of the algorithm: the one implemented in the 
    506. 

  506.   //  Analysis Note
    507. 

  507.   //
    508. 

  508.   if (mc < 0 && mc >-0.75) {  // select value -0.5 that gives the 
    509. 

  509.                               // string parser
    510. 

  510.     
    511. 

  511.     //////// STATISTICAL ERROR CALCULATION /////////////////////////////
    512. 

  512.     double A = (1.0-I)*(FzP-Fz);
    513. 

  513.     double B = I*(FzP+1.0)*(FzP*(c_sig+c_qcd)-(d_sig+d_qcd)) + 
    514. 

  514.       (1+Fz)*(1-I)*((a_sig+a_qcd)-dFzP*(b_sig+b_qcd));
    515. 

  515.     double C = I*(1.+Fz)*(1.+FzP)*((d_sig+d_qcd)*(b_sig+b_qcd) - 
    516. 

  516. 				   (a_sig+a_qcd)*(c_sig+c_qcd));
    517. 

  517.     //
    518. 

  518.     // signal calculation:
    519. 

  519.     double S = Trionym(A,B,C, a_sig+b_sig);
    520. 

  520.     //
    521. 

  521.     double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
    522. 

  522.     double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
    523. 

  523.     double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
    524. 

  524.     double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
    525. 

  525.     //
    526. 

  526.     double Na = a_sig+a_qcd+a_ewk, Nb = b_sig+b_qcd+b_ewk;
    527. 

  527.     double Nc=c_sig+c_qcd+c_ewk,   Nd = d_sig+d_qcd+d_ewk;
    528. 

  528.     double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
    529. 

  529.     if (A != 0) {
    530. 

  530.       
    531. 

  531.       ApI   = -(FzP-Fz);
    532. 

  532.       ApFz  = -(1.0-I);
    533. 

  533.       ApFzP = (1.0-I);
    534. 

  534.       ApNa  = 0.0;
    535. 

  535.       ApNb  = 0.0;
    536. 

  536.       ApNc  = 0.0;
    537. 

  537.       ApNd  = 0.0;
    538. 

  538.       
    539. 

  539.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    540. 

  540.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    541. 

  541.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    542. 

  542.       BpNa  =  (1.0-I)*(1.0+Fz);
    543. 

  543.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    544. 

  544.       BpNc  = I*(FzP+1.0)*Fz;
    545. 

  545.       BpNd  = -I*(FzP+1.0);
    546. 

  546.       
    547. 

  547.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    548. 

  548.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    549. 

  549.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    550. 

  550.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    551. 

  551.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    552. 

  552.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    553. 

  553.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    554. 

  554.       
    555. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  562.     }
    563. 

  563.     else {
    564. 

  564.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    565. 

  565.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    566. 

  566.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    567. 

  567.       BpNa  =  (1.0-I)*(1.0+Fz);
    568. 

  568.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    569. 

  569.       BpNc  = I*(FzP+1.0)*Fz;
    570. 

  570.       BpNd  = -I*(FzP+1.0);
    571. 

  571.       
    572. 

  572.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    573. 

  573.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    574. 

  574.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    575. 

  575.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    576. 

  576.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    577. 

  577.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    578. 

  578.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    579. 

  579.       
    580. 

  580.       SpI   = -CpI/B+C*BpI/(B*B);
    581. 

  581.       SpFz  = -CpFz/B+C*BpFz/(B*B);
    582. 

  582.       SpFzP = -CpFzP/B+C*BpFzP/(B*B);
    583. 

  583.       SpNa  = -CpNa/B+C*BpNa/(B*B);
    584. 

  584.       SpNb  = -CpNb/B+C*BpNb/(B*B);
    585. 

  585.       SpNc  = -CpNc/B+C*BpNc/(B*B);
    586. 

  586.       SpNd  = -CpNd/B+C*BpNd/(B*B);
    587. 

  587.     }
    588. 

  588.     double  DS;
    589. 

  589.     DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
    590. 

  590. 	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
    591. 

  591. 

  592.     ////////////////////////////////////////////////////////////////////////
    593. 

  593.     // SYSTEMATICS CALCULATION /////////////////////////////////////////////
    594. 

  594.     ////////////////////////////////////////////////////////////////////////
    595. 

  595.     // recalculate the basic quantities
    596. 

  596.     double Imc = (a_sig + b_sig) / (a_sig + b_sig + c_sig + d_sig);
    597. 

  597.     double dImc = sqrt(Imc*(1-Imc)/(a_sig + b_sig + c_sig + d_sig));
    598. 

  598.     double Fzmc = a_sig/b_sig;
    599. 

  599.     double e =a_sig/(a_sig + b_sig);
    600. 

  600.     double de = sqrt(e*(1-e)/(a_sig + b_sig));
    601. 

  601.     double alpha = de/(2.*Fzmc-e);
    602. 

  602.     double dFzmc = alpha/(1-alpha);
    603. 

  603.     double FzPmc = d_sig/c_sig;
    604. 

  604.     double ep =d_sig/(c_sig + d_sig);
    605. 

  605.     double dep = sqrt(ep*(1-ep)/(c_sig + d_sig));
    606. 

  606.     double alphap = dep/(2.*FzPmc-ep);
    607. 

  607.     double dFzPmc = alphap/(1-alphap);
    608. 

  608.     //
    609. 

  609.     // calculate the K parameter as it is in MC:
    610. 

  610.     double KMC = (d_qcd/c_qcd)/(a_qcd/b_qcd);
    611. 

  611.     double SMC = a_sig + b_sig;
    612. 

  612.     //
    613. 

  613.     double dfz = Fz -Fzmc;
    614. 

  614.     double di = I - Imc;
    615. 

  615.     double dfzp = FzP - FzPmc;
    616. 

  616.     double fk = fabs(1-KMC);
    617. 

  617.     ////////////////////////////////////////////////////////////////////////
    618. 

  618.     // ewk error: this error has to be inserted by hand
    619. 

  619.     double fm = 1.-ewkerror;
    620. 

  620.     double fp = 1.+ewkerror;
    621. 

  621.     double S_EWK_PLUS = CalcABCD(Imc, Fzmc, FzPmc, KMC, fp, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    622. 

  622.     double S_EWK_MINUS = CalcABCD(Imc, Fzmc, FzPmc, KMC, fm, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    623. 

  623.     // error in K
    624. 

  624.     double S_K= CalcABCD(Imc, Fzmc, FzPmc, 1., 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    625. 

  625.     // error in Fz
    626. 

  626.     double S_FZ= CalcABCD(Imc, Fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    627. 

  627.     // error in FzP
    628. 

  628.     double S_FZP= CalcABCD(Imc, Fzmc, FzP, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    629. 

  629.     // error in I
    630. 

  630.     double S_I = CalcABCD(I, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    631. 

  631.     //
    632. 

  632.     // sanity tets
    633. 

  633.     //cout << "Smc=" << SMC<< ", " << CalcABCD(Imc, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed)
    634. 

  634.     // << endl;
    635. 

  635.     //abort();
    636. 

  636.     //
    637. 

  637.     // ************ plots for the systematics calculation ****************
    638. 

  638.     // ewk plot
    639. 

  639.     int const POINTS = 10;
    640. 

  640.     int const allPOINTS = 2*POINTS;
    641. 

  641.     TGraph g_ewk(allPOINTS);
    642. 

  642.     TGraph g_fz(allPOINTS);
    643. 

  643.     TGraph g_fzp(allPOINTS);
    644. 

  644.     TGraph g_k(allPOINTS);
    645. 

  645.     TGraph g_i(allPOINTS);
    646. 

  646.     //
    647. 

  647.     double ewk_syst_min = EWK_SYST_MIN; // because this is just fraction
    648. 

  648.     double i_syst_min = Imc*(1.-I_SYST_MIN);
    649. 

  649.     double fz_syst_min = Fzmc*(1.-FZ_SYST_MIN);
    650. 

  650.     double fzp_syst_min = FzPmc*(1.-FZP_SYST_MIN);
    651. 

  651.     double k_syst_min = KMC*(1.-K_SYST_MIN);
    652. 

  652.     //
    653. 

  653.     double ewk_syst_max = EWK_SYST_MAX; // because this is just fraction
    654. 

  654.     double i_syst_max = Imc*I_SYST_MAX;
    655. 

  655.     double fz_syst_max = Fzmc*FZ_SYST_MAX;
    656. 

  656.     double fzp_syst_max = FzPmc*FZP_SYST_MAX;
    657. 

  657.     double k_syst_max = KMC*K_SYST_MAX;
    658. 

  658.     //
    659. 

  659.     // negative points
    660. 

  660.     for (int i=0; i<POINTS; ++i) {
    661. 

  661.       double x_ewk = 1.-ewk_syst_min + (ewk_syst_min)*i/POINTS;
    662. 

  662.       double x_fz  = fz_syst_min + (Fzmc-fz_syst_min)*i/POINTS;
    663. 

  663.       double x_fzp = fzp_syst_min + (FzPmc-fzp_syst_min)*i/POINTS;
    664. 

  664.       double x_k   = k_syst_min + (KMC-k_syst_min)*i/POINTS;
    665. 

  665.       double x_i   = i_syst_min + (Imc-i_syst_min)*i/POINTS;
    666. 

  666.       //
    667. 

  667.       double y_ewk= CalcABCD(Imc, Fzmc, FzPmc, KMC, x_ewk, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    668. 

  668.       double y_fz = CalcABCD(Imc, x_fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    669. 

  669.       double y_fzp= CalcABCD(Imc, Fzmc, x_fzp, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    670. 

  670.       double y_k  = CalcABCD(Imc, Fzmc, FzPmc, x_k, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    671. 

  671.       double y_i  = CalcABCD(x_i, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    672. 

  672.       //
    673. 

  673.       g_ewk.SetPoint(i,(x_ewk-1.)*100., 100.*fabs(y_ewk-SMC)/SMC);
    674. 

  674.       g_fz.SetPoint(i,(x_fz-Fzmc)*100./Fzmc, 100.*fabs(y_fz-SMC)/SMC);
    675. 

  675.       g_fzp.SetPoint(i,(x_fzp-FzPmc)*100./FzPmc, 100.*fabs(y_fzp-SMC)/SMC);
    676. 

  676.       g_i.SetPoint(i,(x_i-Imc)*100./Imc, 100.*fabs(y_i-SMC)/SMC);
    677. 

  677.       g_k.SetPoint(i,(x_k-KMC)*100./KMC, 100.*fabs(y_k-SMC)/SMC);
    678. 

  678.     }
    679. 

  679.     //
    680. 

  680.     // positive points
    681. 

  681.     for (int i=0; i<=POINTS; ++i) {
    682. 

  682.       double x_ewk = 1.+ewk_syst_max*i/POINTS; 
    683. 

  683.       double x_fz  = Fzmc+fz_syst_max*i/POINTS;
    684. 

  684.       double x_fzp = FzPmc+fzp_syst_max*i/POINTS;
    685. 

  685.       double x_k   = KMC+k_syst_max*i/POINTS;
    686. 

  686.       double x_i   = Imc+i_syst_max*i/POINTS;
    687. 

  687.       //
    688. 

  688.       double y_ewk= CalcABCD(Imc, Fzmc, FzPmc, KMC, x_ewk, Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    689. 

  689.       double y_fz = CalcABCD(Imc, x_fz, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    690. 

  690.       double y_fzp= CalcABCD(Imc, Fzmc, x_fzp, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    691. 

  691.       double y_k  = CalcABCD(Imc, Fzmc, FzPmc, x_k, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    692. 

  692.       double y_i  = CalcABCD(x_i, Fzmc, FzPmc, KMC, 1., Na, Nb, Nc, Nd, Ea,Eb,Ec,Ed);
    693. 

  693.       //
    694. 

  694.       g_ewk.SetPoint(i+POINTS,(x_ewk-1.)*100., 100.*fabs(y_ewk-SMC)/SMC);
    695. 

  695.       g_fz.SetPoint(i+POINTS,(x_fz-Fzmc)*100./Fzmc, 100.*fabs(y_fz-SMC)/SMC);
    696. 

  696.       g_fzp.SetPoint(i+POINTS,(x_fzp-FzPmc)*100./FzPmc, 100.*fabs(y_fzp-SMC)/SMC);
    697. 

  697.       g_i.SetPoint(i+POINTS,(x_i-Imc)*100./Imc, 100.*fabs(y_i-SMC)/SMC);
    698. 

  698.       g_k.SetPoint(i+POINTS,(x_k-KMC)*100./KMC, 100.*fabs(y_k-SMC)/SMC);
    699. 

  699.     }
    700. 

  700.     TString yaxis("(S-S_{mc})/S_{mc} (%)");
    701. 

  701.     TCanvas c;
    702. 

  702.     g_ewk.SetLineWidth(2);
    703. 

  703.     g_ewk.GetXaxis()->SetTitle("EWK Variation (%)");
    704. 

  704.     g_ewk.GetYaxis()->SetTitle(yaxis);
    705. 

  705.     g_ewk.Draw("AL");
    706. 

  706.     c.Print("ewk_syst_variation.C");
    707. 

  707.     //
    708. 

  708.     g_fz.SetLineWidth(2);
    709. 

  709.     g_fz.GetXaxis()->SetTitle("F_{z} Variation (%)");
    710. 

  710.     g_fz.GetYaxis()->SetTitle(yaxis);
    711. 

  711.     g_fz.Draw("AL");
    712. 

  712.     c.Print("fz_syst_variation.C");
    713. 

  713.     //
    714. 

  714.     g_fzp.SetLineWidth(2);
    715. 

  715.     g_fzp.GetXaxis()->SetTitle("F_{z}' Variation (%)");
    716. 

  716.     g_fzp.GetYaxis()->SetTitle(yaxis);
    717. 

  717.     g_fzp.Draw("AL");
    718. 

  718.     c.Print("fzp_syst_variation.C");
    719. 

  719.     //
    720. 

  720.     g_i.SetLineWidth(2);
    721. 

  721.     g_i.GetXaxis()->SetTitle("I Variation (%)");
    722. 

  722.     g_i.GetYaxis()->SetTitle(yaxis);
    723. 

  723.     g_i.Draw("AL");
    724. 

  724.     c.Print("i_syst_variation.C");
    725. 

  725.     //
    726. 

  726.     g_k.SetLineWidth(2);
    727. 

  727.     g_k.GetXaxis()->SetTitle("K Variation (%)");
    728. 

  728.     g_k.GetYaxis()->SetTitle(yaxis);
    729. 

  729.     g_k.Draw("AL");
    730. 

  730.     c.Print("k_syst_variation.C");
    731. 

  731.     //
    732. 

  732.     // ******************************************************************
    733. 

  733.     //
    734. 

  734.     //
    735. 

  735.     // 
    736. 

  736.     double err_ewk = std::max(fabs(SMC-S_EWK_PLUS),fabs(SMC-S_EWK_MINUS));
    737. 

  737.     double err_fz = fabs(SMC-S_FZ);
    738. 

  738.     double err_fzp = fabs(SMC-S_FZP);
    739. 

  739.     double err_i  = fabs(SMC-S_I);
    740. 

  740.     double err_k = fabs(SMC-S_K);
    741. 

  741.     //
    742. 

  742.     double DS_syst = sqrt(err_ewk*err_ewk + err_fz*err_fz + err_fzp*err_fzp+
    743. 

  743. 			  err_i*err_i + err_k*err_k);
    744. 

  744.     //
    745. 

  745.     cout << "********************************************************" << endl;
    746. 

  746.     cout << "Signal Prediction: " << S << "+-" << DS << "(stat) +-"
    747. 

  747. 	 << DS_syst << "(syst)"  << endl;
    748. 

  748.     cout << "stat error: " << 100.*DS/S <<"%" << endl;
    749. 

  749.     cout << "syt  error: " << 100.*DS_syst/S<< "%"  << endl;
    750. 

  750.     cout << "********************************************************" << endl;
    751. 

  751.     cout << "Parameters used in calculation: " << endl;
    752. 

  752.     cout << "I=  " << I << "+-" << dI << endl;
    753. 

  753.     cout << "Fz= " << Fz << "+-" << dFz << endl;
    754. 

  754.     cout << "FzP=" << FzP << "+-" << dFzP << endl;
    755. 

  755.     cout << "EWK error assumed to be: " << ewkerror << endl;
    756. 

  756.     cout << endl;
    757. 

  757.     cout << "ABCD Regions population:" << endl;
    758. 

  758.     cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
    759. 

  759.          << ", ewk=" << a_ewk << endl;
    760. 

  760.     cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
    761. 

  761.          << ", ewk=" << b_ewk << endl;
    762. 

  762.     cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
    763. 

  763.          << ", ewk=" << c_ewk << endl;
    764. 

  764.     cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
    765. 

  765.          << ", ewk=" << d_ewk << endl;
    766. 

  766.     cout << endl;
    767. 

  767.     cout << "Parameters from MC: " << endl;
    768. 

  768.     cout << "I=  " << Imc << "+-" << dImc << endl;
    769. 

  769.     cout << "Fz= " << Fzmc << "+-" << dFzmc << endl;
    770. 

  770.     cout << "FzP=" << FzPmc << "+-" << dFzPmc << endl;
    771. 

  771.     cout << endl;
    772. 

  772.     cout << "Real value of K=" << KMC << endl;
    773. 

  773.     cout << "Real value of Signal=" << SMC << endl;
    774. 

  774.     cout << endl;
    775. 

  775.     cout << "Difference Measured - MC value (% wrt MC value except K=1): " 
    776. 

  776. 	 << endl;
    777. 

  777.     cout << "Fz : " << dfz  << ", (" << dfz*100./Fzmc << "%)" << endl;
    778. 

  778.     cout << "FzP: " << dfzp << ", (" << dfzp*100./FzPmc << "%)"  << endl;
    779. 

  779.     cout << "I  : " << di   << ", (" << di*100./Imc << "%)"  << endl;
    780. 

  780.     cout << "K  : " << fk   << ", (" << fk*100./1. << "%)"  << endl;
    781. 

  781.     cout << endl;
    782. 

  782.     //
    783. 

  783.     cout << "DETAILS OF THE CALCULATION" << endl;
    784. 

  784.     cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl;
    785. 

  785.     cout << "Statistical Error Summary: " << endl;
    786. 

  786.     cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
    787. 

  787.     cout << "due to FzP= "<< SpFzP*dFzP
    788. 

  788. 	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
    789. 

  789.     cout << "due to  I = "<< SpI*dI
    790. 

  790. 	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
    791. 

  791.     cout << "due to Na = "<< SpNa*sqrt(Na)
    792. 

  792. 	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
    793. 

  793.     cout << "due to Nb = "<< SpNb*sqrt(Nb)
    794. 

  794. 	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
    795. 

  795.     cout << "due to Nc = "<< SpNc*sqrt(Nc)
    796. 

  796. 	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
    797. 

  797.     cout << "due to Nd = "<< SpNd*sqrt(Nd)
    798. 

  798. 	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
    799. 

  799.     cout << "Total Statistical Error: " 
    800. 

  800. 	 << DS << ", (" << DS*100./S << "%)"<< endl;
    801. 

  801.     cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
    802. 

  802.     cout << endl;
    803. 

  803.     cout << "Systematic Error Summary:" << endl;
    804. 

  804.     cout << "due to k   = " << err_k << " ( " << err_k*100./S  << "%)" << endl;
    805. 

  805.     cout << "due to Fz  = " << err_fz << " ( " << err_fz*100./S  << "%)" << endl;
    806. 

  806.     cout << "due to FzP = " << err_fzp << " ( " << err_fzp*100./S  << "%)" << endl;
    807. 

  807.     cout << "due to I   = " << err_i << " ( " << err_i*100./S  << "%)" << endl;
    808. 

  808.     cout << "due to EWK = " << err_ewk << " ( " << err_ewk*100./S  << "%)" << endl;
    809. 

  809. 

  810.     cout << "Syst Error percentages are wrt S prediction, not S mc" << endl;
    811. 

  811.   }
    812. 

  812.   //
    813. 

  813.   //
    814. 

  814.   if (mcOnly < 0 && mcOnly >-0.75) {  // select value -0.5 that gives the
    815. 

  815.                                       // string parser
    816. 

  816.     cout << "=======================================================" << endl;
    817. 

  817.     cout << "Calculating ABCD Result and Stat Error Assuming MC ONLY"  << endl;
    818. 

  818.     cout << "=======================================================" << endl;
    819. 

  819.     cout << "All input parameters that the user have inserted will be "
    820. 

  820. 	 << "ignored and recalculated from MC" << endl;
    821. 

  821.     cout << "This option will not give you systematics estimation" << endl;
    822. 

  822.     // recalculate the basic quantities
    823. 

  823.     I = (a_sig + b_sig) / (a_sig + b_sig + c_sig + d_sig);
    824. 

  824.     dI = sqrt(I*(1-I)/(a_sig + b_sig + c_sig + d_sig));
    825. 

  825.     Fz = a_sig/b_sig;
    826. 

  826.     double e =a_sig/(a_sig + b_sig);
    827. 

  827.     double de = sqrt(e*(1-e)/(a_sig + b_sig));
    828. 

  828.     double alpha = de/(2.*Fz-e);
    829. 

  829.     dFz = alpha/(1-alpha);
    830. 

  830.     FzP = d_sig/c_sig;
    831. 

  831.     double ep =d_sig/(c_sig + d_sig);
    832. 

  832.     double dep = sqrt(ep*(1-ep)/(c_sig + d_sig));
    833. 

  833.     double alphap = dep/(2.*FzP-ep);
    834. 

  834.     dFzP = alphap/(1-alphap);
    835. 

  835.     //
    836. 

  836.     double KMC = (d_qcd/c_qcd)/(a_qcd/b_qcd);
    837. 

  837.     //
    838. 

  838.     // now everything is done from data + input
    839. 

  839.     double A = (1.0-I)*(FzP-Fz);
    840. 

  840.     double B = I*(FzP+1.0)*(FzP*(c_sig+c_qcd)-(d_sig+d_qcd)) + 
    841. 

  841.       (1+Fz)*(1-I)*((a_sig+a_qcd)-dFzP*(b_sig+b_qcd));
    842. 

  842.     double C = I*(1.+Fz)*(1.+FzP)*((d_sig+d_qcd)*(b_sig+b_qcd) - 
    843. 

  843. 				   (a_sig+a_qcd)*(c_sig+c_qcd));
    844. 

  844.     //
    845. 

  845.     // signal calculation:
    846. 

  846.     double S = Trionym(A,B,C, a_sig+b_sig);
    847. 

  847. 

  848.     // the errors now:
    849. 

  849.     // calculate the statistical error now:
    850. 

  850.     double  ApI=0, ApFz=0, ApFzP=0, ApNa=0, ApNb=0, ApNc=0, ApNd=0;
    851. 

  851.     double  BpI=0, BpFz=0, BpFzP=0, BpNa=0, BpNb=0, BpNc=0, BpNd=0;
    852. 

  852.     double  CpI=0, CpFz=0, CpFzP=0, CpNa=0, CpNb=0, CpNc=0, CpNd=0;
    853. 

  853.     double  SpI=0, SpFz=0, SpFzP=0, SpNa=0, SpNb=0, SpNc=0, SpNd=0;
    854. 

  854.     //
    855. 

  855.     double Na = a_sig+a_qcd+a_ewk, Nb = b_sig+b_qcd+b_ewk;
    856. 

  856.     double Nc=c_sig+c_qcd+c_ewk,   Nd = d_sig+d_qcd+d_ewk;
    857. 

  857.     double Ea = a_ewk, Eb = b_ewk, Ec=c_ewk, Ed = d_ewk;
    858. 

  858.     if (A != 0) {
    859. 

  859.       
    860. 

  860.       ApI   = -(FzP-Fz);
    861. 

  861.       ApFz  = -(1.0-I);
    862. 

  862.       ApFzP = (1.0-I);
    863. 

  863.       ApNa  = 0.0;
    864. 

  864.       ApNb  = 0.0;
    865. 

  865.       ApNc  = 0.0;
    866. 

  866.       ApNd  = 0.0;
    867. 

  867.       
    868. 

  868.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    869. 

  869.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    870. 

  870.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    871. 

  871.       BpNa  =  (1.0-I)*(1.0+Fz);
    872. 

  872.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    873. 

  873.       BpNc  = I*(FzP+1.0)*Fz;
    874. 

  874.       BpNd  = -I*(FzP+1.0);
    875. 

  875.       
    876. 

  876.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    877. 

  877.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    878. 

  878.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    879. 

  879.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    880. 

  880.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    881. 

  881.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    882. 

  882.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    883. 

  883.       
    884. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  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. 

  891.     }
    892. 

  892.     else {
    893. 

  893.       BpI   = (FzP+1.0)*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0+Fz)*((Na-Ea)-FzP*(Nb-Eb));
    894. 

  894.       BpFz  = I*(FzP+1.0)*(Nc-Ec)+(1.0-I)*((Na-Ea)-FzP*(Nb-Eb));
    895. 

  895.       BpFzP = I*(Fz*(Nc-Ec)-(Nd-Ed))-(1.0-I)*(1.0+Fz)*(Nb-Eb);
    896. 

  896.       BpNa  =  (1.0-I)*(1.0+Fz);
    897. 

  897.       BpNb  = -(1.0-I)*(1.0+Fz)*FzP;
    898. 

  898.       BpNc  = I*(FzP+1.0)*Fz;
    899. 

  899.       BpNd  = -I*(FzP+1.0);
    900. 

  900.       
    901. 

  901.       CpI   = (1.0+Fz)*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec)); 
    902. 

  902.       CpFz  = I*(1.0+FzP)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    903. 

  903.       CpFzP = I*(1.0+Fz)*((Nd-Ed)*(Nb-Eb)-(Na-Ea)*(Nc-Ec));
    904. 

  904.       CpNa  = -I*(1.0+Fz)*(1.0+FzP)*(Nc-Ec);
    905. 

  905.       CpNb  =  I*(1.0+Fz)*(1.0+FzP)*(Nd-Ed);
    906. 

  906.       CpNc  = -I*(1.0+Fz)*(1.0+FzP)*(Na-Ea);
    907. 

  907.       CpNd  =  I*(1.0+Fz)*(1.0+FzP)*(Nb-Eb);
    908. 

  908.       
    909. 

  909.       SpI   = -CpI/B+C*BpI/(B*B);
    910. 

  910.       SpFz  = -CpFz/B+C*BpFz/(B*B);
    911. 

  911.       SpFzP = -CpFzP/B+C*BpFzP/(B*B);
    912. 

  912.       SpNa  = -CpNa/B+C*BpNa/(B*B);
    913. 

  913.       SpNb  = -CpNb/B+C*BpNb/(B*B);
    914. 

  914.       SpNc  = -CpNc/B+C*BpNc/(B*B);
    915. 

  915.       SpNd  = -CpNd/B+C*BpNd/(B*B);
    916. 

  916.     }
    917. 

  917.     double  DS;
    918. 

  918.     DS = sqrt( SpI*dI*SpI*dI + SpFz*dFz*SpFz*dFz + SpFzP*dFzP*SpFzP*dFzP  +
    919. 

  919. 	       SpNa*SpNa*Na + SpNb*SpNb*Nb + SpNc*SpNc*Nc + SpNd*SpNd*Nd );
    920. 

  920.     // warning: S here denotes the method prediction ..........
    921. 

  921.     cout << "********************************************************" << endl;
    922. 

  922.     cout << "Signal Prediction: " << S << "+-" << DS << "(stat)" << endl;
    923. 

  923.     cout << "********************************************************" << endl;
    924. 

  924.     cout << "Parameters used in calculation: " << endl;
    925. 

  925.     cout << "I=  " << I << "+-" << dI << endl;
    926. 

  926.     cout << "Fz= " << Fz << "+-" << dFz << endl;
    927. 

  927.     cout << "FzP=" << FzP << "+-" << dFzP << endl;
    928. 

  928.     cout << endl;
    929. 

  929.     cout << "ABCD Regions population:" << endl;
    930. 

  930.     cout << "A:  N=" << Na << ", sig=" << a_sig << ", qcd=" << a_qcd
    931. 

  931. 	 << ", ewk=" << a_ewk << endl;
    932. 

  932.     cout << "B:  N=" << Nb << ", sig=" << b_sig << ", qcd=" << b_qcd
    933. 

  933. 	 << ", ewk=" << b_ewk << endl;
    934. 

  934.     cout << "C:  N=" << Nc << ", sig=" << c_sig << ", qcd=" << c_qcd
    935. 

  935. 	 << ", ewk=" << c_ewk << endl;
    936. 

  936.     cout << "D:  N=" << Nd << ", sig=" << d_sig << ", qcd=" << d_qcd
    937. 

  937. 	 << ", ewk=" << d_ewk << endl;
    938. 

  938.     cout << "K value from MC: " << KMC << endl;
    939. 

  939.     cout << endl;
    940. 

  940.     cout << "Statistical Error Summary: " << endl;
    941. 

  941.     cout << "due to Fz = "<< SpFz*dFz<< ", ("<<SpFz*dFz*100./S << "%)"<< endl;
    942. 

  942.     cout << "due to FzP= "<< SpFzP*dFzP
    943. 

  943. 	 << ", ("<<SpFzP*dFzP*100./S << "%)"<< endl; 
    944. 

  944.     cout << "due to  I = "<< SpI*dI
    945. 

  945. 	 << ", ("<<SpI*dI*100./S << "%)"<< endl; 
    946. 

  946.     cout << "due to Na = "<< SpNa*sqrt(Na)
    947. 

  947. 	 << ", ("<< SpNa*sqrt(Na)*100./S << "%)"<< endl; 
    948. 

  948.     cout << "due to Nb = "<< SpNb*sqrt(Nb)
    949. 

  949. 	 << ", ("<< SpNb*sqrt(Nb)*100./S << "%)"<< endl; 
    950. 

  950.     cout << "due to Nc = "<< SpNc*sqrt(Nc)
    951. 

  951. 	 << ", ("<< SpNc*sqrt(Nc)*100./S << "%)"<< endl; 
    952. 

  952.     cout << "due to Nd = "<< SpNd*sqrt(Nd)
    953. 

  953. 	 << ", ("<< SpNd*sqrt(Nd)*100./S << "%)"<< endl; 
    954. 

  954.     cout << "Total Statistical Error: " 
    955. 

  955. 	 << DS << ", (" << DS*100./S << "%)"<< endl;
    956. 

  956.     cout << "Stat Error percentages are wrt S prediction, not S mc" << endl;
    957. 

  957.   }
    958. 

  958. 

  959. 

  960. }
    961. 

  961. 

  962. void plotMaker(TString histoName, TString wzsignal,
    963. 

  963. 	       vector<TString> file, vector<TString> type, 
    964. 

  964. 	       vector<double> weight, TString xtitle)
    965. 

  965. {
    966. 

  966.   gROOT->Reset();
    967. 

  967.   gROOT->ProcessLine(".L tdrstyle.C"); 
    968. 

  968.   gROOT->ProcessLine("setTDRStyle()");
    969. 

  969.   // automatic recognition of histogram dimension
    970. 

  970.   int NBins = 0; double min = 0; double max = -1;
    971. 

  971.   for (int i=0; i< (int) file.size(); ++i) {
    972. 

  972.     if (weight[i]>0) {
    973. 

  973.       TFile f(file[i]);
    974. 

  974.       TH1F *h = (TH1F*) f.Get(histoName);
    975. 

  975.       NBins = h->GetNbinsX();
    976. 

  976.       min = h->GetBinLowEdge(1);
    977. 

  977.       max = h->GetBinLowEdge(NBins+1);
    978. 

  978.       break;
    979. 

  979.     }
    980. 

  980.   }
    981. 

  981.   if (NBins ==0 || (max<min)) {
    982. 

  982.     std::cout << "PlotCombiner::abcd error: Could not find valid histograms in file"
    983. 

  983.               << std::endl;
    984. 

  984.     abort();
    985. 

  985.   }
    986. 

  986.   cout << "Histograms with "<< NBins <<" bins  and range " << min << "-" << max  << endl;
    987. 

  987.   // Wenu Signal .......................................................
    988. 

  988.   TH1F h_wenu("h_wenu", "h_wenu", NBins, min, max);
    989. 

  989.   int fmax = (int) file.size();
    990. 

  990.   for (int i=0; i<fmax; ++i) {
    991. 

  991.     if (type[i] == "sig" && weight[i]>0) {
    992. 

  992.       TFile f(file[i]);
    993. 

  993.       TH1F *h = (TH1F*) f.Get(histoName);
    994. 

  994.       h_wenu.Add(h, weight[i]);
    995. 

  995.     }
    996. 

  996.   }
    997. 

  997.   // Bkgs ..............................................................
    998. 

  998.   //
    999. 

  999.   // QCD light flavor
    1000. 

  1000.   TH1F h_qcd("h_qcd", "h_qcd", NBins, min, max);
    1001. 

  1001.   for (int i=0; i<fmax; ++i) {
    1002. 

  1002.     if (type[i] == "qcd" && weight[i]>0) {
    1003. 

  1003.       TFile f(file[i]);
    1004. 

  1004.       TH1F *h = (TH1F*) f.Get(histoName);
    1005. 

  1005.       h_qcd.Add(h, weight[i]);
    1006. 

  1006.     }
    1007. 

  1007.   }
    1008. 

  1008.   // QCD heavy flavor
    1009. 

  1009.   TH1F h_bce("h_bce", "h_bce", NBins, min, max);
    1010. 

  1010.   for (int i=0; i<fmax; ++i) {
    1011. 

  1011.     if (type[i] == "bce" && weight[i]>0) {
    1012. 

  1012.       TFile f(file[i]);
    1013. 

  1013.       TH1F *h = (TH1F*) f.Get(histoName);
    1014. 

  1014.       h_bce.Add(h, weight[i]);
    1015. 

  1015.     }
    1016. 

  1016.   }
    1017. 

  1017.   // QCD Gjets
    1018. 

  1018.   TH1F h_gj("h_gj", "h_gj", NBins, min, max);
    1019. 

  1019.   for (int i=0; i<fmax; ++i) {
    1020. 

  1020.     if (type[i] == "gje" && weight[i]>0) {
    1021. 

  1021.       TFile f(file[i]);
    1022. 

  1022.       TH1F *h = (TH1F*) f.Get(histoName);
    1023. 

  1023.       h_gj.Add(h, weight[i]);
    1024. 

  1024.     }
    1025. 

  1025.   }
    1026. 

  1026.   // Other EWK bkgs
    1027. 

  1027.   TH1F h_ewk("h_ewk", "h_ewk", NBins, min, max);
    1028. 

  1028.   for (int i=0; i<fmax; ++i) {
    1029. 

  1029.     if (type[i] == "ewk" && weight[i]>0) {
    1030. 

  1030.       TFile f(file[i]);
    1031. 

  1031.       TH1F *h = (TH1F*) f.Get(histoName);
    1032. 

  1032.       h_ewk.Add(h, weight[i]);
    1033. 

  1033.     }
    1034. 

  1034.   }
    1035. 

  1035.   //
    1036. 

  1036.   // ok now decide how to plot them:
    1037. 

  1037.   // first the EWK bkgs
    1038. 

  1038.   h_ewk.SetFillColor(3);
    1039. 

  1039.   //
    1040. 

  1040.   // then the gjets
    1041. 

  1041.   h_gj.Add(&h_ewk);
    1042. 

  1042.   h_gj.SetFillColor(1);
    1043. 

  1043.   // thent the QCD dijets
    1044. 

  1044.   h_bce.Add(&h_qcd);
    1045. 

  1045.   h_bce.Add(&h_gj);
    1046. 

  1046.   h_bce.SetFillColor(2);
    1047. 

  1047.   // and the signal at last
    1048. 

  1048.   TH1F h_tot("h_tot", "h_tot", NBins, min, max);
    1049. 

  1049.   h_tot.Add(&h_bce);
    1050. 

  1050.   h_tot.Add(&h_wenu);
    1051. 

  1051.   h_wenu.SetLineColor(4);  h_wenu.SetLineWidth(2);
    1052. 

  1052.   //
    1053. 

  1053.   TCanvas c;
    1054. 

  1054.   h_tot.GetXaxis()->SetTitle(xtitle);
    1055. 

  1055.   h_tot.Draw("PE");
    1056. 

  1056.   h_bce.Draw("same");
    1057. 

  1057.   h_gj.Draw("same");
    1058. 

  1058.   h_ewk.Draw("same");
    1059. 

  1059.   h_wenu.Draw("same");
    1060. 

  1060. 

  1061.   // the Legend
    1062. 

  1062.   TLegend  leg(0.6,0.65,0.95,0.92);
    1063. 

  1063.   if (wzsignal == "wenu")
    1064. 

  1064.     leg.AddEntry(&h_wenu, "W#rightarrow e#nu","l");
    1065. 

  1065.   else
    1066. 

  1066.     leg.AddEntry(&h_wenu, "Z#rightarrow ee","l");
    1067. 

  1067.   leg.AddEntry(&h_tot, "Signal + Bkg","p");
    1068. 

  1068.   leg.AddEntry(&h_bce, "dijets","f");
    1069. 

  1069.   leg.AddEntry(&h_gj, "#gamma + jets","f");
    1070. 

  1070.   leg.AddEntry(&h_ewk,  "EWK+t#bar t", "f");
    1071. 

  1071.   leg.Draw("same");
    1072. 

  1072. 

  1073.   c.Print("test.png");
    1074. 

  1074. 

  1075. 

  1076. 

  1077. }
    1078. 

  1078. 

  1079. 

  1080. //
    1081. 

  1081. // reads the ABCD string and returns its value
    1082. 

  1082. // value is whatever it exists after the = sign and
    1083. 

  1083. // before the comma or the parenthesis
    1084. 

  1084. //
    1085. 

  1085. // if the string is not contained returns -1
    1086. 

  1086. // if there is no value, but the string is contained -0.5
    1087. 

  1087. // if there is an error in the algorithm return -99 and print error
    1088. 

  1088. // else returns its value
    1089. 

  1089. double searchABCDstring(TString abcdString, TString keyword)
    1090. 

  1090. {
    1091. 

  1091.   int size = keyword.Sizeof();
    1092. 

  1092.   int existsEntry = abcdString.Index(keyword);
    1093. 

  1093.   //
    1094. 

  1094.   if (existsEntry==-1) return -1.;
    1095. 

  1095.   //
    1096. 

  1096.   TString previousVal = abcdString(existsEntry-1);
    1097. 

  1097.   if (!(previousVal == "," || previousVal == " " || 
    1098. 

  1098. 	previousVal == "(" )) return -1.;
    1099. 

  1099.   //
    1100. 

  1100.   TString afterVal = abcdString(existsEntry+size-1);
    1101. 

  1101.   //std::cout << "afterVal=" << afterVal << std::endl;
    1102. 

  1102.   if (afterVal =="," || afterVal==")") return -0.5;
    1103. 

  1103.   else if (afterVal != "=") return -1.;
    1104. 

  1104.   //
    1105. 

  1105.   // now find the comma or the parenthesis after the = sign
    1106. 

  1106.   int comma = abcdString.Index(",",existsEntry);
    1107. 

  1107.   //std::cout << "first comma=" << comma << std::endl;
    1108. 

  1108.   if (comma<0) comma = abcdString.Index(")",existsEntry);
    1109. 

  1109.   if (comma<0) {
    1110. 

  1110.     std::cout << "Error in parcing abcd line, chech syntax " 
    1111. 

  1111. 	      << std::endl;
    1112. 

  1112.     return -99.;
    1113. 

  1113.   }
    1114. 

  1114.   TString svalue=abcdString(existsEntry+size,comma-existsEntry-size);
    1115. 

  1115.   std::cout << "Found ABCD parameter "<< keyword
    1116. 

  1116. 	    << " with value "  << svalue << endl;
    1117. 

  1117.   // convert this to a float
    1118. 

  1118.   double value = svalue.Atof();
    1119. 

  1119.   return value;
    1120. 

  1120. 

  1121. }
    1122. 

  1122. 

  1123. 

  1124. double Trionym(double a, double b, double c, double sum)
    1125. 

  1125. {
    1126. 

  1126.   if (a==0) {
    1127. 

  1127.     return -c/b;
    1128. 

  1128.   }
    1129. 

  1129.   double D2 = b*b - 4.*a*c;
    1130. 

  1130.   //return (-b + sqrt(D2)) / (2.*a);
    1131. 

  1131.   if (D2 > 0) {
    1132. 

  1132.     double s1 = (-b + sqrt(D2)) / (2.*a);
    1133. 

  1133.     double s2 = (-b - sqrt(D2)) / (2.*a);
    1134. 

  1134.     double solution =   fabs(s1-sum)<fabs(s2-sum)?s1:s2;
    1135. 

  1135.     return solution;
    1136. 

  1136.   }
    1137. 

  1137.   else  {
    1138. 

  1138.     return -1.;  
    1139. 

  1139.   }
    1140. 

  1140. }
    1141. 

  1141. 

  1142. //
    1143. 

  1143. // the naming of the variables and the order is in this way for historical
    1144. 

  1144. // reasons
    1145. 

  1145. // your complains for different Nd_ and nd to G.Daskalakis :P
    1146. 

  1146. //
    1147. 

  1147. double CalcABCD
    1148. 

  1148. (double I, double Fz, double FzP, double K, double ewk,
    1149. 

  1149.  double Na_, double Nb_, double Nc_, double Nd_ ,
    1150. 

  1150.  double Ea_, double Eb_, double Ec_, double Ed_)
    1151. 

  1151. {
    1152. 

  1152.   double A, B, C;
    1153. 

  1153.   A = (1.0-I)*(FzP-K*Fz);
    1154. 

  1154.   B = I*(FzP+1.0)*(K*Fz*(Nc_-ewk*Ec_)-(Nd_-ewk*Ed_))+
    1155. 

  1155.     (1.0-I)*(1.0+Fz)*(K*(Na_-ewk*Ea_)-FzP*(Nb_-ewk*Eb_));
    1156. 

  1156.   C = I*(1.0+Fz)*(1.0+FzP)*((Nd_-ewk*Ed_)*(Nb_-ewk*Eb_)-
    1157. 

  1157. 			    K*(Na_-ewk*Ea_)*(Nc_-ewk*Ec_));
    1158. 

  1158.   //
    1159. 

  1159.   return Trionym(A, B, C, Na_+Nb_-Ea_-Eb_);
    1160. 

  1160. 

  1161. }
    1162.