Development/Plotting/iJZB_Analysis.C

#include <iostream>
#include <vector>
#include <sys/stat.h>
#include <getopt.h>
#include <stdio.h>
#include <stdlib.h>

#ifndef Verbosity
#define Verbosity 0
#endif
#ifndef HUSH
#define HUSH 1
#endif

#include "Modules/GeneralToolBox.C"
#include "Modules/SampleClass.C"
#include "Modules/setTDRStyle.C"
#include "Modules/Setup.C"
#include "Modules/Poisson_Calculator.C"
#include "Modules/JSON/JSONSampleLoader.C"
#include "Modules/ActiveSamples.C"
#include "Modules/PeakFinder.C"
#include "Modules/UpperLimitsWithShape.C"
#include "Modules/Plotting_Functions.C"
#include "Modules/LimitCalculation.C"
#include "Modules/ResultModule.C"
#include "Modules/CrossSectionReader.C"
#include "Modules/Systematics.C"
#include "Modules/SugarCoating.C"
#include "Modules/ExclusionPlot.C"
#include "Modules/SUSYScan.C"
#include "Modules/AachenCompatibility.C"
#include "Modules/MetPlotting.C"
#include "Modules/WZStudy.C"
//#include "Modules/FSRStudy.C"
#include "Modules/TTbarAnalysis.C"
#include "Modules/OverlapStudy.C"


#include <TCut.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TColor.h>
#include <TPaveText.h>
#include <TRandom.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TSQLResult.h>

using namespace PlottingSetup;

