1 |
buchmann |
1.1 |
#include <iostream>
|
2 |
|
|
#include <sstream>
|
3 |
|
|
#include <vector>
|
4 |
|
|
#include <stdio.h>
|
5 |
|
|
#include <stdlib.h>
|
6 |
|
|
#include <sys/types.h>
|
7 |
|
|
#include <sys/stat.h>
|
8 |
|
|
|
9 |
|
|
#include <TFile.h>
|
10 |
|
|
#include <TTree.h>
|
11 |
|
|
#include <TCut.h>
|
12 |
|
|
#include <TLegend.h>
|
13 |
|
|
#include <TLatex.h>
|
14 |
|
|
#include <TText.h>
|
15 |
|
|
#include <TGraph.h>
|
16 |
|
|
#include <TH1.h>
|
17 |
buchmann |
1.2 |
#include <TF1.h>
|
18 |
buchmann |
1.1 |
#include <TMath.h>
|
19 |
|
|
#include <TStyle.h>
|
20 |
|
|
#include <TCanvas.h>
|
21 |
|
|
#include <TError.h>
|
22 |
buchmann |
1.3 |
#include <TVirtualPad.h>
|
23 |
buchmann |
1.1 |
#include <TGraphAsymmErrors.h>
|
24 |
|
|
#include <TRandom.h>
|
25 |
|
|
#ifndef Verbosity
|
26 |
|
|
#define Verbosity 0
|
27 |
|
|
#endif
|
28 |
|
|
|
29 |
|
|
/*
|
30 |
|
|
#ifndef SampleClassLoaded
|
31 |
|
|
#include "SampleClass.C"
|
32 |
|
|
#endif
|
33 |
|
|
*/
|
34 |
|
|
#define GeneralToolBoxLoaded
|
35 |
|
|
|
36 |
|
|
using namespace std;
|
37 |
|
|
|
38 |
|
|
bool dopng=false;
|
39 |
|
|
bool doC=false;
|
40 |
|
|
bool doeps=false;
|
41 |
|
|
string basedirectory="";
|
42 |
|
|
|
43 |
|
|
TLegend* make_legend(string title);
|
44 |
|
|
TText* write_title(string title);
|
45 |
|
|
TText* write_title_low(string title);
|
46 |
|
|
|
47 |
|
|
TText* write_text(float xpos,float ypos,string title);
|
48 |
|
|
float computeRatioError(float a, float da, float b, float db);
|
49 |
|
|
float computeProductError(float a, float da, float b, float db);
|
50 |
|
|
TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning);
|
51 |
|
|
void setlumi(float l);
|
52 |
|
|
void CompleteSave(TCanvas *can, string filename, bool feedback);
|
53 |
buchmann |
1.3 |
void CompleteSave(TVirtualPad *can, string filename, bool feedback);
|
54 |
buchmann |
1.1 |
void write_warning(string funcname, string text);
|
55 |
|
|
void write_error(string funcname, string text);
|
56 |
|
|
void write_info(string funcname, string text);
|
57 |
|
|
//-------------------------------------------------------------------------------------
|
58 |
|
|
float lumi;
|
59 |
|
|
template<class A>
|
60 |
|
|
string any2string(const A& a){
|
61 |
|
|
ostringstream out;
|
62 |
|
|
out << a;
|
63 |
|
|
return out.str();
|
64 |
|
|
}
|
65 |
|
|
|
66 |
|
|
void do_png(bool s) { dopng=s;}
|
67 |
|
|
void do_eps(bool s) { doeps=s;}
|
68 |
|
|
void do_C(bool s) { doC=s;}
|
69 |
|
|
|
70 |
|
|
string topdir(string child) {
|
71 |
|
|
string tempdirectory=child;
|
72 |
|
|
if(tempdirectory.substr(tempdirectory.length()-1,1)=="/") tempdirectory=tempdirectory.substr(0,tempdirectory.length());
|
73 |
|
|
//we now have a directory without the trailing slash so we can just look for the last non-slash character :-)
|
74 |
|
|
for(int ichar=tempdirectory.length()-1;ichar>=0;ichar--) {
|
75 |
|
|
if(tempdirectory.substr(ichar,1)=="/") {
|
76 |
|
|
return tempdirectory.substr(0,ichar);
|
77 |
|
|
}
|
78 |
|
|
}
|
79 |
|
|
}
|
80 |
|
|
|
81 |
|
|
void ensure_directory_exists(string thisdirectory) {
|
82 |
|
|
struct stat st;
|
83 |
|
|
if(stat(thisdirectory.c_str(),&st) == 0) {
|
84 |
|
|
if(Verbosity>0) cout << "Directory " << thisdirectory << " exists!" << endl;
|
85 |
|
|
}
|
86 |
|
|
else {
|
87 |
|
|
if(Verbosity>0) cout << "Directory " << thisdirectory << " does not exist. Need to create it!" << endl;
|
88 |
|
|
ensure_directory_exists(topdir(thisdirectory));
|
89 |
|
|
if (mkdir(thisdirectory.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH))
|
90 |
|
|
if(Verbosity>0) cout << "Created the directory " << thisdirectory << endl;
|
91 |
|
|
}
|
92 |
|
|
}
|
93 |
|
|
|
94 |
|
|
void set_directory(string basedir="") {
|
95 |
|
|
if(basedir.substr(0,1)=="/") basedir=basedir.substr(1,basedir.length()-1);
|
96 |
|
|
if(basedir.substr(basedir.length()-1,1)!="/") basedir+="/";
|
97 |
|
|
char currentpath[1024];
|
98 |
|
|
getcwd(currentpath,1024);
|
99 |
|
|
basedirectory=(string)currentpath+"/"+basedir;
|
100 |
|
|
ensure_directory_exists(basedirectory);
|
101 |
|
|
}
|
102 |
|
|
|
103 |
|
|
string extract_directory(string savethis) {
|
104 |
|
|
bool foundslash=false;
|
105 |
|
|
int position=savethis.length();
|
106 |
|
|
while(!foundslash&&position>0) {
|
107 |
|
|
position--;
|
108 |
|
|
if(savethis.substr(position,1)=="/") foundslash=true;
|
109 |
|
|
}
|
110 |
|
|
if(position>0) return savethis.substr(0,position+1);
|
111 |
|
|
else return "";
|
112 |
|
|
}
|
113 |
|
|
|
114 |
|
|
void CompleteSave(TCanvas *can, string filename, bool feedback=false) {
|
115 |
buchmann |
1.3 |
//any change you make here should also be done below in the CompleteSave function for virtual pads
|
116 |
buchmann |
1.1 |
Int_t currlevel=gErrorIgnoreLevel;
|
117 |
|
|
if(!feedback) gErrorIgnoreLevel=1001;
|
118 |
|
|
ensure_directory_exists(extract_directory(basedirectory+filename));
|
119 |
|
|
if(dopng) can->SaveAs((basedirectory+filename+".png").c_str());
|
120 |
|
|
if(doeps) can->SaveAs((basedirectory+filename+".eps").c_str());
|
121 |
|
|
if(doC) can->SaveAs((basedirectory+filename+".C").c_str());
|
122 |
|
|
gErrorIgnoreLevel=currlevel;
|
123 |
buchmann |
1.3 |
cout << "Saved " << filename << " in all requested formats" << endl;
|
124 |
buchmann |
1.1 |
}
|
125 |
buchmann |
1.3 |
|
126 |
|
|
void CompleteSave(TVirtualPad *can, string filename, bool feedback=false) {
|
127 |
|
|
Int_t currlevel=gErrorIgnoreLevel;
|
128 |
|
|
if(!feedback) gErrorIgnoreLevel=1001;
|
129 |
|
|
ensure_directory_exists(extract_directory(basedirectory+filename));
|
130 |
|
|
if(dopng) can->SaveAs((basedirectory+filename+".png").c_str());
|
131 |
|
|
if(doeps) can->SaveAs((basedirectory+filename+".eps").c_str());
|
132 |
|
|
if(doC) can->SaveAs((basedirectory+filename+".C").c_str());
|
133 |
|
|
gErrorIgnoreLevel=currlevel;
|
134 |
|
|
cout << "Saved " << filename << " in all requested formats" << endl;
|
135 |
|
|
}
|
136 |
|
|
|
137 |
buchmann |
1.1 |
|
138 |
|
|
void setlumi(float l) {
|
139 |
|
|
lumi=l;
|
140 |
|
|
}
|
141 |
|
|
|
142 |
|
|
int write_first_line(vector<vector<string> > &entries) {
|
143 |
|
|
if(entries.size()>0) {
|
144 |
|
|
vector<string> firstline = entries[0];
|
145 |
|
|
int ncolumns=firstline.size();
|
146 |
|
|
int ndividers=ncolumns+1;
|
147 |
|
|
int cellwidth=(int)(((float)(60-ndividers))/(ncolumns));
|
148 |
|
|
cout << " |";
|
149 |
|
|
for(int idiv=0;idiv<ncolumns;idiv++) {
|
150 |
|
|
for(int isig=0;isig<cellwidth;isig++) cout << "-";
|
151 |
|
|
cout << "|";
|
152 |
|
|
}
|
153 |
|
|
cout << endl;
|
154 |
|
|
return ncolumns;
|
155 |
|
|
} else {
|
156 |
|
|
return 0;
|
157 |
|
|
}
|
158 |
|
|
}
|
159 |
|
|
|
160 |
|
|
void write_entry(string entry,int width,int iline=0,int ientry=0) {
|
161 |
|
|
int currwidth=entry.size();
|
162 |
|
|
while(currwidth<width) {
|
163 |
|
|
entry=" "+entry;
|
164 |
|
|
if(entry.size()<width) entry=entry+" ";
|
165 |
|
|
currwidth=entry.size();
|
166 |
|
|
}
|
167 |
|
|
bool do_special=false;
|
168 |
|
|
if(iline==1&&ientry==1) { cout << "\033[1;32m" << entry << "\033[0m|";do_special=true;}//observed
|
169 |
|
|
if(iline==1&&ientry==2) { cout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
|
170 |
|
|
if(iline==2&&ientry==1) { cout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
|
171 |
|
|
if(iline==2&&ientry==2) { cout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
|
172 |
|
|
if(!do_special) cout << entry << "|";
|
173 |
|
|
}
|
174 |
|
|
|
175 |
|
|
void make_nice_table(vector<vector <string> > &entries) {
|
176 |
|
|
int ncolumns=write_first_line(entries);
|
177 |
|
|
int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
|
178 |
|
|
for(int iline=0;iline<entries.size();iline++) {
|
179 |
|
|
vector<string> currline = entries[iline];
|
180 |
|
|
cout << " |";
|
181 |
|
|
for(int ientry=0;ientry<currline.size();ientry++) {
|
182 |
|
|
write_entry(currline[ientry],cellwidth);
|
183 |
|
|
}
|
184 |
|
|
cout << endl;
|
185 |
|
|
if(iline==0) write_first_line(entries);
|
186 |
|
|
}
|
187 |
|
|
write_first_line(entries);
|
188 |
|
|
}
|
189 |
|
|
|
190 |
|
|
void make_nice_jzb_table(vector<vector <string> > &entries) {
|
191 |
|
|
int ncolumns=write_first_line(entries);
|
192 |
|
|
int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
|
193 |
|
|
for(int iline=0;iline<entries.size();iline++) {
|
194 |
|
|
vector<string> currline = entries[iline];
|
195 |
|
|
cout << " |";
|
196 |
|
|
for(int ientry=0;ientry<currline.size();ientry++) {
|
197 |
|
|
write_entry(currline[ientry],cellwidth,iline,ientry);
|
198 |
|
|
}
|
199 |
|
|
cout << endl;
|
200 |
|
|
if(iline==0) write_first_line(entries);
|
201 |
|
|
}
|
202 |
|
|
write_first_line(entries);
|
203 |
|
|
}
|
204 |
|
|
|
205 |
|
|
|
206 |
|
|
void write_warning(string funcname, string text) {
|
207 |
|
|
cout << endl << endl;
|
208 |
|
|
cout << "\033[1;33m" << " _ " << endl;
|
209 |
|
|
cout << "\033[1;33m" << " (_) " << endl;
|
210 |
|
|
cout << "\033[1;33m" << "__ ____ _ _ __ _ __ _ _ __ __ _ " << endl;
|
211 |
|
|
cout << "\033[1;33m" << "\\ \\ /\\ / / _` | '__| '_ \\| | '_ \\ / _` |" << endl;
|
212 |
|
|
cout << "\033[1;33m" << " \\ V V / (_| | | | | | | | | | | (_| |" << endl;
|
213 |
|
|
cout << "\033[1;33m" << " \\_/\\_/ \\__,_|_| |_| |_|_|_| |_|\\__, |" << endl;
|
214 |
|
|
cout << "\033[1;33m" << " __/ |" << endl;
|
215 |
|
|
cout << "\033[1;33m" << " |___/ " << endl;
|
216 |
|
|
cout << endl;
|
217 |
|
|
cout << "\033[1;33m [" << funcname << "] " << text << " \033[0m" << endl;
|
218 |
|
|
cout << endl << endl;
|
219 |
|
|
}
|
220 |
|
|
void write_error(string funcname, string text) {
|
221 |
|
|
cout << endl << endl;
|
222 |
|
|
cout << "\033[1;31m ___ _ __ _ __ ___ _ __ " << endl;
|
223 |
|
|
cout << "\033[1;31m / _ \\ __| __/ _ \\| '__|" << endl;
|
224 |
|
|
cout << "\033[1;31m| __/ | | | | (_) | | " << endl;
|
225 |
|
|
cout << "\033[1;31m \\___|_| |_| \\___/|_| " << endl;
|
226 |
|
|
cout << endl;
|
227 |
|
|
cout << "\033[1;31m [" << funcname << "] " << text << " \033[0m" << endl;
|
228 |
|
|
cerr << "\033[1;31m [" << funcname << "] " << text << " \033[0m" << endl;
|
229 |
|
|
cout << endl << endl;
|
230 |
|
|
}
|
231 |
|
|
|
232 |
|
|
void write_info(string funcname, string text) {
|
233 |
|
|
cout << endl << endl;
|
234 |
|
|
cout << "\033[1;34m _____ __ " << endl;
|
235 |
|
|
cout << "\033[1;34m |_ _| / _| " << endl;
|
236 |
|
|
cout << "\033[1;34m | | _ __ | |_ ___ " << endl;
|
237 |
|
|
cout << "\033[1;34m | | | '_ \\| _/ _ \\ " << endl;
|
238 |
|
|
cout << "\033[1;34m _| |_| | | | || (_) | " << endl;
|
239 |
|
|
cout << "\033[1;34m |_____|_| |_|_| \\___/ " << endl;
|
240 |
|
|
cout << endl;
|
241 |
|
|
cout << "\033[1;34m [" << funcname << "] " << text << " \033[0m" << endl;
|
242 |
|
|
cerr << "\033[1;34m [" << funcname << "] " << text << " \033[0m" << endl;
|
243 |
|
|
cout << endl << endl;
|
244 |
|
|
}
|
245 |
|
|
|
246 |
|
|
TText* write_text(float xpos,float ypos,string title)
|
247 |
|
|
{
|
248 |
|
|
TLatex* titlebox = new TLatex (xpos,ypos,title.c_str());
|
249 |
|
|
titlebox->SetNDC(true);
|
250 |
|
|
titlebox->SetTextFont(42);
|
251 |
|
|
titlebox->SetTextSize(0.04);
|
252 |
|
|
titlebox->SetTextAlign(21);
|
253 |
|
|
return titlebox;
|
254 |
|
|
}
|
255 |
|
|
|
256 |
|
|
TText* write_title(string title)
|
257 |
|
|
{
|
258 |
|
|
TText* titlebox = write_text(0.5,0.96,title);
|
259 |
|
|
return titlebox;
|
260 |
|
|
}
|
261 |
|
|
|
262 |
|
|
TText* write_title_low(string title)
|
263 |
|
|
{
|
264 |
|
|
TText* titlebox = write_text(0.5,0.94,title);
|
265 |
|
|
return titlebox;
|
266 |
|
|
}
|
267 |
|
|
|
268 |
|
|
TLegend* make_legend(string title="")
|
269 |
|
|
{
|
270 |
|
|
// TLegend *leg = new TLegend(0.65,0.65,0.89,0.89);
|
271 |
|
|
gStyle->SetTextFont(42);
|
272 |
|
|
TLegend *leg = new TLegend(0.65,0.65,0.89,0.77);
|
273 |
|
|
if(title!="") leg->SetHeader(title.c_str());
|
274 |
|
|
leg->SetTextFont(42);
|
275 |
|
|
leg->SetFillColor(kWhite);
|
276 |
|
|
leg->SetBorderSize(0);
|
277 |
|
|
leg->SetLineColor(kWhite);
|
278 |
|
|
TText *writeline1 = write_text(0.77,0.87,"CMS Preliminary 2011");
|
279 |
|
|
stringstream lumitext;
|
280 |
|
|
lumitext<<"#sqrt{s}=7, L="<<lumi<<" pb^{-1}";
|
281 |
|
|
TText *writeline2 = write_text(0.77,0.83,lumitext.str().c_str());
|
282 |
|
|
writeline1->SetTextSize(0.03);
|
283 |
|
|
writeline2->SetTextSize(0.03);
|
284 |
|
|
writeline1->Draw();
|
285 |
|
|
writeline2->Draw();
|
286 |
|
|
return leg;
|
287 |
|
|
}
|
288 |
|
|
|
289 |
|
|
TGraph* make_nice_ratio(int nbins,float binning[],TH1F* histo)
|
290 |
|
|
{
|
291 |
|
|
float errorsquared[nbins];
|
292 |
|
|
float errors[nbins];
|
293 |
|
|
float bincontent[nbins];
|
294 |
|
|
for (int i=0;i<nbins;i++) {
|
295 |
|
|
errorsquared[i]=0;
|
296 |
|
|
bincontent[i]=0;
|
297 |
|
|
errors[i]=0;
|
298 |
|
|
}
|
299 |
|
|
float currlimit=binning[0];
|
300 |
|
|
int currtoplim=1;
|
301 |
|
|
for(int ibin=1;ibin<=histo->GetNbinsX();ibin++)
|
302 |
|
|
{
|
303 |
|
|
if(binning[currtoplim]<histo->GetBinCenter(ibin)) currtoplim++;
|
304 |
|
|
cout << "Bin i=" << ibin << " with bin center " << histo->GetBinCenter(ibin) << " contains " << histo->GetBinContent(ibin) << " is within " << binning[currtoplim-1] << " and " << binning[currtoplim] << endl;
|
305 |
|
|
|
306 |
|
|
}
|
307 |
|
|
|
308 |
|
|
return 0;
|
309 |
|
|
}
|
310 |
|
|
|
311 |
|
|
float statErrorN(float x){return x - 0.5*TMath::ChisquareQuantile(0.3173/2,2*x);}
|
312 |
|
|
float statErrorP(float x){return 0.5*TMath::ChisquareQuantile(1-0.3173/2,2*(x+1))-x;}
|
313 |
|
|
float lowLimit(float a, float x){return 0.5*TMath::ChisquareQuantile(a,2*x);}
|
314 |
|
|
float highLimit(float a,float x){return 0.5*TMath::ChisquareQuantile(1-a,2*(x+1));}
|
315 |
|
|
|
316 |
|
|
float computeRatioError(float a, float da, float b, float db)
|
317 |
|
|
{
|
318 |
|
|
float val=0.;
|
319 |
|
|
float errorSquare = (a/b)*(a/b)*( (da/a)*(da/a) + (db/b)*(db/b));
|
320 |
|
|
val = TMath::Sqrt(errorSquare);
|
321 |
|
|
return val;
|
322 |
|
|
|
323 |
|
|
}
|
324 |
|
|
float computeProductError(float a, float da, float b, float db)
|
325 |
|
|
{
|
326 |
|
|
float val=0.;
|
327 |
|
|
float errorSquare = (a*b)*(a*b)*( (da/a)*(da/a) + (db/b)*(db/b));
|
328 |
|
|
val = TMath::Sqrt(errorSquare);
|
329 |
|
|
return val;
|
330 |
|
|
}
|
331 |
|
|
|
332 |
|
|
TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning)
|
333 |
|
|
{
|
334 |
|
|
int absJZBbinsNumber = binning.size()-1;
|
335 |
|
|
TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
|
336 |
|
|
|
337 |
|
|
for(unsigned int i=0;i<absJZBbinsNumber;i++)
|
338 |
|
|
{
|
339 |
|
|
float xCenter=h1->GetBinCenter(i+1);
|
340 |
|
|
float xWidth=(h1->GetBinWidth(i+1))*0.5;
|
341 |
|
|
float nominatorError = h1->GetBinError(i+1);
|
342 |
|
|
float nominator=h1->GetBinContent(i+1);
|
343 |
|
|
float denominatorError=h2->GetBinError(i+1);
|
344 |
|
|
float denominator=h2->GetBinContent(i+1);
|
345 |
|
|
float errorN = 0;
|
346 |
|
|
float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
|
347 |
|
|
if(id==1) // (is data)
|
348 |
|
|
{
|
349 |
|
|
errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
|
350 |
|
|
errorN = errorP; // symmetrize using statErrorP
|
351 |
|
|
} else {
|
352 |
|
|
errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
|
353 |
|
|
errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
|
354 |
|
|
}
|
355 |
|
|
if(denominator!=0) {
|
356 |
|
|
graph->SetPoint(i, xCenter, nominator/denominator);
|
357 |
|
|
graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
|
358 |
|
|
}
|
359 |
|
|
else {
|
360 |
|
|
graph->SetPoint(i, xCenter, -999);
|
361 |
|
|
graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
|
362 |
|
|
}
|
363 |
|
|
}
|
364 |
|
|
return graph;
|
365 |
|
|
}
|
366 |
|
|
|
367 |
|
|
string print_range(float cent, float down, float up) {//note that up&down can be flipped, we don't care, but the central value needs to come 1st!
|
368 |
|
|
float uperr=0,downerr=0;
|
369 |
|
|
if(down>up&&down>cent) uperr=down-cent;
|
370 |
|
|
if(up>down&&up>cent) uperr=up-cent;
|
371 |
|
|
if(down<cent&&down<up) downerr=cent-down;
|
372 |
|
|
if(up<cent&&up<down) downerr=cent-up;
|
373 |
|
|
if(cent>up&¢>down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE LARGER THAN BOTH UP&DOWN!");
|
374 |
|
|
if(cent<up&¢<down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE SMALLER THAN BOTH UP&DOWN!");
|
375 |
|
|
stringstream result;
|
376 |
|
|
result << cent << " + " << uperr << " - " << downerr;
|
377 |
|
|
return result.str();
|
378 |
|
|
}
|
379 |
|
|
|
380 |
|
|
void bubbleSort ( int arr [ ], int size, int order [ ]) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
|
381 |
|
|
{
|
382 |
|
|
int last = size - 2;
|
383 |
|
|
int isChanged = 1;
|
384 |
|
|
|
385 |
|
|
while ( last >= 0 && isChanged )
|
386 |
|
|
{
|
387 |
|
|
isChanged = 0;
|
388 |
|
|
for ( int k = 0; k <= last; k++ )
|
389 |
|
|
if ( arr[k] > arr[k+1] )
|
390 |
|
|
{
|
391 |
|
|
swap ( arr[k], arr[k+1] );
|
392 |
|
|
isChanged = 1;
|
393 |
|
|
int bkp=order[k];
|
394 |
|
|
order[k]=order[k+1];
|
395 |
|
|
order[k+1]=bkp;
|
396 |
|
|
}
|
397 |
|
|
last--;
|
398 |
|
|
}
|
399 |
|
|
}
|
400 |
|
|
|
401 |
|
|
void swapvec(vector<float> &vec,int j, int k) {
|
402 |
|
|
float bkp=vec[j];
|
403 |
|
|
vec[j]=vec[k];
|
404 |
|
|
vec[k]=bkp;
|
405 |
|
|
}
|
406 |
|
|
|
407 |
|
|
void bubbleSort ( vector<float> &arr , vector<int> &order) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
|
408 |
|
|
{
|
409 |
|
|
int last = arr.size() - 2;
|
410 |
|
|
int isChanged = 1;
|
411 |
|
|
|
412 |
|
|
while ( last >= 0 && isChanged )
|
413 |
|
|
{
|
414 |
|
|
isChanged = 0;
|
415 |
|
|
for ( int k = 0; k <= last; k++ )
|
416 |
|
|
if ( arr[k] > arr[k+1] )
|
417 |
|
|
{
|
418 |
|
|
swapvec (arr,k,k+1);
|
419 |
|
|
isChanged = 1;
|
420 |
|
|
int bkp=order[k];
|
421 |
|
|
order[k]=order[k+1];
|
422 |
|
|
order[k+1]=bkp;
|
423 |
|
|
}
|
424 |
|
|
last--;
|
425 |
|
|
}
|
426 |
|
|
}
|
427 |
|
|
|
428 |
|
|
int numerichistoname=0;
|
429 |
|
|
string GetNumericHistoName() {
|
430 |
|
|
stringstream b;
|
431 |
|
|
b << "h_" << numerichistoname;
|
432 |
|
|
numerichistoname++;
|
433 |
|
|
return b.str();
|
434 |
|
|
}
|
435 |
|
|
|
436 |
|
|
//********************** BELOW : CUT INTERPRETATION **************************//
|
437 |
|
|
void splitupcut(string incut, vector<string> &partvector)
|
438 |
|
|
{
|
439 |
|
|
//idea: go thru the string called incut; if a parantheses is opened, then the cut cannot be split up until the parantheses is closed.
|
440 |
|
|
//ok anyway screw the parantheses.
|
441 |
|
|
int paranthesis_open=0;
|
442 |
|
|
int substr_start=0;
|
443 |
|
|
string currchar="";
|
444 |
|
|
for (int ichar=0;ichar<incut.length();ichar++)
|
445 |
|
|
{
|
446 |
|
|
currchar=incut.substr(ichar,1);
|
447 |
|
|
// if(currchar=="(") paranthesis_open++;
|
448 |
|
|
// if(currchar==")") paranthesis_open--;
|
449 |
|
|
if(currchar=="&"&&incut.substr(ichar+1,1)=="&"&¶nthesis_open==0) {
|
450 |
|
|
partvector.push_back(incut.substr(substr_start,ichar-substr_start));
|
451 |
|
|
substr_start=ichar+2;
|
452 |
|
|
}
|
453 |
|
|
}
|
454 |
|
|
partvector.push_back(incut.substr(substr_start,incut.length()-substr_start));
|
455 |
|
|
if(Verbosity>1) {
|
456 |
|
|
cout << "[ splitupcut() ] : The cut vector now contains the following elements: "<< endl;
|
457 |
|
|
for (int ipart=0;ipart<partvector.size();ipart++)
|
458 |
|
|
{
|
459 |
|
|
cout << " - " << partvector[ipart] << endl;
|
460 |
|
|
}
|
461 |
|
|
}
|
462 |
|
|
}
|
463 |
|
|
|
464 |
|
|
int atleastvalue(string expression, int &morethanlessthan) // takes in an expression such as ">2" or ">=3" and returns e.g. 3 (in both examples)
|
465 |
|
|
{
|
466 |
|
|
int retval=0;
|
467 |
|
|
if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
|
468 |
|
|
// cout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
|
469 |
|
|
morethanlessthan=1;
|
470 |
|
|
return atoi(expression.substr(2,1).c_str());
|
471 |
|
|
}
|
472 |
|
|
if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
|
473 |
|
|
// cout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
|
474 |
|
|
morethanlessthan=0;
|
475 |
|
|
return atoi(expression.substr(1,1).c_str());
|
476 |
|
|
}
|
477 |
|
|
if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
|
478 |
|
|
// cout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
|
479 |
|
|
morethanlessthan=-1;
|
480 |
|
|
return 1+atoi(expression.substr(1,1).c_str());
|
481 |
|
|
}
|
482 |
|
|
if(expression.substr(0,1)==">") {
|
483 |
|
|
// cout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
|
484 |
|
|
morethanlessthan=1;
|
485 |
|
|
return 1+atoi(expression.substr(2,1).c_str());
|
486 |
|
|
}
|
487 |
|
|
if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
|
488 |
|
|
// cout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
|
489 |
|
|
morethanlessthan=-1;
|
490 |
|
|
return 1+atoi(expression.substr(2,1).c_str());
|
491 |
|
|
}
|
492 |
|
|
}
|
493 |
|
|
|
494 |
|
|
int do_jet_cut(string incut, int *nJets) {
|
495 |
|
|
string expression=(incut.substr(12,incut.size()-12));
|
496 |
|
|
cout << "Going to analyze the jet cut : " << expression << " with 0,1 being " << expression.substr(0,1) << " and 1,1 being " << expression.substr(1,1) << endl;
|
497 |
|
|
if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
|
498 |
|
|
int nJet=atoi(expression.substr(2,1).c_str());
|
499 |
|
|
for(int i=nJet+1;i<20;i++) nJets[i]=0;
|
500 |
|
|
cout << "Is of type <=" << endl;
|
501 |
|
|
return 0;
|
502 |
|
|
}
|
503 |
|
|
if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
|
504 |
|
|
int nJet=atoi(expression.substr(2,1).c_str());
|
505 |
|
|
for(int i=0;i<20&&i!=nJet;i++) nJets[i]=0;
|
506 |
|
|
cout << "Is of type ==" << endl;
|
507 |
|
|
return 0;
|
508 |
|
|
}
|
509 |
|
|
if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
|
510 |
|
|
int nJet=atoi(expression.substr(2,1).c_str());
|
511 |
|
|
for(int i=0;i<nJet&&i!=nJet;i++) nJets[i]=0;
|
512 |
|
|
cout << "Is of type >=" << endl;
|
513 |
|
|
return 0;
|
514 |
|
|
}
|
515 |
|
|
if(expression.substr(0,1)=="<") {
|
516 |
|
|
int nJet=atoi(expression.substr(1,1).c_str());
|
517 |
|
|
for(int i=nJet;i<20;i++) nJets[i]=0;
|
518 |
|
|
cout << "Is of type <" << endl;
|
519 |
|
|
return 0;
|
520 |
|
|
}
|
521 |
|
|
if(expression.substr(0,1)==">") {
|
522 |
|
|
int nJet=atoi(expression.substr(1,1).c_str());
|
523 |
|
|
for(int i=0;i<nJet+1&&i!=nJet;i++) nJets[i]=0;
|
524 |
|
|
cout << "Is of type >" << endl;
|
525 |
|
|
return 0;
|
526 |
|
|
}
|
527 |
|
|
}
|
528 |
|
|
|
529 |
|
|
string interpret_cut(string incut, bool &isJetCut, int *permittednJets)
|
530 |
|
|
{
|
531 |
|
|
// isJetCut=false;nJets=-1;
|
532 |
|
|
if(incut=="()") return "";
|
533 |
|
|
while(incut.substr(0,1)=="(") incut=incut.substr(1,incut.length()-1);
|
534 |
|
|
while(incut.length()>0&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(0,incut.length()-1);
|
535 |
|
|
// if(incut.substr(0,1)=="("&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(1,incut.length()-2); //this is to make (cut) to cut.
|
536 |
|
|
|
537 |
|
|
if(Verbosity>0) {
|
538 |
|
|
cout << "Now interpreting cut " << incut << endl;
|
539 |
|
|
}
|
540 |
|
|
/*
|
541 |
|
|
if(incut=="ch1*ch2<0") return "OS";
|
542 |
|
|
if(incut=="id1==id2") return "SF";
|
543 |
|
|
if(incut=="id1!=id2") return "OF";
|
544 |
|
|
*/
|
545 |
|
|
if(incut=="ch1*ch2<0") return "";
|
546 |
|
|
if(incut=="id1==id2") return "";
|
547 |
|
|
if(incut=="id1!=id2") return "";
|
548 |
|
|
if(incut=="mll>2") return "";
|
549 |
|
|
|
550 |
|
|
if(incut=="mll>0") return ""; // my typical "fake cut"
|
551 |
|
|
|
552 |
|
|
if(incut=="passed_triggers||!is_data") return "Triggers";
|
553 |
|
|
if(incut=="pfjzb[0]>-998") return "";
|
554 |
|
|
|
555 |
|
|
|
556 |
|
|
if(incut=="id1==0") return "ee";
|
557 |
|
|
if(incut=="id1==1") return "#mu#mu";
|
558 |
|
|
if(incut=="abs(mll-91.2)<20") return "|m_{l^{+}l^{-}}-m_{Z}|<20";
|
559 |
|
|
if(incut=="pfJetGoodID[0]") return "";
|
560 |
|
|
if(incut=="pfJetGoodID[1]") return "";
|
561 |
|
|
if((int)incut.find("pfJetGoodNum")>-1) {
|
562 |
|
|
//do_jet_cut(incut,permittednJets);
|
563 |
|
|
stringstream result;
|
564 |
|
|
result << "nJets" << incut.substr(12,incut.size()-12);
|
565 |
|
|
/* cout << "Dealing with a jet cut: " << incut << endl;
|
566 |
|
|
stringstream result;
|
567 |
|
|
result << "nJets" << incut.substr(12,incut.size()-12);
|
568 |
|
|
isJetCut=true;
|
569 |
|
|
if(exactjetcut(incut,nJets))
|
570 |
|
|
// nJets=atleastvalue((incut.substr(12,incut.size()-12)),morethanlessthan);
|
571 |
|
|
return result.str();*/
|
572 |
|
|
return result.str();
|
573 |
|
|
}
|
574 |
|
|
return incut;
|
575 |
|
|
}
|
576 |
|
|
|
577 |
|
|
string interpret_nJet_range(int *nJets) {
|
578 |
|
|
for (int i=0;i<20;i++) cout << i << " : " << nJets[i] << endl;
|
579 |
|
|
return "hello";
|
580 |
|
|
}
|
581 |
|
|
|
582 |
|
|
string interpret_cuts(vector<string> &cutparts)
|
583 |
|
|
{
|
584 |
|
|
stringstream nicecut;
|
585 |
|
|
int nJets;
|
586 |
|
|
bool isJetCut;
|
587 |
|
|
int finalJetCut=-1;
|
588 |
|
|
int permittednJets[20];
|
589 |
|
|
for(int ijet=0;ijet<20;ijet++) permittednJets[ijet]=1;
|
590 |
|
|
int morethanlessthan=0;//-1: less than, 0: exactly, 1: more than
|
591 |
|
|
for(int icut=0;icut<cutparts.size();icut++)
|
592 |
|
|
{
|
593 |
|
|
if(icut==0) nicecut<<interpret_cut(cutparts[icut],isJetCut,permittednJets);
|
594 |
|
|
else {
|
595 |
|
|
string nice_this_cut = interpret_cut(cutparts[icut],isJetCut,permittednJets);//blublu
|
596 |
|
|
if(nice_this_cut.length()>0&&nicecut.str().length()>0) {
|
597 |
|
|
if(!isJetCut) nicecut<<"&&"<<nice_this_cut;
|
598 |
|
|
else {
|
599 |
|
|
if(nJets>finalJetCut) finalJetCut=nJets;
|
600 |
|
|
}
|
601 |
|
|
}
|
602 |
|
|
if(nice_this_cut.length()>0&&nicecut.str().length()==0) {
|
603 |
|
|
if(!isJetCut) {
|
604 |
|
|
nicecut<<nice_this_cut;
|
605 |
|
|
}
|
606 |
|
|
else {
|
607 |
|
|
if(nJets>finalJetCut) finalJetCut=nJets;
|
608 |
|
|
}
|
609 |
|
|
}
|
610 |
|
|
}
|
611 |
|
|
}
|
612 |
|
|
if(finalJetCut>-1) {
|
613 |
|
|
if(nicecut.str().length()==0) {
|
614 |
|
|
nicecut << "nJets#geq" << finalJetCut;
|
615 |
|
|
}
|
616 |
|
|
else
|
617 |
|
|
{
|
618 |
|
|
nicecut << "&&nJets#geq " << finalJetCut;
|
619 |
|
|
}
|
620 |
|
|
}
|
621 |
|
|
|
622 |
|
|
// cout << "The nJet allowed range is given by: " << interpret_nJet_range(permittednJets) << endl;
|
623 |
|
|
|
624 |
|
|
return nicecut.str();
|
625 |
|
|
}
|
626 |
|
|
|
627 |
|
|
string decipher_cut(TCut originalcut,TCut ignorethispart)
|
628 |
|
|
{
|
629 |
|
|
string incut=(const char*)originalcut;
|
630 |
|
|
string ignore=(const char*)ignorethispart;
|
631 |
|
|
|
632 |
|
|
if(ignore.length()>0 && incut.find(ignore)!=string::npos) incut=incut.replace(incut.find(ignore),ignore.length(),"");
|
633 |
|
|
|
634 |
|
|
vector<string>cutparts;
|
635 |
|
|
splitupcut(incut,cutparts);
|
636 |
|
|
string write_cut=interpret_cuts(cutparts);
|
637 |
|
|
return write_cut;
|
638 |
|
|
}
|
639 |
|
|
|
640 |
|
|
//********************** ABOVE : CUT INTERPRETATION **************************//
|
641 |
buchmann |
1.2 |
|
642 |
|
|
Double_t GausRandom(Double_t mu, Double_t sigma) {
|
643 |
|
|
return gRandom->Gaus(mu,sigma);// real deal
|
644 |
|
|
//return mu;//debugging : no smearing.
|
645 |
|
|
}
|
646 |
|
|
|
647 |
buchmann |
1.3 |
int functionalhistocounter=0;
|
648 |
buchmann |
1.2 |
TH1F * makehistofromfunction(TF1 *f1,TH1F *model) {
|
649 |
|
|
TH1F *histo = (TH1F*)model->Clone();
|
650 |
buchmann |
1.3 |
functionalhistocounter++;
|
651 |
|
|
stringstream histoname;
|
652 |
|
|
histoname << "histo_based_on_function_" << f1->GetName() << "__"<<functionalhistocounter;
|
653 |
|
|
histo->SetTitle(histoname.str().c_str());
|
654 |
|
|
histo->SetName(histoname.str().c_str());
|
655 |
buchmann |
1.2 |
int nbins=histo->GetNbinsX();
|
656 |
|
|
float low=histo->GetBinLowEdge(1);
|
657 |
|
|
float hi=histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX());
|
658 |
|
|
|
659 |
|
|
for(int i=0;i<=nbins;i++) {
|
660 |
buchmann |
1.3 |
histo->SetBinContent(i,f1->Eval(histo->GetBinCenter(i)));
|
661 |
|
|
histo->SetBinError(i,TMath::Sqrt(f1->Eval(histo->GetBinCenter(i))));
|
662 |
buchmann |
1.2 |
}
|
663 |
|
|
|
664 |
|
|
return histo;
|
665 |
|
|
} |