ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/iJZB_Analysis.C
Revision: 1.3
Committed: Wed May 22 07:50:52 2013 UTC (11 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +17 -16 lines
Error occurred while calculating annotation data.
Log Message:
Added calculation of global significance

File Contents

# Content
1 #include <iostream>
2 #include <vector>
3 #include <sys/stat.h>
4 #include <getopt.h>
5 #include <stdio.h>
6 #include <stdlib.h>
7
8 #ifndef Verbosity
9 #define Verbosity 0
10 #endif
11 #ifndef HUSH
12 #define HUSH 1
13 #endif
14
15 #include "Modules/GeneralToolBox.C"
16 #include "Modules/SampleClass.C"
17 #include "Modules/setTDRStyle.C"
18 #include "Modules/Setup.C"
19 #include "Modules/Poisson_Calculator.C"
20 #include "Modules/JSON/JSONSampleLoader.C"
21 #include "Modules/ActiveSamples.C"
22 #include "Modules/PeakFinder.C"
23 #include "Modules/UpperLimitsWithShape.C"
24 #include "Modules/Plotting_Functions.C"
25 #include "Modules/LimitCalculation.C"
26 #include "Modules/ResultModule.C"
27 #include "Modules/CrossSectionReader.C"
28 #include "Modules/Systematics.C"
29 #include "Modules/SugarCoating.C"
30 #include "Modules/ExclusionPlot.C"
31 #include "Modules/SUSYScan.C"
32 #include "Modules/AachenCompatibility.C"
33 #include "Modules/MetPlotting.C"
34 #include "Modules/WZStudy.C"
35 //#include "Modules/FSRStudy.C"
36 #include "Modules/TTbarAnalysis.C"
37 #include "Modules/OverlapStudy.C"
38
39
40 #include <TCut.h>
41 #include <TROOT.h>
42 #include <TCanvas.h>
43 #include <TMath.h>
44 #include <TColor.h>
45 #include <TPaveText.h>
46 #include <TRandom.h>
47 #include <TH1.h>
48 #include <TH2.h>
49 #include <TF1.h>
50 #include <TSQLResult.h>
51
52 using namespace PlottingSetup;
53
54 void usage(int passed=0 ) {
55 std::cout << "USAGE : " << std::endl;
56 std::cout << "You can use different options when running this program : " << std::endl;
57 std::cout << std::endl;
58 std::cout << "\033[1;34m all\033[0m \t\t All processes of the standard workflow" << std::endl;
59 std::cout << "\033[1;34m peak\033[0m \t\t Peak finding" << std::endl;
60 std::cout << "\033[1;34m pileup\033[0m\t\t Pileup correction determination" << std::endl;
61 std::cout << "\033[1;34m kin\033[0m \t\t Kinematic plots" << std::endl;
62 std::cout << "\033[1;34m pfkin\033[0m \t Kinematic plots for PF variables" << std::endl;
63 std::cout << "\033[1;34m region\033[0m \t Region comparison plots" << std::endl;
64 std::cout << "\033[1;34m lepton\033[0m \t Lepton comparison plots" << std::endl;
65 std::cout << "\033[1;34m jzb\033[0m \t\t JZB plots" << std::endl;
66 std::cout << "\033[1;34m pred\033[0m \t\t Prediction plots" << std::endl;
67 std::cout << "\033[1;34m signalbg\033[0m \t Signal vs. background plots (shape comparison)" << std::endl;
68 std::cout << "\033[1;34m ttbar\033[0m \t TTbar sideband comparison" << std::endl;
69 std::cout << "\033[1;34m ttbaranalysis\033[0m \t TTbar analysis. To prove we can." << std::endl;
70 std::cout << "\033[1;34m ttbartest\033[0m \t TTbar prediction tests (peak, response)" << std::endl;
71 std::cout << "\033[1;34m zjets\033[0m \t ZJets plots" << std::endl;
72 std::cout << "\033[1;34m diboson\033[0m \t Diboson plots" << std::endl;
73 std::cout << "\033[1;34m results\033[0m \t Compute results (JZB>X)" << std::endl;
74 std::cout << "\033[1;34m syst\033[0m \t\t Compute systematics" << std::endl;
75 std::cout << "\033[1;34m eff\033[0m \t\t Compute all (partial) efficiencies" << std::endl;
76 std::cout << "\033[1;34m shapes\033[0m \t Prepare ROOT files for computing limits using shapes" << std::endl;
77 std::cout << "\033[1;34m counting\033[0m \t Establish limits with counting experiment" << std::endl;
78 std::cout << "\033[1;34m shapelimit\033[0m \t Establish limits with shapes" << std::endl;
79 std::cout << "\033[1;34m qcd\033[0m \t\t Show that QCD is negligible" << endl;
80 std::cout << "\033[1;34m response\033[0m \t Compute response correction" << std::endl;
81 std::cout << "\033[1;34m pickup\033[0m \t Pick up interesting events" << std::endl;
82 std::cout << "\033[1;34m yields\033[0m \t Compute MC yields for this lumi" << std::endl;
83 std::cout << "\033[1;34m sideassess\033[0m\t Assessment of the mass sidebands definition in the data" << std::endl;
84 std::cout << "\033[1;34m jes\033[0m \t\t Make Jet Energy Scale (JES) plot" << std::endl;
85 std::cout << "\033[1;34m nosidebands\033[0m \t\t Switch off sidebands for cJZB" << std::endl;
86 std::cout << "\033[1;34m effcurve\033[0m \t\t Make JZB efficiency plot" << std::endl;
87 std::cout << "\033[1;34m save\033[0m \t\t Save results in a template file (for reuse when establishing limits for scans)" << std::endl;
88 std::cout << "\033[1;34m png\033[0m \t\t Save all plots as pngs (also available: --pdf, --root, --eps, --C) standard: png&pdf" << std::endl;
89 std::cout << "\033[1;34m offpeak\033[0m \t Sets offpeak to true (i.e. carry out off-peak analysis!)" << std::endl;
90 std::cout << "\033[1;34m onpeak\033[0m \t Sets onpeak to true (i.e. carry out standard on-peak analysis!)" << std::endl;
91 std::cout << "\033[1;34m btagged (b)\033[0m \t Sets DoBTag to true (i.e. carry out b tagged analysis)" << std::endl;
92 std::cout << "\033[1;34m jzbcuts\033[0m \t Sets the jzb cuts (please provide them in the format 50,100,150,200 and so on (comma separated)" << endl;
93 std::cout << "\033[1;34m ptsanity\033[0m \t Checks whether all samples have the correct pt range (additional canvas for you to check)" << endl;
94 std::cout << "\033[1;34m mlls\033[0m \t Draw dilepton invariant mass for all signals" << endl;
95 std::cout << "\033[1;34m rmue\033[0m \t Compute rmu" << std::endl;
96 std::cout << "\033[1;34m overlap\033[0m \t Study overlap between iJZB and METJZB" << endl;
97 // std::cout << "\033[1;34m edge\033[0m \t Be on the edge of development and try out the new edge fitter :-)" << endl;
98 std::cout << "\033[1;34m wz \033[0m \t WZ" << std::endl;
99 std::cout << "\033[1;34m ittbar\033[0m\t ttbar closure test for iJZB" << std::endl;
100 // std::cout << "\033[1;34m fsr \033[0m \t Carry out all FSR studies (you can switch off the ones you don't want in FSRStudy.C)" << endl;
101 std::cout << std::endl;
102 std::cout << "\033[1;34m ijzb \033[om\t This is iJZB. Rerun peak positions and everything" << std::endl;
103 std::cout << "\033[1;34m dir (d)\033[0m\t Directory where all plots will be saved" << std::endl;
104 std::cout << "\033[1;34m public \033[0m\t Switch to be activated when using the algorithm in public (no branding etc.)" << std::endl;
105 std::cout << std::endl;
106 std::cout << std::endl;
107 std::cout << "\033[1;34m met\033[0m \t\t Compute obs/pred for met cut vs jzb cut (in data, lm4 and lm8)" << std::endl;
108 std::cout << "\033[1;34m metplots\033[0m \t\t Produce the MET plots (i.e. mll plots with MET cut etc.)" << std::endl;
109 std::cout << "\033[1;34m metvsjzb\033[0m \t\t Produce MET:JZB plots" << std::endl;
110 std::cout << "\033[1;34m metpred\033[0m \t\t Experimental JZB prediction for MET signal region" << std::endl;
111 std::cout << "\033[1;34m aachen\033[0m \t\t Produce Aachen root file (please provide source root file)" << std::endl;
112 std::cout << std::endl;
113 std::cout << std::endl;
114 std::cout << "\033[1;34m paper\033[0m \t\t Switch paper mode on (no \" Preliminary \" on plots)" << std::endl;
115 std::cout << "\033[1;34m approved\033[0m \t\t Only approved plots" << std::endl;
116 exit(-1);
117 }
118
119 int main (int argc, char ** argv)
120 {
121
122
123 PlottingSetup::IsIJZBAnalysis=true;
124 int do_all = false;/// DONE
125 int do_global_significance = false; /// DONE
126 int do_peak_finding = false; /// DONE
127 int do_pileup = false;
128 int do_kinematic_variables = false; ///DONE
129 int do_kinematic_PF_variables = false; ///DONE
130 int do_region_comparison = false; /// Comparison between signal and control regions
131 int do_lepton_comparison = false; ///DONE
132 int do_jzb_plots = false; /// DONE
133 int do_pred = false; /// DONE
134 bool overlay_signal = false; /// Overlay LM4 signal to predictions
135 int do_signal_bg_comparison_plot = false; /// DONE
136 int do_ttbar_comparison = false; /// DONE
137 int do_ttbaranalysis = false;
138 int do_overlap = false;
139 int do_ttbartest = false; /// DONE
140 int do_jzb_efficiency_curve = false; /// DONE
141 int do_zjets_comparison = false; /// DONE
142 int do_diboson_plots = false;
143 // int calculate_pred_and_observed = false; /// DONE -- now rerouted to the new way of computing results
144 // int calculate_yields = false; ///DONE /// superfluous (replaced by get_new_results)
145 int get_new_results = false; ///DONE
146 bool dopoisson = false;//should we calculate stat err with poisson? ///DONE
147 bool verbose = false;//prints out a lot more information ... ///DONE
148 //--------------------------------------------
149 // Systematics and limits
150 int do_compute_systematics=false; /// DONE (might want to add something related to the peak?)
151 bool requireZ=true; // should we require a Z for MC efficiency?
152 int do_compute_efficiency=false; /// DONE (might want to add something related to the peak?)
153 int do_prepare_limits_using_shapes=false; /// ALMOST DONE **************************************
154 int do_compute_upper_limits_from_counting_experiment=false; /// DONE
155 bool doobserved=false; // want to get the observed limit as well?
156 int doquick=1; // this will cause the number of predicted/observed events to hurry up (0 = do mc and all results, 1 = all results but no mc, 2 = don't save results in library and don't compute mc stuff)
157 int do_compute_upper_limits_from_shapes=false;
158
159 //--------------------------------------------
160 // More not-so-standard stuff
161 int do_model_scan=false; // this can take days - use the grid version!
162 int do_zjet_ttbar_shapes=false; /// DONE
163 int do_test=false;/// DONE (just any test you wish)
164 int do_pick_up_events=false; ///DONE
165
166 //--------------------------------------------
167 // To be run manually !
168 int do_response_correction=false; ///DONE
169 //use this to find out the correction factor; originally this was done automatically but now you should do this separately and update Modules/Setup.C
170
171 int do_sideband_assessment=false;
172 int do_qcd=false;
173
174 int do_save_template=false;
175 int do_mc_yields=false; // compute MC yields
176 int do_metjzb=false;
177
178 int savepdf=true;
179 int saveC=true;
180 int saveRoot=true;
181 int savepng=true;
182 int saveeps=false;
183
184 int dooffpeak=true;
185 int doonpeak=false;
186 int dobtag=false;
187
188 int do_JES=false;
189 int do_ptsanity=false;
190 int do_mlls=false;
191 int do_rmue=false;
192 // int do_edge=false;
193 int do_wz=false;
194 int do_ittbar=false;
195 int do_ijzb=false;
196 int do_fsr=false;
197 int do_met_vs_jzb=false;
198 int do_aachen=false;
199 int do_metplots=false;
200 int do_metpred=false;
201 string AachenFileName="";
202
203 std::string directory="";
204 int option_iterator;
205 int option_counter=0;
206 bool moreoptions=true;
207
208 int SwitchOffSidebands=0;
209
210 string jzbcuts_string="";
211 float local_signif=0.0;
212
213 while(moreoptions) {
214 static struct option long_options[] =
215 {
216 /* These options set a flag. */
217 {"all", no_argument, &do_all, 1},
218 {"peak", no_argument, &do_peak_finding, 1},
219 {"pileup", no_argument, &do_pileup, 1},
220 {"test", no_argument, &do_test, 1},
221 {"kin", no_argument, &do_kinematic_variables, 1},
222 {"pfkin", no_argument, &do_kinematic_PF_variables, 1},
223 {"region", no_argument, &do_region_comparison, 1},
224 {"lepton", no_argument, &do_lepton_comparison, 1},
225 {"jzb", no_argument, &do_jzb_plots, 1},
226 {"pred", no_argument, &do_pred, 1},
227 {"signalbg", no_argument, &do_signal_bg_comparison_plot, 1},
228 {"ttbar", no_argument, &do_ttbar_comparison, 1},
229 {"ttbaranalysis", no_argument, &do_ttbaranalysis, 1},
230 {"ttbartest",no_argument, &do_ttbartest, 1},
231 {"overlap" ,no_argument, &do_overlap, 1},
232 {"zjets", no_argument, &do_zjets_comparison, 1},
233 {"diboson", no_argument, &do_diboson_plots, 1},
234 {"results", no_argument, &get_new_results, 1},
235 {"syst", no_argument, &do_compute_systematics, 1},
236 {"eff", no_argument, &do_compute_efficiency, 1},
237 {"shapes", no_argument, &do_prepare_limits_using_shapes, 1},
238 {"counting", no_argument, &do_compute_upper_limits_from_counting_experiment, 1},
239 {"shapelimit",no_argument, &do_compute_upper_limits_from_shapes, 1},
240 {"ptsanity", no_argument, &do_ptsanity, 1},
241 {"response", no_argument, &do_response_correction, 1},
242 {"pickup", no_argument, &do_pick_up_events, 1},
243 {"yields", no_argument, &do_mc_yields, 1},
244 {"sideassess",no_argument, &do_sideband_assessment, 1},
245 {"metvsjzb", no_argument, &do_met_vs_jzb, 1},
246 {"qcd", no_argument, &do_qcd, 1},
247 {"save", no_argument, &do_save_template,1},
248 {"png", no_argument, &savepng,1},
249 {"eps", no_argument, &saveeps,1},
250 {"pdf", no_argument, &savepdf,1},
251 {"root", no_argument, &saveRoot,1},
252 {"C", no_argument, &saveC,1},
253 {"offpeak", no_argument, &dooffpeak,1},
254 {"onpeak", no_argument, &doonpeak,1},
255 {"btagged", no_argument, &dobtag,1},
256 {"met", no_argument, &do_metjzb,1},
257 {"jes", no_argument, &do_JES,1},
258 {"effcurve", no_argument, &do_jzb_efficiency_curve,1},
259 {"public", no_argument, &PlottingSetup::publicmode,1},
260 {"paper", no_argument, &PlottingSetup::PaperMode,1},
261 {"approved", no_argument, &PlottingSetup::Approved,1},
262 {"nosidebands",no_argument, &SwitchOffSidebands,1},
263 {"mlls", no_argument, &do_mlls,1},
264 {"rmue", no_argument, &do_rmue,1},
265 // {"edge", no_argument, &do_edge,1},
266 {"wz", no_argument, &do_wz,1},
267 {"ijzb", no_argument, &do_ijzb,1},
268 {"ittbar", no_argument, &do_ittbar,1},
269 {"fsr", no_argument, &do_fsr,1},
270 {"metplots", no_argument, &do_metplots,1},
271 {"metpred", no_argument, &do_metpred,1},
272 /* The following options store values.*/
273 {"jzbcuts", required_argument, 0, 'j'},
274 {"dir", required_argument, 0, 'd'},
275 {"aachen", required_argument, 0, 'c'},
276 {"gsig", required_argument, 0, 'g'},
277 {0, 0, 0, 0}
278 };
279 int option_index = 0;
280 option_iterator = getopt_long(argc, argv, "abc:d:g:j:",long_options, &option_index);
281 if(option_iterator == -1) moreoptions=false;
282 else {
283 option_counter++;
284 switch (option_iterator)
285 {
286 case 0:
287 if (long_options[option_index].flag != 0)
288 break;
289 printf ("option %s", long_options[option_index].name);
290 if (optarg)
291 printf (" with arg %s", optarg);
292 printf ("\n");
293 break;
294 case 'a':
295 do_all=true;
296 break;
297 case 'b':
298 dobtag=true;
299 break;
300 case 'c':
301 do_aachen=true;
302 AachenFileName=(std::string)optarg;
303 cout << "Received file to convert to Aachen format: " << AachenFileName << endl;
304 break;
305 case 'd':
306 directory=(std::string)optarg;
307 std::cout<<"Option directory was passed with argument " << optarg << std::endl;
308 break;
309 case 'j':
310 jzbcuts_string=(std::string)optarg;
311 std::cout<<"JZB cuts were manually defined:" << optarg << std::endl;
312 break;
313 case 'g':
314 do_global_significance=true;
315 local_signif=atof(optarg);
316 break;
317 case '?':
318 usage(option_iterator);
319 break;
320 default:
321 usage(option_iterator);
322 }
323 }
324 }
325
326 if(directory!="") PlottingSetup::directoryname=directory;
327 if(dooffpeak) PlottingSetup::RestrictToMassPeak=false;
328 if(doonpeak) PlottingSetup::RestrictToMassPeak=true;
329 if(dobtag) PlottingSetup::DoBTag=true;
330 if(option_counter==0) usage();
331 if(SwitchOffSidebands) PlottingSetup::UseSidebandsForcJZB=false;
332
333 ///----------------------------------- BELOW THIS LINE: NO MORE OPTIONS BUT ACTUAL FUNCTION CALLS! ---------------------------------------------------------
334 gROOT->SetStyle("Plain");
335 bool do_fat_line=false; // if you want to have HistLineWidth=1 and FuncWidth=1 as it was before instead of 2
336 setTDRStyle(do_fat_line);
337 gStyle->SetTextFont(42);
338 bool showList=true;
339 set_directory(directoryname);//Indicate the directory name where you'd like to save the output files in Setup.C
340
341 if(do_wz) {
342 IsWZAnalysis=true;
343 IsJZBAnalysis=false;
344 }
345
346
347 write_info(__FUNCTION__,"Making sure that we are only considering category A !!");
348
349 string backup_basicqualitycut = (const char*) basicqualitycut;
350 string backup_essentialcut = (const char*) essentialcut;
351 string backup_basiccut = (const char*) basiccut;
352 string backup_leptoncut = (const char*) leptoncut;
353 string backup_essential = (const char*) essential;
354
355 string Sbasicqualitycut = backup_basicqualitycut;
356 Sbasicqualitycut = ReplaceAll(Sbasicqualitycut,")<2.4",")<1.4");
357 basicqualitycut=TCut(Sbasicqualitycut.c_str());
358
359 string Sbasiccut = backup_basiccut;
360 Sbasiccut = ReplaceAll(Sbasiccut,")<2.4",")<1.4");
361 basiccut=TCut(Sbasiccut.c_str());
362
363 string Sessentialcut = backup_essentialcut;
364 Sessentialcut = ReplaceAll(Sessentialcut,")<2.4",")<1.4");
365 essentialcut=TCut(Sessentialcut.c_str());
366
367 string Sleptoncut = backup_leptoncut;
368 Sleptoncut = ReplaceAll(Sleptoncut,")<2.4",")<1.4");
369 leptoncut=TCut(Sleptoncut.c_str());
370
371 string Sessential = backup_essential;
372 Sessential = ReplaceAll(Sessential,")<2.4",")<1.4");
373 essential = TCut(Sessential.c_str());
374
375
376 if(do_global_significance) {
377 compute_global_significance(local_signif);
378 return 0;
379 }
380
381 define_samples(showList,allsamples,signalsamples,scansample,raresample,systsamples,qcdsamples,comparesamples);
382 setlumi(luminosity);
383 setessentialcut(essential);//this sets the essential cut; this one is used in the draw command so it is AUTOMATICALLY applied everywhere. IMPORTANT: Do NOT store weights here!
384 stringstream resultsummary;
385
386 write_analysis_type(PlottingSetup::RestrictToMassPeak,PlottingSetup::DoBTag);
387 do_png(savepng);
388 do_pdf(savepdf);
389 do_eps(saveeps);
390 do_C(saveC);
391 do_root(saveRoot);
392
393 write_info(__FUNCTION__,"Imposing standard JZB definition now");
394 /*
395 jzbvariabledata="((jzb[1]+0.0595776*pt-0.279212*numVtx)- 0.375633 )";
396 jzbvariablemc="((jzb[1]+0.0348755*pt-0.279212*numVtx)- 0.224927 )";
397 */
398 /*
399 jzbvariabledata="(jzb[1]+0.0595776*pt)";
400 jzbvariablemc="(jzb[1]+0.0348755*pt)";
401 */
402 jzbvariabledata="(jzb[1]+0.0595776*pt- 4.18673)";
403 jzbvariablemc="(jzb[1]+0.0348755*pt- 3.8028)";
404
405 /* write_warning(__FUNCTION__,"WATCH OUT, EXPERIMENTAL PT CORRECTION");
406 jzbvariabledata="jzb[1]+(1-(9.57810e-01 + -3.53527e+00/pt + 1.85206e+02/pow(pt,2)))*pt";
407 jzbvariablemc ="jzb[1]+(1-(9.60547e-01 + -1.16558e+00/pt + 2.12197e+02/pow(pt,2)))*pt";
408
409 jzbvariabledata="(jzb[1]+(1-(9.57810e-01 + -3.53527e+00/pt + 1.85206e+02/pow(pt,2)))*pt- 1.92518 )";
410 jzbvariablemc ="(jzb[1]+(1-(9.60547e-01 + -1.16558e+00/pt + 2.12197e+02/pow(pt,2)))*pt+ 0.0553847 )";
411 */
412
413 /*write_warning(__FUNCTION__,"No response correction");
414 jzbvariabledata="(jzb[1]+ 0.91084)";
415 jzbvariablemc="(jzb[1]- 0.694762)";*/
416
417
418
419 global_ratio_binning.push_back(0);
420 global_ratio_binning.push_back(5);
421 global_ratio_binning.push_back(10);
422 global_ratio_binning.push_back(20);
423 global_ratio_binning.push_back(50);
424 global_ratio_binning.push_back(100);
425 global_ratio_binning.push_back(150);
426 global_ratio_binning.push_back(200);
427 global_ratio_binning.push_back(PlottingSetup::jzbHigh);
428 //global_ratio_binning.push_back(500);
429
430 //these are the JZB cuts defining our search regions
431 vector<float> jzb_cut; //starting where, please?
432 if(jzbcuts_string=="") {
433 // jzb_cut.push_back(50);
434 // jzb_cut.push_back(75);
435 jzb_cut.push_back(80);
436 jzb_cut.push_back(100);
437 // jzb_cut.push_back(125);
438 jzb_cut.push_back(150);
439 // jzb_cut.push_back(175);
440 // jzb_cut.push_back(200);
441 // jzb_cut.push_back(225);
442 // jzb_cut.push_back(250);
443 // jzb_cut.push_back(275);
444 // jzb_cut.push_back(300);
445 } else {
446 manually_set_jzb_cuts(jzb_cut,jzbcuts_string);
447 if(jzb_cut.size()==0) {
448 write_error(__FUNCTION__,"There are no jzb cuts!");
449 return -1;
450 }
451 }
452
453 todo();
454
455
456 //**** part 1 : peak finding
457 float MCPeak=0,MCPeakError=0,DataPeak=0,DataPeakError=0,MCSigma=10,DataSigma=10;
458 method=Kostasmethod;//Kostasmethod;//dogaus3sigma;// options: dogaus,doKM,dogaus2sigma,dogaus3sigma
459 stringstream datajzb;
460 stringstream mcjzb;
461
462
463 if(do_pileup) FindPeakPileUpCorrection();
464
465 if ( do_peak_finding ) {
466 //if(do_peak_finding||do_pick_up_events||PlottingSetup::Approved||do_mlls||do_zjets_comparison||do_sideband_assessment||do_metjzb||do_zjet_ttbar_shapes||do_ttbartest||do_diboson_plots||do_region_comparison||do_jzb_plots||do_lepton_comparison||do_signal_bg_comparison_plot||do_pred||do_mc_yields||do_save_template||do_compute_upper_limits_from_counting_experiment||do_compute_upper_limits_from_shapes||do_compute_systematics||do_model_scan||do_prepare_limits_using_shapes||do_all||get_new_results||do_kinematic_variables||do_kinematic_PF_variables||do_aachen) find_peaks(MCPeak,MCPeakError, DataPeak, DataPeakError,resultsummary,true,datajzb,mcjzb);
467 find_peaks(MCPeak,MCPeakError, DataPeak, DataPeakError,resultsummary,true,datajzb,mcjzb);
468
469 if(datajzb.str().size()<1) datajzb<<"("<<jzbvariabledata<<")";
470 if(mcjzb.str().size()<1) mcjzb<<"("<<jzbvariablemc<<")";
471
472 } else {
473 write_warning(__FUNCTION__,"NOT DOING ANY PEAK FINDING ATM");
474 // datajzb << "((jzb[1]+0.059979*pt)- 4.18335 )";
475 // mcjzb << "((jzb[1]+0.034665*pt)- 3.58273 )";
476 // datajzb << "((jzb[1]+0.0595776*pt)- 4.17885 )";
477 // mcjzb << "((jzb[1]+0.0348755*pt)- 3.91542 )";
478 dout << "Not only not doing peak finding, but actually just copying global variable" << endl; datajzb << jzbvariabledata; mcjzb << jzbvariablemc;
479 }
480
481
482 if(do_pred || get_new_results) {
483
484 /* datajzb.str("");
485 mcjzb.str("");
486 jzbvariabledata="(((id1==0&&id1==id2)*(jzb[1]+0.0647272*pt))+((id1==1&&id1==id2)*(jzb[1]+0.053997*pt))+((id1!=id2)*(jzb[1]+0.0596296*pt)))";
487 jzbvariablemc="(((id1==0&&id1==id2)*(jzb[1]+0.0363773*pt))+((id1==1&&id1==id2)*(jzb[1]+0.0337406*pt))+((id1!=id2)*(jzb[1]+0.0350181*pt)))";
488
489 find_peaks(MCPeak,MCPeakError, DataPeak, DataPeakError,resultsummary,true,datajzb,mcjzb);
490
491 if(datajzb.str().size()<1) datajzb<<"("<<jzbvariabledata<<")";
492 if(mcjzb.str().size()<1) mcjzb<<"("<<jzbvariablemc<<")";
493 */
494 }
495
496
497 dout << "With peak correction, we get : " << endl;
498 dout << " Data : " << datajzb.str() << endl;
499 dout << " MC : " << mcjzb.str() << endl;
500
501 //**** part 2 : kinematic plots
502 if(do_kinematic_variables||(do_all&&!PlottingSetup::Approved)) do_kinematic_plots(mcjzb.str(),datajzb.str());
503 if(do_kinematic_PF_variables) do_kinematic_PF_plots(mcjzb.str(),datajzb.str());
504
505 //**** part 3: Leptonic comparison (ee vs mm, eemm vs emu)
506 if(do_lepton_comparison||(do_all&&!PlottingSetup::Approved)) lepton_comparison_plots();
507
508
509 //**** part 3b: comparison between control regions (SFZPJZBPOS,SFZPJZBNEG,...)
510 if (do_region_comparison||(do_all&&!PlottingSetup::Approved)) region_comparison_plots(mcjzb.str(),datajzb.str(),jzb_cut);
511
512 //**** part 4: JZB plots (OSOF, OSSF)
513 if(do_jzb_plots||do_all) jzb_plots(mcjzb.str(),datajzb.str(),global_ratio_binning);
514
515 //**** part 5 : Prediction plots
516 if(do_pred||do_all) do_prediction_plots(mcjzb.str(),datajzb.str(),DataSigma,MCSigma,overlay_signal);
517
518 //**** part 6: Ratio plots
519 //included in prediction plots.
520
521 //**** part 7: Some systematics plots
522 if(do_signal_bg_comparison_plot||(do_all&&!PlottingSetup::Approved)) signal_bg_comparison();
523 if(do_ttbar_comparison||(do_all&&!PlottingSetup::Approved)) ttbar_sidebands_comparison(mcjzb.str(),global_ratio_binning);
524 if(do_zjets_comparison||(do_all&&!PlottingSetup::Approved)) zjets_prediction_comparison(mcjzb.str());
525
526 if(do_diboson_plots) diboson_plots(mcjzb.str(),datajzb.str(),global_ratio_binning);
527
528 if(do_JES||(do_all&&!PlottingSetup::Approved)) {
529 make_JES_plot(cutmass&&cutOSSF&&basiccut,"_ossf");
530 make_JES_plot(cutmass&&cutOSOF&&basiccut,"_osof");
531 }
532
533 if (do_jzb_efficiency_curve) plot_jzb_sel_eff(mcjzb.str(),signalsamples,jzb_cut);
534
535 //**** part 8: observed and predicted!
536 if(do_all||get_new_results) doquick=0;
537 if(do_save_template) doquick=2;
538 // if(do_edge&&doquick>1) doquick=1;//we need at least doquick=1 (i.e. with storage) for edge fitting
539
540 if(do_all||do_save_template||do_compute_upper_limits_from_counting_experiment||do_model_scan||get_new_results) {
541
542 get_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,verbose,dopoisson,doquick);
543 // Also exclusive bins for low JZB
544 get_exclusive_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,0.,30.,verbose,dopoisson,doquick);
545 get_exclusive_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,30.,50.,verbose,dopoisson,doquick);
546 if(!PlottingSetup::RestrictToMassPeak) {
547 dout << " ***************************** " << endl;
548 dout << " Now also going to compute results in 20<mll<70 GeV range !" << endl;
549 TCut cutMassbkp=cutmass;
550 cutmass=TCut(cutmass&&TCut("mll>20&&mll<70"));
551 get_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,verbose,dopoisson,doquick);
552 cutmass=cutMassbkp;
553 }
554 }
555
556
557 vector<float>jzb_shape_limit_bins;
558 jzb_shape_limit_bins.push_back(50);
559 jzb_shape_limit_bins.push_back(65.05);
560 jzb_shape_limit_bins.push_back(73.4);
561 jzb_shape_limit_bins.push_back(88.7);
562 jzb_shape_limit_bins.push_back(115.25);
563 jzb_shape_limit_bins.push_back(200);
564 // jzb_shape_limit_bins.push_back(250);
565 if(do_prepare_limits_using_shapes) prepare_limits(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_shape_limit_bins);
566
567 /* if(do_edge) {
568 // write_info(__FUNCTION__,"Doing edge fit for data");
569 // DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,data,cutmass&&cutnJets);
570 write_info(__FUNCTION__,"Doing edge fit for MC");
571 DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,mc,cutmass&&cutnJets);
572 // write_info(__FUNCTION__,"Doing edge fit for MC with signal");
573 // DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,mcwithsignal,cutmass&&cutnJets);
574 }*/
575
576 //------------------------ end of standard functions
577
578 // if(calculate_yields||do_all) calculate_all_yields(mcjzb.str(),jzb_cut); // now outdated - get_results now replaces this.
579
580
581 vector<vector<float> > all_systematics;
582 vector<vector<float> > all_efficiency;
583
584 if(do_compute_systematics||do_compute_upper_limits_from_counting_experiment||do_all) all_systematics=compute_systematics(mcjzb.str(),MCPeakError,alwaysflip,datajzb.str(),signalsamples,jzb_cut,requireZ);
585 // if(do_compute_efficiency) compute_efficiency(mcjzb.str(),allsamples, 50, -200, 300);
586 if(do_compute_efficiency) write_warning(__FUNCTION__,"efficiency computation deactivated");
587
588 if(do_compute_upper_limits_from_counting_experiment||do_all) compute_upper_limits_from_counting_experiment(all_systematics,jzb_cut,mcjzb.str(),doobserved,alwaysflip);
589
590 //------------------------ end of analysis parts - below: special applications
591
592 if(do_model_scan) scan_SUSY_parameter_space(mcjzb.str(),datajzb.str(),jzb_cut,requireZ,MCPeakError,DataPeakError);
593
594 if(do_zjet_ttbar_shapes) draw_ttbar_and_zjets_shape(mcjzb.str(),datajzb.str());
595
596 if(do_response_correction) {
597 find_correction_factors(jzbvariabledata,jzbvariablemc);
598 dout << "Please update Modules/Setup.C to reflect the following values (round them first ... )" <<endl;
599 dout << "Corrected JZB variable definition " << endl;
600 cout << " Data: " << jzbvariabledata << endl;
601 cout << " MC: " << jzbvariablemc << endl;
602 dout << endl;
603 dout << "If you're feeling lazy, copy & paste this : " << endl;
604 dout << " string jzbvariabledata=\""<<jzbvariabledata<<"\";"<<endl;
605 dout << " string jzbvariablemc=\"" <<jzbvariablemc<<"\";"<<endl;
606 }
607
608 if(do_pick_up_events) {
609
610 TCut essentialcutbkp = essentialcut;
611 essentialcut=TCut("1.0");
612
613
614 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>20)))&&(((met[4]>100)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>20))))&&(pfJetGoodNum40>=3))&&mll>20&&mll<70&&(abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>20)&&pt1>20&&pt2>20)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","ETH_SF__LOOSE.txt");
615
616 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)))&&(((met[4]>150)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10))))&&(pfJetGoodNum40>=2&&pfTightHT>100))&&mll>20&&mll<70&&(abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)&&pt1>20&&pt2>10)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","Aachen_SF__LOOSE.txt");
617
618 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)))&&(((met[4]>100)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20))))&&(pfJetGoodNum40>=3))&&mll>20&&mll<70&&(abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)&&pt1>20&&pt2>20)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","ETH_SF__TIGHT.txt");
619
620 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>10)))&&(((met[4]>150)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>10))))&&(pfJetGoodNum40>=2&&pfTightHT>100))&&mll>20&&mll<70&&(abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>10)&&pt1>20&&pt2>10)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","Aachen_SF__TIGHT.txt");
621
622
623
624
625 /*
626 cout << "LOW MASS, SF, ETH " << endl;
627 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)))&&(((met[4]>100)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20))))&&(pfJetGoodNum40>=3))&&mll>20&&mll<70&&(abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)&&pt1>20&&pt2>20)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","ETH_SF.txt");
628
629 cout << "LOW MASS, OF, ETH " << endl;
630 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)))&&(((met[4]>100)&&((id1!=id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20))))&&(pfJetGoodNum40>=3))&&mll>20&&mll<70&&(abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>30 && pt2>20)&&pt1>20&&pt2>20)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","ETH_OF.txt");
631
632 cout << "LOW MASS, SF, Aachen " << endl;
633 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)))&&(((met[4]>150)&&((id1==id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10))))&&(pfJetGoodNum40>=2&&pfTightHT>100))&&mll>20&&mll<70&&(abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)&&pt1>20&&pt2>10)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","Aachen_SF.txt");
634
635 cout << "LOW MASS, OF, Aachen " << endl;
636 pick_up_events("((((((passed_triggers ||!is_data)))&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)))&&(((met[4]>150)&&((id1!=id2)&&(ch1*ch2<0)))&&((((pfJetGoodNum40>=2&&pfJetGoodID[0]!=0)&&(pfJetGoodNum40>=2&&pfJetGoodID[1]!=0))&&((mll>2)&&((abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10))))&&(pfJetGoodNum40>=2&&pfTightHT>100))&&mll>20&&mll<70&&(abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>30 && pt2>10)&&pt1>20&&pt2>10)))*((weight*(weight<1000)*(is_data+(!is_data)*((id1==id2)*0.95+(id1!=id2)*0.94))))","Aachen_OF.txt");
637 essentialcut=essentialcutbkp;
638 /*
639 dout << "Observed: " << endl;
640 pick_up_events((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
641 dout << "Predicted (JZB<-100) OSSF" << endl;
642 pick_up_events((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
643 dout << "Predicted (emu, JZB>100) OSOF" << endl;
644 pick_up_events((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
645 dout << "Predicted (emu, JZB<-100) OSOF" << endl;
646 pick_up_events((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
647 dout << "Predicted (SB SF, JZB>100) OSSF" << endl;
648 pick_up_events((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
649 dout << "Predicted (SB SF, JZB<-100) OSSF" << endl;
650 pick_up_events((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
651 dout << "Predicted (SB emu, JZB>100) OSOF" << endl;
652 pick_up_events((const char*)(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
653 dout << "Predicted (SB emu, JZB<-100) OSOF" << endl;
654 pick_up_events((const char*)(cutOSOF&&cutnJets&&basiccut&&sidebandcut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));*/
655 }
656
657 if(do_save_template||do_compute_upper_limits_from_shapes) {
658 // write_error(__FUNCTION__,"Currently commented out shape storage");
659 save_template(mcjzb.str(),datajzb.str(),jzb_cut,MCPeakError,DataPeakError,jzb_shape_limit_bins);
660 // PrepareDataShapes(mcjzb.str(), datajzb.str(), jzb_shape_limit_bins, DataPeakError);
661 // PrepareEdgeShapes(mcjzb.str(), datajzb.str(), jzb_shape_limit_bins, DataPeakError);
662 }
663
664 if(do_compute_upper_limits_from_shapes) compute_upper_limits_from_shapes(mcjzb.str(),datajzb.str(),jzb_cut,requireZ, MCPeakError, DataPeakError);
665
666 if(do_sideband_assessment) sideband_assessment(datajzb.str(),30.0,50.0);
667
668 if(do_mc_yields) {
669 TCut MassCuts[3] = {TCut("mll>20"),Restrmasscut,TCut("mll>20&&mll<70")};
670 string sMassCut[3] = {"NoMassCut" ,"ZWindow" , "LowMassRegion" };
671 for(int i=0;i<3;i++) {
672 cutmass=MassCuts[i];
673 cout << "_____________________________________________" << endl;
674 cout << "Considering mass region : " << sMassCut[i] << endl;
675 compute_Improved_MC_yields(mcjzb.str(),jzb_cut);
676 }
677 }
678
679 if(do_metjzb) met_jzb_cut(datajzb.str(),mcjzb.str(),jzb_cut);
680
681 if(do_qcd||do_all) qcd_plots(datajzb.str(),mcjzb.str(),jzb_cut);
682
683 if(do_mlls) do_mlls_plot(mcjzb.str());
684
685 if(do_ptsanity) check_ptsanity();
686
687 if(do_metpred) ExperimentalMetPrediction(false,true);
688
689 if(do_met_vs_jzb) met_vs_jzb_plots(datajzb.str(),mcjzb.str());
690
691 if(do_aachen) {
692 if(Contains(AachenFileName,"Data")) CastETHtoAachenNtuple(datajzb.str(),AachenFileName);
693 else CastETHtoAachenNtuple(mcjzb.str(),AachenFileName);
694 }
695
696 if(do_metplots) {
697 DoMetPlots(datajzb.str(),mcjzb.str());
698 }
699
700 if(do_wz) WZstudy();
701
702 if(do_rmue) ComputeRMuE();
703
704 // if(do_fsr) FSRstudy();
705
706 if(do_ttbartest) ttbar_tests(mcjzb.str());
707
708 if(do_test) test();
709
710 if(do_ttbaranalysis) ttbar_analysis();
711
712 if(do_overlap) overlapanalysis(datajzb.str(),mcjzb.str());
713
714 if(do_ittbar) make_ijzb_ttbar_closure_test(datajzb.str(),mcjzb.str());
715 return 0;
716 }
717