void usage(int passed=0 ) {
  std::cout << "USAGE : " << std::endl;
  std::cout << "You can use different options when running this program : " << std::endl;
  std::cout << std::endl;
  std::cout << "\033[1;34m all\033[0m  \t\t  All processes of the standard workflow" << std::endl;
  std::cout << "\033[1;34m peak\033[0m  \t\t  Peak finding" << std::endl;
  std::cout << "\033[1;34m pileup\033[0m\t\t  Pileup correction determination" << std::endl;
  std::cout << "\033[1;34m kin\033[0m  \t\t  Kinematic plots" << std::endl;
  std::cout << "\033[1;34m pfkin\033[0m  \t  Kinematic plots for PF variables" << std::endl;
  std::cout << "\033[1;34m region\033[0m  \t  Region comparison plots" << std::endl;
  std::cout << "\033[1;34m lepton\033[0m  \t  Lepton comparison plots" << std::endl;
  std::cout << "\033[1;34m jzb\033[0m  \t\t  JZB plots" << std::endl;
  std::cout << "\033[1;34m pred\033[0m  \t\t  Prediction plots" << std::endl;
  std::cout << "\033[1;34m signalbg\033[0m  \t  Signal vs. background plots (shape comparison)" << std::endl;
  std::cout << "\033[1;34m ttbar\033[0m  \t  TTbar sideband comparison" << std::endl;
  std::cout << "\033[1;34m ttbaranalysis\033[0m  \t  TTbar analysis. To prove we can." << std::endl;
  std::cout << "\033[1;34m ttbartest\033[0m  \t  TTbar prediction tests (peak, response)" << std::endl;
  std::cout << "\033[1;34m zjets\033[0m  \t  ZJets plots" << std::endl;
  std::cout << "\033[1;34m diboson\033[0m  \t  Diboson plots" << std::endl;
  std::cout << "\033[1;34m results\033[0m  \t  Compute results (JZB>X)" << std::endl;
  std::cout << "\033[1;34m syst\033[0m  \t\t  Compute systematics" << std::endl;
  std::cout << "\033[1;34m eff\033[0m  \t\t  Compute all (partial) efficiencies" << std::endl;
  std::cout << "\033[1;34m shapes\033[0m  \t  Prepare ROOT files for computing limits using shapes" << std::endl;
  std::cout << "\033[1;34m counting\033[0m  \t  Establish limits with counting experiment" << std::endl;
  std::cout << "\033[1;34m shapelimit\033[0m  \t  Establish limits with shapes" << std::endl;
  std::cout << "\033[1;34m qcd\033[0m     \t\t  Show that QCD is negligible" << endl;
  std::cout << "\033[1;34m response\033[0m  \t  Compute response correction" << std::endl;
  std::cout << "\033[1;34m pickup\033[0m  \t  Pick up interesting events" << std::endl;
  std::cout << "\033[1;34m yields\033[0m  \t  Compute MC yields for this lumi" << std::endl;
  std::cout << "\033[1;34m sideassess\033[0m\t  Assessment of the mass sidebands definition in the data" << std::endl;
  std::cout << "\033[1;34m jes\033[0m \t\t  Make Jet Energy Scale (JES) plot" << std::endl;
  std::cout << "\033[1;34m nosidebands\033[0m \t\t  Switch off sidebands for cJZB" << std::endl;
  std::cout << "\033[1;34m effcurve\033[0m \t\t  Make JZB efficiency plot" << std::endl;
  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;
  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;
  std::cout << "\033[1;34m offpeak\033[0m   \t  Sets offpeak to true (i.e. carry out off-peak analysis!)" << std::endl;
  std::cout << "\033[1;34m onpeak\033[0m   \t  Sets onpeak to true (i.e. carry out standard on-peak analysis!)" << std::endl;
  std::cout << "\033[1;34m btagged (b)\033[0m  \t  Sets DoBTag to true (i.e. carry out b tagged analysis)" << std::endl;
  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;
  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;
  std::cout << "\033[1;34m mlls\033[0m  \t  Draw dilepton invariant mass for all signals" << endl;
  std::cout << "\033[1;34m rmue\033[0m  \t  Compute rmu" << std::endl;
  std::cout << "\033[1;34m overlap\033[0m  \t  Study overlap between iJZB and METJZB" << endl;
//  std::cout << "\033[1;34m edge\033[0m  \t  Be on the edge of development and try out the new edge fitter :-)" << endl;
  std::cout << "\033[1;34m wz  \033[0m  \t  WZ" << std::endl;
  std::cout << "\033[1;34m ittbar\033[0m\t  ttbar closure test for iJZB" << std::endl;
//  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;
  std::cout << std::endl;
  std::cout << "\033[1;34m ijzb   \033[om\t  This is iJZB. Rerun peak positions and everything" << std::endl;
  std::cout << "\033[1;34m dir (d)\033[0m\t  Directory where all plots will be saved" << std::endl;
  std::cout << "\033[1;34m public \033[0m\t  Switch to be activated when using the algorithm in public (no branding etc.)" << std::endl;
  std::cout << std::endl;
  std::cout << std::endl;
  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;
  std::cout << "\033[1;34m metplots\033[0m \t\t  Produce the MET plots (i.e. mll plots with MET cut etc.)" << std::endl;
  std::cout << "\033[1;34m metvsjzb\033[0m \t\t Produce MET:JZB plots" << std::endl;
  std::cout << "\033[1;34m metpred\033[0m \t\t  Experimental JZB prediction for MET signal region" << std::endl;
  std::cout << "\033[1;34m aachen\033[0m  \t\t  Produce Aachen root file (please provide source root file)" << std::endl;
  std::cout << std::endl;
  std::cout << std::endl;
  std::cout << "\033[1;34m paper\033[0m \t\t  Switch paper mode on (no \" Preliminary \" on plots)" << std::endl;
  std::cout << "\033[1;34m approved\033[0m \t\t  Only approved plots" << std::endl;
  exit(-1);
}

