TopJetCombination/bin/TtSemiLRJetComb_SoverSplusBLoop.cpp

0001 #include <iostream>
0002 #include <cassert>
0003 #include <TROOT.h>
0004 #include <TSystem.h>
0005 #include <TFile.h>
0006 #include <TTree.h>
0007 #include <TBranch.h>
0008 #include <TCanvas.h>
0009 #include <TH1.h>
0010 #include <TH2.h>
0011 #include <TGraph.h>
0012 #include <TF1.h>
0013 #include <TFormula.h>
0014 #include <TStyle.h>
0015 #include <TKey.h>
0016 #include <vector>
0017 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0018 #include "AnalysisDataFormats/TopObjects/interface/TtSemiEvtSolution.h"
0019 #include "TopQuarkAnalysis/TopTools/interface/LRHelpFunctions.h"
0020
0021 ///////////////////////
0022 // Constants         //
0023 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0024
0025 //input files
0026 const int nrFiles = 51;
0027 const TString path = "dcap://maite.iihe.ac.be:/pnfs/iihe/becms/heyninck/TtSemiMuEvents_TopRex_Juni/TtSemiMuEvents_";
0028
0029 //matching variables
0030 const bool useSpaceAngle = true;
0031 const double SumAlphaCut = 0.7;
0032
0033 //observable histogram variables (include all defined observables!!!)
0034 const int nrJetCombObs = 7;
0035 const int JetCombObs[nrJetCombObs] = {1, 2, 3, 4, 5, 6, 7};
0036 const int nrJetCombHistBins = 50;
0037 const double JetCombObsMin[nrJetCombObs] = {0, 0, 0, 0, 0, -15, -6};
0038 const double JetCombObsMax[nrJetCombObs] = {3, 3.5, 5, 5, 5, 70, 0};
0039
0040 //observable fit functions
0041 const char *JetCombObsFits[nrJetCombObs] = {
0042     "[0]/(1 + 1/exp([1]*([2] - x)))",                //obs1
0043     "[0]/(1 + 1/exp([1]*([2] - x)))",                //obs2
0044     "gaus",                                          //obs3
0045     "gaus",                                          //obs4
0046     "gaus",                                          //obs5
0047     "([0]+[3]*abs(x)/x)*(1-exp([1]*(abs(x)-[2])))",  //obs6
0048     "[0]/(1 + 1/exp([1]*([2] - x)))"                 //obs7
0049 };
0050
0051 //output files ps/root
0052 const TString JetCombOutfileName = "../data/TtSemiLRJetCombAllObs.root";
0053 const TString JetCombPSfile = "../data/TtSemiLRJetCombAllObs.ps";
0054
0055 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0056
0057 //
0058 // Global variables
0059 //
0060 LRHelpFunctions *myLRhelper;
0061 void doEventloop();
0062 std::vector<int> obsNrs;
0063 std::vector<double> obsMin, obsMax;
0064 std::vector<const char *> obsFits;
0065
0066 //
0067 // Main analysis
0068 //
0069
0070 int main() {
0071   gSystem->Load("libFWCoreFWLite");
0072   FWLiteEnabler::enable();
0073
0074   // define all histograms & fit functions
0075   //to replace with something more elegant
0076   for (int j = 0; j < nrJetCombObs; j++) {
0077     obsNrs.push_back(JetCombObs[j]);
0078     obsMin.push_back(JetCombObsMin[j]);
0079     obsMax.push_back(JetCombObsMax[j]);
0080     obsFits.push_back(JetCombObsFits[j]);
0081   }
0082   myLRhelper = new LRHelpFunctions(obsNrs, nrJetCombHistBins, obsMin, obsMax, obsFits);
0083
0084   std::vector<double> parsFobs6;
0085   parsFobs6.push_back(0.8);
0086   parsFobs6.push_back(-0.1);
0087   parsFobs6.push_back(-0.8);
0088   parsFobs6.push_back(0.2);
0089   myLRhelper->setObsFitParameters(6, parsFobs6);
0090
0091   // fill signal and background contributions to S and B histograms
0092   doEventloop();
0093
0094   // normalize the S and B histograms to construct the pdf's
0095   myLRhelper->normalizeSandBhists();
0096
0097   // produce and fit the S/S+N histograms
0098   myLRhelper->makeAndFitSoverSplusBHists();
0099
0100   // store histograms and fits in root-file
0101   myLRhelper->storeToROOTfile(JetCombOutfileName);
0102
0103   // make some control plots and put them in a .ps file
0104   myLRhelper->storeControlPlots(JetCombPSfile);
0105 }
0106
0107 //
0108 // Loop over the events (with the definition of what is considered signal and background)
0109 //
0110
0111 void doEventloop() {
0112   std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP ****" << std::endl;
0113   int okEvents = 0, totNrEv = 0;
0114   for (int nr = 1; nr <= nrFiles; nr++) {
0115     TString ft = path;
0116     ft += nr - 1;
0117     ft += ".root";
0118     if (!gSystem->AccessPathName(ft)) {
0119       TFile *file = TFile::Open(ft);
0120       TTree *events = dynamic_cast<TTree *>(file->Get("Events"));
0121       assert(events != nullptr);
0122       TBranch *solsbranch = events->GetBranch("TtSemiEvtSolutions_solutions__TtEventReco.obj");
0123       assert(solsbranch != nullptr);
0124       std::vector<TtSemiEvtSolution> sols;
0125       solsbranch->SetAddress(&sols);
0126
0127       //loop over all events in a file
0128       for (int ev = 0; ev < events->GetEntries(); ++ev) {
0129         ++totNrEv;
0130         if ((double)((totNrEv * 1.) / 1000.) == (double)(totNrEv / 1000))
0131           std::cout << "  Processing event " << totNrEv << std::endl;
0132         solsbranch->GetEntry(ev);
0133         if (sols.size() == 12) {
0134           // check if good matching solution exists
0135           bool trueSolExists = false;
0136           for (int s = 0; s < 12; s++) {
0137             if (sols[s].getMCBestSumAngles() < SumAlphaCut)
0138               trueSolExists = true;
0139           }
0140           if (trueSolExists) {
0141             for (int s = 0; s < 12; s++) {
0142               // get observable values
0143               std::vector<double> obsVals;
0144               for (int j = 0; j < nrJetCombObs; j++) {
0145                 if (myLRhelper->obsFitIncluded((unsigned int)obsNrs[j]))
0146                   obsVals.push_back(sols[s].getLRJetCombObsVal((unsigned int)obsNrs[j]));
0147               }
0148               // Fill the observables for each jet combination if a good matching exists
0149               // signal: best matching solution
0150               // background: all other solutions
0151               if (sols[s].getMCBestJetComb() == s) {
0152                 myLRhelper->fillToSignalHists(obsVals);
0153                 ++okEvents;
0154               } else {
0155                 myLRhelper->fillToBackgroundHists(obsVals);
0156               }
0157             }
0158           }
0159         }
0160       }
0161       file->Close();
0162     } else {
0163       std::cout << ft << " doesn't exist" << std::endl;
0164     }
0165   }
0166   std::cout << std::endl << "***********************  STATISTICS  *************************" << std::endl;
0167   std::cout << " Probability that a correct jet combination exists:" << std::endl;
0168   std::cout << " (fraction events with ";
0169   if (useSpaceAngle)
0170     std::cout << "min SumAngle_jp < ";
0171   if (!useSpaceAngle)
0172     std::cout << "min DR_jp < ";
0173   std::cout << SumAlphaCut << " )" << std::endl;
0174   std::cout << std::endl << "                 " << (100. * okEvents) / (1. * totNrEv) << " %" << std::endl;
0175   std::cout << std::endl << "******************************************************************" << std::endl;
0176 }