TopEventSelection/bin/TtSemiLRSignalSel_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 signal_nrDir = 5;
0027 const int signal_nrFiles[signal_nrDir] = {45, 40, 30, 15, 17};
0028 const TString signal_path[signal_nrDir] = {
0029     "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt0j/Alpgen_tt0j_TtSemiMuEvents_",
0030     "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt1j/Alpgen_tt1j_TtSemiMuEvents_",
0031     "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt2j/Alpgen_tt2j_TtSemiMuEvents_",
0032     "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt3j/Alpgen_tt3j_TtSemiMuEvents_",
0033     "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt4j/Alpgen_tt4j_TtSemiMuEvents_"};
0034 const int signal_NrEv[signal_nrDir] = {76250, 68000, 40000, 16000, 24425};
0035
0036 const int bckgd_nrDir = 8;
0037 const int bckgd_nrFiles[bckgd_nrDir] = {45, 40, 30, 15, 17, 10, 6, 2};  //{9,4,2};//{35,30,20,10,10};
0038 const TString bckgd_path[bckgd_nrDir] = {"dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt0j/Alpgen_tt0j_TtOtherEvents_",
0039                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt1j/Alpgen_tt1j_TtOtherEvents_",
0040                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt2j/Alpgen_tt2j_TtOtherEvents_",
0041                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt3j/Alpgen_tt3j_TtOtherEvents_",
0042                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/tt4j/Alpgen_tt4j_TtOtherEvents_",
0043                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/W4j/Alpgen_W4j_TtSemiMuEvents_",
0044                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/W5j/Alpgen_W5j_TtSemiMuEvents_",
0045                                          "dcap:///pnfs/iihe/becms/ghammad/TQAF136Final/W6j/Alpgen_W6j_TtSemiMuEvents_"};
0046
0047 const int bckgd_NrEv[bckgd_nrDir] = {
0048     76250, 68000, 40000, 16000, 24425, 69700, 18000, 12500};  //{69700,18000,12500};//{152500,136000,80000,32000,49000};
0049
0050 //observable histogram variables
0051 const int nrSignalSelObs = 7;
0052 const int SignalSelObs[nrSignalSelObs] = {1, 3, 6, 12, 15, 16, 31};
0053 const TString SignalSelInputFileName = "./TtSemiLRSignalSelAllObs.root";
0054
0055 //likelihood histogram variables
0056 const int nrSignalSelLRtotBins = 35;
0057 const double SignalSelLRtotMin = -10;
0058 const double SignalSelLRtotMax = 8;
0059 const char *SignalSelLRtotFitFunction = "[0]/(1 + 1/exp([1]*([2] - x)))+[3]";
0060
0061 //output files ps/root
0062 const TString SignalSelOutfileName = "./TtSemiLRSignalSelSelObsAndPurity.root";
0063 const TString SignalSelPSfile = "./TtSemiLRSignalSelSelObsAndPurity.ps";
0064
0065 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0066
0067 //
0068 // Global variables
0069 //
0070 LRHelpFunctions *myLRhelper;
0071 void doEventloop();
0072 std::vector<int> obsNrs;
0073
0074 //
0075 // Main analysis
0076 //
0077
0078 int main() {
0079   gSystem->Load("libFWCoreFWLite");
0080   FWLiteEnabler::enable();
0081
0082   // define all histograms & fit functions
0083   //to replace with something more elegant
0084   for (int j = 0; j < nrSignalSelObs; j++) {
0085     obsNrs.push_back(SignalSelObs[j]);
0086   }
0087   myLRhelper =
0088       new LRHelpFunctions(nrSignalSelLRtotBins, SignalSelLRtotMin, SignalSelLRtotMax, SignalSelLRtotFitFunction);
0089
0090   // read in S/S+N fits of observables to use
0091   myLRhelper->readObsHistsAndFits(SignalSelInputFileName, obsNrs, false);
0092
0093   // fill calculated LR value for each signal or background contributions
0094   doEventloop();
0095
0096   // make Purity vs logLR and Purity vs. Efficiency plot
0097   myLRhelper->makeAndFitPurityHists();
0098
0099   // store histograms and fits in root-file
0100   myLRhelper->storeToROOTfile(SignalSelOutfileName);
0101
0102   // make some control plots and put them in a .ps file
0103   myLRhelper->storeControlPlots(SignalSelPSfile);
0104 }
0105
0106 //
0107 // Loop over the events (with the definition of what is considered signal and background)
0108 //
0109
0110 void doEventloop() {
0111   std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP FOR SIGNAL ****" << std::endl;
0112
0113   /********************************************** for the signal **********************************************/
0114
0115   int okEvents = 0, totNrEv = 0;
0116   for (int nrDir = 0; nrDir < signal_nrDir; nrDir++) {
0117     std::cout << " Signal : " << signal_path[nrDir] << std::endl;
0118
0119     int Signal_totNrEv = 0, Signal_okEvents = 0;
0120     for (int nr = 1; nr <= signal_nrFiles[nrDir]; nr++) {
0121       TString signal_ft = signal_path[nrDir];
0122       signal_ft += nr;
0123       signal_ft += ".root";
0124       if (!gSystem->AccessPathName(signal_ft)) {
0125         TFile *signal_file = TFile::Open(signal_ft);
0126         TTree *signal_events = dynamic_cast<TTree *>(signal_file->Get("Events"));
0127         assert(signal_events != nullptr);
0128
0129         TBranch *signal_solsbranch = signal_events->GetBranch("TtSemiEvtSolutions_solutions__TEST.obj");
0130         assert(signal_solsbranch != nullptr);
0131         std::vector<TtSemiEvtSolution> signal_sols;
0132         signal_solsbranch->SetAddress(&signal_sols);
0133
0134         //loop over all events in a file
0135         for (int ev = 0; ev < signal_events->GetEntries(); ++ev) {
0136           if (Signal_totNrEv > signal_NrEv[nrDir] && signal_NrEv[nrDir] != -1)
0137             continue;
0138           ++Signal_totNrEv;
0139           ++totNrEv;
0140           if ((double)((totNrEv * 1.) / 5000.) == (double)(totNrEv / 5000))
0141             std::cout << "  Processing signal event " << totNrEv << std::endl;
0142           signal_solsbranch->GetEntry(ev);
0143           if (signal_sols.size() == 12) {
0144             // get observable values
0145             std::vector<double> signal_obsVals;
0146             for (int j = 0; j < nrSignalSelObs; j++) {
0147               if (myLRhelper->obsFitIncluded(obsNrs[j]))
0148                 signal_obsVals.push_back(signal_sols[0].getLRSignalEvtObsVal(obsNrs[j]));
0149             }
0150
0151             double logLR = myLRhelper->calcLRval(signal_obsVals);
0152             myLRhelper->fillLRSignalHist(logLR);
0153             ++Signal_okEvents;
0154             ++okEvents;
0155           }
0156         }
0157         signal_file->Close();
0158       } else {
0159         std::cout << signal_ft << " doesn't exist" << std::endl;
0160       }
0161     }
0162     std::cout << std::endl
0163               << "********************  STATISTICS FOR SIGNAL " << signal_path[nrDir] << " ***********************"
0164               << std::endl;
0165     std::cout << std::endl << " Nb of processed events  :" << (Signal_totNrEv) << std::endl;
0166     std::cout << std::endl << " Nb of events filled in the histo :" << (Signal_okEvents) << std::endl;
0167     std::cout << std::endl << "******************************************************************" << std::endl;
0168   }
0169   std::cout << std::endl << "********************  STATISTICS FOR SIGNAL ***********************" << std::endl;
0170   std::cout << std::endl << " Nb of processed events  :" << (totNrEv) << std::endl;
0171   std::cout << std::endl << " Nb of events filled in the histo :" << (okEvents) << std::endl;
0172   std::cout << std::endl << "******************************************************************" << std::endl;
0173
0174   std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP FOR BCKGD ****" << std::endl;
0175
0176   /********************************************** for the background **********************************************/
0177
0178   okEvents = 0, totNrEv = 0;
0179   for (int nrDir = 0; nrDir < bckgd_nrDir; nrDir++) {
0180     std::cout << " Background : " << bckgd_path[nrDir] << std::endl;
0181
0182     int Bckgd_totNrEv = 0, Bckgd_okEvents = 0;
0183     for (int nr = 1; nr <= bckgd_nrFiles[nrDir]; nr++) {
0184       TString bckgd_ft = bckgd_path[nrDir];
0185       bckgd_ft += nr;
0186       bckgd_ft += ".root";
0187       if (!gSystem->AccessPathName(bckgd_ft)) {
0188         TFile *bckgd_file = TFile::Open(bckgd_ft);
0189         TTree *bckgd_events = dynamic_cast<TTree *>(bckgd_file->Get("Events"));
0190         assert(bckgd_events != nullptr);
0191
0192         TBranch *bckgd_solsbranch = bckgd_events->GetBranch("TtSemiEvtSolutions_solutions__TEST.obj");
0193         assert(bckgd_solsbranch != nullptr);
0194         std::vector<TtSemiEvtSolution> bckgd_sols;
0195         bckgd_solsbranch->SetAddress(&bckgd_sols);
0196
0197         //loop over all events in a file
0198         for (int ev = 0; ev < bckgd_events->GetEntries(); ++ev) {
0199           if (Bckgd_totNrEv > bckgd_NrEv[nrDir])
0200             continue;
0201           ++Bckgd_totNrEv;
0202           ++totNrEv;
0203           if ((double)((totNrEv * 1.) / 1000.) == (double)(totNrEv / 1000))
0204             std::cout << "  Processing bckgd event " << totNrEv << std::endl;
0205           bckgd_solsbranch->GetEntry(ev);
0206           if (bckgd_sols.size() == 12) {
0207             // get observable values
0208             std::vector<double> bckgd_obsVals;
0209             for (int j = 0; j < nrSignalSelObs; j++) {
0210               if (myLRhelper->obsFitIncluded(obsNrs[j]))
0211                 bckgd_obsVals.push_back(bckgd_sols[0].getLRSignalEvtObsVal(obsNrs[j]));
0212             }
0213
0214             double logLR = myLRhelper->calcLRval(bckgd_obsVals);
0215             myLRhelper->fillLRBackgroundHist(logLR);
0216             ++okEvents;
0217             ++Bckgd_okEvents;
0218           }
0219         }
0220         bckgd_file->Close();
0221       } else {
0222         std::cout << bckgd_ft << " doesn't exist" << std::endl;
0223       }
0224     }
0225     std::cout << std::endl
0226               << "********************  STATISTICS FOR BCKGD " << bckgd_path[nrDir] << " ***********************"
0227               << std::endl;
0228     std::cout << std::endl << " Nb of processed events  :" << (Bckgd_totNrEv) << std::endl;
0229     std::cout << std::endl << " Nb of events filled in the histo :" << (Bckgd_okEvents) << std::endl;
0230     std::cout << std::endl << "******************************************************************" << std::endl;
0231   }
0232   std::cout << std::endl << "********************  STATISTICS FOR BCKGD ***********************" << std::endl;
0233   std::cout << std::endl << " Nb of processed events  :" << (totNrEv) << std::endl;
0234   std::cout << std::endl << " Nb of events filled in the histo :" << (okEvents) << std::endl;
0235   std::cout << std::endl << "******************************************************************" << std::endl;
0236 }