int main (int argc, char ** argv)
{
  
  
  PlottingSetup::IsIJZBAnalysis=true;
  int do_all           = false;/// DONE
  int do_global_significance       = false; /// DONE
  int do_peak_finding              = false; /// DONE
  int do_pileup                    = false;
  int do_kinematic_variables       = false; ///DONE
  int do_kinematic_PF_variables    = false; ///DONE
  int do_region_comparison         = false; /// Comparison between signal and control regions
  int do_lepton_comparison         = false; ///DONE
  int do_jzb_plots                 = false; /// DONE
  int do_pred                      = false; /// DONE
        bool overlay_signal               = false; /// Overlay LM4 signal to predictions
  int do_signal_bg_comparison_plot = false; /// DONE
  int do_ttbar_comparison          = false; /// DONE
  int do_ttbaranalysis             = false;
  int do_overlap                   = false;
  int do_ttbartest                 = false; /// DONE
  int do_jzb_efficiency_curve      = false; /// DONE
  int do_zjets_comparison          = false; /// DONE
  int do_diboson_plots             = false;
//  int calculate_pred_and_observed  = false; /// DONE -- now rerouted to the new way of computing results
//  int calculate_yields             = false; ///DONE /// superfluous (replaced by get_new_results)
  int get_new_results              = false; ///DONE
        bool dopoisson              = false;//should we calculate stat err with poisson? ///DONE
        bool verbose                = false;//prints out a lot more information ... ///DONE
  //--------------------------------------------
  // Systematics and limits
  int do_compute_systematics=false; /// DONE (might want to add something related to the peak?)
        bool requireZ=true; // should we require a Z for MC efficiency?
  int do_compute_efficiency=false; /// DONE (might want to add something related to the peak?)
  int do_prepare_limits_using_shapes=false; /// ALMOST DONE **************************************
  int do_compute_upper_limits_from_counting_experiment=false; /// DONE
        bool doobserved=false; // want to get the observed limit as well?
        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)
  int do_compute_upper_limits_from_shapes=false;
  
  //--------------------------------------------
  // More not-so-standard stuff
  int do_model_scan=false; // this can take days - use the grid version!
  int do_zjet_ttbar_shapes=false; /// DONE
  int do_test=false;/// DONE    (just any test you wish)
  int do_pick_up_events=false; ///DONE
  
  //--------------------------------------------
  // To be run manually !
  int do_response_correction=false; ///DONE 
               //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
  
  int do_sideband_assessment=false;
  int do_qcd=false;

  int do_save_template=false;
  int do_mc_yields=false; // compute MC yields
  int do_metjzb=false;
  
  int savepdf=true;
  int saveC=true;
  int saveRoot=true;
  int savepng=true;
  int saveeps=false;

  int dooffpeak=true;
  int doonpeak=false;
  int dobtag=false;

  int do_JES=false;
  int do_ptsanity=false;
  int do_mlls=false;
  int do_rmue=false;
