diff --git a/Install.py b/Install.py new file mode 100644 index 0000000..f63c11f --- /dev/null +++ b/Install.py @@ -0,0 +1,47 @@ +#!/usr/bin/python + +import argparse +import subprocess + +def Install() : + """ + Setup the environment. + """ + print '---------------------' + print 'Setup the environment' + + command = 'make -j4' + print '>', command + subprocess.call(command,shell=True,cwd='../') + + command = 'make install' + print '>', command + subprocess.call(command,shell=True,cwd='../') + + command = '/afs/cern.ch/project/eos/installation/0.3.15/bin/eos.select -b fuse mount ~/eos' + print '>', command + subprocess.call(command,shell=True) + + command = 'chmod +x run.sh' + print '>', command + subprocess.call(command,shell=True,cwd='../') + + command = './run.sh bash' + print '>', command + # subprocess.Popen(command,shell=True,cwd='./KeplerDev_v3r0/') + subprocess.call([command,'bash'],shell=True,cwd='../') + + print '---------------------' + + return + + + +############### +# +# Main function +# + +if __name__ == "__main__" : + + Install() diff --git a/TbUT/options/Makefile b/TbUT/options/Makefile new file mode 100755 index 0000000..ae6e18c --- /dev/null +++ b/TbUT/options/Makefile @@ -0,0 +1,31 @@ +# +CC=$(CXX) +glib_cflags=$(shell pkg-config --cflags glib-2.0 gio-2.0) +glib_libs=$(shell pkg-config --libs glib-2.0 gio-2.0) + +ROOTC=$(shell root-config --cflags) +ROOTL=$(shell root-config --libs) +OPT=-g -fno-inline #-std=c++11 +CppFLAGS=$(OPT) -I. $(glib_cflags) +CXXFLAGS=-fPIC $(CppFLAGS) + + +UNAME_S := $(shell uname -s) +ifeq ($(UNAME_S),Linux) + SO=so + SO_FLAGS=-shared + CXXFLAGS += -D LINUX +endif +ifeq ($(UNAME_S),Darwin) + SO=dylib + SO_FLAGS=-dynamiclib -undefined dynamic_lookup -install_name @rpath/$@ + CXXFLAGS += -D OSX +endif + +all: TbUTLandauFits + +TbUTLandauFits: TbUTLandauFits.C + c++ -I$(OPT) $(CXXFLAGS) $(ROOTC) -o $@ $^ $(LDLIBS) $(ROOTL) -lRooFit -lRooFitCore -lMinuit -lRooStats $(gliblibs) + +clean: + rm -f *.o TbUTLandauFits.C diff --git a/TbUT/options/TbUTLandauFits.C b/TbUT/options/TbUTLandauFits.C new file mode 100644 index 0000000..0e6d58a --- /dev/null +++ b/TbUT/options/TbUTLandauFits.C @@ -0,0 +1,395 @@ +//----------------------------------------------------------------------- +// +// Convoluted Landau and Gaussian Fitting Function +// (using ROOT's Landau and Gauss functions) +// +// Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at) +// Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and +// Markus Friedl (Markus.Friedl@cern.ch) +// +// to execute this example, do: +// root > .x langaus.C +// or +// root > .x langaus.C++ +// +//----------------------------------------------------------------------- + +#include "../../Tools/Lib.C" +#include "../../Tools/lhcbStyle.C" +#include "../../Tools/Style.C" +#include "../../Tools/FindAndReplaceString.C" + +Double_t langaufun(Double_t *x, Double_t *par) { + + //Fit parameters: + //par[0]=Width (scale) parameter of Landau density + //par[1]=Most Probable (MP, location) parameter of Landau density + //par[2]=Total area (integral -inf to inf, normalization constant) + //par[3]=Width (sigma) of convoluted Gaussian function + // + //In the Landau distribution (represented by the CERNLIB approximation), + //the maximum is located at x=-0.22278298 with the location parameter=0. + //This shift is corrected within this function, so that the actual + //maximum is identical to the MP parameter. + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // number of convolution steps + Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx; + Double_t mpc; + Double_t fland; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp-xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for(i=1.0; i<=np/2; i++) { + xx = xlow + (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + + xx = xupp - (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + } + + return (par[2] * step * sum * invsq2pi / par[3]); +} + + + +TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) +{ + // Once again, here are the Landau * Gaussian parameters: + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + // + // Variables for langaufit call: + // his histogram to fit + // fitrange[2] lo and hi boundaries of fit range + // startvalues[4] reasonable start values for the fit + // parlimitslo[4] lower parameter limits + // parlimitshi[4] upper parameter limits + // fitparams[4] returns the final fit parameters + // fiterrors[4] returns the final fit errors + // ChiSqr returns the chi square + // NDF returns ndf + + Int_t i; + Char_t FunName[100]; + + sprintf(FunName,"Fitfcn_%s",his->GetName()); + + TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); + if (ffitold) delete ffitold; + + TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); + ffit->SetParameters(startvalues); + ffit->SetParNames("Width","MP","Area","GSigma"); + + for (i=0; i<4; i++) { + ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); + } + + his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot + + ffit->GetParameters(fitparams); // obtain fit parameters + for (i=0; i<4; i++) { + fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors + } + ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 + NDF[0] = ffit->GetNDF(); // obtain ndf + + return (ffit); // return fit function + +} + + +Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) { + + // Seaches for the location (x value) at the maximum of the + // Landau-Gaussian convolute and its full width at half-maximum. + // + // The search is probably not very efficient, but it's a first try. + + Double_t p,x,fy,fxr,fxl; + Double_t step; + Double_t l,lold; + Int_t i = 0; + Int_t MAXCALLS = 10000; + + + // Search for maximum + + p = params[1] - 0.1 * params[0]; + step = 0.05 * params[0]; + lold = -2.0; + l = -1.0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = langaufun(&x,params); + + if (l < lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-1); + + maxx = x; + + fy = l/2; + + + // Search for right x location of fy + + p = maxx + params[0]; + step = params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-2); + + fxr = x; + + + // Search for left x location of fy + + p = maxx - 0.5 * params[0]; + step = -params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-3); + + + fxl = x; + + FWHM = fxr - fxl; + return (0); +} + +int N = 512; +// string filename = "~flionett/work/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/output/Run_Bias_Scan-B6-A-299-8431_Tuple.root"; +// string filename = "~flionett/work/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/output/Plots/AnalysisOutput_A6_300_1_8431.root"; +// int channel = -1; + +void TbUTLandauFits(string filename, int channel); + +int main(int argc, char *argv[]) +{ + getLHCbStyle(); + PersonalStyle(); + + if ((argc ==2) && (string(argv[1]) == "--info")) + { + cout << "**************************************************" << endl; + cout << "Put comments here." << endl; + cout << "**************************************************" << endl; + return 0; + } + else if (argc < 3) + { + cout << "**************************************************" << endl; + cout << "Error! Arguments missing..." << endl; + cout << "Please use the following format:" << endl; + cout << "./TbUTLandauFits [1] [2]" << endl; + cout << "with:" << endl; + cout << "[1] = filename (full path);" << endl; + cout << "[2] = channel." << endl; + cout << "Type ./TbUTLandauFits --info for more information." << endl; + cout << "**************************************************" << endl; + return 0; + } + else if (argc == 3) + { + cout << "**************************************************" << endl; + cout << "Filename: " << argv[1] << endl; + cout << "Channel: " << argv[2] << endl; + cout << "**************************************************" << endl; + TbUTLandauFits(string(argv[1]),atoi(argv[2])); + return 0; + } + else + { + cout << "**************************************************" << endl; + cout << "Error! Too many arguments provided..." << endl; + cout << "**************************************************" << endl; + return 0; + } +} + +// The variable is the filename with the complete path. +void TbUTLandauFits(string filename, int channel) { + TString inFilename = filename.substr(0,filename.find(".root")); // No file extension. + cout << inFilename << endl; + TString outFilename = filename.substr(0,filename.find(".root"))+"_LandauFits"; // No file extension. + cout << outFilename << endl; + + TString pathToInput = filename.substr(0,filename.find_last_of("/")); + cout << pathToInput << endl; + TString pathToOutput = filename.substr(0,filename.find_last_of("/"))+"/LandauFits"; + cout << pathToOutput << endl; + TString pathToFigures = filename.substr(0,filename.find_last_of("/"))+"/LandauFits/Figures"; + + // Create output folders if they do not exist. + TString pathToOutputMake = "mkdir -p "+pathToOutput; + cout << "Create folder " << pathToOutput << " where output will be saved" << endl; + system(string(pathToOutputMake).c_str()); + TString path_to_make_figures = "mkdir -p "+pathToFigures; + system(string(path_to_make_figures).c_str()); + cout << "Setting figures to: " << path_to_make_figures << endl; + + TFile *fIn = new TFile(inFilename+".root","READ"); + TFile *fOut = new TFile(outFilename+".root","RECREATE"); + + TH1F *hS = 0; + TH1F *hS1ch[N]; + + if (channel == -1) { + TList *list = new TList; + for (int i=0; iGet("hlandau_"+TString(ch)); + InitHist(hS1ch[i],Form("Signal distribution on channel %d",i),"ADC counts",""); + list->Add(hS1ch[i]); + } + hS=(TH1F *)hS1ch[0]->Clone("hlandau_0"); + InitHist(hS,"Signal distribution","ADC counts",""); + hS->Reset(); + hS->Merge(list); + /* + TTree *t = (TTree *)fIn->Get("TbUT/Clusters"); + vector clustersCharge(10); + + t->SetBranchAddress("clustersCharge",&(&clustersCharge[10])); + if (fIn->GetListOfKeys()->Contains("TbUT/Clusters/clustersCharge")) { + + + hS = (TH1F *)fIn->Get("TbUT/Clusters/clustersCharge"); + } + else { + cout << "ERROR! No input object found." << endl; + return ; + } + */ + } + else { + // Convert the integer to string. + stringstream s; + s << channel; + string ch(s.str()); + + hS = (TH1F *)fIn->Get("hlandau_"+TString(ch)); + InitHist(hS,"Signal distribution","ADC counts",""); + if (hS->GetEntries() == 0) { + cout << "ERROR! The selected histogram is empty." << endl; + return; + } + } + + // Fitting S histo + printf("Fitting...\n"); + + // Setting fit range and start values + Double_t fr[2]; + Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4]; + fr[0]=0.4*hS->GetMean(); + fr[1]=3.0*hS->GetMean(); + + pllo[0]=0.; pllo[1]=100.0; pllo[2]=1.; pllo[3]=0.01; + plhi[0]=50.; plhi[1]=350.0; plhi[2]=100000000.; plhi[3]=50.; + sv[0]=30.; sv[1]=250.; sv[2]=500.; sv[3]=3.; + + Double_t chisqr; + Int_t ndf; + TF1 *fitsnr = langaufit(hS,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf); + + Double_t SPeak, SFWHM; + langaupro(fp,SPeak,SFWHM); + + printf("Fitting done\nPlotting results...\n"); + + // Global style settings + gStyle->SetOptStat(1111); + gStyle->SetOptFit(111); + gStyle->SetLabelSize(0.03,"x"); + gStyle->SetLabelSize(0.03,"y"); + + // hS->GetXaxis()->SetRange(0,70); + // hS->Draw(); + // fitsnr->Draw("lsame"); + + TCanvas *cS = new TCanvas("cS","",800,600); + cS->cd(); + fitsnr->SetLineColor(kRed); + DrawHistFunc(cS,hS,fitsnr,pathToFigures); + + fOut->cd(); + hS->Write(); + fIn->Close(); + fOut->Close(); + + return; +} diff --git a/TbUT/options/TbUTPedestal.py b/TbUT/options/TbUTPedestal.py index b15557c..fcce402 100644 --- a/TbUT/options/TbUTPedestal.py +++ b/TbUT/options/TbUTPedestal.py @@ -2,7 +2,7 @@ sys.path.append( '../Tb/TbUT/options/python' ) from TbUTPedestalRunner import TbUTPedestalRunner app=TbUTPedestalRunner() -app.inputData="/afs/cern.ch/work/i/ibezshyi/public/TESTBEAM/testsoft/KeplerDev_v3r0/Tb/Kepler/eos/lhcb/testbeam/ut/OfficialData/July2015/BoardA6/RawData/Pedestal-B6-A-209-0.dat" +app.inputData="/afs/cern.ch/work/f/flionett/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/Tb/Kepler/eos/lhcb/testbeam/ut/OfficialData/July2015/BoardA6/RawData/Pedestal-B6-A-209-0.dat" app.isAType=True app.eventMax=20000 app.pedestalOutputData ='$KEPLERROOT/../../../output/Pedestal-BoardA6.dat' diff --git a/TbUT/options/TbUTRun.py b/TbUT/options/TbUTRun.py index 208a374..d5a84ae 100644 --- a/TbUT/options/TbUTRun.py +++ b/TbUT/options/TbUTRun.py @@ -2,7 +2,7 @@ sys.path.append( '../Tb/TbUT/options/python' ) from TbUTClusterizator import TbUTClusterizator app=TbUTClusterizator() -app.inputData="/afs/cern.ch/work/i/ibezshyi/public/TESTBEAM/testsoft/KeplerDev_v3r0/Tb/Kepler/eos/lhcb/testbeam/ut/OfficialData/July2015/BoardA6/RawData/Run_Bias_Scan-B6-A-299-8431.dat" +app.inputData="/afs/cern.ch/work/f/flionett/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/Tb/Kepler/eos/lhcb/testbeam/ut/OfficialData/July2015/BoardA6/RawData/Run_Bias_Scan-B6-A-299-8431.dat" app.isAType=True app.sensorType = 'PType' app.eventMax = 100000 diff --git a/TbUT/scripts/ClusterWithTrackAna.C b/TbUT/scripts/ClusterWithTrackAna.C index b190d65..8638dc6 100644 --- a/TbUT/scripts/ClusterWithTrackAna.C +++ b/TbUT/scripts/ClusterWithTrackAna.C @@ -8,6 +8,8 @@ #include #include +#include "../../Tools/Landau.C" + void addGraphics(TH1 *h, int iCol = 1, TString XTitle="", TString YTitle="") { h->SetXTitle(XTitle); @@ -240,6 +242,12 @@ + + + + + + TH1F* hcAll = new TH1F("hcAll","Cluster charge",100,0.0,1000.0); TH1F* hcTrk = new TH1F("hcTrk","Cluster charge",100,0.0,1000.0); TH1F* hcTrkCorr = new TH1F("hcTrkCorr","Cluster charge",100,0.0,1000.0); @@ -781,7 +789,42 @@ hcTrk->Draw("same"); hcTrk1->Draw("same"); hcTrk2->Draw("same"); - + + + // Fitting S histo + printf("Fitting...\n"); + + // Setting fit range and start values + Double_t fr[2]; + Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4]; + fr[0]=0.4*hcTrk->GetMean(); + fr[1]=3.0*hcTrk->GetMean(); + + pllo[0]=0.; pllo[1]=100.0; pllo[2]=1.; pllo[3]=0.01; + plhi[0]=50.; plhi[1]=350.0; plhi[2]=100000000.; plhi[3]=50.; + sv[0]=30.; sv[1]=250.; sv[2]=500.; sv[3]=3.; + + Double_t chisqr; + Int_t ndf; + TF1 *fitsnr = langaufit(hcTrk,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf); + + Double_t SPeak, SFWHM; + langaupro(fp,SPeak,SFWHM); + + printf("Fitting done\nPlotting results...\n"); + + // Global style settings + gStyle->SetOptStat(1111); + gStyle->SetOptFit(111); + gStyle->SetLabelSize(0.03,"x"); + gStyle->SetLabelSize(0.03,"y"); + + fitsnr->SetLineColor(kRed); + fitsnr->Draw(); + + + + TLegend* legend2 = new TLegend(0.15,0.70,0.94,0.89); legend2->SetFillStyle(0); legend2->SetBorderSize(0); diff --git a/Tools/AddNewBranch.py b/Tools/AddNewBranch.py new file mode 100644 index 0000000..968c162 --- /dev/null +++ b/Tools/AddNewBranch.py @@ -0,0 +1,78 @@ +""" +Author: Federica Lionett +Date: March 2nd, 2016 + +Description: + Add a branch called branchname and filled with branchvalue to a tree identified by (filename,treename) and save the output in a new file. + +How to run it: + python AddNewBranch.py --filename [filename] --treename [treename] --branchname [branchname] --branchvalue [branchvalue] +""" + +import ROOT +import argparse +from array import array + +def AddNewBranch(filename,treename,branchname,branchvalue) : + print 'Adding a new branch...' + print 'Input file: ', filename + print 'Branch name: ', branchname + print 'Branch value: ', branchvalue + + fIn = ROOT.TFile(filename,'READ') + + # Check that the tree exists. + if not fIn.GetListOfKeys().Contains(treename) and not fIn.GetListOfKeys().Contains(treename.rpartition('/')[2]) : + print 'ERROR! Object not found in the file.' + return + + tIn = fIn.Get(treename) + + fOut = ROOT.TFile(filename.replace('.root','')+'_With'+branchname+'Label'+'.root','RECREATE') + + fOut.cd() + + tOut = tIn.CloneTree() + + new = array('i',[0]) + # new = int(branchvalue) + new[0] = int(branchvalue) + print new + bNew = tOut.Branch(branchname,new,branchname+'/I') + + for iEvent in range(tOut.GetEntries()) : + tOut.GetEntry(iEvent) + bNew.Fill() + + # Write output. + fOut.Write() + + # Close files. + fIn.Close(); + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Add a branch to a tree.') + + parser.add_argument('--filename',required=True,help='filename') + parser.add_argument('--treename',required=True,help='treename') + parser.add_argument('--branchname',required=True,help='branchname') + parser.add_argument('--branchvalue',required=True,help='branchvalue') + + args = parser.parse_args() + + # Parameters and configuration. + filename = args.filename + treename = args.treename + branchname = args.branchname + branchvalue = args.branchvalue + + AddNewBranch(filename,treename,branchname,branchvalue) + diff --git a/Tools/CreateReducedTree.py b/Tools/CreateReducedTree.py new file mode 100644 index 0000000..0357a45 --- /dev/null +++ b/Tools/CreateReducedTree.py @@ -0,0 +1,72 @@ +""" +Author: Federica Lionett +Date: February 5th, 2016 + +Description: + Create a reduced tree containing only a subset of the branches of the original tree. + +How to run it: + python CreateReducedTree.py --filename [filename] --treename [treename] --branches [list of (branches,type) pairs] +""" + +import ROOT +import argparse + +def CreateReducedTree(filename,treename,branches) : + fIn = ROOT.TFile(filename,'READ') + + # Check that the tree exists. + if not fIn.GetListOfKeys().Contains(treename) : + print 'ERROR! Object not found.' + return + + tIn = fIn.Get(treename) + + # Check that the branches exist. + for branch in branches : + if not tIn.GetListOfBranches().Contains(branch) : + print 'ERROR! Branch not found.' + return + + print 'Creating a reduced tree...' + print 'File: ', filename + print 'Tree: ', treename + print 'Branches: ', branches + + tIn.SetBranchStatus('*',0) + for branch in branches : + tIn.SetBranchStatus(branch,1) + + fOut = ROOT.TFile(filename.replace('.root','_ReducedTree.root'),'RECREATE') + tOut = tIn.CloneTree() + + # Write output. + fOut.cd() + fOut.Write() + + # Close files. + fIn.Close(); + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Create a reduced tree.') + + parser.add_argument('--filename',required=True,help='filename') + parser.add_argument('--treename',required=True,help='treename') + parser.add_argument('--branches',required=True,nargs='+',help='branches') # Each branch should be separated from the others by a space. + + args = parser.parse_args() + + # Parameters and configuration. + filename = args.filename + treename = args.treename + branches = args.branches + + CreateReducedTree(filename,treename,branches) diff --git a/Tools/FindAndReplaceString.C b/Tools/FindAndReplaceString.C new file mode 100644 index 0000000..ca1a8fd --- /dev/null +++ b/Tools/FindAndReplaceString.C @@ -0,0 +1,36 @@ +//************************************************ +// Author: Federica Lionetto +// Created on: 28/09/2015 +//************************************************ + +/* +Find and replace a string. +*/ + +// Header guard. +#ifndef __FINDANDREPLACESTRING_C_INCLUDED__ +#define __FINDANDREPLACESTRING_C_INCLUDED__ + +#include "Lib.C" + +//************************************************ +// +// Declarations. +// + +void ReplaceStringInPlace(std::string& subject, const std::string& search, const std::string& replace); +//************************************************ +// +// Definitions. +// + +void ReplaceStringInPlace(std::string& subject, const std::string& search, const std::string& replace) { + size_t pos = 0; + while((pos = subject.find(search, pos)) != std::string::npos) { + subject.replace(pos, search.length(), replace); + pos += replace.length(); + } + return; +} + +#endif diff --git a/Tools/ImportModules.py b/Tools/ImportModules.py new file mode 100644 index 0000000..a6b7d13 --- /dev/null +++ b/Tools/ImportModules.py @@ -0,0 +1,105 @@ +""" +Author: Federica Lionetto +Date: November 6th, 2015 + +Description: +Import modules used for training. + +How to run it: +Include it in other python scripts. +""" + +import ROOT as ROOT + +import os +import time + +import sys +sys.path.insert(0,'../CreateNtuples') +from CreateNtuples_ToolsFile import * + +import subprocess + +import yaml +import argparse + +import collections + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from pylab import * + +import numpy as np +import root_numpy as rnp +import pandas as pd +import root_pandas as rpd + +import pandas.core.common as com +from pandas.core.index import Index + +from rep.utils import train_test_split +from rep.utils import calc_feature_correlation_matrix +from sklearn.metrics import roc_auc_score +from sklearn.metrics import roc_curve +from sklearn.metrics import auc +from sklearn.metrics import zero_one_loss + +# from rep.data.storage import LabeledDataStorage +# from rep.report import ClassificationReport + +from sklearn.tree import DecisionTreeClassifier + +from sklearn.grid_search import GridSearchCV +from sklearn.cross_validation import train_test_split +from sklearn.metrics import classification_report + +from rep.estimators import TMVAClassifier +from rep.estimators import SklearnClassifier +from rep.estimators.xgboost import XGBoostClassifier +from rep.estimators.pybrain import PyBrainClassifier +from rep.estimators.neurolab import NeurolabClassifier +# from rep.estimators.theanets import TheanetsClassifier + +# Classifiers for uniformity +from hep_ml.uboost import uBoostClassifier +from hep_ml.gradientboosting import UGradientBoostingClassifier, KnnAdaLossFunction, KnnFlatnessLossFunction + +from sklearn.ensemble import GradientBoostingClassifier +from sklearn.ensemble import ExtraTreesClassifier + +from hep_ml.commonutils import train_test_split +from hep_ml import uboost, gradientboosting as ugb, losses +from hep_ml.metrics import KnnBasedSDE, KnnBasedCvM + +from rep.report.metrics import RocAuc, LogLoss, OptimalMetric, ams + +from rep.metaml.gridsearch import RandomParameterOptimizer, AnnealingParameterOptimizer, SubgridParameterOptimizer, RegressionParameterOptimizer +from rep.metaml import FoldingClassifier +from rep.metaml import FoldingScorer +from rep.metaml import GridOptimalSearchCV +from rep.metaml import ClassifiersFactory + +from copy import deepcopy + +# Import modules for recursive feature elimination. +from sklearn.svm import SVC +from sklearn.cross_validation import StratifiedKFold +from sklearn.cross_validation import LeaveOneOut +from sklearn.feature_selection import RFECV +from sklearn.feature_selection import RFE +from sklearn.datasets import make_classification + +from sklearn.gaussian_process import GaussianProcess + +from sklearn.decomposition import PCA + +import itertools +from itertools import combinations +from itertools import chain + +import fnmatch + +from array import array + +import cPickle as pickle diff --git a/Tools/Landau.C b/Tools/Landau.C new file mode 100644 index 0000000..0500e36 --- /dev/null +++ b/Tools/Landau.C @@ -0,0 +1,234 @@ +//************************************************ +// Author: Federica Lionetto +// Created on: 07/05/2016 +//************************************************ + +// Header guard. +#ifndef __LANDAU_C_INCLUDED__ +#define __LANDAU_C_INCLUDED__ + +//----------------------------------------------------------------------- +// +// Convoluted Landau and Gaussian Fitting Function +// (using ROOT's Landau and Gauss functions) +// +// Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at) +// Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and +// Markus Friedl (Markus.Friedl@cern.ch) +// +// to execute this example, do: +// root > .x langaus.C +// or +// root > .x langaus.C++ +// +//----------------------------------------------------------------------- + +Double_t langaufun(Double_t *x, Double_t *par) { + + //Fit parameters: + //par[0]=Width (scale) parameter of Landau density + //par[1]=Most Probable (MP, location) parameter of Landau density + //par[2]=Total area (integral -inf to inf, normalization constant) + //par[3]=Width (sigma) of convoluted Gaussian function + // + //In the Landau distribution (represented by the CERNLIB approximation), + //the maximum is located at x=-0.22278298 with the location parameter=0. + //This shift is corrected within this function, so that the actual + //maximum is identical to the MP parameter. + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // number of convolution steps + Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx; + Double_t mpc; + Double_t fland; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp-xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for(i=1.0; i<=np/2; i++) { + xx = xlow + (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + + xx = xupp - (i-.5) * step; + fland = TMath::Landau(xx,mpc,par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0],xx,par[3]); + } + + return (par[2] * step * sum * invsq2pi / par[3]); +} + + + +TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) +{ + // Once again, here are the Landau * Gaussian parameters: + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + // + // Variables for langaufit call: + // his histogram to fit + // fitrange[2] lo and hi boundaries of fit range + // startvalues[4] reasonable start values for the fit + // parlimitslo[4] lower parameter limits + // parlimitshi[4] upper parameter limits + // fitparams[4] returns the final fit parameters + // fiterrors[4] returns the final fit errors + // ChiSqr returns the chi square + // NDF returns ndf + + Int_t i; + Char_t FunName[100]; + + sprintf(FunName,"Fitfcn_%s",his->GetName()); + + TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); + if (ffitold) delete ffitold; + + TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); + ffit->SetParameters(startvalues); + ffit->SetParNames("Width","MP","Area","GSigma"); + + for (i=0; i<4; i++) { + ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); + } + + his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot + + ffit->GetParameters(fitparams); // obtain fit parameters + for (i=0; i<4; i++) { + fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors + } + ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2 + NDF[0] = ffit->GetNDF(); // obtain ndf + + return (ffit); // return fit function + +} + + + +Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) { + + // Seaches for the location (x value) at the maximum of the + // Landau-Gaussian convolute and its full width at half-maximum. + // + // The search is probably not very efficient, but it's a first try. + + Double_t p,x,fy,fxr,fxl; + Double_t step; + Double_t l,lold; + Int_t i = 0; + Int_t MAXCALLS = 10000; + + + // Search for maximum + + p = params[1] - 0.1 * params[0]; + step = 0.05 * params[0]; + lold = -2.0; + l = -1.0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = langaufun(&x,params); + + if (l < lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-1); + + maxx = x; + + fy = l/2; + + + // Search for right x location of fy + + p = maxx + params[0]; + step = params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-2); + + fxr = x; + + + // Search for left x location of fy + + p = maxx - 0.5 * params[0]; + step = -params[0]; + lold = -2.0; + l = -1e300; + i = 0; + + while ( (l != lold) && (i < MAXCALLS) ) { + i++; + + lold = l; + x = p + step; + l = TMath::Abs(langaufun(&x,params) - fy); + + if (l > lold) + step = -step/10; + + p += step; + } + + if (i == MAXCALLS) + return (-3); + + + fxl = x; + + FWHM = fxr - fxl; + return (0); +} + +#endif diff --git a/Tools/Lib.C b/Tools/Lib.C new file mode 100644 index 0000000..6e03bc6 --- /dev/null +++ b/Tools/Lib.C @@ -0,0 +1,80 @@ +//************************************************ +// Author: Federica Lionetto +// Created on: 06/27/2014 +//************************************************ + +/* +List of useful libraries. +*/ + +// Header guard. +#ifndef __LIB_C_INCLUDED__ +#define __LIB_C_INCLUDED__ + +#include "TAxis.h" +#include "TBenchmark.h" +#include "TCanvas.h" +#include "TChain.h" +#include "TFile.h" +#include "TF1.h" +#include "TGaxis.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TLine.h" +#include "TLorentzVector.h" +#include "TMath.h" +#include "TMultiGraph.h" +#include "TNtuple.h" +#include "TPad.h" +#include "TPaveStats.h" +#include "TPaveText.h" +#include "TProfile.h" +#include "TRandom.h" +#include "TRandom3.h" +#include "TROOT.h" +#include "TString.h" +#include "TStyle.h" +#include "TSystem.h" +#include "TTree.h" +#include "TVirtualPad.h" + +#include +#include + +// #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// #include "RooAbsPdf.h" +#include "RooAddPdf.h" +#include "RooCBShape.h" +#include "RooConstVar.h" +#include "RooCurve.h" +// #include "RooDataSet.h" +#include "RooExponential.h" +#include "RooExtendPdf.h" +// #include "RooFitResult.h" +#include "RooGaussian.h" +#include "RooHist.h" +#include "RooPlot.h" +#include "RooRealVar.h" +#include "RooStats/SPlot.h" + +using namespace std; +using namespace RooFit; +using namespace RooStats; + +#endif diff --git a/Tools/Merge3Trees.py b/Tools/Merge3Trees.py new file mode 100644 index 0000000..d7fca25 --- /dev/null +++ b/Tools/Merge3Trees.py @@ -0,0 +1,94 @@ +""" +Author: Federica Lionett +Date: March 2nd, 2016 + +Description: + Merge three trees identified by (filename1,treename1), (filename2,treename2), and (filename3,treename3) and save the output in (filename,treename). + +How to run it: + python MergeTrees.py --filename1 [filename1] --treename1 [treename1] --filename2 [filename2] --treename2 [treename2] --filename3 [filename3] --treename3 [treename3] --filename [filename] --treename [treename] +""" + +import ROOT +import argparse + +def Merge3Trees(filename1,treename1,filename2,treename2,filename3,treename3,filename,treename) : + fIn1 = ROOT.TFile(filename1,'READ') + fIn2 = ROOT.TFile(filename2,'READ') + fIn3 = ROOT.TFile(filename3,'READ') + + # Check that the tree exists. + if not fIn1.GetListOfKeys().Contains(treename1) : + print 'ERROR! Object not found in the first file.' + return + if not fIn2.GetListOfKeys().Contains(treename2) : + print 'ERROR! Object not found in the second file.' + return + if not fIn3.GetListOfKeys().Contains(treename3) : + print 'ERROR! Object not found in the third file.' + return + + tIn1 = fIn1.Get(treename1) + tIn2 = fIn2.Get(treename2) + tIn3 = fIn3.Get(treename3) + tList = ROOT.TList() + tList.Add(tIn1) + tList.Add(tIn2) + tList.Add(tIn3) + + print 'Merging the three trees...' + print 'First file: ', filename1 + print 'First tree: ', treename1 + print 'Second file: ', filename2 + print 'Second tree: ', treename2 + print 'Third file: ', filename3 + print 'Third tree: ', treename3 + + fOut = ROOT.TFile(filename,'RECREATE') + + fOut.cd() + + tOut = ROOT.TTree.MergeTrees(tList) + tOut.SetName(treename) + + # Write output. + fOut.Write() + + # Close files. + fIn1.Close(); + fIn2.Close(); + fIn3.Close(); + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Merge three trees.') + + parser.add_argument('--filename1',required=True,help='filename1') + parser.add_argument('--treename1',required=True,help='treename1') + parser.add_argument('--filename2',required=True,help='filename2') + parser.add_argument('--treename2',required=True,help='treename2') + parser.add_argument('--filename3',required=True,help='filename3') + parser.add_argument('--treename3',required=True,help='treename3') + parser.add_argument('--filename',required=True,help='filename') + parser.add_argument('--treename',required=True,help='treename') + + args = parser.parse_args() + + # Parameters and configuration. + filename1 = args.filename1 + treename1 = args.treename1 + filename2 = args.filename2 + treename2 = args.treename2 + filename3 = args.filename3 + treename3 = args.treename3 + filename = args.filename + treename = args.treename + + MergeTrees(filename1,treename1,filename2,treename2,filename3,treename3,filename,treename) diff --git a/Tools/Merge4Trees.py b/Tools/Merge4Trees.py new file mode 100644 index 0000000..8368aee --- /dev/null +++ b/Tools/Merge4Trees.py @@ -0,0 +1,107 @@ +""" +Author: Federica Lionett +Date: March 2nd, 2016 + +Description: + Merge four trees identified by (filename1,treename1), (filename2,treename2), (filename3,treename3), and (filename4,treename4) and save the output in (filename,treename). + +How to run it: + python MergeTrees.py --filename1 [filename1] --treename1 [treename1] --filename2 [filename2] --treename2 [treename2] --filename3 [filename3] --treename3 [treename3] --filename4 [filename4] --treename4 [treename4] --filename [filename] --treename [treename] +""" + +import ROOT +import argparse + +def Merge4Trees(filename1,treename1,filename2,treename2,filename3,treename3,filename4,treename4,filename,treename) : + fIn1 = ROOT.TFile(filename1,'READ') + fIn2 = ROOT.TFile(filename2,'READ') + fIn3 = ROOT.TFile(filename3,'READ') + fIn4 = ROOT.TFile(filename4,'READ') + + # Check that the tree exists. + if not fIn1.GetListOfKeys().Contains(treename1) : + print 'ERROR! Object not found in the first file.' + return + if not fIn2.GetListOfKeys().Contains(treename2) : + print 'ERROR! Object not found in the second file.' + return + if not fIn3.GetListOfKeys().Contains(treename3) : + print 'ERROR! Object not found in the third file.' + return + if not fIn4.GetListOfKeys().Contains(treename4) : + print 'ERROR! Object not found in the fourth file.' + return + + tIn1 = fIn1.Get(treename1) + tIn2 = fIn2.Get(treename2) + tIn3 = fIn3.Get(treename3) + tIn4 = fIn4.Get(treename4) + tList = ROOT.TList() + tList.Add(tIn1) + tList.Add(tIn2) + tList.Add(tIn3) + tList.Add(tIn4) + + print 'Merging the four trees...' + print 'First file: ', filename1 + print 'First tree: ', treename1 + print 'Second file: ', filename2 + print 'Second tree: ', treename2 + print 'Third file: ', filename3 + print 'Third tree: ', treename3 + print 'Fourth file: ', filename4 + print 'Fourth tree: ', treename4 + + fOut = ROOT.TFile(filename,'RECREATE') + + fOut.cd() + + tOut = ROOT.TTree.MergeTrees(tList) + tOut.SetName(treename) + + # Write output. + fOut.Write() + + # Close files. + fIn1.Close(); + fIn2.Close(); + fIn3.Close(); + fIn4.Close(); + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Merge four trees.') + + parser.add_argument('--filename1',required=True,help='filename1') + parser.add_argument('--treename1',required=True,help='treename1') + parser.add_argument('--filename2',required=True,help='filename2') + parser.add_argument('--treename2',required=True,help='treename2') + parser.add_argument('--filename3',required=True,help='filename3') + parser.add_argument('--treename3',required=True,help='treename3') + parser.add_argument('--filename4',required=True,help='filename4') + parser.add_argument('--treename4',required=True,help='treename4') + parser.add_argument('--filename',required=True,help='filename') + parser.add_argument('--treename',required=True,help='treename') + + args = parser.parse_args() + + # Parameters and configuration. + filename1 = args.filename1 + treename1 = args.treename1 + filename2 = args.filename2 + treename2 = args.treename2 + filename3 = args.filename3 + treename3 = args.treename3 + filename4 = args.filename4 + treename4 = args.treename4 + filename = args.filename + treename = args.treename + + MergeTrees(filename1,treename1,filename2,treename2,filename3,treename3,filename4,treename4,filename,treename) diff --git a/Tools/MergeBranchesWithTree.py b/Tools/MergeBranchesWithTree.py new file mode 100644 index 0000000..23176e8 --- /dev/null +++ b/Tools/MergeBranchesWithTree.py @@ -0,0 +1,114 @@ +""" +Author: Federica Lionetto +Date: April 19th, 2016 + +Description: + Create a new tree, consisting in the original tree, plus a branch from another tree. + +How to run it: + python MergeBranchesWithTree.py --filename1 [filename1] --treename1 [treename1] --filename2 [filename2] --treename2 [treename2] --branches [list of (branches,type) pairs] +""" + +import ROOT +import argparse +from array import array +import numpy as np + +def MergeBranchesWithTree(filename1,treename1,filename2,treename2,branches,branchTypes) : + fIn1 = ROOT.TFile(filename1,'READ') + fIn2 = ROOT.TFile(filename2,'READ') + + # Check that the tree exists. + if not fIn1.GetListOfKeys().Contains(treename1) : + print 'ERROR! Object not found in the first file.' + return + if not fIn2.GetListOfKeys().Contains(treename2) : + print 'ERROR! Object not found in the second file.' + return + + tIn1 = fIn1.Get(treename1) + tIn2 = fIn2.Get(treename2) + + # Check that the branches exist. + for branch in branches : + if not tIn2.GetListOfBranches().Contains(branch) : + print 'ERROR! Branch not found.' + return + + filename = filename1.replace('.root','_Extended.root') + treename = treename1 + + print 'Creating a new tree...' + print 'File: ', filename + print 'Tree: ', treename + print 'Branches: ', branches + + + fOut = ROOT.TFile(filename,'RECREATE') + tOut = tIn1.CloneTree() + + nEntries1 = tIn1.GetEntries() + nEntries2 = tIn2.GetEntries() + if not nEntries1 == nEntries2 : + print 'ERROR! Different number of entries in the two trees.' + return + + nEntries = nEntries2 + print 'Number of entries: ', nEntries + + branchInList = [] + branchOutList = [] + + for iBranch, branch in enumerate(branches) : + if not tIn1.GetListOfBranches().Contains(branch) : + print 'Branch', branch + branchInList.append(array(branchTypes[iBranch],[0])) + tIn2.SetBranchAddress(branch,branchInList[-1]) + + branchOutList.append(tOut.Branch(branch,branchInList[-1],branch+'/'+branchTypes[iBranch].upper())) + + for iEvent in range(nEntries) : + tIn2.GetEntry(iEvent) + for iBranch in range(len(branchInList)) : + if iEvent < 5 : + print branchInList[iBranch][0] + print branchInList[iBranch] + branchOutList[iBranch].Fill() + + # Write output. + fOut.cd() + fOut.Write() + + # Close files. + fIn1.Close() + fIn2.Close() + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Create a new tree, consisting in the original tree, plus some branches from another tree.') + + parser.add_argument('--filename1',required=True,help='filename1') + parser.add_argument('--treename1',required=True,help='treename1') + parser.add_argument('--filename2',required=True,help='filename2') + parser.add_argument('--treename2',required=True,help='treename2') + parser.add_argument('--branches',required=True,nargs='+',help='branches') # Each branch should be separated from the others by a space. + parser.add_argument('--branchTypes',required=True,nargs='+',help='branch types') # Each branch type should be separated from the others by a space. + + args = parser.parse_args() + + # Parameters and configuration. + filename1 = args.filename1 + treename1 = args.treename1 + filename2 = args.filename2 + treename2 = args.treename2 + branches = args.branches + branchTypes = args.branchTypes + + MergeBranchesWithTree(filename1,treename1,filename2,treename2,branches,branchTypes) diff --git a/Tools/MergeTrees.py b/Tools/MergeTrees.py new file mode 100644 index 0000000..7774242 --- /dev/null +++ b/Tools/MergeTrees.py @@ -0,0 +1,81 @@ +""" +Author: Federica Lionett +Date: March 2nd, 2016 + +Description: + Merge two trees identified by (filename1,treename1) and (filename2,treename2) and save the output in (filename,treename). + +How to run it: + python MergeTrees.py --filename1 [filename1] --treename1 [treename1] --filename2 [filename2] --treename2 [treename2] --filename [filename] --treename [treename] +""" + +import ROOT +import argparse + +def MergeTrees(filename1,treename1,filename2,treename2,filename,treename) : + fIn1 = ROOT.TFile(filename1,'READ') + fIn2 = ROOT.TFile(filename2,'READ') + + # Check that the tree exists. + if not fIn1.GetListOfKeys().Contains(treename1) : + print 'ERROR! Object not found in the first file.' + return + if not fIn2.GetListOfKeys().Contains(treename2) : + print 'ERROR! Object not found in the second file.' + return + + tIn1 = fIn1.Get(treename1) + tIn2 = fIn2.Get(treename2) + tList = ROOT.TList() + tList.Add(tIn1) + tList.Add(tIn2) + + print 'Merging the two trees...' + print 'First file: ', filename1 + print 'First tree: ', treename1 + print 'Second file: ', filename2 + print 'Second tree: ', treename2 + + fOut = ROOT.TFile(filename,'RECREATE') + + fOut.cd() + + tOut = ROOT.TTree.MergeTrees(tList) + tOut.SetName(treename) + + # Write output. + fOut.Write() + + # Close files. + fIn1.Close(); + fIn2.Close(); + fOut.Close() + + return + +############### +# +# Main function +# + +if __name__ == "__main__" : + parser = argparse.ArgumentParser(description='Merge two trees.') + + parser.add_argument('--filename1',required=True,help='filename1') + parser.add_argument('--treename1',required=True,help='treename1') + parser.add_argument('--filename2',required=True,help='filename2') + parser.add_argument('--treename2',required=True,help='treename2') + parser.add_argument('--filename',required=True,help='filename') + parser.add_argument('--treename',required=True,help='treename') + + args = parser.parse_args() + + # Parameters and configuration. + filename1 = args.filename1 + treename1 = args.treename1 + filename2 = args.filename2 + treename2 = args.treename2 + filename = args.filename + treename = args.treename + + MergeTrees(filename1,treename1,filename2,treename2,filename,treename) diff --git a/Tools/Style.C b/Tools/Style.C new file mode 100644 index 0000000..3acae52 --- /dev/null +++ b/Tools/Style.C @@ -0,0 +1,947 @@ +//************************************************ +// Author: Federica Lionetto +// Created on: 06/27/2014 +//************************************************ + +/* + List of useful functions. + */ + +// Header guard. +#ifndef __STYLE_C_INCLUDED__ +#define __STYLE_C_INCLUDED__ + +#include "Lib.C" +#include "lhcbStyle.C" + +//************************************************ +// +// Declarations. +// + +void dca(); + +void PersonalStyle(); + +void InitHist(TH1 *hist, TString title = "", TString x = "", TString y = ""); + +void InitHist2(TH2 *hist, TString title = "", TString x = "", TString y = "", TString z = ""); + +void DrawHist(TCanvas *canvas, TH1 *hist, TString option = "", TString folder = ""); + +void DrawHistCompare(TCanvas *canvas, TH1 *hist1, TString option1, TH1 *hist2, TString option2, TLegend *leg, TString folder = ""); + +void DrawHistCompare3(TCanvas *canvas, TH1 *hist1, TString option1, TH1 *hist2, TString option2, TH1 *hist3, TString option3, TLegend *leg, TString folder = ""); + +void DrawHistFunc(TCanvas *canvas, TH1 *hist, TF1 *func, TString folder = ""); + +void DrawHist2(TCanvas *canvas, TH2 *hist, TString option = "", TString folder = ""); + +void DrawHistStats(TCanvas *canvas, TH1 *hist, TPaveStats *stats, TString option = "", TString folder = ""); + +void DrawHistStats2(TCanvas *canvas, TH2 *hist, TPaveStats *stats, TString option = "", TString folder = ""); + +void InitGraph(TGraph *graph, TString title = "", TString x = "", TString y = ""); + +void InitGraphErrors(TGraphErrors *graph, TString title = "", TString x = "", TString y = ""); + +void DrawGraph(TCanvas *canvas, TGraph *graph, TString option = "APC", TString folder = ""); + +void DrawGraphErrors(TCanvas *canvas, TGraphErrors *graph, TString option = "APC", TString folder = ""); + +void DrawGraphFunc(TCanvas *canvas, TGraph *graph, TF1 *func, TString option = "APC", TString folder = ""); + +void DrawGraphFunc2(TCanvas *canvas, TGraph *graph, TF1 *func1, TF1 *func2, TString option = "APC", TString folder = ""); + +void DrawGraphErrorsFunc(TCanvas *canvas, TGraphErrors *graph, TF1 *func, TString option = "APC", TString folder = ""); + +void DrawGraphErrorsFunc2(TCanvas *canvas, TGraphErrors *graph, TF1 *func1, TF1 *func2, TString option = "APC", TString folder = ""); + +void DrawGraphCompare(TCanvas *canvas, TMultiGraph *graph, TLegend *leg, TString title = "", TString x = "", TString y = "", TString option = "APE", TString folder = ""); + +void DrawFrame(TCanvas *canvas, RooPlot *frame, TLegend *leg, bool logy = false, bool res = false, TString folder = ""); + +TLegend *CreateLegend2(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TString option = "w", Double_t x1 = 0.64, Double_t y1 = 0.59, Double_t x2 = 0.94, Double_t y2 = 0.89); + +TLegend *CreateLegend3(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TString option = "w", Double_t x1 = 0.64, Double_t y1 = 0.59, Double_t x2 = 0.94, Double_t y2 = 0.89); + +TLegend *CreateLegend4(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TObject *obj4, TString lab4, TString option4, TString option = "w", Double_t x1 = 0.64, Double_t y1 = 0.59, Double_t x2 = 0.94, Double_t y2 = 0.89); + +TLegend *CreateLegend5(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TObject *obj4, TString lab4, TString option4, TObject *obj5, TString lab5, TString option5, TString option = "w", Double_t x1 = 0.64, Double_t y1 = 0.59, Double_t x2 = 0.94, Double_t y2 = 0.89); + +TPaveStats *CreateStats(TH1 *hist, TString option = ""); + +double Residual(Double_t datum, Double_t pdf); + +void AddLinesToResidualHist(const TH1F *hist); + +TH1F *CreateResidualHist(const RooHist *rhist, const RooCurve *curve); + +//************************************************ +// +// Definitions. +// + +// Thanks to Angelo Di Canto. +void dca() +{ + TStyle *dcastyle = new TStyle("dcastyle","The real perfect style"); + + // canvas + dcastyle->SetCanvasColor(0); + dcastyle->SetCanvasBorderSize(10); + dcastyle->SetCanvasBorderMode(0); + dcastyle->SetCanvasDefW(600); + dcastyle->SetCanvasDefH(600); + + // pads + dcastyle->SetPadColor(0); + dcastyle->SetPadBorderSize(10); + dcastyle->SetPadBorderMode(0); + dcastyle->SetPadLeftMargin(.18); + dcastyle->SetPadRightMargin(.05); + dcastyle->SetPadBottomMargin(.12); + dcastyle->SetPadTopMargin(.1); + //dcastyle->SetPadGridX(1); + //dcastyle->SetPadGridY(1); + dcastyle->SetPadTickX(0); + dcastyle->SetPadTickY(0); + + // frame + dcastyle->SetFrameBorderMode(0); + dcastyle->SetFrameBorderSize(10); + dcastyle->SetFrameFillStyle(0); + dcastyle->SetFrameFillColor(0); + dcastyle->SetFrameLineColor(1); + dcastyle->SetFrameLineStyle(0); + dcastyle->SetFrameLineWidth(1); + + // histogram + dcastyle->SetHistFillColor(0); + dcastyle->SetHistFillStyle(1001);// solid + dcastyle->SetHistLineColor(1); + dcastyle->SetHistLineStyle(0); + dcastyle->SetHistLineWidth(1.5); +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,16,0) + dcastyle->SetHistMinimumZero(); +#endif + dcastyle->SetOptStat(0); + dcastyle->SetOptFit(0); + dcastyle->SetStatColor(0); + dcastyle->SetStatBorderSize(1); + dcastyle->SetStatFontSize(.05); + + // graph + dcastyle->SetEndErrorSize(0); + dcastyle->SetErrorX(0.5); + + // marker + dcastyle->SetMarkerStyle(20); + dcastyle->SetMarkerColor(kBlack); + dcastyle->SetMarkerSize(0.6); + + // title + dcastyle->SetOptTitle(1); +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,16,0) + dcastyle->SetTitleAlign(33); + dcastyle->SetTitleX(.95); +#else + dcastyle->SetTitleX(.15); +#endif + dcastyle->SetTitleFillColor(0); + dcastyle->SetTitleBorderSize(0); + dcastyle->SetTitleStyle(0); + + // axes + dcastyle->SetNdivisions(505,"X"); + dcastyle->SetNdivisions(505,"Y"); + + dcastyle->SetTitleSize(.05,"X");//.055 + dcastyle->SetTitleOffset(1.,"X");//1.2,0.9 + dcastyle->SetLabelOffset(0.003,"X"); + dcastyle->SetLabelSize(.05,"X"); + dcastyle->SetLabelFont(42,"X"); + + dcastyle->SetTitleSize(.05,"Y");//.055 + dcastyle->SetTitleOffset(1.8,"Y"); + dcastyle->SetLabelOffset(0.008,"Y"); + dcastyle->SetLabelSize(.05,"Y"); + dcastyle->SetLabelFont(42,"Y"); + + dcastyle->SetStripDecimals(kFALSE); + + dcastyle->SetTitleSize(0.05,"Z"); + dcastyle->SetTitleOffset(1.800,"Z"); + dcastyle->SetLabelOffset(0.008,"Z"); + dcastyle->SetLabelSize(0.05,"Z"); + dcastyle->SetLabelFont(42,"Z"); + + // fonts + dcastyle->SetTextSize(.05);//.055 + dcastyle->SetTextFont(42); + dcastyle->SetStatFont(42); + dcastyle->SetTitleFont(42,""); + dcastyle->SetTitleFont(42,"Z"); + dcastyle->SetTitleFont(42,"X"); + dcastyle->SetTitleFont(42,"Y"); + + // function + dcastyle->SetFuncColor(kBlue); + dcastyle->SetFuncStyle(0); + dcastyle->SetFuncWidth(1); + + // legend +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,16,0) + dcastyle->SetLegendBorderSize(1); +#endif + + // palette + dcastyle->SetPalette(1); + dcastyle->SetNumberContours(20); + + // set dcastyle as current style + gROOT->SetStyle("dcastyle"); + + gSystem->ProcessEvents(); +} + +// PersonalStyle modifies the official LHCb style. +void PersonalStyle() { + TStyle *style = getLHCbStyle(); + // style->SetOptStat("emruo"); + style->SetTitleOffset(1.47,"Y"); + style->SetStatY(0.88); + // TGaxis::SetMaxDigits(3); + style->SetPadTopMargin(0.0735); + style->SetPadLeftMargin(0.2); + style->SetPadRightMargin(0.1); + + // Set the line style + style->SetLineStyleString(kDashed, "40 40"); + style->SetLineStyleString(kDotted, "4 12"); + style->SetLineStyleString(kDashDotted, "20 20 4 20"); + + return; +} + +// InitHist defines some attributes of the "hist" histogram: the title "title" and the x-axis and y-axis labels "x" and "y". +void InitHist(TH1 *hist, TString title, TString x, TString y) { + hist->SetTitle(title); + hist->GetXaxis()->SetTitle(x); + if (y == "") + hist->GetYaxis()->SetTitle(Form("Entries / %g",hist->GetBinWidth(1))); + else + hist->GetYaxis()->SetTitle(y); + + return; +} + +// InitHist2 defines some attributes of the "hist" histogram: the title "title" and the x-axis, y-axis, and z-axis labels "x", "y", and "z". +void InitHist2(TH2 *hist, TString title, TString x, TString y, TString z) { + hist->SetTitle(title); + hist->GetXaxis()->SetTitle(x); + hist->GetYaxis()->SetTitle(y); + if (z == "") + hist->GetZaxis()->SetTitle(Form("Entries / (%g*%g)",hist->GetXaxis()->GetBinWidth(1),hist->GetYaxis()->GetBinWidth(1))); + else + hist->GetZaxis()->SetTitle(z); + hist->SetLabelOffset(0.005,"Z"); + hist->SetTitleOffset(0.65,"Z"); + + return; +} + +// DrawHist draws the "hist" histogram in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHist(TCanvas *canvas, TH1 *hist, TString option, TString folder) { + TString path; + /* + if (folder == "") + path = "Figures/"; + else + path = folder + "/" + "Figures/"; + */ + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist->Draw(option); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHistCompare draws the "hist1" and "hist2" histograms in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHistCompare(TCanvas *canvas, TH1 *hist1, TString option1, TH1 *hist2, TString option2, TLegend *leg, TString folder) { + TString path; + /* + if (folder == "") + path = "Figures/"; + else + path = folder + "/" + "Figures/"; + */ + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist1->SetMarkerColor(kViolet); + hist2->SetMarkerColor(kTeal-5); + hist1->SetLineColor(kViolet); + hist2->SetLineColor(kTeal-5); + hist1->SetMarkerStyle(20); + hist2->SetMarkerStyle(21); + if (hist1->GetMaximum() > hist2->GetMaximum()) { + hist1->Draw(option1); + hist2->Draw(option2+"SAME"); + } + else { + hist2->Draw(option2); + hist1->Draw(option1+"SAME"); + } + leg->Draw(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHistCompare3 draws the "hist1", "hist2", and "hist3" histograms in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHistCompare3(TCanvas *canvas, TH1 *hist1, TString option1, TH1 *hist2, TString option2, TH1 *hist3, TString option3, TLegend *leg, TString folder) { + TString path; + /* + if (folder == "") + path = "Figures/"; + else + path = folder + "/" + "Figures/"; + */ + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist1->SetMarkerColor(kViolet); + hist2->SetMarkerColor(kTeal-5); + hist3->SetMarkerColor(kBlue); + hist1->SetLineColor(kViolet); + hist2->SetLineColor(kTeal-5); + hist3->SetLineColor(kBlue); + hist1->SetMarkerStyle(20); + hist2->SetMarkerStyle(21); + hist3->SetMarkerStyle(22); + if (hist1->GetMaximum() > hist2->GetMaximum()) { + if (hist1->GetMaximum() > hist3->GetMaximum()) + { + hist1->Draw(option1); + hist2->Draw(option2+"SAME"); + hist3->Draw(option3+"SAME"); + } + else + { + hist3->Draw(option3); + hist1->Draw(option1+"SAME"); + hist2->Draw(option2+"SAME"); + } + } + else + { + if (hist2->GetMaximum() > hist3->GetMaximum()) + { + hist2->Draw(option2); + hist1->Draw(option1+"SAME"); + hist3->Draw(option3+"SAME"); + } + else + { + hist3->Draw(option3); + hist1->Draw(option1+"SAME"); + hist2->Draw(option2+"SAME"); + } + } + leg->Draw(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHistFunc draws the "hist" histogram and the "func" function in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHistFunc(TCanvas *canvas, TH1 *hist, TF1 *func, TString folder){ + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist->Draw("E"); + func->Draw("SAME"); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHist2 draws the "hist" histogram in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHist2(TCanvas *canvas, TH2 *hist, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist->Draw(option); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHistStats draws the "hist" histogram and the "stats" stats in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHistStats(TCanvas *canvas, TH1 *hist, TPaveStats *stats, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist->Draw(option); + stats->Draw("SAME"); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawHistStats2 draws the "hist" histogram and the "stats" stats in the "canvas" canvas, by default with the "" option, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawHistStats2(TCanvas *canvas, TH2 *hist, TPaveStats *stats, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + hist->Draw(option); + stats->Draw("SAME"); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// InitGraph defines some attributes of the "graph" graph: the title "title" and the x-axis and y-axis labels "x" and "y". +void InitGraph(TGraph *graph, TString title, TString x, TString y) { + graph->SetTitle(title); + graph->GetXaxis()->SetTitle(x); + graph->GetYaxis()->SetTitle(y); + graph->SetFillColor(kWhite); + + return; +} + +// InitGraphErrors defines some attributes of the "graph" graph: the title "title" and the x-axis and y-axis labels "x" and "y". +void InitGraphErrors(TGraphErrors *graph, TString title, TString x, TString y) { + graph->SetTitle(title); + graph->GetXaxis()->SetTitle(x); + graph->GetYaxis()->SetTitle(y); + graph->SetFillColor(kWhite); + + return; +} + +// DrawGraph draws the "graph" graph in the "canvas" canvas, by default with the "APC" option, adds the "leg" legend, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +// "A" = with axes +// "P" = with points +// "C" = with smooth curves +void DrawGraph(TCanvas *canvas, TGraph *graph, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraph draws the "graph" graph in the "canvas" canvas, by default with the "APC" option, adds the "leg" legend, and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +// "A" = with axes +// "P" = with points +// "C" = with smooth curves +void DrawGraphErrors(TCanvas *canvas, TGraphErrors *graph, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraphFunc draws the "graph" graph and the "func" function in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawGraphFunc(TCanvas *canvas, TGraph *graph, TF1 *func, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + func->Draw("SAME"); + func->SetLineWidth(2); + func->SetLineStyle(1); + func->SetLineColor(kMagenta); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraphFunc2 draws the "graph" graph and the "func1" and "func2" functions in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawGraphFunc2(TCanvas *canvas, TGraph *graph, TF1 *func1, TF1 *func2, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + func1->Draw("SAME"); + func1->SetLineWidth(2); + func1->SetLineStyle(1); + func1->SetLineColor(kOrange); + func2->Draw("SAME"); + func2->SetLineWidth(2); + func2->SetLineStyle(1); + func2->SetLineColor(kMagenta); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraphErrorsFunc draws the "graph" graph and the "func" function in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawGraphErrorsFunc(TCanvas *canvas, TGraphErrors *graph, TF1 *func, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + func->Draw("SAME"); + func->SetLineWidth(2); + func->SetLineStyle(1); + func->SetLineColor(kMagenta); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraphErrorsFunc2 draws the "graph" graph and the "func1" and "func2" functions in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawGraphErrorsFunc2(TCanvas *canvas, TGraphErrors *graph, TF1 *func1, TF1 *func2, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option); + graph->SetFillColor(38); + graph->SetMarkerSize(1); + graph->SetMarkerStyle(20); + graph->SetMarkerColor(1); + graph->SetLineWidth(2); + graph->SetLineStyle(1); + graph->SetLineColor(2); + func1->Draw("SAME"); + func1->SetLineWidth(2); + func1->SetLineStyle(1); + func1->SetLineColor(kOrange); + func2->Draw("SAME"); + func2->SetLineWidth(2); + func2->SetLineStyle(1); + func2->SetLineColor(kMagenta); + gPad->Modified(); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawGraphCompare draws the "graph" graph in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file named "canvas.pdf". +void DrawGraphCompare(TCanvas *canvas, TMultiGraph *graph, TLegend *leg, TString title, TString x, TString y, TString option, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + graph->Draw(option.Data()); + leg->Draw(); + // graph->SetFillColor(38); + // graph->SetMarkerSize(1); + // graph->SetMarkerStyle(20); + // graph->SetMarkerColor(1); + // graph->SetLineWidth(2); + // graph->SetLineStyle(1); + // graph->SetLineColor(2); + graph->GetYaxis()->SetTitleOffset(0.8); + gPad->Modified(); + graph->SetTitle(title); + graph->GetXaxis()->SetTitle(x); + graph->GetYaxis()->SetTitle(y); + canvas->Update(); + canvas->Write(); + name = path + "/" + name + ".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// DrawFrame draws the "frame" frame in the "canvas" canvas and saves the result in the "folder" folder, in a pdf file names "canvas.pdf". +void DrawFrame(TCanvas *canvas, RooPlot *frame, TLegend *leg, bool logy, bool res, TString folder) { + TString path; + path = folder; + TString name = canvas->GetName(); + canvas->cd(); + canvas->SetTickx(1); + canvas->SetTicky(1); + + if (res == false) { + if (logy == true) + canvas->SetLogy(); + frame->Draw(); + leg->Draw(); + } + + else if (res == true) { + RooHist *histogram = frame->getHist("histogram"); + RooCurve *curve = frame->getCurve("curve"); + TH1F *hRes = CreateResidualHist(histogram,curve); + hRes->GetXaxis()->SetLabelSize((1./(1.+ .35)) * hRes->GetXaxis()->GetLabelSize()); + hRes->GetYaxis()->SetLabelSize((1./(1.+ .35)) * hRes->GetYaxis()->GetLabelSize()); + + canvas->Divide(1,2,.1,.1); + TPad *pHist = (TPad *)canvas->cd(1); + TPad *pRes = (TPad *)canvas->cd(2); + + pHist->SetPad(0.,0.25,1.,1.); + if (logy == true) + pHist->SetLogy(); + + pRes->SetPad(0.,0.,1.,0.25); + pRes->SetBottomMargin(0.3); + pRes->SetTopMargin(0.1); + + pHist->cd(); + frame->Draw(); + leg->Draw(); + pRes->cd(); + hRes->Draw(); + AddLinesToResidualHist(hRes); + hRes->Draw("SAME"); + } + + canvas->Update(); + canvas->Write(); + name = path+"/"+name+".pdf"; + canvas->SaveAs(name); + canvas->Close(); + + return; +} + +// CreateLegend2 creates and returns the legend for "obj1" (with "lab1" label and "option1" option) and "obj2" (with "lab2" label and "option2" option) objects. It allows to specify the drawing options ("lpfw" by default) and the posizion and size inside the canvas. +// l = lines +// p = points +// f = filled +// w = white +TLegend *CreateLegend2(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TString option, Double_t x1, Double_t y1, Double_t x2, Double_t y2) { + TLegend *leg = new TLegend(x1,y1,x2,y2); + leg->AddEntry(obj1,lab1,option1.Data()); + leg->AddEntry(obj2,lab2,option2.Data()); + leg->SetTextSize(0.05); + if(option.Contains("w")) + leg->SetFillColor(kWhite); + + return leg; +} + +// CreateLegend3 creates and returns the legend for 3 objects. +// l = lines +// p = points +// f = filled +// w = white +TLegend *CreateLegend3(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TString option, Double_t x1, Double_t y1, Double_t x2, Double_t y2) { + TLegend *leg = new TLegend(x1,y1,x2,y2); + leg->AddEntry(obj1,lab1,option1.Data()); + leg->AddEntry(obj2,lab2,option2.Data()); + leg->AddEntry(obj3,lab3,option3.Data()); + leg->SetTextSize(0.05); + if(option.Contains("w")) + leg->SetFillColor(kWhite); + + return leg; +} + +// CreateLegend4 creates and returns the legend for 4 objects. +// l = lines +// p = points +// f = filled +// w = white +TLegend *CreateLegend4(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TObject *obj4, TString lab4, TString option4, TString option, Double_t x1, Double_t y1, Double_t x2, Double_t y2) { + TLegend *leg = new TLegend(x1,y1,x2,y2); + leg->AddEntry(obj1,lab1,option1.Data()); + leg->AddEntry(obj2,lab2,option2.Data()); + leg->AddEntry(obj3,lab3,option3.Data()); + leg->AddEntry(obj4,lab4,option4.Data()); + leg->SetTextSize(0.05); + if(option.Contains("w")) + leg->SetFillColor(kWhite); + + return leg; +} + +// CreateLegend5 creates and returns the legend for 5 objects. +// l = lines +// p = points +// f = filled +// w = white +TLegend *CreateLegend5(TObject *obj1, TString lab1, TString option1, TObject *obj2, TString lab2, TString option2, TObject *obj3, TString lab3, TString option3, TObject *obj4, TString lab4, TString option4, TObject *obj5, TString lab5, TString option5, TString option, Double_t x1, Double_t y1, Double_t x2, Double_t y2) { + TLegend *leg = new TLegend(x1,y1,x2,y2); + leg->AddEntry(obj1,lab1,option1.Data()); + leg->AddEntry(obj2,lab2,option2.Data()); + leg->AddEntry(obj3,lab3,option3.Data()); + leg->AddEntry(obj4,lab4,option4.Data()); + leg->AddEntry(obj5,lab5,option5.Data()); + leg->SetTextSize(0.05); + if(option.Contains("w")) + leg->SetFillColor(kWhite); + + return leg; +} + +// CreateStats creates and returns the stats for "hist" histogram. It allows to specify the drawing options ("" by default). +// I = integral +// S = skewness +// K = kurtosis +// U = underflow +// O = overflow +TPaveStats *CreateStats(TH1 *hist, TString option) { + TPaveStats *stats = new TPaveStats(0.629396,0.657343,0.932953,0.879860,"brNDC"); + stats->SetName("stats"); + stats->SetBorderSize(1); + stats->SetFillColor(0); + stats->SetTextAlign(12); + TText *text = stats->AddText(Form("Entries = %g",hist->GetEntries())); + if(option.Contains("I")) + text = stats->AddText(Form("Integral = %g",hist->Integral())); + Int_t N = hist->GetDimension(); + if (N==1) { + text = stats->AddText(Form("Mean = %g",hist->GetMean())); + text = stats->AddText(Form("RMS = %g",hist->GetRMS())); + if(option.Contains("S")) + text = stats->AddText(Form("Skewness = %g",hist->GetSkewness())); + if(option.Contains("K")) + text = stats->AddText(Form("Kurtosis = %g",hist->GetKurtosis())); + if(option.Contains("U")) + text = stats->AddText(Form("Underflow = %f",hist->GetBinContent(0))); + if(option.Contains("O")) { + Int_t over = hist->GetNbinsX()+1; + text = stats->AddText(Form("Overflow = %f",hist->GetBinContent(over))); + } + } else { + text = stats->AddText(Form("Mean x = %g",hist->GetMean(1))); + text = stats->AddText(Form("Mean y = %g",hist->GetMean(2))); + text = stats->AddText(Form("RMS x = %g",hist->GetRMS(1))); + text = stats->AddText(Form("RMS y = %g",hist->GetRMS(2))); + if(option.Contains("S")) { + text = stats->AddText(Form("Skewness x = %g",hist->GetSkewness(1))); + text = stats->AddText(Form("Skewness y = %g",hist->GetSkewness(2))); + } + if(option.Contains("K")) { + text = stats->AddText(Form("Kurtosis x = %g",hist->GetKurtosis(1))); + text = stats->AddText(Form("Kurtosis y = %g",hist->GetKurtosis(2))); + } + if(option.Contains("U") || option.Contains("O")) { + Int_t nbinsX = hist->GetNbinsX(); + Int_t nbinsY = hist->GetNbinsY(); + Double_t bl = hist->GetBinContent(0,0); + Double_t bc = 0; for(int i=1; i<=nbinsX; i++) bc += hist->GetBinContent(i,0); + Double_t br = hist->GetBinContent(nbinsX+1,0); + Double_t cl = 0; for(int i=1; i<=nbinsY; i++) cl += hist->GetBinContent(0,i); + Double_t cc = hist->GetEffectiveEntries(); + Double_t cr = 0; for(int i=1; i<=nbinsX; i++) cr += hist->GetBinContent(nbinsX+1,i); + Double_t tl = hist->GetBinContent(0,nbinsY+1); + Double_t tc = 0; for(int i=1; i<=nbinsX; i++) tc += hist->GetBinContent(i,nbinsY+1); + Double_t tr = hist->GetBinContent(nbinsX+1,nbinsY+1); + text = stats->AddText(Form("%g| %g| %g",tl,tc,tr)); + text = stats->AddText(Form("%g| %g| %g",cl,cc,cr)); + text = stats->AddText(Form("%g| %g| %g",bl,bc,br)); + } + } + + return stats; +} + +// Residual evaluates the residual. +double Residual(Double_t datum, Double_t pdf) +{ + double chi2 = 0.; + + if ( pdf > 0 ) + chi2 += 2. * ( pdf - datum ); + + if ( datum > 0 && pdf > 0 ) + chi2 += 2. * datum * log( datum / pdf ); + + return ( ( datum >= pdf ) ? sqrt( chi2 ) : -sqrt( chi2 ) ); +} + +// AddLinesToResidualHist adds the lines to the histogram of the residuals. +void AddLinesToResidualHist(const TH1F *hist) +{ + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + + TLine *midLine = new TLine( xMin, 0., xMax, 0. ); + TLine *uppLine = new TLine( xMin, 3., xMax, 3. ); + TLine *lowLine = new TLine( xMin, -3., xMax, -3. ); + + uppLine->SetLineColor( kGray ); + lowLine->SetLineColor( kGray ); + + uppLine->SetLineStyle(7); + lowLine->SetLineStyle(7); + + midLine->Draw( "same" ); + uppLine->Draw( "same" ); + lowLine->Draw( "same" ); + + return; +} + +// CreateResidualHist returns the histogram of the residuals. +TH1F *CreateResidualHist(const RooHist *rhist, const RooCurve *curve) +{ + double r = 0.2; + double sr = 1. / r; + + // Grab info from the histogram. + int n = rhist->GetN(); + double* x = rhist->GetX(); + double* y = rhist->GetY(); + + // Create residual histogram. + double xMin = x[ 0 ]; + double xMax = x[ n - 1 ]; + TH1F* residuals_temp = new TH1F( "r", "", n, xMin, xMax ); + + double datum = 0.; + double pdf = 0.; + + // Fill the histogram. + if ( curve ) + for ( int bin = 0; bin < n; bin++ ) + { + datum = y[ bin ]; + pdf = curve->Eval( x[ bin ] ); + residuals_temp->SetBinContent( bin + 1, Residual( datum, pdf ) ); + residuals_temp->SetBinError ( bin + 1, 1. ); + } + + residuals_temp->SetMinimum ( -5. ); + residuals_temp->SetMaximum ( 5. ); + residuals_temp->SetStats ( false ); + residuals_temp->SetMarkerStyle( 8 ); + residuals_temp->SetMarkerSize ( .8 ); + + TAxis* xAxis = residuals_temp->GetXaxis(); + xAxis->SetTickLength ( sr * xAxis->GetTickLength() ); + xAxis->SetLabelSize ( sr * xAxis->GetLabelSize() ); + xAxis->SetTitleSize ( sr * xAxis->GetTitleSize() ); + xAxis->SetLabelOffset( sr * xAxis->GetLabelOffset() ); + + TAxis* yAxis = residuals_temp->GetYaxis(); + yAxis->SetNdivisions ( 504 ); + yAxis->SetLabelSize ( sr * yAxis->GetLabelSize() ); + + return residuals_temp; +} + +#endif diff --git a/Tools/Style.py b/Tools/Style.py new file mode 100644 index 0000000..0e43047 --- /dev/null +++ b/Tools/Style.py @@ -0,0 +1,531 @@ +""" +Author: Federica Lionetto +Date: November 16th, 2015 + +Description: +Define some useful functions to be used to accomplish common tasks with ROOT objects. +- LHCbStyle() +- InitCanvas(canvas,top,bottom,left,right) +- InitHist(hist,title,x,y) +- InitHist2(hist,title,x,y,z) +- InitGraph(graph,title,x,y) +- SetStyleObj(obj,markerColor,lineColor,fillColor,fillStyle,markerSize,markerStyle,lineWidth,lineStyle) +- AdjustOffsets(hist,axis,labelOffset,titleOffset) +- DrawObj(canvas,obj,leg,option,folder) +- DrawOjbCompare(canvas,objList,optionDict,leg,folder) +- CreateLegend(objList,labelDict,optionDict,x1,y1,x2,y2) +- CreateStats(hist,option,x1,y1,x2,y2) +- TestStyle() + +How to run it: +Include it in the other python scripts. +""" + +from array import array + +import numpy as np + +import ROOT + +def ColorList(n) : + default_list = [ROOT.kViolet, + ROOT.kTeal-5, + ROOT.kRed, + ROOT.kViolet+2, + ROOT.kAzure+1, + ROOT.kAzure-4, + ROOT.kOrange, + ROOT.kPink+9, + ROOT.kOrange+8, + ROOT.kGreen] + result = [] + if (n <= len(default_list)) : + result = default_list[:n] + else : + mod = n%len(default_list) + div = n/len(default_list) + for i in range(div) : + result = result + default_list + result = result + default_list[:mod] + print result + return result + +def LineStyleList(n) : + default_list = [1, + 2, + 3, + 4, + 5] + result = [] + for i in range(n) : + result.append(default_list[(i/10)%5]) + print result + return result + +def LHCbStyle() : + ''' + black = 1 + red = 2 + green = 3 + blue = 4 + yellow = 5 + magenta = 6 + cyan = 7 + purple = 9 + ''' + + lhcbFont = 132 + lhcbWidth = 2 + lhcbTSize = 0.06 + + # Fine. + ROOT.gROOT.SetStyle('Plain') + lhcbStyle = ROOT.TStyle('lhcbStyle','LHCb plots style') + + # Fine. + lhcbStyle.SetFillColor(0) + lhcbStyle.SetFillStyle(1001) + lhcbStyle.SetFrameFillColor(0) + lhcbStyle.SetFrameBorderMode(0) + lhcbStyle.SetPadBorderMode(0) + lhcbStyle.SetPadColor(0) + lhcbStyle.SetCanvasBorderMode(0) + lhcbStyle.SetCanvasColor(0) + lhcbStyle.SetStatColor(0) + lhcbStyle.SetLegendBorderSize(0) + lhcbStyle.SetLegendFont(132) + + # Fine. + # If you want the usual gradient palette (blue -> red) + lhcbStyle.SetPalette(1) + # If you want colors that correspond to gray scale in black and white + # colors=array('i',[0,5,7,3,6,2,4,1]) + # lhcbStyle.SetPalette(len(colors),colors) + + # Fine. + lhcbStyle.SetPaperSize(20,26) + lhcbStyle.SetPadTopMargin(0.05) + lhcbStyle.SetPadRightMargin(0.05) + lhcbStyle.SetPadBottomMargin(0.16) + lhcbStyle.SetPadLeftMargin(0.14) + + # Fine. + lhcbStyle.SetTextFont(lhcbFont) + lhcbStyle.SetTextSize(lhcbTSize) + lhcbStyle.SetLabelFont(lhcbFont,"x") + lhcbStyle.SetLabelFont(lhcbFont,"y") + lhcbStyle.SetLabelFont(lhcbFont,"z") + lhcbStyle.SetLabelSize(lhcbTSize,"x") + lhcbStyle.SetLabelSize(lhcbTSize,"y") + lhcbStyle.SetLabelSize(lhcbTSize,"z") + lhcbStyle.SetTitleFont(lhcbFont) + lhcbStyle.SetTitleFont(lhcbFont,"x") + lhcbStyle.SetTitleFont(lhcbFont,"y") + lhcbStyle.SetTitleFont(lhcbFont,"z") + lhcbStyle.SetTitleSize(1.2*lhcbTSize,"x") + lhcbStyle.SetTitleSize(1.2*lhcbTSize,"y") + lhcbStyle.SetTitleSize(1.2*lhcbTSize,"z") + + # Fine. + lhcbStyle.SetLineWidth(lhcbWidth) + lhcbStyle.SetFrameLineWidth(lhcbWidth) + lhcbStyle.SetHistLineWidth(lhcbWidth) + lhcbStyle.SetFuncWidth(lhcbWidth) + lhcbStyle.SetGridWidth(lhcbWidth) + lhcbStyle.SetLineStyleString(2,"[12 12]") + lhcbStyle.SetMarkerStyle(20) + lhcbStyle.SetMarkerSize(1.0) + + # Fine. + lhcbStyle.SetLabelOffset(0.010,"X") + lhcbStyle.SetLabelOffset(0.010,"Y") + + # Fine. + lhcbStyle.SetOptStat(0) + lhcbStyle.SetStatFormat("6.3g") + lhcbStyle.SetOptTitle(0) + lhcbStyle.SetOptFit(0) + lhcbStyle.SetTitleOffset(0.95,"X") + lhcbStyle.SetTitleOffset(0.95,"Y") + lhcbStyle.SetTitleOffset(1.2,"Z") + lhcbStyle.SetTitleFillColor(0) + lhcbStyle.SetTitleStyle(0) + lhcbStyle.SetTitleBorderSize(0) + lhcbStyle.SetTitleFont(lhcbFont,"title") + lhcbStyle.SetTitleX(0.0) + lhcbStyle.SetTitleY(1.0) + lhcbStyle.SetTitleW(1.0) + lhcbStyle.SetTitleH(0.05) + + # Fine. + lhcbStyle.SetStatBorderSize(0) + lhcbStyle.SetStatFont(lhcbFont) + lhcbStyle.SetStatFontSize(0.05) + lhcbStyle.SetStatX(0.9) + lhcbStyle.SetStatY(0.9) + lhcbStyle.SetStatW(0.25) + lhcbStyle.SetStatH(0.15) + + # Fine. + lhcbStyle.SetPadTickX(1) + lhcbStyle.SetPadTickY(1) + + # Fine. + lhcbStyle.SetNdivisions(505,"x") + lhcbStyle.SetNdivisions(510,"y") + + # Fine. + lhcbStyle.SetPaintTextFormat("5.2f") + + # Fine. + ROOT.gROOT.SetStyle("lhcbStyle") + # ROOT.gStyle = ROOT.gROOT.GetGlobal('gStyle',1) + ROOT.gROOT.ForceStyle() + + # Fine. + lhcbName = ROOT.TPaveText(ROOT.gStyle.GetPadLeftMargin() + 0.05,0.87 - ROOT.gStyle.GetPadTopMargin(),ROOT.gStyle.GetPadLeftMargin() + 0.20,0.95 - ROOT.gStyle.GetPadTopMargin(),'BRNDC') + lhcbName.AddText("LHCb") + lhcbName.SetFillColor(0) + lhcbName.SetTextAlign(12) + lhcbName.SetBorderSize(0) + + # Fine. + lhcbLabel = ROOT.TText() + lhcbLabel.SetTextFont(lhcbFont) + lhcbLabel.SetTextColor(1) + lhcbLabel.SetTextSize(lhcbTSize) + lhcbLabel.SetTextAlign(12) + + # Fine. + lhcbLatex = ROOT.TLatex() + lhcbLatex.SetTextFont(lhcbFont) + lhcbLatex.SetTextColor(1) + lhcbLatex.SetTextSize(lhcbTSize) + lhcbLatex.SetTextAlign(12) + + print '-------------------------' + print 'Set LHCb Style - Feb 2012' + print '-------------------------' + + return lhcbStyle + +def InitCanvas(canvas,top=0.15,bottom=None,left=None,right=None) : + """ + Initialise canvas with top margin , bottom margin , left margin and right margin . + """ + if top is not None : + canvas.SetTopMargin(top) + if bottom is not None : + canvas.SetBottomMargin(bottom) + if left is not None : + canvas.SetLeftMargin(left) + # Set the right margin to 0.19 for COLZ plots. + if right is not None : + canvas.SetRightMargin(right) + return + +def InitHist(hist,title='',x='',y='') : + """ + Initialise 1D histogram with title , x axis title <x>, and y axis title <y>. + If <y> is not given, a y axis title based on the number of entries per bin is set. + """ + hist.SetTitle(title) + hist.GetXaxis().SetTitle(x) + if (y=='') : + hist.GetYaxis().SetTitle('Entries / {0}'.format(round(hist.GetBinWidth(1),2))) + else : + hist.GetYaxis().SetTitle(y) + return + +def InitHist2(hist,title='',x='',y='',z='') : + """ + Initialise 2D histogram <hist> with title <title>, x axis title <x>, y axis title <y>, and z axis title <z>. + If <z> is not given, a z axis title based on the number of entries per bin is set. + """ + hist.SetTitle(title) + hist.GetXaxis().SetTitle(x) + hist.GetYaxis().SetTitle(y) + if (z=='') : + hist.GetZaxis().SetTitle('Entries / ({0}*{1})'.format(round(hist.GetXaxis().GetBinWidth(1),2),round(hist.GetYaxis().GetBinWidth(1),2))) + else : + hist.GetZaxis().SetTitle(z) + return + +def InitGraph(graph,title='',x='',y='') : + """ + Initialise graph <graph> with title <title>, x axis title <x>, and y axis title <y>. + """ + graph.SetTitle(title) + graph.GetXaxis().SetTitle(x) + graph.GetYaxis().SetTitle(y) + return + +def SetStyleObj(obj,markerColor=ROOT.kBlack,lineColor=ROOT.kMagenta,fillColor=None,fillStyle=1001,markerSize=1,markerStyle=20,lineWidth=2,lineStyle=1) : + """ + Set the style of object <obj>: + - marker color + - line color + - fill color + - marker size + - marker style + - line width + - line style + A default value exists for the last four. + """ + obj.SetMarkerColor(markerColor) + obj.SetLineColor(lineColor) + obj.SetFillStyle(fillStyle) + if fillColor is not None : + obj.SetFillColor(fillColor) + obj.SetMarkerSize(markerSize) + obj.SetMarkerStyle(markerStyle) + obj.SetLineWidth(lineWidth) + obj.SetLineStyle(lineStyle) + return + +def AdjustOffsets(hist,axis,labelOffset=0.005,titleOffset=0.85) : + """ + Adjust the offsets of histogram <hist> along the axis <axis>, by setting the label offset <labelOffset> and the title offset <titleOffset>. + The argument <axis> can be 'X', 'Y', or 'Z'. + The default values are those used for a TH2 object one wants to draw with the COLZ option. + """ + hist.SetLabelOffset(labelOffset,axis) + hist.SetTitleOffset(titleOffset,axis) + return + +def DrawObj(canvas,obj,leg=None,option='',folder='./') : + """ + Draw object <obj> in canvas <canvas>, together with legend <leg> (if provided), with option <option> and save the result in a pdf file having the same name as the canvas, located in the folder <folder>. + """ + name=canvas.GetName() + canvas.cd() + obj.Draw(option) + if ('A' in option) : + canvas.SetTickx(1) + canvas.SetTicky(1) + if (leg is not None) : + leg.Draw() + canvas.Update() + canvas.Write() + canvas.Print(folder+name+'.pdf') + canvas.Close() + return + +def DrawObjCompare(canvas,objList,optionDict,leg,folder='./') : + """ + Draw list of objects <objList> in canvas <canvas>, together with legend <leg>, with option <optionDict> and save the result in a pdf file having the same name as the canvas, located in the folder <folder>. + This allows to draw histograms, graphs, and functions in the same canvas. + However, since graphs do not have a 'GetMaximum()' method, the object to be drawn first is searched for among histograms and functions only (whose names start with 'h' and 'f', respectively). The remaining objects are drawn after. + """ + # print 'DrawObjCompare, optionDict: ', optionDict + + name=canvas.GetName() + canvas.cd() + + containsTMultiGraph = False + + maxList = [] + for obj in objList : + if (isinstance(obj,ROOT.TH1F) or isinstance(obj,ROOT.TF1)) : + maxList.append(obj.GetMaximum()) + if (isinstance(obj,ROOT.TMultiGraph)) : + containsTMultiGraph = True + + # print 'DrawObjCompare, maxList: ', maxList + + # print 'DrawObjCompare, containsTMultiGraph: ', containsTMultiGraph + + if containsTMultiGraph : + for obj in objList : + if (isinstance(obj,ROOT.TMultiGraph)) : + obj.Draw(optionDict[obj.GetName()]) + + if (len(maxList) != 0): + firstDraw = int(maxList.index(max(maxList))) + # print 'DrawObjCompare, firstDraw: ', firstDraw + + objList[firstDraw].Draw(optionDict[objList[firstDraw].GetName()]) + # print 'DrawObjCompare, objList: ', objList + for obj in objList : + if (obj != objList[firstDraw] and not isinstance(obj,ROOT.TMultiGraph)) : + obj.Draw(optionDict[obj.GetName()]+'SAME') + # print 'DrawObjCompare, obj: ', obj + + # print 'DrawObjCompare, optionDict: ', optionDict + + if leg is not None : + leg.Draw() + canvas.Update() + canvas.Write() + canvas.SaveAs(folder+name+'.pdf') + canvas.Close() + return + +def CreateLegend(objList,labelDict,optionDict,x1=0.64,y1=0.59,x2=0.94,y2=0.89) : + """ + Create legend <leg> starting from list of objects <objList>, dictionary of labels <labelDict>, and dictionary of options <optionDict>. + The legend is displayed in the position identified by <x1>, <y1>, <x2>, and <y2>. + If the option 'w' is specified, the fill color is set to white. + The text size is set to '0.05'. + """ + leg = ROOT.TLegend(x1,y1,x2,y2) + leg.SetFillColor(ROOT.kWhite) + for obj in objList : + label = labelDict[obj.GetName()] + option = optionDict[obj.GetName()] + leg.AddEntry(obj,label,option) + leg.SetTextSize(0.05) + return leg + +def CreateStats(hist,option,x1=0.64,y1=0.59,x2=0.94,y2=0.89) : + """ + Create stats <stats> of histogram <hist> depending on option <option>. + The stats are displayed in the position identified by <x1>, <y1>, <x2>, and <y2>. + """ + stats = ROOT.TPaveStats(x1,y1,x2,y2,'brNDC') + stats.SetName('stats') + stats.SetBorderSize(1) + stats.SetFillColor(0) + stats.SetTextAlign(12) + text = stats.AddText('Entries = {0}'.format(hist.GetEntries())) + if ('I' in option) : + text = stats.AddText('Integral = {0}'.format(hist.Integral())) + N = hist.GetDimension() + if (N == 1) : + text = stats.AddText('Mean = {0}'.format(hist.GetMean())); + text = stats.AddText('RMS = {0}'.format(hist.GetRMS())); + if ('S' in option) : + text = stats.AddText('Skewness = {0}'.format(hist.GetSkewness())) + if ('K' in option) : + text = stats.AddText('Kurtosis = {0}'.format(hist.GetKurtosis())) + if ('U' in option) : + text = stats.AddText('Underflow = {0}'.format(hist.GetBinContent(0))) + if ('O' in option) : + over = hist.GetNbinsX()+1 + text = stats.AddText('Overflow = {0}'.format(hist.GetBinContent(over))) + else : + text = stats.AddText('Mean x = {0}'.format(hist.GetMean(1))) + text = stats.AddText('Mean y = {0}'.format(hist.GetMean(2))) + text = stats.AddText('RMS x = {0}'.format(hist.GetRMS(1))) + text = stats.AddText('RMS y = {0}'.format(hist.GetRMS(2))) + if ('S' in option) : + text = stats.AddText('Skewness x = {0}'.format(hist.GetSkewness(1))) + text = stats.AddText('Skewness y = {0}'.format(hist.GetSkewness(2))) + if ('K' in option) : + text = stats.AddText('Kurtosis x = {0}'.format(hist.GetKurtosis(1))) + text = stats.AddText('Kurtosis y = {0}'.format(hist.GetKurtosis(2))) + if ('U' in option or 'O' in option) : + nbinsX = hist.GetNbinsX() + nbinsY = hist.GetNbinsY() + bl = hist.GetBinContent(0,0) + bc = 0 + for i in range(1,nbinsX+1) : + bc += hist.GetBinContent(i,0) + br = hist.GetBinContent(nbinsX+1,0) + cl = 0 + for i in range(1,nbinsY+1) : + cl += hist.GetBinContent(0,i) + cc = hist.GetEffectiveEntries() + cr = 0 + for i in range(1,nbinsX+1) : + cr += hist.GetBinContent(nbinsX+1,i) + tl = hist.GetBinContent(0,nbinsY+1) + tc = 0 + for i in range(1,nbinsX+1) : + tc += hist.GetBinContent(i,nbinsY+1) + tr = hist.GetBinContent(nbinsX+1,nbinsY+1) + text = stats.AddText('{0}| {1}| {2}'.format(tl,tc,tr)) + text = stats.AddText('{0}| {1}| {2}'.format(cl,cc,cr)) + text = stats.AddText('{0}| {1}| {2}'.format(bl,bc,br)) + return stats + +def TestStyle() : + """ + Test the functions defined above in a quick and automatic way. + """ + lhcbStyle = LHCbStyle() + + canvasDrawObj = ROOT.TCanvas('canvasDrawObj','canvasDrawObj',400,300) + canvasDrawObjCompare = ROOT.TCanvas('canvasDrawObjCompare','canvasDrawObjCompare',400,300) + canvasDrawObjCompare2 = ROOT.TCanvas('canvasDrawObjCompare2','canvasDrawObjCompare2',400,300) + + hist1 = ROOT.TH1F('hist1','hist1',100,0,10) + hist2 = ROOT.TH1F('hist2','hist2',100,0,10) + hist3 = ROOT.TH1F('hist3','hist3',100,0,10) + + hist2D = ROOT.TH2F('hist2D','hist2D',100,0,10,10,0,2) + + hist1.Fill(2,10) + hist2.Fill(3,2) + hist3.Fill(5,1) + + hist2D.Fill(2,1) + hist2D.Fill(4,0) + hist2D.Fill(5,2,10) + + ax = np.arange(0.,10.,2.) + ay = ax*5. + asize = len(ax) + print ax + print ay + print asize + graph = ROOT.TGraph(asize,ax,ay) + graph.SetName('graph') + + # Test InitHist. + InitHist(hist1,title='title for hist 1',x='x axis',y='y axis') + InitHist(hist2,title='title for hist 2',x='x axis') + InitHist(hist3,title='') + + # Test InitHist2. + InitHist2(hist2D) + + # Test InitGraph. + InitGraph(graph) + + markerColor=ROOT.kMagenta + lineColor=ROOT.kMagenta + fillColor=ROOT.kMagenta + # Test SetStyleObj. + SetStyleObj(hist1,markerColor,lineColor,fillColor) + SetStyleObj(graph,markerColor,lineColor,fillColor) + + # Test AdjustOffsets. + AdjustOffsets(hist2D,'Z') + + objList=[hist1,hist2,hist3] + colorDict={hist1.GetName():ROOT.kRed,hist2.GetName():ROOT.kBlue,hist3.GetName():ROOT.kGreen} + optionDrawDict={hist1.GetName():'E',hist2.GetName():'',hist3.GetName():'E'} + labelDict={hist1.GetName():'first',hist2.GetName():'second',hist3.GetName():'third'} + optionLegDict={hist1.GetName():'l',hist2.GetName():'p',hist3.GetName():'fw'} + + objList2=[hist1,graph] + colorDict2={hist1.GetName():ROOT.kRed,graph.GetName():ROOT.kBlue} + optionDrawDict2={hist1.GetName():'E',graph.GetName():'APC'} + labelDict2={hist1.GetName():'first',graph.GetName():'second'} + optionLegDict2={hist1.GetName():'l',graph.GetName():'pw'} + + print objList2 + print colorDict2 + print optionDrawDict2 + print labelDict2 + print optionLegDict2 + + # Test DrawObj. + DrawObj(canvasDrawObj,hist1) + + # Test CreateLegend. + leg = CreateLegend(objList,labelDict,optionLegDict,0.5,0.5,0.8,0.7) + leg2 = CreateLegend(objList2,labelDict2,optionLegDict2,0.6,0.2,0.9,0.5) + + # Test DrawObjCompare. + DrawObjCompare(canvasDrawObjCompare,objList,colorDict,optionDrawDict,leg) + DrawObjCompare(canvasDrawObjCompare2,objList2,colorDict2,optionDrawDict2,leg2) + + option = 'IUO' + # Test CreateStats. + CreateStats(hist1,option) + + return + +LHCbStyle() diff --git a/Tools/lhcbStyle.C b/Tools/lhcbStyle.C new file mode 100644 index 0000000..eac9d57 --- /dev/null +++ b/Tools/lhcbStyle.C @@ -0,0 +1,206 @@ +// Header guard. +#ifndef __LHCBSTYLE_C_INCLUDED__ +#define __LHCBSTYLE_C_INCLUDED__ + +#include "Lib.C" + +// all users - please change the name of this file to lhcbStyle.C +// Commits to lhcbdocs svn of .C files are not allowed +TStyle *getLHCbStyle(); + +TStyle *getLHCbStyle() +{ + + // define names for colours + Int_t black = 1; + Int_t red = 2; + Int_t green = 3; + Int_t blue = 4; + Int_t yellow = 5; + Int_t magenta= 6; + Int_t cyan = 7; + Int_t purple = 9; + + + //////////////////////////////////////////////////////////////////// + // PURPOSE: + // + // This macro defines a standard style for (black-and-white) + // "publication quality" LHCb ROOT plots. + // + // USAGE: + // + // Include the lines + // gROOT->ProcessLine(".L lhcbstyle.C"); + // lhcbStyle(); + // at the beginning of your root macro. + // + // Example usage is given in myPlot.C + // + // COMMENTS: + // + // Font: + // + // The font is chosen to be 132, this is Times New Roman (like the text of + // your document) with precision 2. + // + // "Landscape histograms": + // + // The style here is designed for more or less square plots. + // For longer histograms, or canvas with many pads, adjustements are needed. + // For instance, for a canvas with 1x5 histograms: + // TCanvas* c1 = new TCanvas("c1", "L0 muons", 600, 800); + // c1->Divide(1,5); + // Adaptions like the following will be needed: + // gStyle->SetTickLength(0.05,"x"); + // gStyle->SetTickLength(0.01,"y"); + // gStyle->SetLabelSize(0.15,"x"); + // gStyle->SetLabelSize(0.1,"y"); + // gStyle->SetStatW(0.15); + // gStyle->SetStatH(0.5); + // + // Authors: Thomas Schietinger, Andrew Powell, Chris Parkes, Niels Tuning + // Maintained by Editorial board member (currently Niels) + /////////////////////////////////////////////////////////////////// + + // Use times new roman, precision 2 + Int_t lhcbFont = 132; // Old LHCb style: 62; + // Line thickness + Double_t lhcbWidth = 2.00; // Old LHCb style: 3.00; + // Text size + Double_t lhcbTSize = 0.06; + + // use plain black on white colors + gROOT->SetStyle("Plain"); + TStyle *lhcbStyle= new TStyle("lhcbStyle","LHCb plots style"); + + //lhcbStyle->SetErrorX(0); // don't suppress the error bar along X + + lhcbStyle->SetFillColor(1); + lhcbStyle->SetFillStyle(1001); // solid + lhcbStyle->SetFrameFillColor(0); + lhcbStyle->SetFrameBorderMode(0); + lhcbStyle->SetPadBorderMode(0); + lhcbStyle->SetPadColor(0); + lhcbStyle->SetCanvasBorderMode(0); + lhcbStyle->SetCanvasColor(0); + lhcbStyle->SetStatColor(0); + lhcbStyle->SetLegendBorderSize(0); + lhcbStyle->SetLegendFont(132); + + // If you want the usual gradient palette (blue -> red) + lhcbStyle->SetPalette(1); + // If you want colors that correspond to gray scale in black and white: + int colors[8] = {0,5,7,3,6,2,4,1}; + lhcbStyle->SetPalette(8,colors); + + // set the paper & margin sizes + lhcbStyle->SetPaperSize(20,26); + lhcbStyle->SetPadTopMargin(0.05); + lhcbStyle->SetPadRightMargin(0.05); // increase for colz plots + lhcbStyle->SetPadBottomMargin(0.16); + lhcbStyle->SetPadLeftMargin(0.14); + + // use large fonts + lhcbStyle->SetTextFont(lhcbFont); + lhcbStyle->SetTextSize(lhcbTSize); + lhcbStyle->SetLabelFont(lhcbFont,"x"); + lhcbStyle->SetLabelFont(lhcbFont,"y"); + lhcbStyle->SetLabelFont(lhcbFont,"z"); + lhcbStyle->SetLabelSize(lhcbTSize,"x"); + lhcbStyle->SetLabelSize(lhcbTSize,"y"); + lhcbStyle->SetLabelSize(lhcbTSize,"z"); + lhcbStyle->SetTitleFont(lhcbFont); + lhcbStyle->SetTitleFont(lhcbFont,"x"); + lhcbStyle->SetTitleFont(lhcbFont,"y"); + lhcbStyle->SetTitleFont(lhcbFont,"z"); + lhcbStyle->SetTitleSize(1.2*lhcbTSize,"x"); + lhcbStyle->SetTitleSize(1.2*lhcbTSize,"y"); + lhcbStyle->SetTitleSize(1.2*lhcbTSize,"z"); + + // use medium bold lines and thick markers + lhcbStyle->SetLineWidth(lhcbWidth); + lhcbStyle->SetFrameLineWidth(lhcbWidth); + lhcbStyle->SetHistLineWidth(lhcbWidth); + lhcbStyle->SetFuncWidth(lhcbWidth); + lhcbStyle->SetGridWidth(lhcbWidth); + lhcbStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes + lhcbStyle->SetMarkerStyle(20); + lhcbStyle->SetMarkerSize(1.0); + + // label offsets + lhcbStyle->SetLabelOffset(0.010,"X"); + lhcbStyle->SetLabelOffset(0.010,"Y"); + + // by default, do not display histogram decorations: + lhcbStyle->SetOptStat(0); + // lhcbStyle->SetOptStat("emr"); // show only nent -e , mean - m , rms -r + // full opts at http://root.cern.ch/root/html/TStyle.html#TStyle:SetOptStat + lhcbStyle->SetStatFormat("6.3g"); // specified as c printf options + lhcbStyle->SetOptTitle(0); + lhcbStyle->SetOptFit(0); + //lhcbStyle->SetOptFit(1011); // order is probability, Chi2, errors, parameters + //titles + lhcbStyle->SetTitleOffset(0.95,"X"); + lhcbStyle->SetTitleOffset(0.95,"Y"); + lhcbStyle->SetTitleOffset(1.2,"Z"); + lhcbStyle->SetTitleFillColor(0); + lhcbStyle->SetTitleStyle(0); + lhcbStyle->SetTitleBorderSize(0); + lhcbStyle->SetTitleFont(lhcbFont,"title"); + lhcbStyle->SetTitleX(0.0); + lhcbStyle->SetTitleY(1.0); + lhcbStyle->SetTitleW(1.0); + lhcbStyle->SetTitleH(0.05); + + // look of the statistics box: + lhcbStyle->SetStatBorderSize(0); + lhcbStyle->SetStatFont(lhcbFont); + lhcbStyle->SetStatFontSize(0.05); + lhcbStyle->SetStatX(0.9); + lhcbStyle->SetStatY(0.9); + lhcbStyle->SetStatW(0.25); + lhcbStyle->SetStatH(0.15); + + // put tick marks on top and RHS of plots + lhcbStyle->SetPadTickX(1); + lhcbStyle->SetPadTickY(1); + + // histogram divisions: only 5 in x to avoid label overlaps + lhcbStyle->SetNdivisions(505,"x"); + lhcbStyle->SetNdivisions(510,"y"); + + gROOT->SetStyle("lhcbStyle"); + gROOT->ForceStyle(); + + // add LHCb label + TPaveText *lhcbName = new TPaveText(gStyle->GetPadLeftMargin() + 0.05, + 0.87 - gStyle->GetPadTopMargin(), + gStyle->GetPadLeftMargin() + 0.20, + 0.95 - gStyle->GetPadTopMargin(), + "BRNDC"); + lhcbName->AddText("LHCb"); + lhcbName->SetFillColor(0); + lhcbName->SetTextAlign(12); + lhcbName->SetBorderSize(0); + + TText *lhcbLabel = new TText(); + lhcbLabel->SetTextFont(lhcbFont); + lhcbLabel->SetTextColor(1); + lhcbLabel->SetTextSize(lhcbTSize); + lhcbLabel->SetTextAlign(12); + + TLatex *lhcbLatex = new TLatex(); + lhcbLatex->SetTextFont(lhcbFont); + lhcbLatex->SetTextColor(1); + lhcbLatex->SetTextSize(lhcbTSize); + lhcbLatex->SetTextAlign(12); + + cout << "-------------------------" << endl; + cout << "Set LHCb Style - Feb 2012" << endl; + cout << "-------------------------" << endl; + + return lhcbStyle; +} + +#endif