File indexing completed on 2024-04-06 12:31:15
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
0023
0024
0025
0026 const int signal_nrDir = 5;
0027 const int signal_nrFiles[signal_nrDir] = {30, 25, 20, 15, 10};
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] = {30, 25, 20, 15, 10, 10, 10, 5};
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};
0049
0050
0051 const int nrSignalSelObs = 38;
0052 const int SignalSelObs[nrSignalSelObs] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
0053 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38};
0054 const int nrSignalSelHistBins = 35;
0055 const double SignalSelObsMin[nrSignalSelObs] = {10, 0, 0, 0, 0, 0.5, 0, -5, -2, -10, 0, 0, 0, 0.15, 0, 0, 0, 0, 0,
0056 0, 0.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.55, 0, 0, 0, 0, 0, 0, 0};
0057 const double SignalSelObsMax[nrSignalSelObs] = {200, 200, 1000, 300, 0.5, 1, 1600, 5, 20, 15, 50, 1, 1,
0058 1, 1, 0.3, 0.9, 0.3, 0.35, 1000, 500, 0.8, 70, 2, 1, 1,
0059 0.8, 1, 0.65, 0.75, 1, 0.7, 0.65, 0.6, 1, 0.35, 1, 0.4};
0060
0061
0062
0063 TFormula gauss("gauss", "gaus");
0064 TFormula symgauss("symgauss", "[0]*(exp(-0.5*(x/[1])**2))");
0065 TFormula dblgauss("dblgauss", "[0]*(exp(-0.5*((x-[1])/[2])**2)+exp(-0.5*((x+[3])/[4])**2))");
0066 TFormula symdblgauss("symdblgauss", "[0]*(exp(-0.5*((x-[1])/[2])**2)+exp(-0.5*((x+[1])/[2])**2))");
0067 TFormula sigm("sigm", "[0]/(1 + 1/exp([1]*([2] - x)))");
0068 TFormula sigmc("sigmc", "[0]/(1 + 1/exp([1]*([2] - x)))+[3]");
0069 TFormula dblsigm("dblsigm", "[0]/(1 + 1/exp([1]**2*([2] - x)))/(1 + 1/exp([3]**2*(x - [4])))");
0070 TFormula symdblsigm("symdblsigm", "[0]/(1 + 1/exp([1]**2*([2] - x)))/(1 + 1/exp([1]**2*([2] + x)))");
0071
0072 const char *SignalSelObsFits[nrSignalSelObs] = {
0073
0074 "[0]/(1 + 1/exp([1]*([2] - x)))",
0075 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0076 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol4",
0077 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0078 "[0]/(1 + 1/exp([1]*([2] - x)))",
0079 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]+pol3(4)+pol5(8)",
0080 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol4",
0081 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol3+[4]",
0082 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol3+[4]",
0083 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0084 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0085 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0086 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]",
0087 "[0]/(1 + 1/exp([1]*([2] - x)))",
0088 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol3",
0089 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol4",
0090 "pol9",
0091 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol4",
0092 "pol6",
0093 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*exp([4]*([5] - x))",
0094 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*pol4",
0095 "[0]/(1 + 1/exp([1]*([2] - x)))",
0096 "[0]/(1 + 1/exp([1]*([2] - x)))",
0097 "[0]/(1 + 1/exp([1]*([2] - x)))",
0098 "[0]/(1 + 1/exp([1]*([2] - x)))",
0099 "[0]/(1 + 1/exp([1]*([2] - x)))",
0100 "pol6",
0101 "[0]/(1 + 1/exp([1]*([2] - x)))",
0102 "[0]/(1 + 1/exp([1]*([2] - x)))+[3]*exp([4]*([5] - x))",
0103 "[0]/(1 + 1/exp([1]*([2] - x)))",
0104 "[0]/(1 + 1/exp([1]*([2] - x)))",
0105 "pol6",
0106 "[0]/(1 + 1/exp([1]*([2] - x)))",
0107 "[0]/(1 + 1/exp([1]*([2] - x)))",
0108 "pol3+[3]*exp([4]*([5] - x))",
0109 "[0]/(1 + 1/exp([1]*([2] - x)))",
0110 "pol3+[3]*exp([4]*([5] - x))",
0111 "pol5+[3]*exp([4]*([5] - x))",
0112
0113
0114
0115
0116
0117 };
0118
0119
0120 const TString SignalSelOutfileName = "./TtSemiLRSignalSelAllObs.root";
0121 const TString SignalSelPSfile = "./TtSemiLRSignalSelAllObs.ps";
0122
0123
0124
0125
0126
0127
0128 LRHelpFunctions *myLRhelper;
0129 void doEventloop();
0130 std::vector<int> obsNrs;
0131 std::vector<double> obsMin, obsMax;
0132 std::vector<const char *> obsFits;
0133
0134 bool MuonIso = true;
0135
0136
0137
0138
0139
0140 int main() {
0141 gSystem->Load("libFWCoreFWLite");
0142 FWLiteEnabler::enable();
0143
0144
0145
0146 for (int j = 0; j < nrSignalSelObs; j++) {
0147 obsNrs.push_back(SignalSelObs[j]);
0148 obsMin.push_back(SignalSelObsMin[j]);
0149 obsMax.push_back(SignalSelObsMax[j]);
0150 obsFits.push_back(SignalSelObsFits[j]);
0151 }
0152 myLRhelper = new LRHelpFunctions(obsNrs, nrSignalSelHistBins, obsMin, obsMax, obsFits);
0153
0154
0155 std::vector<double> parsFobs1;
0156 parsFobs1.push_back(20);
0157 parsFobs1.push_back(0.04);
0158 parsFobs1.push_back(21);
0159 parsFobs1.push_back(0.04);
0160 myLRhelper->setObsFitParameters(1, parsFobs1);
0161
0162 std::vector<double> parsFobs2;
0163 parsFobs2.push_back(0.495);
0164 parsFobs2.push_back(-0.148);
0165 parsFobs2.push_back(60.33);
0166 parsFobs2.push_back(0.396);
0167 myLRhelper->setObsFitParameters(2, parsFobs2);
0168
0169 std::vector<double> parsFobs3;
0170 parsFobs3.push_back(7.60418e-01);
0171 parsFobs3.push_back(-3.31635e-02);
0172 parsFobs3.push_back(1.57387e+02);
0173 parsFobs3.push_back(-1.23931e-08);
0174 parsFobs3.push_back(-1.90918e-04);
0175 myLRhelper->setObsFitParameters(3, parsFobs3);
0176
0177 std::vector<double> parsFobs4;
0178 parsFobs4.push_back(1.087);
0179 parsFobs4.push_back(-0.1978);
0180 parsFobs4.push_back(22.803);
0181 parsFobs4.push_back(-0.126);
0182 myLRhelper->setObsFitParameters(4, parsFobs4);
0183
0184 std::vector<double> parsFobs7;
0185 parsFobs7.push_back(0.606878);
0186 parsFobs7.push_back(-1.52796e-02);
0187 parsFobs7.push_back(2.50574e+02);
0188 parsFobs7.push_back(-4.46936e-10);
0189 parsFobs7.push_back(-1.48804e-04);
0190 myLRhelper->setObsFitParameters(7, parsFobs7);
0191
0192 std::vector<double> parsFobs8;
0193 parsFobs8.push_back(3.30611e-01);
0194 parsFobs8.push_back(-8.34406e+00);
0195 parsFobs8.push_back(1.04307e+00);
0196 parsFobs8.push_back(-1.75190e-03);
0197 parsFobs8.push_back(5.66972e-01);
0198 myLRhelper->setObsFitParameters(8, parsFobs8);
0199
0200 std::vector<double> parsFobs9;
0201 parsFobs9.push_back(6.37793e-01);
0202 parsFobs9.push_back(-1.71768e+00);
0203 parsFobs9.push_back(1.88952e+00);
0204 parsFobs9.push_back(-1.03833e-03);
0205 parsFobs9.push_back(3.30284e-01);
0206 myLRhelper->setObsFitParameters(9, parsFobs9);
0207
0208 std::vector<double> parsFobs10;
0209 parsFobs10.push_back(0.618);
0210 parsFobs10.push_back(-1.579);
0211 parsFobs10.push_back(-0.10812);
0212 parsFobs10.push_back(0.342);
0213 myLRhelper->setObsFitParameters(10, parsFobs10);
0214
0215 std::vector<double> parsFobs11;
0216 parsFobs11.push_back(0.7624);
0217 parsFobs11.push_back(-0.64975);
0218 parsFobs11.push_back(3.1225);
0219 parsFobs11.push_back(0.218675);
0220 myLRhelper->setObsFitParameters(11, parsFobs11);
0221
0222 std::vector<double> parsFobs12;
0223 parsFobs12.push_back(1.57736e-01);
0224 parsFobs12.push_back(-2.01467e+01);
0225 parsFobs12.push_back(5.97867e-01);
0226 parsFobs12.push_back(3.81101e-01);
0227 myLRhelper->setObsFitParameters(12, parsFobs12);
0228
0229 std::vector<double> parsFobs13;
0230 parsFobs13.push_back(1.57736e-01);
0231 parsFobs13.push_back(-2.01467e+01);
0232 parsFobs13.push_back(5.97867e-01);
0233 parsFobs13.push_back(3.81101e-01);
0234 myLRhelper->setObsFitParameters(13, parsFobs13);
0235
0236 std::vector<double> parsFobs15;
0237 parsFobs15.push_back(0.6672);
0238 parsFobs15.push_back(-9.3022);
0239 parsFobs15.push_back(0.03384);
0240 parsFobs15.push_back(0.00014967);
0241 parsFobs15.push_back(-4315.96);
0242 myLRhelper->setObsFitParameters(15, parsFobs15);
0243
0244 std::vector<double> parsFobs16;
0245 parsFobs16.push_back(0.56855);
0246 parsFobs16.push_back(-165.768);
0247 parsFobs16.push_back(0.0021543);
0248 parsFobs16.push_back(0.0148839);
0249 parsFobs16.push_back(4391.8);
0250 myLRhelper->setObsFitParameters(16, parsFobs16);
0251
0252 std::vector<double> parsFobs17;
0253 parsFobs17.push_back(0.45862);
0254 parsFobs17.push_back(-42.3119);
0255 parsFobs17.push_back(0.0024431);
0256 parsFobs17.push_back(-0.0082168);
0257 parsFobs17.push_back(-41.3239);
0258 myLRhelper->setObsFitParameters(17, parsFobs17);
0259
0260 std::vector<double> parsFobs18;
0261 parsFobs18.push_back(0.57713);
0262 parsFobs18.push_back(-88.4547);
0263 parsFobs18.push_back(-0.0079014);
0264 parsFobs18.push_back(-0.025394);
0265 parsFobs18.push_back(4512.33);
0266 myLRhelper->setObsFitParameters(18, parsFobs18);
0267
0268 std::vector<double> parsFobs33;
0269 parsFobs33.push_back(5.99882e-01);
0270 parsFobs33.push_back(-1.33575e+01);
0271 parsFobs33.push_back(1.24161e-01);
0272 myLRhelper->setObsFitParameters(33, parsFobs33);
0273
0274 std::vector<double> parsFobs35;
0275 parsFobs35.push_back(2.49026e-01);
0276 parsFobs35.push_back(1.08819e+00);
0277 parsFobs35.push_back(-7.26373e-01);
0278 parsFobs35.push_back(1.26367e-07);
0279 parsFobs35.push_back(5.51754e+02);
0280 parsFobs35.push_back(3.94562e-02);
0281 myLRhelper->setObsFitParameters(35, parsFobs35);
0282
0283 std::vector<double> parsFobs37;
0284 parsFobs37.push_back(1.43676e-01);
0285 parsFobs37.push_back(2.44475e+00);
0286 parsFobs37.push_back(-4.56374e+00);
0287 parsFobs37.push_back(3.01449e+00);
0288 parsFobs37.push_back(4.65671e+01);
0289 parsFobs37.push_back(-4.40296e-02);
0290 myLRhelper->setObsFitParameters(37, parsFobs37);
0291
0292
0293 doEventloop();
0294
0295
0296
0297 myLRhelper->normalizeSandBhists();
0298
0299
0300 myLRhelper->makeAndFitSoverSplusBHists();
0301
0302
0303 myLRhelper->storeToROOTfile(SignalSelOutfileName);
0304
0305
0306 myLRhelper->storeControlPlots(SignalSelPSfile);
0307 }
0308
0309
0310
0311
0312
0313 void doEventloop() {
0314 std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP FOR SIGNAL ****" << std::endl;
0315
0316
0317
0318 int okEvents = 0, totNrEv = 0;
0319 for (int nrDir = 0; nrDir < signal_nrDir; nrDir++) {
0320 std::cout << " Signal : " << signal_path[nrDir] << std::endl;
0321
0322 int Signal_totNrEv = 0, Signal_okEvents = 0;
0323 for (int nr = 1; nr <= signal_nrFiles[nrDir]; nr++) {
0324 TString signal_ft = signal_path[nrDir];
0325 signal_ft += nr;
0326 signal_ft += ".root";
0327 if (!gSystem->AccessPathName(signal_ft)) {
0328 TFile *signal_file = TFile::Open(signal_ft);
0329 TTree *signal_events = dynamic_cast<TTree *>(signal_file->Get("Events"));
0330 assert(signal_events != nullptr);
0331
0332 TBranch *signal_solsbranch = signal_events->GetBranch("TtSemiEvtSolutions_solutions__TEST.obj");
0333 assert(signal_solsbranch != nullptr);
0334 std::vector<TtSemiEvtSolution> signal_sols;
0335
0336
0337
0338 for (int ev = 0; ev < signal_events->GetEntries(); ++ev) {
0339 if (Signal_totNrEv > signal_NrEv[nrDir] && signal_NrEv[nrDir] != -1)
0340 continue;
0341 ++Signal_totNrEv;
0342 ++totNrEv;
0343 if ((double)((totNrEv * 1.) / 5000.) == (double)(totNrEv / 5000))
0344 std::cout << " Processing signal event " << totNrEv << std::endl;
0345 signal_solsbranch->SetAddress(&signal_sols);
0346 signal_solsbranch->GetEntry(ev);
0347 signal_events->GetEntry(ev, 0);
0348 if (signal_sols.size() == 12) {
0349
0350 std::vector<double> signal_obsVals;
0351 for (int j = 0; j < nrSignalSelObs; j++) {
0352 if (myLRhelper->obsFitIncluded(obsNrs[j]))
0353 signal_obsVals.push_back(signal_sols[0].getLRSignalEvtObsVal(obsNrs[j]));
0354 }
0355
0356
0357
0358 myLRhelper->fillToSignalHists(signal_obsVals);
0359 ++Signal_okEvents;
0360 ++okEvents;
0361 }
0362 }
0363 signal_file->Close();
0364
0365 } else {
0366 std::cout << signal_ft << " doesn't exist" << std::endl;
0367 }
0368 }
0369 std::cout << std::endl
0370 << "******************** STATISTICS FOR SIGNAL " << signal_path[nrDir] << " ***********************"
0371 << std::endl;
0372 std::cout << std::endl << " Nb of requested events :" << (signal_NrEv[nrDir]) << std::endl;
0373 std::cout << std::endl << " Nb of processed events :" << (Signal_totNrEv) << std::endl;
0374 std::cout << std::endl << " Nb of events filled in the histo :" << (Signal_okEvents) << std::endl;
0375 std::cout << std::endl << "******************************************************************" << std::endl;
0376 }
0377 std::cout << std::endl << "******************** STATISTICS FOR SIGNAL ***********************" << std::endl;
0378 std::cout << std::endl << " Nb of processed events :" << (totNrEv) << std::endl;
0379 std::cout << std::endl << " Nb of events filled in the histo :" << (okEvents) << std::endl;
0380 std::cout << std::endl << "******************************************************************" << std::endl;
0381
0382 std::cout << std::endl << std::endl << "**** STARTING EVENT LOOP FOR BCKGD ****" << std::endl;
0383
0384
0385
0386 okEvents = 0, totNrEv = 0;
0387 for (int nrDir = 0; nrDir < bckgd_nrDir; nrDir++) {
0388 std::cout << " Background : " << bckgd_path[nrDir] << std::endl;
0389
0390 int Bckgd_totNrEv = 0, Bckgd_okEvents = 0;
0391 for (int nr = 1; nr <= bckgd_nrFiles[nrDir]; nr++) {
0392 TString bckgd_ft = bckgd_path[nrDir];
0393 bckgd_ft += nr;
0394 bckgd_ft += ".root";
0395 if (!gSystem->AccessPathName(bckgd_ft)) {
0396 TFile *bckgd_file = TFile::Open(bckgd_ft);
0397 TTree *bckgd_events = dynamic_cast<TTree *>(bckgd_file->Get("Events"));
0398 assert(bckgd_events != nullptr);
0399
0400 TBranch *bckgd_solsbranch = bckgd_events->GetBranch("TtSemiEvtSolutions_solutions__TEST.obj");
0401 assert(bckgd_solsbranch != nullptr);
0402 std::vector<TtSemiEvtSolution> bckgd_sols;
0403
0404
0405
0406 for (int ev = 0; ev < bckgd_events->GetEntries(); ++ev) {
0407 if (Bckgd_totNrEv > bckgd_NrEv[nrDir] && bckgd_NrEv[nrDir] != -1)
0408 continue;
0409 ++Bckgd_totNrEv;
0410 ++totNrEv;
0411 if ((double)((totNrEv * 1.) / 5000.) == (double)(totNrEv / 5000))
0412 std::cout << " Processing bckgd event " << totNrEv << std::endl;
0413 bckgd_solsbranch->SetAddress(&bckgd_sols);
0414 bckgd_solsbranch->GetEntry(ev);
0415 bckgd_events->GetEntry(ev, 0);
0416 if (bckgd_sols.size() == 12) {
0417
0418 std::vector<double> bckgd_obsVals;
0419 for (int j = 0; j < nrSignalSelObs; j++) {
0420 if (myLRhelper->obsFitIncluded(obsNrs[j]))
0421 bckgd_obsVals.push_back(bckgd_sols[0].getLRSignalEvtObsVal(obsNrs[j]));
0422 }
0423
0424
0425 myLRhelper->fillToBackgroundHists(bckgd_obsVals);
0426 ++okEvents;
0427 ++Bckgd_okEvents;
0428 }
0429 }
0430 bckgd_file->Close();
0431
0432 } else {
0433 std::cout << bckgd_ft << " doesn't exist" << std::endl;
0434 }
0435 }
0436 std::cout << std::endl
0437 << "******************** STATISTICS FOR BCKGD " << bckgd_path[nrDir] << " ***********************"
0438 << std::endl;
0439 std::cout << std::endl << " Nb of requested events :" << (bckgd_NrEv[nrDir]) << std::endl;
0440 std::cout << std::endl << " Nb of processed events :" << (Bckgd_totNrEv) << std::endl;
0441 std::cout << std::endl << " Nb of events filled in the histo :" << (Bckgd_okEvents) << std::endl;
0442 std::cout << std::endl << "******************************************************************" << std::endl;
0443 }
0444 std::cout << std::endl << "******************** STATISTICS FOR BCKGD ***********************" << std::endl;
0445 std::cout << std::endl << " Nb of processed events :" << (totNrEv) << std::endl;
0446 std::cout << std::endl << " Nb of events filled in the histo :" << (okEvents) << std::endl;
0447 std::cout << std::endl << "******************************************************************" << std::endl;
0448 }