//  int do_edge=false;
  int do_wz=false;
  int do_ittbar=false;
  int do_ijzb=false;
  int do_fsr=false;
  int do_met_vs_jzb=false;
  int do_aachen=false;
  int do_metplots=false;
  int do_metpred=false;
  string AachenFileName="";

  std::string directory="";
  int option_iterator;
  int option_counter=0;
  bool moreoptions=true;
  
  int SwitchOffSidebands=0;

  string jzbcuts_string="";
  float local_signif=0.0;

  while(moreoptions) {
    static struct option long_options[] =
         {
           /* These options set a flag. */
           {"all",      no_argument,       &do_all, 1},
           {"peak",     no_argument,       &do_peak_finding, 1},
               {"pileup",   no_argument,       &do_pileup, 1},
               {"test",      no_argument,       &do_test, 1},
           {"kin",      no_argument,       &do_kinematic_variables, 1},
           {"pfkin",    no_argument,       &do_kinematic_PF_variables, 1},
           {"region",   no_argument,       &do_region_comparison, 1},
           {"lepton",   no_argument,       &do_lepton_comparison, 1},
           {"jzb",      no_argument,       &do_jzb_plots, 1},
           {"pred",     no_argument,       &do_pred, 1},
           {"signalbg", no_argument,       &do_signal_bg_comparison_plot, 1},
           {"ttbar",    no_argument,       &do_ttbar_comparison, 1},
           {"ttbaranalysis",    no_argument,       &do_ttbaranalysis, 1},
           {"ttbartest",no_argument,       &do_ttbartest, 1},
               {"overlap"  ,no_argument,       &do_overlap, 1},
           {"zjets",    no_argument,       &do_zjets_comparison, 1},
               {"diboson",  no_argument,       &do_diboson_plots, 1},
           {"results",  no_argument,       &get_new_results, 1},
           {"syst",     no_argument,       &do_compute_systematics, 1},
           {"eff",      no_argument,       &do_compute_efficiency, 1},
           {"shapes",   no_argument,       &do_prepare_limits_using_shapes, 1},
           {"counting", no_argument,       &do_compute_upper_limits_from_counting_experiment, 1},
           {"shapelimit",no_argument,      &do_compute_upper_limits_from_shapes, 1},
           {"ptsanity", no_argument,       &do_ptsanity, 1},
           {"response", no_argument,       &do_response_correction, 1},
           {"pickup",   no_argument,       &do_pick_up_events, 1},
           {"yields",   no_argument,       &do_mc_yields, 1},
           {"sideassess",no_argument,      &do_sideband_assessment, 1},
               {"metvsjzb", no_argument,       &do_met_vs_jzb, 1},
           {"qcd",      no_argument,       &do_qcd, 1},
           {"save",     no_argument,       &do_save_template,1},
           {"png",      no_argument,       &savepng,1},
           {"eps",      no_argument,       &saveeps,1},
           {"pdf",      no_argument,       &savepdf,1},
           {"root",     no_argument,       &saveRoot,1},
           {"C",        no_argument,       &saveC,1},
           {"offpeak",  no_argument,       &dooffpeak,1},
           {"onpeak",   no_argument,       &doonpeak,1},
           {"btagged",  no_argument,       &dobtag,1},
           {"met",      no_argument,       &do_metjzb,1},
           {"jes",      no_argument,       &do_JES,1},
           {"effcurve", no_argument,       &do_jzb_efficiency_curve,1},
           {"public",   no_argument,       &PlottingSetup::publicmode,1},
           {"paper",    no_argument,       &PlottingSetup::PaperMode,1},
           {"approved", no_argument,       &PlottingSetup::Approved,1},
           {"nosidebands",no_argument,     &SwitchOffSidebands,1},
           {"mlls",     no_argument,       &do_mlls,1},
           {"rmue",     no_argument,       &do_rmue,1},
//         {"edge",     no_argument,       &do_edge,1},
           {"wz",       no_argument,       &do_wz,1},
           {"ijzb",     no_argument,       &do_ijzb,1},
           {"ittbar",   no_argument,       &do_ittbar,1},
           {"fsr",      no_argument,       &do_fsr,1},
           {"metplots", no_argument,       &do_metplots,1},
           {"metpred",  no_argument,       &do_metpred,1},
           /* The following options store values.*/
           {"jzbcuts", required_argument, 0, 'j'},
           {"dir",     required_argument, 0, 'd'},
               {"aachen",  required_argument, 0, 'c'},
           {"gsig",    required_argument,  0, 'g'},
               {0, 0, 0, 0}
         };
    int option_index = 0;
    option_iterator = getopt_long(argc, argv, "abc:d:g:j:",long_options, &option_index);
    if(option_iterator == -1) moreoptions=false;
    else {
      option_counter++;
      switch (option_iterator)
             {
             case 0:
               if (long_options[option_index].flag != 0)
                 break;
               printf ("option %s", long_options[option_index].name);
               if (optarg)
                 printf (" with arg %s", optarg);
               printf ("\n");
               break;
             case 'a':
               do_all=true;
               break;
             case 'b':
               dobtag=true;
               break;
             case 'c':
               do_aachen=true;
               AachenFileName=(std::string)optarg;
               cout << "Received file to convert to Aachen format: " << AachenFileName << endl;
               break;
             case 'd':
               directory=(std::string)optarg;
               std::cout<<"Option directory was passed with argument " << optarg << std::endl;
               break;
             case 'j':
               jzbcuts_string=(std::string)optarg;
               std::cout<<"JZB cuts were manually defined:" << optarg << std::endl;
               break;
             case 'g':
                     do_global_significance=true;
                     local_signif=atof(optarg);
                     break;
             case '?':
               usage(option_iterator);
               break;
             default:
               usage(option_iterator);
             }
    }
  }
  
  if(directory!="") PlottingSetup::directoryname=directory;
  if(dooffpeak) PlottingSetup::RestrictToMassPeak=false;
  if(doonpeak) PlottingSetup::RestrictToMassPeak=true;
  if(dobtag) PlottingSetup::DoBTag=true;
  if(option_counter==0) usage();
  if(SwitchOffSidebands) PlottingSetup::UseSidebandsForcJZB=false;
  
  ///----------------------------------- BELOW THIS LINE: NO MORE OPTIONS BUT ACTUAL FUNCTION CALLS! ---------------------------------------------------------
  gROOT->SetStyle("Plain");
  bool do_fat_line=false; // if you want to have HistLineWidth=1 and FuncWidth=1 as it was before instead of 2
  setTDRStyle(do_fat_line);
  gStyle->SetTextFont(42);
  bool showList=true;
  set_directory(directoryname);//Indicate the directory name where you'd like to save the output files in Setup.C

  if(do_wz) {
    IsWZAnalysis=true;
    IsJZBAnalysis=false;
  }

  
  write_info(__FUNCTION__,"Making sure that we are only considering category A !!");
  
  string backup_basicqualitycut = (const char*) basicqualitycut;
  string backup_essentialcut    = (const char*) essentialcut;
  string backup_basiccut        = (const char*) basiccut;
  string backup_leptoncut       = (const char*) leptoncut;
  string backup_essential       = (const char*) essential;
  
  string Sbasicqualitycut = backup_basicqualitycut;
  Sbasicqualitycut = ReplaceAll(Sbasicqualitycut,")<2.4",")<1.4");
  basicqualitycut=TCut(Sbasicqualitycut.c_str());
  
  string Sbasiccut = backup_basiccut;
  Sbasiccut = ReplaceAll(Sbasiccut,")<2.4",")<1.4");
  basiccut=TCut(Sbasiccut.c_str());
  
  string Sessentialcut = backup_essentialcut;
  Sessentialcut = ReplaceAll(Sessentialcut,")<2.4",")<1.4");
  essentialcut=TCut(Sessentialcut.c_str());
  
  string Sleptoncut = backup_leptoncut;
  Sleptoncut = ReplaceAll(Sleptoncut,")<2.4",")<1.4");
  leptoncut=TCut(Sleptoncut.c_str());
  
  string Sessential = backup_essential;
  Sessential        = ReplaceAll(Sessential,")<2.4",")<1.4");
  essential         = TCut(Sessential.c_str());
  
  
  if(do_global_significance) {
        compute_global_significance(local_signif);
        return 0;
  }

  define_samples(showList,allsamples,signalsamples,scansample,raresample,systsamples,qcdsamples,comparesamples);
  setlumi(luminosity);
  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!
  stringstream resultsummary;

  write_analysis_type(PlottingSetup::RestrictToMassPeak,PlottingSetup::DoBTag);
  do_png(savepng);
  do_pdf(savepdf);
  do_eps(saveeps);
  do_C(saveC);
  do_root(saveRoot);

  write_info(__FUNCTION__,"Imposing standard JZB definition now");
  /*
  jzbvariabledata="((jzb[1]+0.0595776*pt-0.279212*numVtx)- 0.375633 )";
  jzbvariablemc="((jzb[1]+0.0348755*pt-0.279212*numVtx)- 0.224927 )";
*/
/*
   jzbvariabledata="(jzb[1]+0.0595776*pt)";
   jzbvariablemc="(jzb[1]+0.0348755*pt)";
   */
   jzbvariabledata="(jzb[1]+0.0595776*pt- 4.18673)";
   jzbvariablemc="(jzb[1]+0.0348755*pt- 3.8028)";
   
