Development/Plotting/Selective_Plot_Generator.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 <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 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" << 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" << 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 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)
{
  int do_all           = false;/// DONE
  int do_peak_finding              = false; /// DONE
  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_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_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="";

  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},
           {"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},
           {"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},
           {"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'},
               {0, 0, 0, 0}
         };
    int option_index = 0;
    option_iterator = getopt_long(argc, argv, "abc:d: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 '?':
               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
  set_treename("events");//you can set the treename here to be used; options are "events" (for reco) for "PFevents" (for particle flow)
//  define_samples(showList,allsamples);

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

    
  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);

  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;
        }
  }

  if ( !PlottingSetup::openBox ) {
     std::cout << std::endl <<  "\033[1;35m KEEP COOL: WE'RE NOT OPENING THE BOX YET!\033[0m" << std::endl << std::endl;
  } else {
     std::cout << std::endl <<  "\033[1;31m WARNING: WE'RE LOOKING INTO THE SIGNAL REGION!\033[0m" << std::endl << std::endl;
  }
  

  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;

  
//    DoPeakVsPU();
  
  if ( 0 ) {
    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);

    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 )";

  }

  
  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);
  }

  
  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");
    
     
   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) compute_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,false);
  
  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();
  
  return 0;  
}

Revision:	1.42
Committed:	Thu Apr 4 13:06:19 2013 UTC (12 years, 1 month ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.41:	+15 -5 lines
Log Message:	Added quick ttbar analysis to reproduce AN-12-375 (dileptonic ttbar cross section measurement)