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