/*   write_warning(__FUNCTION__,"WATCH OUT, EXPERIMENTAL PT CORRECTION");
   jzbvariabledata="jzb[1]+(1-(9.57810e-01 + -3.53527e+00/pt + 1.85206e+02/pow(pt,2)))*pt";
   jzbvariablemc  ="jzb[1]+(1-(9.60547e-01 + -1.16558e+00/pt + 2.12197e+02/pow(pt,2)))*pt";
   
   jzbvariabledata="(jzb[1]+(1-(9.57810e-01 + -3.53527e+00/pt + 1.85206e+02/pow(pt,2)))*pt- 1.92518 )";
   jzbvariablemc  ="(jzb[1]+(1-(9.60547e-01 + -1.16558e+00/pt + 2.12197e+02/pow(pt,2)))*pt+ 0.0553847 )";
*/
   
/*write_warning(__FUNCTION__,"No response correction");
   jzbvariabledata="(jzb[1]+ 0.91084)";
   jzbvariablemc="(jzb[1]- 0.694762)";*/


  global_ratio_binning.push_back(0);
  global_ratio_binning.push_back(5);
  global_ratio_binning.push_back(10);
  global_ratio_binning.push_back(20);
  global_ratio_binning.push_back(50);
  global_ratio_binning.push_back(100);
  global_ratio_binning.push_back(150);
  global_ratio_binning.push_back(200);
  global_ratio_binning.push_back(PlottingSetup::jzbHigh);
  //global_ratio_binning.push_back(500);

  //these are the JZB cuts defining our search regions
  vector<float> jzb_cut; //starting where, please?
  if(jzbcuts_string=="") {
//      jzb_cut.push_back(50);
        //  jzb_cut.push_back(75);
        jzb_cut.push_back(80);
        jzb_cut.push_back(100);
        //  jzb_cut.push_back(125);
        jzb_cut.push_back(150);
        //  jzb_cut.push_back(175);
//      jzb_cut.push_back(200); 
        //  jzb_cut.push_back(225);
//      jzb_cut.push_back(250);
        //  jzb_cut.push_back(275);
//      jzb_cut.push_back(300);
  } else {
        manually_set_jzb_cuts(jzb_cut,jzbcuts_string);
        if(jzb_cut.size()==0) {
                write_error(__FUNCTION__,"There are no jzb cuts!");
                return -1;
        }
  }

  todo();
  

  //**** part 1 : peak finding
  float MCPeak=0,MCPeakError=0,DataPeak=0,DataPeakError=0,MCSigma=10,DataSigma=10;
  method=Kostasmethod;//Kostasmethod;//dogaus3sigma;// options: dogaus,doKM,dogaus2sigma,dogaus3sigma
  stringstream datajzb;
  stringstream mcjzb;

  
  if(do_pileup) FindPeakPileUpCorrection();
  
  if ( do_peak_finding ) {
    //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);
  find_peaks(MCPeak,MCPeakError, DataPeak, DataPeakError,resultsummary,true,datajzb,mcjzb);

    if(datajzb.str().size()<1) datajzb<<"("<<jzbvariabledata<<")";
    if(mcjzb.str().size()<1) mcjzb<<"("<<jzbvariablemc<<")";
    
  } else {
    write_warning(__FUNCTION__,"NOT DOING ANY PEAK FINDING ATM");
//    datajzb << "((jzb[1]+0.059979*pt)- 4.18335 )";
//    mcjzb << "((jzb[1]+0.034665*pt)- 3.58273 )";
//    datajzb << "((jzb[1]+0.0595776*pt)- 4.17885 )";
//    mcjzb   << "((jzb[1]+0.0348755*pt)- 3.91542 )";
    dout << "Not only not doing peak finding, but actually just copying global variable" << endl; datajzb << jzbvariabledata; mcjzb << jzbvariablemc; 
  }

  
  if(do_pred || get_new_results) {
    
/*    datajzb.str("");
    mcjzb.str("");
    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)))";
    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)))";
    
    find_peaks(MCPeak,MCPeakError, DataPeak, DataPeakError,resultsummary,true,datajzb,mcjzb);

    if(datajzb.str().size()<1) datajzb<<"("<<jzbvariabledata<<")";
    if(mcjzb.str().size()<1) mcjzb<<"("<<jzbvariablemc<<")";
  */  
  }
  
  
  dout << "With peak correction, we get : " << endl;
  dout << "    Data : " << datajzb.str() << endl;
  dout << "    MC : " << mcjzb.str() << endl;

  //**** part 2 : kinematic plots
   if(do_kinematic_variables||(do_all&&!PlottingSetup::Approved)) do_kinematic_plots(mcjzb.str(),datajzb.str());
  if(do_kinematic_PF_variables) do_kinematic_PF_plots(mcjzb.str(),datajzb.str());
  
  //**** part 3: Leptonic comparison (ee vs mm, eemm vs emu)
  if(do_lepton_comparison||(do_all&&!PlottingSetup::Approved)) lepton_comparison_plots();


  //**** part 3b: comparison between control regions (SFZPJZBPOS,SFZPJZBNEG,...)
  if (do_region_comparison||(do_all&&!PlottingSetup::Approved)) region_comparison_plots(mcjzb.str(),datajzb.str(),jzb_cut);

  //**** part 4: JZB plots (OSOF, OSSF) 
  if(do_jzb_plots||do_all) jzb_plots(mcjzb.str(),datajzb.str(),global_ratio_binning);

  //**** part 5 : Prediction plots
  if(do_pred||do_all) do_prediction_plots(mcjzb.str(),datajzb.str(),DataSigma,MCSigma,overlay_signal);
  
  //**** part 6: Ratio plots
   //included in prediction plots.
  
  //**** part 7: Some systematics plots 
  if(do_signal_bg_comparison_plot||(do_all&&!PlottingSetup::Approved)) signal_bg_comparison();
  if(do_ttbar_comparison||(do_all&&!PlottingSetup::Approved)) ttbar_sidebands_comparison(mcjzb.str(),global_ratio_binning);
  if(do_zjets_comparison||(do_all&&!PlottingSetup::Approved)) zjets_prediction_comparison(mcjzb.str());
  
  if(do_diboson_plots) diboson_plots(mcjzb.str(),datajzb.str(),global_ratio_binning);

  if(do_JES||(do_all&&!PlottingSetup::Approved)) {
    make_JES_plot(cutmass&&cutOSSF&&basiccut,"_ossf");
    make_JES_plot(cutmass&&cutOSOF&&basiccut,"_osof");
  }

  if (do_jzb_efficiency_curve) plot_jzb_sel_eff(mcjzb.str(),signalsamples,jzb_cut);
  
  //**** part 8: observed and predicted!
  if(do_all||get_new_results) doquick=0;
  if(do_save_template) doquick=2;
