TopJetCombination/bin/TtSemiLRJetComb_PurityLoop.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 double SumAlphaCut = 0.7;
0031
0032 //select which observables to use
0033 const int nrJetCombObs = 7;
0034 const int JetCombObs[nrJetCombObs] = {1, 2, 3, 4, 5, 6, 7};
0035 const TString JetCombInputFileName = "../data/TtSemiLRJetCombAllObs.root";
0036
0037 //likelihood histogram variables
0038 const int nrJetCombLRtotBins = 30;
0039 const double JetCombLRtotMin = -5;
0040 const double JetCombLRtotMax = 7;
0041 const char *JetCombLRtotFitFunction = "[0]/(1 + 1/exp([1]*([2] - x)))";
0042
0043 //output files ps/root
0044 const TString JetCombOutFileName = "../data/TtSemiLRJetCombSelObsAndPurity.root";
0045 const TString JetCombPSfile = "../data/TtSemiLRJetCombSelObsAndPurity.ps";
0046
0047 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0048
0049 //
0050 // Global variables
0051 //
0052 LRHelpFunctions *myLRhelper;
0053 void doEventloop();
0054 std::vector<int> obsNrs;
0055
0056 //
0057 // Main analysis
0058 //
0059
0060 int main() {
0061   gSystem->Load("libFWCoreFWLite");
0062   FWLiteEnabler::enable();
0063
0064   // define all histograms & fit functions
0065   //to replace with something more elegant
0066   for (int j = 0; j < nrJetCombObs; j++) {
0067     obsNrs.push_back(JetCombObs[j]);
0068   }
0069   myLRhelper = new LRHelpFunctions(nrJetCombLRtotBins, JetCombLRtotMin, JetCombLRtotMax, JetCombLRtotFitFunction);
0070
0071   // read in S/S+N fits of observables to use
0072   myLRhelper->readObsHistsAndFits(JetCombInputFileName, obsNrs, false);
0073
0074   // fill calculated LR value for each signal or background contributions
0075   doEventloop();
0076
0077   // make Purity vs logLR and Purity vs. Efficiency plot
0078   myLRhelper->makeAndFitPurityHists();
0079
0080   // store histograms and fits in root-file
0081   myLRhelper->storeToROOTfile(JetCombOutFileName);
0082
0083   // make some control plots and put them in a .ps file
0084   myLRhelper->storeControlPlots(JetCombPSfile);
0085 }
0086
0087 //
0088 // Loop over the events (with the definition of what is considered signal and background)
0089 //
0090
0091 void doEventloop() {
0092   std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP ****" << std::endl;
0093   int totNrEv = 0;
0094   for (int nr = 1; nr <= nrFiles; nr++) {
0095     TString ft = path;
0096     ft += nr - 1;
0097     ft += ".root";
0098     if (!gSystem->AccessPathName(ft)) {
0099       TFile *file = TFile::Open(ft);
0100       TTree *events = dynamic_cast<TTree *>(file->Get("Events"));
0101       assert(events != nullptr);
0102       TBranch *solsbranch = events->GetBranch("TtSemiEvtSolutions_solutions__TtEventReco.obj");
0103       //TBranch * solsbranch = events->GetBranch( "TtSemiEvtSolutions_solutions__CommonBranchSel.obj" );
0104       assert(solsbranch != nullptr);
0105       std::vector<TtSemiEvtSolution> sols;
0106       solsbranch->SetAddress(&sols);
0107
0108       //loop over all events in a file
0109       for (int ev = 0; ev < events->GetEntries(); ++ev) {
0110         ++totNrEv;
0111         if ((double)((totNrEv * 1.) / 1000.) == (double)(totNrEv / 1000))
0112           std::cout << "  Processing event " << totNrEv << std::endl;
0113         solsbranch->GetEntry(ev);
0114         if (sols.size() == 12) {
0115           // check if good matching solution exists
0116           bool trueSolExists = false;
0117           for (int s = 0; s < 12; s++) {
0118             if (sols[s].getMCBestSumAngles() < SumAlphaCut)
0119               trueSolExists = true;
0120           }
0121           if (trueSolExists) {
0122             double maxLogLRVal = -999.;
0123             int maxLogLRSol = -999;
0124             //loop over solutions
0125             for (int s = 0; s < 12; s++) {
0126               // get observable values
0127               std::vector<double> obsVals;
0128               for (int j = 0; j < nrJetCombObs; j++) {
0129                 if (myLRhelper->obsFitIncluded(obsNrs[j]))
0130                   obsVals.push_back(sols[s].getLRJetCombObsVal(obsNrs[j]));
0131               }
0132               double logLR = myLRhelper->calcLRval(obsVals);
0133               if (logLR > maxLogLRVal) {
0134                 maxLogLRVal = logLR;
0135                 maxLogLRSol = s;
0136               };
0137             }
0138             if (sols[maxLogLRSol].getMCBestSumAngles() < SumAlphaCut &&
0139                 sols[maxLogLRSol].getMCBestJetComb() == maxLogLRSol) {
0140               myLRhelper->fillLRSignalHist(maxLogLRVal);
0141               //std::cout << "mxLR " << maxLogLRVal << std::endl;
0142             } else {
0143               myLRhelper->fillLRBackgroundHist(maxLogLRVal);
0144               //std::cout << "mxLR (bg) " << maxLogLRVal << std::endl;
0145             }
0146           }
0147         }
0148       }
0149       file->Close();
0150     } else {
0151       std::cout << ft << " doesn't exist" << std::endl;
0152     }
0153   }
0154 }