RecoB/src/Tools.cc

0001 #include "DQMOffline/RecoB/interface/Tools.h"
0002
0003 #include "TROOT.h"
0004 #include "TSystem.h"
0005 #include "TStyle.h"
0006 #include <cmath>
0007
0008 using namespace std;
0009 using namespace RecoBTag;
0010
0011 //
0012 //
0013 // TOOLS
0014 //
0015 //
0016
0017 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0018 double RecoBTag::HistoBinWidth(const TH1F* theHisto, const int& iBin) {
0019   const int& nBins = theHisto->GetSize();  // includes underflow/overflow
0020   // return 0.0 , if invalid bin
0021   if (iBin < 0 || iBin >= nBins)
0022     return 0.0;
0023   // return first binwidth, if underflow bin
0024   if (iBin == 0)
0025     return theHisto->GetBinWidth(1);
0026   // return last real binwidth, if overflow bin
0027   if (iBin == nBins - 1)
0028     return theHisto->GetBinWidth(nBins - 2);
0029   // return binwidth from histo, if within range
0030   return theHisto->GetBinWidth(iBin);
0031 }
0032 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0033
0034 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0035 double RecoBTag::IntegrateHistogram(const TH1F* theHisto) {
0036   // include underflow and overflow: assign binwidth of first/last bin to them!!
0037   // integral = sum ( entry_i * binwidth_i )
0038   //
0039   double histoIntegral = 0.0;
0040   const int& nBins = theHisto->GetSize();
0041   //
0042   // loop over bins:
0043   // bin 0       : underflow
0044   // bin nBins-1 : overflow
0045   for (int iBin = 0; iBin != nBins; ++iBin) {
0046     const double& binWidth = HistoBinWidth(theHisto, iBin);
0047     histoIntegral += (*theHisto)[iBin] * binWidth;
0048   }
0049   //
0050   return histoIntegral;
0051 }
0052 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0053
0054 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0055 void RecoBTag::HistoToNormalizedArrays(const TH1F* theHisto,
0056                                        TArrayF& theNormalizedArray,
0057                                        TArrayF& theLeftOfBinArray,
0058                                        TArrayF& theBinWidthArray) {
0059   const int& nBins = theHisto->GetSize();
0060
0061   // check that all arrays/histo have the same size
0062   if (nBins == theNormalizedArray.GetSize() && nBins == theLeftOfBinArray.GetSize() &&
0063       nBins == theBinWidthArray.GetSize()) {
0064     const double& histoIntegral = IntegrateHistogram(theHisto);
0065
0066     for (int iBin = 0; iBin != nBins; ++iBin) {
0067       theNormalizedArray[iBin] = (*theHisto)[iBin] / histoIntegral;
0068       theLeftOfBinArray[iBin] = theHisto->GetBinLowEdge(iBin);
0069       theBinWidthArray[iBin] = HistoBinWidth(theHisto, iBin);
0070     }
0071
0072   } else {
0073     cout << "============>>>>>>>>>>>>>>>>" << endl
0074          << "============>>>>>>>>>>>>>>>>" << endl
0075          << "============>>>>>>>>>>>>>>>>" << endl
0076          << "============>>>>>>>>>>>>>>>>" << endl
0077          << "============>>>>>>>>>>>>>>>> HistoToNormalizedArrays failed: not equal sizes of all arrays!!" << endl
0078          << "============>>>>>>>>>>>>>>>>" << endl
0079          << "============>>>>>>>>>>>>>>>>" << endl
0080          << "============>>>>>>>>>>>>>>>>" << endl
0081          << "============>>>>>>>>>>>>>>>>" << endl;
0082   }
0083 }
0084 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0085
0086 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0087 double RecoBTag::IntegrateArray(const TArrayF& theArray, const TArrayF& theBinWidth) {
0088   double arrayIntegral = 0.0;
0089   const int& nBins = theArray.GetSize();
0090   //
0091   for (int iBin = 0; iBin != nBins; ++iBin) {
0092     arrayIntegral += theArray[iBin] * theBinWidth[iBin];
0093   }
0094   //
0095   return arrayIntegral;
0096 }
0097 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0098
0099 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0100 void RecoBTag::PrintCanvasHistos(TCanvas* canvas,
0101                                  const std::string& psFile,
0102                                  const std::string& epsFile,
0103                                  const std::string& gifFile) {
0104   //
0105   //
0106   // to create gif in 'batch mode' (non-interactive) see
0107   // http://root.cern.ch/cgi-bin/print_hit_bold.pl/root/roottalk/roottalk00/0402.html?gifbatch#first_hit
0108   //
0109   // ROOT 4 can do it!!??
0110   //
0111   // if string = "" don't print to corresponding file
0112   //
0113   if (!psFile.empty())
0114     canvas->Print(psFile.c_str());
0115   if (!epsFile.empty())
0116     canvas->Print(epsFile.c_str(), "eps");
0117   // if in batch: use a converter tool
0118   const std::string& rootVersion(gROOT->GetVersion());
0119   const bool& rootCanGif = rootVersion.find('4') == 0 || rootVersion.find('5') == 0;
0120   if (!gifFile.empty()) {
0121     if (!(gROOT->IsBatch()) || rootCanGif) {  // to find out if running in batch mode
0122       cout << "--> Print directly gif!" << endl;
0123       canvas->Print(gifFile.c_str(), "gif");
0124     } else {
0125       if (!epsFile.empty()) {  // eps file must have been created before
0126         cout << "--> Print gif via scripts!" << endl;
0127         const std::string& executeString1 = "pstopnm -ppm -xborder 0 -yborder 0 -portrait " + epsFile;
0128         gSystem->Exec(executeString1.c_str());
0129         const std::string& ppmFile = epsFile + "001.ppm";
0130         const std::string& executeString2 = "ppmtogif " + ppmFile + " > " + gifFile;
0131         gSystem->Exec(executeString2.c_str());
0132         const std::string& executeString3 = "rm " + ppmFile;
0133         gSystem->Exec(executeString3.c_str());  // delete the intermediate file
0134       }
0135     }
0136   }
0137 }
0138 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0139
0140 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0141 TObjArray RecoBTag::getHistArray(TFile* histoFile, const std::string& baseName) {
0142   //
0143   // return the TObjArray built from the basename
0144   //
0145   //
0146   TObjArray histos(3);  // reserve 3
0147   //
0148   const std::string nameB(baseName + "B");
0149   const std::string nameC(baseName + "C");
0150   const std::string nameDUSG(baseName + "DUSG");
0151   //
0152   histos.Add((TH1F*)histoFile->Get(nameB.c_str()));
0153   histos.Add((TH1F*)histoFile->Get(nameC.c_str()));
0154   histos.Add((TH1F*)histoFile->Get(nameDUSG.c_str()));
0155   //
0156   return histos;
0157 }
0158 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0159
0160 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0161 std::string RecoBTag::flavour(const int& flav) {
0162   switch (flav) {
0163     case 1:
0164       return "d";
0165     case 2:
0166       return "u";
0167     case 3:
0168       return "s";
0169     case 4:
0170       return "c";
0171     case 5:
0172       return "b";
0173     case 21:
0174       return "g";
0175     default:
0176       return "";
0177   }
0178 }
0179 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0180
0181 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0182 bool RecoBTag::flavourIsD(const int& flav) { return flav == 1; }
0183 bool RecoBTag::flavourIsU(const int& flav) { return flav == 2; }
0184 bool RecoBTag::flavourIsS(const int& flav) { return flav == 3; }
0185 bool RecoBTag::flavourIsC(const int& flav) { return flav == 4; }
0186 bool RecoBTag::flavourIsB(const int& flav) { return flav == 5; }
0187 bool RecoBTag::flavourIsG(const int& flav) { return flav == 21; }
0188
0189 bool RecoBTag::flavourIsDUS(const int& flav) { return (flavourIsD(flav) || flavourIsU(flav) || flavourIsS(flav)); }
0190 bool RecoBTag::flavourIsDUSG(const int& flav) { return (flavourIsDUS(flav) || flavourIsG(flav)); }
0191
0192 bool RecoBTag::flavourIsNI(const int& flav) {
0193   return !(flavourIsD(flav) || flavourIsU(flav) || flavourIsS(flav) || flavourIsC(flav) || flavourIsB(flav) ||
0194            flavourIsG(flav));
0195 }
0196 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0197
0198 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0199 int RecoBTag::checkCreateDirectory(const std::string& directory) {
0200   cout << "====>>>> ToolsC:checkCreateDirectory() : " << endl;
0201   int exists = gSystem->Exec(("ls -d " + directory).c_str());
0202   // create it if it doesn't exist
0203   if (exists != 0) {
0204     cout << "====>>>> ToolsC:checkCreateDirectory() : The directory does not exist : " << directory << endl;
0205     cout << "====>>>> ToolsC:checkCreateDirectory() : I'll try to create it" << endl;
0206     const int& create = gSystem->Exec(("mkdir " + directory).c_str());
0207     if (create != 0) {
0208       cout << "====>>>> ToolsC:checkCreateDirectory() : Creation of directory failed : " << directory << endl
0209            << "====>>>> ToolsC:checkCreateDirectory() : Please check your write permissions!" << endl;
0210     } else {
0211       cout << "====>>>> ToolsC:checkCreateDirectory() : Creation of directory successful!" << endl;
0212       // check again if it exists now
0213       cout << "====>>>> ToolsC:checkCreateDirectory() : " << endl;
0214       exists = gSystem->Exec(("ls -d " + directory).c_str());
0215       if (exists != 0)
0216         cout << "ToolsC:checkCreateDirectory() : However, it still doesn't exist!?" << endl;
0217     }
0218   }
0219   cout << endl;
0220   return exists;
0221 }
0222 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0223
0224 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0225 int RecoBTag::findBinClosestYValue(const TH1F* histo, const float& yVal, const float& yLow, const float& yHigh) {
0226   //
0227   // Find the bin in a 1-dim. histogram which has its y-value closest to
0228   // the given value yVal where the value yVal has to be in the range yLow < yVal < yHigh.
0229   // If it is outside this range the corresponding bin number is returned as negative value.
0230   // Currently, there is no protection if there are many bins with the same value!
0231   // The user has to take care to interpret the output correctly.
0232   //
0233
0234   // init
0235   const int& nBins = histo->GetNbinsX() - 2;  // -2 because we don't include under/overflow alos in this loop
0236   int iBinClosestInit = 0;
0237   // init start value properly: must avoid that the real one is not filled
0238   float yClosestInit;
0239   //
0240   const float& maxInHisto = histo->GetMaximum();
0241   const float& minInHisto = histo->GetMinimum();
0242   //
0243   // if yVal is smaller than max -> take any value well above the maximum
0244   if (yVal <= maxInHisto) {
0245     yClosestInit = maxInHisto + 1;
0246   } else {
0247     // if yVal is greater than max value -> take a value < minimum
0248     yClosestInit = minInHisto - 1.0;
0249   }
0250
0251   int iBinClosest = iBinClosestInit;
0252   float yClosest = yClosestInit;
0253
0254   // loop over bins of histogram
0255   for (int iBin = 1; iBin <= nBins; ++iBin) {
0256     const float& yBin = histo->GetBinContent(iBin);
0257     if (fabs(yBin - yVal) < fabs(yClosest - yVal)) {
0258       yClosest = yBin;
0259       iBinClosest = iBin;
0260     }
0261   }
0262
0263   // check if in interval
0264   if (yClosest < yLow || yClosest > yHigh) {
0265     iBinClosest *= -1;
0266   }
0267
0268   // check that not the initialization bin (would mean that init value was the closest)
0269   if (iBinClosest == iBinClosestInit) {
0270     cout << "====>>>> ToolsC=>findBinClosestYValue() : WARNING: returned bin is the initialization bin!!" << endl;
0271   }
0272
0273   return iBinClosest;
0274 }
0275 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0276
0277 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0278 vector<int> RecoBTag::findBinClosestYValueAtFixedZ(const TH2F* histoY,
0279                                                    const float& yVal,
0280                                                    const float& yLow,
0281                                                    const float& yHigh,
0282                                                    const TH2F* histoZ,
0283                                                    const vector<double>& zVal) {
0284   //
0285   // Find the bin in a 2-dim. histogram which has its y-value closest to
0286   // the given value yVal where the value yVal has to be in the range yLow < yVal < yHigh.
0287   // If it is outside this range the corresponding bin number is returned as negative value.
0288   // The bin should also correspond to a value of z=zVal within the same precision as yVal.
0289   // Currently, there is no protection if there are many bins with the same value!
0290   // The user has to take care to interpret the output correctly.
0291   //
0292
0293   // init
0294   const int& nBinsX = histoY->GetNbinsX() - 2;  // -2 because we don't include under/overflow alos in this loop
0295   const int& nBinsY = histoY->GetNbinsY() - 2;  // -2 because we don't include under/overflow alos in this loop
0296   int iBinClosestInit = 0;
0297   // init start value properly: must avoid that the real one is not filled
0298   vector<float> yClosestInit(zVal.size());
0299   //
0300   const float& maxInHisto = histoY->GetMaximum();
0301   const float& minInHisto = histoY->GetMinimum();
0302   //
0303   // if yVal is smaller than max -> take any value well above the maximum
0304   for (unsigned int i = 0; i < yClosestInit.size(); i++) {
0305     if (yVal <= maxInHisto) {
0306       yClosestInit[i] = maxInHisto + 1;
0307     } else {
0308       // if yVal is greater than max value -> take a value < minimum
0309       yClosestInit[i] = minInHisto - 1.0;
0310     }
0311   }
0312
0313   vector<int> iBinClosest(zVal.size(), iBinClosestInit);
0314   vector<float> yClosest(yClosestInit);
0315
0316   // loop over bins of histogram
0317   for (int iBinX = 1; iBinX <= nBinsX; ++iBinX) {
0318     for (int iBinY = 1; iBinY <= nBinsY; ++iBinY) {
0319       const float& yBin = histoY->GetBinContent(iBinX, iBinY);
0320       for (unsigned int i = 0; i < zVal.size(); i++) {
0321         if (fabs(yBin - yVal) < fabs(yClosest[i] - yVal)) {
0322           const float& zLow = zVal[i] - (yVal - yLow);
0323           const float& zHigh = zVal[i] + (yHigh - yVal);
0324           const float& zBin = histoZ->GetBinContent(iBinX, iBinY);
0325           if (zBin < zLow || zBin > zHigh)
0326             continue;
0327           yClosest[i] = yBin;
0328           iBinClosest[i] = histoY->GetBin(iBinX, iBinY);
0329         }
0330       }
0331     }
0332   }
0333   // check if in interval
0334   for (unsigned int i = 0; i < yClosest.size(); i++) {
0335     if (yClosest[i] < yLow || yClosest[i] > yHigh)
0336       iBinClosest[i] *= -1;
0337   }
0338
0339   return iBinClosest;
0340 }
0341 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0342
0343 TStyle* RecoBTag::setTDRStyle() {
0344   TStyle* tdrStyle = new TStyle("tdrStyle", "Style for P-TDR");
0345
0346   // For the canvas:
0347   tdrStyle->SetCanvasBorderMode(0);
0348   tdrStyle->SetCanvasColor(kWhite);
0349   tdrStyle->SetCanvasDefH(600);  //Height of canvas
0350   tdrStyle->SetCanvasDefW(600);  //Width of canvas
0351   tdrStyle->SetCanvasDefX(0);    //POsition on screen
0352   tdrStyle->SetCanvasDefY(0);
0353
0354   // For the Pad:
0355   tdrStyle->SetPadBorderMode(0);
0356   // tdrStyle->SetPadBorderSize(Width_t size = 1);
0357   tdrStyle->SetPadColor(kWhite);
0358   tdrStyle->SetPadGridX(false);
0359   tdrStyle->SetPadGridY(false);
0360   tdrStyle->SetGridColor(0);
0361   tdrStyle->SetGridStyle(3);
0362   tdrStyle->SetGridWidth(1);
0363
0364   // For the frame:
0365   tdrStyle->SetFrameBorderMode(0);
0366   tdrStyle->SetFrameBorderSize(1);
0367   tdrStyle->SetFrameFillColor(0);
0368   tdrStyle->SetFrameFillStyle(0);
0369   tdrStyle->SetFrameLineColor(1);
0370   tdrStyle->SetFrameLineStyle(1);
0371   tdrStyle->SetFrameLineWidth(1);
0372
0373   // For the histo:
0374   // tdrStyle->SetHistFillColor(1);
0375   // tdrStyle->SetHistFillStyle(0);
0376   tdrStyle->SetHistLineColor(1);
0377   tdrStyle->SetHistLineStyle(0);
0378   tdrStyle->SetHistLineWidth(1);
0379   // tdrStyle->SetLegoInnerR(Float_t rad = 0.5);
0380   // tdrStyle->SetNumberContours(Int_t number = 20);
0381
0382   tdrStyle->SetEndErrorSize(15);
0383   //   tdrStyle->SetErrorMarker(20);
0384   tdrStyle->SetErrorX(1);
0385
0386   tdrStyle->SetMarkerStyle(21);
0387   tdrStyle->SetMarkerSize(1.);
0388
0389   //For the fit/function:
0390   tdrStyle->SetOptFit(0);
0391   tdrStyle->SetFitFormat("5.4g");
0392   tdrStyle->SetFuncColor(2);
0393   tdrStyle->SetFuncStyle(1);
0394   tdrStyle->SetFuncWidth(1);
0395
0396   //For the date:
0397   tdrStyle->SetOptDate(0);
0398   // tdrStyle->SetDateX(Float_t x = 0.01);
0399   // tdrStyle->SetDateY(Float_t y = 0.01);
0400
0401   // For the statistics box:
0402   tdrStyle->SetOptFile(1111);
0403   tdrStyle->SetOptStat(0);  // To display the mean and RMS:   SetOptStat("mr");
0404   tdrStyle->SetStatColor(kWhite);
0405   tdrStyle->SetStatFont(42);
0406   tdrStyle->SetStatFontSize(0.025);
0407   tdrStyle->SetStatTextColor(1);
0408   tdrStyle->SetStatFormat("6.4g");
0409   tdrStyle->SetStatBorderSize(1);
0410   tdrStyle->SetStatH(0.2);
0411   tdrStyle->SetStatW(0.15);
0412   // tdrStyle->SetStatStyle(Style_t style = 1001);
0413   // tdrStyle->SetStatX(Float_t x = 0);
0414   // tdrStyle->SetStatY(Float_t y = 0);
0415
0416   // Margins:
0417   tdrStyle->SetPadTopMargin(0.05);
0418   tdrStyle->SetPadBottomMargin(0.13);
0419   tdrStyle->SetPadLeftMargin(0.16);
0420   tdrStyle->SetPadRightMargin(0.02);
0421
0422   // For the Global title:
0423
0424   tdrStyle->SetOptTitle(0);
0425   tdrStyle->SetTitleW(0.8);  // Set the width of the title box
0426
0427   tdrStyle->SetTitleFont(42);
0428   tdrStyle->SetTitleColor(1);
0429   tdrStyle->SetTitleTextColor(1);
0430   tdrStyle->SetTitleFillColor(10);
0431   tdrStyle->SetTitleFontSize(0.05);
0432   // tdrStyle->SetTitleH(0); // Set the height of the title box
0433   // tdrStyle->SetTitleX(0); // Set the position of the title box
0434   // tdrStyle->SetTitleY(0.985); // Set the position of the title box
0435   // tdrStyle->SetTitleStyle(Style_t style = 1001);
0436   // tdrStyle->SetTitleBorderSize(2);
0437
0438   // For the axis titles:
0439
0440   tdrStyle->SetTitleColor(1, "XYZ");
0441   tdrStyle->SetTitleFont(42, "XYZ");
0442   tdrStyle->SetTitleSize(0.06, "XYZ");
0443   // tdrStyle->SetTitleXSize(Float_t size = 0.02); // Another way to set the size?
0444   // tdrStyle->SetTitleYSize(Float_t size = 0.02);
0445   tdrStyle->SetTitleXOffset(0.75);
0446   tdrStyle->SetTitleYOffset(0.75);
0447   // tdrStyle->SetTitleOffset(1.1, "Y"); // Another way to set the Offset
0448
0449   // For the axis labels:
0450
0451   tdrStyle->SetLabelColor(1, "XYZ");
0452   tdrStyle->SetLabelFont(42, "XYZ");
0453   tdrStyle->SetLabelOffset(0.007, "XYZ");
0454   tdrStyle->SetLabelSize(0.05, "XYZ");
0455
0456   // For the axis:
0457
0458   tdrStyle->SetAxisColor(1, "XYZ");
0459   tdrStyle->SetStripDecimals(kTRUE);
0460   tdrStyle->SetTickLength(0.03, "XYZ");
0461   tdrStyle->SetNdivisions(510, "XYZ");
0462   tdrStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
0463   tdrStyle->SetPadTickY(1);
0464
0465   // Change for log plots:
0466   tdrStyle->SetOptLogx(0);
0467   tdrStyle->SetOptLogy(0);
0468   tdrStyle->SetOptLogz(0);
0469
0470   // Postscript options:
0471   tdrStyle->SetPaperSize(21., 28.);
0472   //  tdrStyle->SetPaperSize(20.,20.);
0473   // tdrStyle->SetLineScalePS(Float_t scale = 3);
0474   // tdrStyle->SetLineStyleString(Int_t i, const char* text);
0475   // tdrStyle->SetHeaderPS(const char* header);
0476   // tdrStyle->SetTitlePS(const char* pstitle);
0477
0478   // tdrStyle->SetBarOffset(Float_t baroff = 0.5);
0479   // tdrStyle->SetBarWidth(Float_t barwidth = 0.5);
0480   // tdrStyle->SetPaintTextFormat(const char* format = "g");
0481   // tdrStyle->SetPalette(Int_t ncolors = 0, Int_t* colors = 0);
0482   // tdrStyle->SetTimeOffset(Double_t toffset);
0483   // tdrStyle->SetHistMinimumZero(kTRUE);
0484
0485   tdrStyle->cd();
0486   return tdrStyle;
0487 }
0488 // tdrGrid: Turns the grid lines on (true) or off (false)
0489
0490 void RecoBTag::tdrGrid(const bool& gridOn) {
0491   TStyle* tdrStyle = setTDRStyle();
0492   tdrStyle->SetPadGridX(gridOn);
0493   tdrStyle->SetPadGridY(gridOn);
0494   tdrStyle->cd();
0495 }
0496
0497 // fixOverlay: Redraws the axis
0498
0499 void RecoBTag::fixOverlay() { gPad->RedrawAxis(); }
0500
0501 string RecoBTag::itos(const int& i)  // convert int to string
0502 {
0503   ostringstream s;
0504   s << i;
0505   return s.str();
0506 }