//  if(do_edge&&doquick>1) doquick=1;//we need at least doquick=1 (i.e. with storage) for edge fitting

  if(do_all||do_save_template||do_compute_upper_limits_from_counting_experiment||do_model_scan||get_new_results) {
    
     get_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,verbose,dopoisson,doquick);
     // Also exclusive bins for low JZB
     get_exclusive_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,0.,30.,verbose,dopoisson,doquick);
     get_exclusive_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,30.,50.,verbose,dopoisson,doquick);
    if(!PlottingSetup::RestrictToMassPeak) {
      dout << " ***************************** " << endl;
      dout << "   Now also going to compute results in 20<mll<70 GeV range !" << endl;
      TCut cutMassbkp=cutmass;
      cutmass=TCut(cutmass&&TCut("mll>20&&mll<70"));
      get_result(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,verbose,dopoisson,doquick);
      cutmass=cutMassbkp;
    }
  }

  
  vector<float>jzb_shape_limit_bins;
  jzb_shape_limit_bins.push_back(50);
  jzb_shape_limit_bins.push_back(65.05);
  jzb_shape_limit_bins.push_back(73.4);
  jzb_shape_limit_bins.push_back(88.7);
  jzb_shape_limit_bins.push_back(115.25);
  jzb_shape_limit_bins.push_back(200);
