TopTools/src/LRHelpFunctions.cc

0001 #include "TopQuarkAnalysis/TopTools/interface/LRHelpFunctions.h"
0002 #include "TopQuarkAnalysis/TopTools/test/tdrstyle.C"
0003
0004 // constructors
0005 LRHelpFunctions::LRHelpFunctions() {
0006   constructPurity = false;
0007   setTDRStyle();
0008   gStyle->SetCanvasDefW(900);
0009 }
0010
0011 LRHelpFunctions::LRHelpFunctions(const std::vector<int>& obsNr,
0012                                  int nrBins,
0013                                  const std::vector<double>& obsMin,
0014                                  const std::vector<double>& obsMax,
0015                                  const std::vector<const char*>& functions) {
0016   obsNumbers = obsNr;
0017   constructPurity = false;
0018   setTDRStyle();
0019   gStyle->SetCanvasDefW(900);
0020   for (size_t o = 0; o < obsNr.size(); o++) {
0021     // create Signal, background and s/(s+b) histogram
0022     TString htS = "Obs";
0023     htS += obsNr[o];
0024     htS += "_S";
0025     TString htB = "Obs";
0026     htB += obsNr[o];
0027     htB += "_B";
0028     TString htSB = "Obs";
0029     htSB += obsNr[o];
0030     htSB += "_SoverSplusB";
0031     hObsS.push_back(new TH1F(htS, htS, nrBins, obsMin[o], obsMax[o]));
0032     hObsB.push_back(new TH1F(htB, htB, nrBins, obsMin[o], obsMax[o]));
0033     hObsSoverSplusB.push_back(new TH1F(htSB, htSB, nrBins, obsMin[o], obsMax[o]));
0034
0035     // create the correlation 2D plots among the observables (for signal)
0036     for (size_t o2 = o + 1; o2 < obsNr.size(); o2++) {
0037       TString hcorr = "Corr_Obs";
0038       hcorr += obsNr[o];
0039       hcorr += "_Obs";
0040       hcorr += obsNr[o2];
0041       hObsCorr.push_back(new TH2F(hcorr, hcorr, nrBins, obsMin[o], obsMax[o], nrBins, obsMin[o2], obsMax[o2]));
0042     }
0043
0044     // create fit functions
0045     TString ftSB = "F_Obs";
0046     ftSB += obsNr[o];
0047     ftSB += "_SoverSplusB";
0048     fObsSoverSplusB.push_back(
0049         new TF1(ftSB, functions[o], hObsS[o]->GetXaxis()->GetXmin(), hObsS[o]->GetXaxis()->GetXmax()));
0050   }
0051 }
0052
0053 void LRHelpFunctions::recreateFitFct(const std::vector<int>& obsNr, const std::vector<const char*>& functions) {
0054   if (!fObsSoverSplusB.empty())
0055     fObsSoverSplusB.clear();
0056   for (size_t o = 0; o < obsNr.size(); o++) {
0057     // create fit functions
0058     TString ftSB = "F_Obs";
0059     ftSB += obsNr[o];
0060     ftSB += "_SoverSplusB";
0061     fObsSoverSplusB.push_back(
0062         new TF1(ftSB, functions[o], hObsS[o]->GetXaxis()->GetXmin(), hObsS[o]->GetXaxis()->GetXmax()));
0063   }
0064 }
0065
0066 LRHelpFunctions::LRHelpFunctions(int nrLRbins, double LRmin, double LRmax, const char* LRfunction) {
0067   constructPurity = true;
0068   setTDRStyle();
0069   gStyle->SetCanvasDefW(900);
0070
0071   // create LR histograms
0072   hLRtotS = new TH1F("hLRtotS", "hLRtotS", nrLRbins, LRmin, LRmax);
0073   hLRtotS->GetXaxis()->SetTitle("Combined LR");
0074   hLRtotB = new TH1F("hLRtotB", "hLRtotB", nrLRbins, LRmin, LRmax);
0075   hLRtotB->GetXaxis()->SetTitle("Combined LR");
0076   hLRtotSoverSplusB = new TH1F("hLRtotSoverSplusB", "hLRtotSoverSplusB", nrLRbins, LRmin, LRmax);
0077
0078   // create LR fit function
0079   fLRtotSoverSplusB = new TF1("fLRtotSoverSplusB", LRfunction, LRmin, LRmax);
0080 }
0081
0082 // destructor
0083 LRHelpFunctions::~LRHelpFunctions() {}
0084
0085 // member function to initialize the LR hists and fits
0086 void LRHelpFunctions::initLRHistsAndFits(int nrLRbins, double LRmin, double LRmax, const char* LRfunction) {
0087   constructPurity = true;
0088   // create LR histograms
0089   hLRtotS = new TH1F("hLRtotS", "hLRtotS", nrLRbins, LRmin, LRmax);
0090   hLRtotB = new TH1F("hLRtotB", "hLRtotB", nrLRbins, LRmin, LRmax);
0091   hLRtotSoverSplusB = new TH1F("hLRtotSoverSplusB", "hLRtotSoverSplusB", nrLRbins, LRmin, LRmax);
0092   // create LR fit function
0093   fLRtotSoverSplusB = new TF1("fLRtotSoverSplusB", LRfunction, LRmin, LRmax);
0094 }
0095
0096 // member function to set initial values to the observable fit function
0097 void LRHelpFunctions::setObsFitParameters(int obs, const std::vector<double>& fitPars) {
0098   for (size_t fit = 0; fit < fObsSoverSplusB.size(); fit++) {
0099     TString fn = "_Obs";
0100     fn += obs;
0101     if (((TString)fObsSoverSplusB[fit]->GetName()).Contains(fn)) {
0102       for (size_t p = 0; p < fitPars.size(); p++) {
0103         fObsSoverSplusB[fit]->SetParameter(p, fitPars[p]);
0104       }
0105     }
0106   }
0107 }
0108
0109 // member function to add observable values to the signal histograms
0110 void LRHelpFunctions::fillToSignalHists(const std::vector<double>& obsVals, double weight) {
0111   int hIndex = 0;
0112   for (size_t o = 0; o < obsVals.size(); o++) {
0113     hObsS[o]->Fill(obsVals[o], weight);
0114     for (size_t o2 = o + 1; o2 < obsVals.size(); o2++) {
0115       hObsCorr[hIndex]->Fill(obsVals[o], obsVals[o2], weight);
0116       ++hIndex;
0117     }
0118   }
0119 }
0120
0121 // member function to add observable values to the signal histograms
0122 void LRHelpFunctions::fillToSignalHists(int obsNbr, double obsVal, double weight) {
0123   TString obs = "Obs";
0124   obs += obsNbr;
0125   obs += "_";
0126   for (size_t f = 0; f < hObsS.size(); f++) {
0127     if (((TString)(hObsS[f]->GetName())).Contains(obs)) {
0128       hObsS[f]->Fill(obsVal, weight);
0129       return;
0130     }
0131   }
0132 }
0133
0134 // member function to add observable values to the signal histograms
0135 void LRHelpFunctions::fillToSignalCorrelation(int obsNbr1, double obsVal1, int obsNbr2, double obsVal2, double weight) {
0136   TString hcorr = "Corr_Obs";
0137   hcorr += std::min(obsNbr1, obsNbr2);
0138   hcorr += "_Obs";
0139   hcorr += std::max(obsNbr1, obsNbr2);
0140   for (size_t i = 0; i < hObsCorr.size(); i++) {
0141     if (((TString)(hObsCorr[i]->GetName())) == hcorr) {
0142       if (obsNbr1 < obsNbr2) {
0143         hObsCorr[i]->Fill(obsVal1, obsVal2, weight);
0144       } else {
0145         hObsCorr[i]->Fill(obsVal2, obsVal1, weight);
0146       }
0147       return;
0148     }
0149   }
0150 }
0151
0152 // member function to add observable values to the background histograms
0153 void LRHelpFunctions::fillToBackgroundHists(const std::vector<double>& obsVals, double weight) {
0154   for (size_t o = 0; o < obsVals.size(); o++)
0155     hObsB[o]->Fill(obsVals[o], weight);
0156 }
0157
0158 // member function to add observable values to the signal histograms
0159 void LRHelpFunctions::fillToBackgroundHists(int obsNbr, double obsVal, double weight) {
0160   TString obs = "Obs";
0161   obs += obsNbr;
0162   obs += "_";
0163   for (size_t f = 0; f < hObsB.size(); f++) {
0164     if (((TString)(hObsB[f]->GetName())).Contains(obs)) {
0165       hObsB[f]->Fill(obsVal, weight);
0166       return;
0167     }
0168   }
0169 }
0170
0171 // member function to normalize the S and B spectra
0172 void LRHelpFunctions::normalizeSandBhists() {
0173   for (size_t o = 0; o < hObsS.size(); o++) {
0174     // count entries in each histo. Do it this way instead of GetEntries method,
0175     // since the latter does not account for the weights.
0176     double nrSignEntries = 0, nrBackEntries = 0;
0177     for (int i = 0; i <= hObsS[o]->GetNbinsX() + 1; i++) {
0178       nrSignEntries += hObsS[o]->GetBinContent(i);
0179       nrBackEntries += hObsB[o]->GetBinContent(i);
0180     }
0181     for (int b = 0; b <= hObsS[o]->GetNbinsX() + 1; b++) {
0182       hObsS[o]->SetBinContent(b, hObsS[o]->GetBinContent(b) / (nrSignEntries));
0183       hObsB[o]->SetBinContent(b, hObsB[o]->GetBinContent(b) / (nrBackEntries));
0184       hObsS[o]->SetBinError(b, hObsS[o]->GetBinError(b) / (nrSignEntries));
0185       hObsB[o]->SetBinError(b, hObsB[o]->GetBinError(b) / (nrBackEntries));
0186     }
0187     //std::cout<<"Integral for obs"<<o<<" S: "<<hObsS[o]->Integral(0,10000)<<" & obs"<<o<<" B: "<<hObsB[o]->Integral(0,10000)<<std::endl;
0188   }
0189 }
0190
0191 // member function to produce and fit the S/(S+B) histograms
0192 void LRHelpFunctions::makeAndFitSoverSplusBHists() {
0193   for (size_t o = 0; o < hObsS.size(); o++) {
0194     for (int b = 0; b <= hObsS[o]->GetNbinsX() + 1; b++) {
0195       if ((hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b)) > 0) {
0196         hObsSoverSplusB[o]->SetBinContent(
0197             b, hObsS[o]->GetBinContent(b) / (hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b)));
0198         double error = sqrt(pow(hObsS[o]->GetBinError(b) * hObsB[o]->GetBinContent(b) /
0199                                     pow(hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b), 2),
0200                                 2) +
0201                             pow(hObsB[o]->GetBinError(b) * hObsS[o]->GetBinContent(b) /
0202                                     pow(hObsS[o]->GetBinContent(b) + hObsB[o]->GetBinContent(b), 2),
0203                                 2));
0204         hObsSoverSplusB[o]->SetBinError(b, error);
0205       }
0206     }
0207     hObsSoverSplusB[o]->Fit(fObsSoverSplusB[o]->GetName(), "RQ");
0208   }
0209 }
0210
0211 // member function to read the observable hists & fits from a root-file
0212 void LRHelpFunctions::readObsHistsAndFits(const TString& fileName,
0213                                           const std::vector<int>& observables,
0214                                           bool readLRplots) {
0215   hObsS.clear();
0216   hObsB.clear();
0217   hObsSoverSplusB.clear();
0218   fObsSoverSplusB.clear();
0219   TFile* fitFile = new TFile(fileName, "READ");
0220   if (observables[0] == -1) {
0221     std::cout << " ... will read hists and fit for all available observables in file " << fileName << std::endl;
0222     TList* list = fitFile->GetListOfKeys();
0223     TIter next(list);
0224     TKey* el;
0225     while ((el = (TKey*)next())) {
0226       TString keyName = el->GetName();
0227       if (keyName.Contains("F_") && keyName.Contains("_SoverSplusB")) {
0228         fObsSoverSplusB.push_back(new TF1(*((TF1*)el->ReadObj())));
0229       } else if (keyName.Contains("_SoverSplusB")) {
0230         hObsSoverSplusB.push_back(new TH1F(*((TH1F*)el->ReadObj())));
0231       } else if (keyName.Contains("_S")) {
0232         hObsS.push_back(new TH1F(*((TH1F*)el->ReadObj())));
0233       } else if (keyName.Contains("_B")) {
0234         hObsB.push_back(new TH1F(*((TH1F*)el->ReadObj())));
0235       } else if (keyName.Contains("Corr")) {
0236         hObsCorr.push_back(new TH2F(*((TH2F*)el->ReadObj())));
0237       }
0238     }
0239   } else {
0240     obsNumbers = observables;
0241     for (unsigned int obs = 0; obs < observables.size(); obs++) {
0242       std::cout << "  ... will read hists and fit for obs " << observables[obs] << " from file " << fileName
0243                 << std::endl;
0244       TString hStitle = "Obs";
0245       hStitle += observables[obs];
0246       hStitle += "_S";
0247       hObsS.push_back(new TH1F(*((TH1F*)fitFile->GetKey(hStitle)->ReadObj())));
0248       TString hBtitle = "Obs";
0249       hBtitle += observables[obs];
0250       hBtitle += "_B";
0251       hObsB.push_back(new TH1F(*((TH1F*)fitFile->GetKey(hBtitle)->ReadObj())));
0252       TString hSBtitle = "Obs";
0253       hSBtitle += observables[obs];
0254       hSBtitle += "_SoverSplusB";
0255       TString fSBtitle = "F_";
0256       fSBtitle += hSBtitle;
0257       hObsSoverSplusB.push_back(new TH1F(*((TH1F*)fitFile->GetKey(hSBtitle)->ReadObj())));
0258       fObsSoverSplusB.push_back(new TF1(*((TF1*)fitFile->GetKey(fSBtitle)->ReadObj())));
0259       for (unsigned int obs2 = obs + 1; obs2 < observables.size(); obs2++) {
0260         TString hCorrtitle = "Corr_Obs";
0261         hCorrtitle += observables[obs];
0262         hCorrtitle += "_Obs";
0263         hCorrtitle += observables[obs2];
0264         hObsCorr.push_back(new TH2F(*((TH2F*)fitFile->GetKey(hCorrtitle)->ReadObj())));
0265       }
0266     }
0267   }
0268
0269   if (readLRplots) {
0270     constructPurity = true;
0271     std::cout << "  ... will LR s and B histos from file " << fileName << std::endl;
0272     hLRtotS = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotS")->ReadObj()));
0273     hLRtotB = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotB")->ReadObj()));
0274     hLRtotSoverSplusB = new TH1F(*((TH1F*)fitFile->GetKey("hLRtotSoverSplusB")->ReadObj()));
0275     fLRtotSoverSplusB = new TF1(*((TF1*)fitFile->GetKey("fLRtotSoverSplusB")->ReadObj()));
0276   }
0277 }
0278
0279 // member function to store all observable plots and fits to a ROOT-file
0280 void LRHelpFunctions::storeToROOTfile(const TString& fname) {
0281   TFile fOut(fname, "RECREATE");
0282   fOut.cd();
0283   for (size_t o = 0; o < hObsS.size(); o++) {
0284     hObsS[o]->Write();
0285     hObsB[o]->Write();
0286     hObsSoverSplusB[o]->Write();
0287     fObsSoverSplusB[o]->Write();
0288   }
0289   int hIndex = 0;
0290   for (size_t o = 0; o < hObsS.size(); o++) {
0291     for (size_t o2 = o + 1; o2 < hObsS.size(); o2++) {
0292       hObsCorr[hIndex]->Write();
0293       ++hIndex;
0294     }
0295   }
0296   if (constructPurity) {
0297     hLRtotS->Write();
0298     hLRtotB->Write();
0299     hLRtotSoverSplusB->Write();
0300     fLRtotSoverSplusB->Write();
0301     hEffvsPur->Write();
0302     hLRValvsPur->Write();
0303     hLRValvsEff->Write();
0304   }
0305   fOut.cd();
0306   fOut.Write();
0307   fOut.Close();
0308 }
0309
0310 // member function to make some simple control plots and store them in a ps-file
0311 void LRHelpFunctions::storeControlPlots(const TString& fname) {
0312   TCanvas c("dummy", "", 1);
0313   c.Print(fname + "[", "landscape");
0314   for (size_t o = 0; o < hObsS.size(); o++) {
0315     TCanvas c2("c2", "", 1);
0316     c2.Divide(2, 1);
0317     c2.cd(1);
0318     hObsS[o]->SetLineColor(2);
0319     if (hObsB[o]->GetMaximum() > hObsS[o]->GetMaximum()) {
0320       hObsB[o]->Draw("hist");
0321       hObsS[o]->Draw("histsame");
0322     } else {
0323       hObsS[o]->Draw("hist");
0324       hObsB[o]->Draw("histsame");
0325     }
0326     c2.cd(2);
0327     hObsSoverSplusB[o]->Draw();
0328     fObsSoverSplusB[o]->Draw("same");
0329     c2.Print(fname, "Landscape");
0330   }
0331
0332   int hIndex = 0;
0333   for (size_t o = 0; o < hObsS.size(); o++) {
0334     for (size_t o2 = o + 1; o2 < hObsS.size(); o2++) {
0335       TCanvas cc("cc", "", 1);
0336       hObsCorr[hIndex]->Draw("box");
0337       TPaveText pt(0.5, 0.87, 0.98, 0.93, "blNDC");
0338       pt.SetFillStyle(1);
0339       pt.SetFillColor(0);
0340       pt.SetBorderSize(0);
0341       TString tcorr = "Corr. of ";
0342       tcorr += (int)(100. * hObsCorr[hIndex]->GetCorrelationFactor());
0343       tcorr += " %";
0344       //TText *text = pt.AddText(tcorr);
0345       pt.Draw("same");
0346       ++hIndex;
0347       cc.Print(fname, "Landscape");
0348     }
0349   }
0350
0351   if (constructPurity) {
0352     TCanvas c3("c3", "", 1);
0353     c3.Divide(2, 1);
0354     c3.cd(1);
0355     hLRtotB->Draw();
0356     hLRtotS->SetLineColor(2);
0357     hLRtotS->Draw("same");
0358     c3.cd(2);
0359     hLRtotSoverSplusB->Draw();
0360     c3.Print(fname, "Landscape");
0361
0362     TCanvas c4("c4", "", 1);
0363     hEffvsPur->Draw("AL*");
0364     c4.Print(fname, "Landscape");
0365
0366     TCanvas c5("c5", "", 1);
0367     hLRValvsPur->Draw("AL*");
0368     c5.Print(fname, "Landscape");
0369
0370     TCanvas c6("c6", "", 1);
0371     hLRValvsEff->Draw("AL*");
0372     c6.Print(fname, "Landscape");
0373   }
0374   c.Print(fname + "]", "landscape");
0375 }
0376
0377 // member function to calculate the LR value, using the S/N definition
0378 double LRHelpFunctions::calcLRval(const std::vector<double>& vals) {
0379   double logLR = 0.;
0380   for (size_t o = 0; o < fObsSoverSplusB.size(); o++) {
0381     double SoverSplusN = fObsSoverSplusB[o]->Eval(vals[o]);
0382     double SoverN = 1. / ((1. / SoverSplusN) - 1.);
0383     logLR += log(SoverN);
0384   }
0385   return logLR;
0386 }
0387
0388 // member function to calculate the LR value, using the definition that was
0389 // used in the P-TDR: S/(S+N)
0390
0391 double LRHelpFunctions::calcPtdrLRval(const std::vector<double>& vals, bool useCorrelation) {
0392   double logLR = 1.;
0393   for (size_t o = 0; o < fObsSoverSplusB.size(); o++) {
0394     double SoverSplusN = fObsSoverSplusB[o]->Eval(vals[o]);
0395     if (SoverSplusN < 0.0001)
0396       SoverSplusN = 0.0001;
0397     if (useCorrelation) {
0398       double corr = 0;
0399       for (size_t j = 0; j < fObsSoverSplusB.size(); j++) {
0400         if (o == j)
0401           ++corr;
0402         else {
0403           TString hcorr = "Corr_Obs";
0404           hcorr += std::min(obsNumbers[o], obsNumbers[j]);
0405           hcorr += "_Obs";
0406           hcorr += std::max(obsNumbers[o], obsNumbers[j]);
0407           for (size_t i = 0; i < hObsCorr.size(); i++) {
0408             if (((TString)(hObsCorr[i]->GetName())) == hcorr)
0409               corr += fabs(hObsCorr[i]->GetCorrelationFactor());
0410           }
0411         }
0412       }
0413       logLR *= pow(SoverSplusN / fObsSoverSplusB[o]->GetMaximum(), 1. / corr);
0414     } else
0415       logLR *= SoverSplusN / fObsSoverSplusB[o]->GetMaximum();
0416   }
0417   //std::cout <<logLR<<std::endl;
0418   return logLR;
0419 }
0420
0421 // member function to fill a signal contribution to the LR histogram
0422 void LRHelpFunctions::fillLRSignalHist(double val, double weight) {
0423   //  std::cout<< "@@@===> LRHelpf Signal LRval = "<< val<<std::endl;
0424   hLRtotS->Fill(val, weight);
0425 }
0426
0427 // member function to fill a background contribution to the LR histogram
0428 void LRHelpFunctions::fillLRBackgroundHist(double val, double weight) {
0429   // std::cout<< "@@@===> LRHelpf Backgroud LRval = "<< val<<std::endl;
0430   hLRtotB->Fill(val, weight);
0431 }
0432
0433 // member function to make and fit the purity vs. LRval histogram, and the purity vs. efficiency graph
0434 void LRHelpFunctions::makeAndFitPurityHists() {
0435   for (int b = 0; b <= hLRtotS->GetNbinsX(); b++) {
0436     float Sint = hLRtotS->Integral(b, hLRtotS->GetNbinsX() + 1);
0437     float Bint = hLRtotB->Integral(b, hLRtotB->GetNbinsX() + 1);
0438     if (Sint + Bint > 0) {
0439       hLRtotSoverSplusB->SetBinContent(b, 1. * Sint / (Sint + Bint));
0440       hLRtotSoverSplusB->SetBinError(
0441           b, sqrt((pow(Sint * sqrt(Bint), 2) + pow(Bint * sqrt(Sint), 2))) / pow((Sint + Bint), 2));
0442     }
0443   }
0444
0445   hLRtotS->GetXaxis()->SetTitle("Combined LR ratio");
0446   hLRtotB->GetXaxis()->SetTitle("Combined LR ratio");
0447   hLRtotSoverSplusB->GetXaxis()->SetTitle("Cut on Combined LR");
0448   hLRtotSoverSplusB->GetYaxis()->SetTitle("Purity");
0449
0450   hLRtotSoverSplusB->Fit(fLRtotSoverSplusB->GetName(), "RQ");
0451   double totSignal = hLRtotS->Integral(0, hLRtotS->GetNbinsX() + 1);
0452   double Eff[200], Pur[200], LRVal[200];
0453   if (hLRtotS->GetNbinsX() > 200) {
0454     std::cout << "Number of bins of LR histograms can not execeed 200!" << std::endl;
0455     return;
0456   }
0457   for (int cut = 0; (cut <= hLRtotS->GetNbinsX()) && (cut < 200); cut++) {
0458     double LRcutVal = hLRtotS->GetBinLowEdge(cut);
0459     Eff[cut] = hLRtotS->Integral(cut, hLRtotS->GetNbinsX() + 1) / totSignal;
0460     Pur[cut] = fLRtotSoverSplusB->Eval(LRcutVal);
0461     LRVal[cut] = LRcutVal;
0462   }
0463   hEffvsPur = new TGraph(hLRtotS->GetNbinsX(), Eff, Pur);
0464   hEffvsPur->SetName("hEffvsPur");
0465   hEffvsPur->SetTitle("");
0466   hEffvsPur->GetXaxis()->SetTitle((TString)("Efficiency of cut on log combined LR"));
0467   hEffvsPur->GetYaxis()->SetTitle((TString)("Purity"));
0468   hEffvsPur->GetYaxis()->SetRangeUser(0, 1.1);
0469
0470   hLRValvsPur = new TGraph(hLRtotS->GetNbinsX(), LRVal, Pur);
0471   hLRValvsPur->SetName("hLRValvsPur");
0472   hLRValvsPur->SetTitle("");
0473   hLRValvsPur->GetXaxis()->SetTitle((TString)("Cut on the log combined LR value"));
0474   hLRValvsPur->GetYaxis()->SetTitle((TString)("Purity"));
0475   hLRValvsPur->GetYaxis()->SetRangeUser(0, 1.1);
0476
0477   hLRValvsEff = new TGraph(hLRtotS->GetNbinsX(), LRVal, Eff);
0478   hLRValvsEff->SetName("hLRValvsEff");
0479   hLRValvsEff->SetTitle("");
0480   hLRValvsEff->GetXaxis()->SetTitle((TString)("Cut on the log combined LR value"));
0481   hLRValvsEff->GetYaxis()->SetTitle((TString)("Efficiency of cut on log combined LR"));
0482   hLRValvsEff->GetYaxis()->SetRangeUser(0, 1.1);
0483 }
0484
0485 // member function to calculate the probability for signal given a logLR value
0486 double LRHelpFunctions::calcProb(double logLR) { return fLRtotSoverSplusB->Eval(logLR); }
0487
0488 // member to check if a certain S/S+B fit function is read from the root-file
0489 bool LRHelpFunctions::obsFitIncluded(int o) {
0490   bool included = false;
0491   TString obs = "_Obs";
0492   obs += o;
0493   obs += "_";
0494   for (size_t f = 0; f < fObsSoverSplusB.size(); f++) {
0495     if (((TString)(fObsSoverSplusB[f]->GetName())).Contains(obs))
0496       included = true;
0497   }
0498   return included;
0499 }
0500
0501 void LRHelpFunctions::setXlabels(const std::vector<std::string>& xLabels) {
0502   if (hObsS.size() != xLabels.size()) {
0503     std::cout << "LRHelpFunctions::setXlabels: Number of labels (" << xLabels.size()
0504               << ") does not match number of obervables(" << hObsS.size() << ").\n";
0505     return;
0506   }
0507   for (size_t i = 0; i < hObsS.size(); ++i) {
0508     hObsS[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
0509     hObsB[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
0510     hObsSoverSplusB[i]->GetXaxis()->SetTitle(TString(xLabels[i]));
0511     hObsSoverSplusB[i]->GetYaxis()->SetTitle(TString("S/(S+B)"));
0512   }
0513 }
0514
0515 void LRHelpFunctions::setYlabels(const std::vector<std::string>& yLabels) {
0516   if (hObsS.size() != yLabels.size()) {
0517     std::cout << "LRHelpFunctions::setYlabels: Number of labels (" << yLabels.size()
0518               << ") does not match number of obervables(" << hObsS.size() << ").\n";
0519     return;
0520   }
0521   for (size_t i = 0; i < hObsS.size(); ++i) {
0522     hObsS[i]->GetYaxis()->SetTitle(TString(yLabels[i]));
0523     hObsB[i]->GetYaxis()->SetTitle(TString(yLabels[i]));
0524   }
0525 }
0526
0527 void LRHelpFunctions::singlePlot(const TString& fname, int obsNbr, const TString& extension) {
0528   if (!obsFitIncluded(obsNbr))
0529     return;
0530
0531   TStyle* tdrStyle = gROOT->GetStyle("tdrStyle");
0532   tdrStyle->SetOptFit(0);
0533   tdrStyle->SetOptStat(0);
0534   tdrStyle->SetOptTitle(0);
0535   //   tdrStyle->SetPadTopMargin(0.01);
0536   //   tdrStyle->SetPadBottomMargin(0.01);
0537   //   tdrStyle->SetPadLeftMargin(0.01);
0538   //   tdrStyle->SetPadRightMargin(0.01);
0539
0540   TCanvas c2("c2", "", 600, 300);
0541   c2.SetTopMargin(0.01);
0542   c2.SetBottomMargin(0.01);
0543   c2.SetLeftMargin(0.01);
0544   c2.SetRightMargin(0.01);
0545   std::cout << fname << std::endl;
0546   c2.Divide(2, 1);
0547
0548   TString obs = "Obs";
0549   obs += obsNbr;
0550   obs += "_";
0551   for (size_t o = 0; o < hObsB.size(); ++o) {
0552     if (((TString)(hObsB[o]->GetName())).Contains(obs)) {
0553       c2.cd(1);
0554       hObsS[o]->SetLineColor(2);
0555       if (hObsB[o]->GetMaximum() > hObsS[o]->GetMaximum()) {
0556         hObsB[o]->Draw("hist");
0557         hObsS[o]->Draw("histsame");
0558       } else {
0559         hObsS[o]->Draw("hist");
0560         hObsB[o]->Draw("histsame");
0561       }
0562       c2.cd(2);
0563
0564       hObsSoverSplusB[o]->Draw();
0565       fObsSoverSplusB[o]->Draw("same");
0566     }
0567   }
0568   std::cout << fname + "." + extension << std::endl;
0569   c2.Print(fname + "." + extension);
0570 }
0571
0572 void LRHelpFunctions::purityPlot(const TString& fname, const TString& extension) {
0573   TStyle* tdrStyle = gROOT->GetStyle("tdrStyle");
0574   tdrStyle->SetOptFit(0);
0575   tdrStyle->SetOptStat(0);
0576   tdrStyle->SetOptTitle(0);
0577   //   tdrStyle->SetPadTopMargin(0.01);
0578   //   tdrStyle->SetPadBottomMargin(0.01);
0579   //   tdrStyle->SetPadLeftMargin(0.01);
0580   //   tdrStyle->SetPadRightMargin(0.01);
0581
0582   TCanvas c2("c2", "", 600, 300);
0583   c2.SetTopMargin(0.01);
0584   c2.SetBottomMargin(0.01);
0585   c2.SetLeftMargin(0.01);
0586   c2.SetRightMargin(0.01);
0587   std::cout << fname << std::endl;
0588   c2.Divide(2, 1);
0589
0590   c2.cd(1);
0591   hLRValvsPur->Draw("AP");
0592   c2.cd(2);
0593   hEffvsPur->Draw("AP");
0594   c2.Print(fname + "Purity." + extension);
0595
0596   hLRtotS->GetXaxis()->SetNdivisions(505);
0597   hLRtotB->GetXaxis()->SetNdivisions(505);
0598   TCanvas c3("c2", "", 300, 300);
0599   //   c3.SetTopMargin(0.01);
0600   //   c3.SetBottomMargin(0.01);
0601   //   c3.SetLeftMargin(0.01);
0602   //   c3.SetRightMargin(0.01);
0603   hLRtotS->GetXaxis()->SetTitle("Combined LR");
0604   hLRtotB->GetXaxis()->SetTitle("Combined LR");
0605
0606   hLRtotB->Draw();
0607   hLRtotS->SetLineColor(2);
0608   hLRtotS->Draw("same");
0609   c3.Print(fname + "Dist." + extension);
0610
0611   TCanvas c5("c3", "", 900, 600);
0612   c5.Divide(2, 2);
0613   c5.cd(1);
0614   hLRtotB->Draw();
0615   hLRtotS->SetLineColor(2);
0616   hLRtotS->Draw("same");
0617   c5.cd(3);
0618   hLRValvsEff->Draw("AP");
0619   c5.cd(2);
0620   hEffvsPur->Draw("AP");
0621   c5.cd(4);
0622   hLRtotSoverSplusB->Draw();
0623   c5.Print(fname + "all." + extension);
0624 }