//  jzb_shape_limit_bins.push_back(250);
  if(do_prepare_limits_using_shapes) prepare_limits(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_shape_limit_bins);
  
/*  if(do_edge) {
//    write_info(__FUNCTION__,"Doing edge fit for data");
//    DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,data,cutmass&&cutnJets);
    write_info(__FUNCTION__,"Doing edge fit for MC");
    DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,mc,cutmass&&cutnJets);
//    write_info(__FUNCTION__,"Doing edge fit for MC with signal");
//    DoEdgeFit(mcjzb.str(),datajzb.str(),DataPeakError,MCPeakError,jzb_cut,mcwithsignal,cutmass&&cutnJets);
  }*/
  
  //------------------------ end of standard functions
  
//  if(calculate_yields||do_all) calculate_all_yields(mcjzb.str(),jzb_cut); // now outdated - get_results now replaces this.
  
  
  vector<vector<float> > all_systematics;
  vector<vector<float> > all_efficiency;

  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);
//  if(do_compute_efficiency) compute_efficiency(mcjzb.str(),allsamples, 50, -200, 300);
  if(do_compute_efficiency) write_warning(__FUNCTION__,"efficiency computation deactivated");

  if(do_compute_upper_limits_from_counting_experiment||do_all) compute_upper_limits_from_counting_experiment(all_systematics,jzb_cut,mcjzb.str(),doobserved,alwaysflip);  

  //------------------------ end of analysis parts - below: special applications
  
  if(do_model_scan) scan_SUSY_parameter_space(mcjzb.str(),datajzb.str(),jzb_cut,requireZ,MCPeakError,DataPeakError);
  
  if(do_zjet_ttbar_shapes) draw_ttbar_and_zjets_shape(mcjzb.str(),datajzb.str());
  
  if(do_response_correction) {
    find_correction_factors(jzbvariabledata,jzbvariablemc);
    dout << "Please update Modules/Setup.C to reflect the following values (round them first ... )" <<endl;
    dout << "Corrected JZB variable definition " << endl;
    cout << "  Data: " << jzbvariabledata << endl;
    cout << "  MC:   " << jzbvariablemc << endl;
    dout << endl;
    dout << "If you're feeling lazy, copy & paste this : " << endl;
    dout << "  string jzbvariabledata=\""<<jzbvariabledata<<"\";"<<endl;
    dout << "  string jzbvariablemc=\""  <<jzbvariablemc<<"\";"<<endl;
  }
  
   if(do_pick_up_events) {
     
     TCut essentialcutbkp = essentialcut;
     essentialcut=TCut("1.0");
    
     
        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");
   
   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");
   
   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");
   
   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");
   

     /*
   cout << "LOW MASS, SF, ETH " << endl;
   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");
   
   cout << "LOW MASS, OF, ETH " << endl;
   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");
   
   cout << "LOW MASS, SF, Aachen " << endl;
   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");
   
   cout << "LOW MASS, OF, Aachen " << endl;
   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");
   essentialcut=essentialcutbkp;
   /*     
    dout << "Observed: " << endl;
    pick_up_events((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
    dout << "Predicted (JZB<-100) OSSF" << endl;
    pick_up_events((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
    dout << "Predicted (emu, JZB>100) OSOF" << endl;
    pick_up_events((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
    dout << "Predicted (emu, JZB<-100) OSOF" << endl;
    pick_up_events((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
    dout << "Predicted (SB SF, JZB>100) OSSF" << endl;
    pick_up_events((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
    dout << "Predicted (SB SF, JZB<-100) OSSF" << endl;
    pick_up_events((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));
    dout << "Predicted (SB emu, JZB>100) OSOF" << endl;
    pick_up_events((const char*)(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&((string)"(("+datajzb.str()+")"+">100)").c_str()));
    dout << "Predicted (SB emu, JZB<-100) OSOF" << endl;
    pick_up_events((const char*)(cutOSOF&&cutnJets&&basiccut&&sidebandcut&&((string)"(("+datajzb.str()+")"+"<-100)").c_str()));*/
  }

  if(do_save_template||do_compute_upper_limits_from_shapes) {
//    write_error(__FUNCTION__,"Currently commented out shape storage");
    save_template(mcjzb.str(),datajzb.str(),jzb_cut,MCPeakError,DataPeakError,jzb_shape_limit_bins);
//    PrepareDataShapes(mcjzb.str(), datajzb.str(), jzb_shape_limit_bins, DataPeakError);
//    PrepareEdgeShapes(mcjzb.str(), datajzb.str(), jzb_shape_limit_bins, DataPeakError);
  }
  
  if(do_compute_upper_limits_from_shapes) compute_upper_limits_from_shapes(mcjzb.str(),datajzb.str(),jzb_cut,requireZ, MCPeakError, DataPeakError);

  if(do_sideband_assessment) sideband_assessment(datajzb.str(),30.0,50.0);
  
  if(do_mc_yields) {
    TCut MassCuts[3]   = {TCut("mll>20"),Restrmasscut,TCut("mll>20&&mll<70")};
    string sMassCut[3] = {"NoMassCut"   ,"ZWindow"   , "LowMassRegion"      };
    for(int i=0;i<3;i++) {
      cutmass=MassCuts[i];
      cout << "_____________________________________________" << endl;
      cout << "Considering mass region : " << sMassCut[i] << endl;
      compute_Improved_MC_yields(mcjzb.str(),jzb_cut);
    }
  }
  
  if(do_metjzb) met_jzb_cut(datajzb.str(),mcjzb.str(),jzb_cut);
  
  if(do_qcd||do_all) qcd_plots(datajzb.str(),mcjzb.str(),jzb_cut);
  
  if(do_mlls) do_mlls_plot(mcjzb.str());
  
  if(do_ptsanity) check_ptsanity();
    
  if(do_metpred) ExperimentalMetPrediction(false,true);
  
  if(do_met_vs_jzb) met_vs_jzb_plots(datajzb.str(),mcjzb.str());
  
  if(do_aachen) {
    if(Contains(AachenFileName,"Data")) CastETHtoAachenNtuple(datajzb.str(),AachenFileName);
    else CastETHtoAachenNtuple(mcjzb.str(),AachenFileName);
  }
  
  if(do_metplots) {
    DoMetPlots(datajzb.str(),mcjzb.str());
  }
  
  if(do_wz) WZstudy();
    
  if(do_rmue) ComputeRMuE();
  
//  if(do_fsr) FSRstudy();
  
  if(do_ttbartest) ttbar_tests(mcjzb.str());
  
  if(do_test) test();
  
  if(do_ttbaranalysis) ttbar_analysis();
  
  if(do_overlap) overlapanalysis(datajzb.str(),mcjzb.str());
  
  if(do_ittbar) make_ijzb_ttbar_closure_test(datajzb.str(),mcjzb.str());
  return 0;  
}

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
Log Message:	Added calculation of global significance