ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/SampleClass.C
Revision: 1.20
Committed: Thu Sep 29 07:37:27 2011 UTC (13 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: Honeypot, cbaf_2p1ifb
Changes since 1.19: +12 -1 lines
Log Message:
Added method to remove sample from sample collection; this can be very useful for scans to reduce memory consumption

File Contents

# Content
1 #include <iostream>
2 #include <sstream>
3 #include <vector>
4 #include <utility>
5
6 #include <TFile.h>
7 #include <TTree.h>
8 #include <TH1.h>
9 #include <TH2.h>
10 #include <TCut.h>
11 #include <THStack.h>
12 #include <TColor.h>
13 #include <TCanvas.h>
14 #include <TError.h>
15 #include <TText.h>
16 #include <TLegend.h>
17 #include <TError.h>
18 #include <TTreeFormula.h>
19
20 #define SampleClassLoaded
21
22 #ifndef Verbosity
23 #define Verbosity 0
24 #endif
25 #ifndef HUSH
26 #define HUSH 1
27 #endif
28 #ifndef GeneralToolBoxLoaded
29 #include "GeneralToolBox.C"
30 #endif
31 using namespace std;
32
33 bool doesROOTFileExist(string filename);
34 bool addoverunderflowbins=false;
35 string treename="PFevents";
36
37 class sample
38 {
39 public:
40 string filename;
41 string samplename;
42 long Nentries;
43 float xs;//cross section
44 float weight;
45 bool is_data;
46 bool is_signal;
47 bool is_active;
48 int groupindex;
49 Color_t samplecolor;
50
51 TFile *tfile;
52 TTree *events;
53
54 sample(string m_filename,string m_samplename,long m_Nentries, float m_xs,bool is_data, bool m_is_signal, int m_groupindex, Color_t color);
55 void closeFile();
56 };
57
58 string write_mc_or_data(bool is_data)
59 {
60 if(is_data) return "data";
61 return "MC";
62 }
63
64 sample::sample(string m_filename, string m_samplename, long m_Nentries, float m_xs,bool m_is_data, bool m_is_signal, int m_groupindex, Color_t mycolor)
65 {
66 this->filename=m_filename;
67 this->samplename=m_samplename;
68 this->Nentries=m_Nentries;
69 this->xs=m_xs;
70 this->is_data=m_is_data;
71 this->is_signal=m_is_signal;
72 this->weight=(xs/(float)Nentries);
73 this->groupindex=m_groupindex;
74 this->is_active=true;
75 this->samplecolor=mycolor;
76 if(!doesROOTFileExist(this->filename)) {
77 stringstream message;
78 message << "The " << write_mc_or_data(is_data) << " sample " << this->samplename << " is invalid because the associated file path, " << this->filename << " is incorrect.";
79 write_error("Sample::Sample",message.str());
80 this->is_active=false;
81 }
82
83 if(doesROOTFileExist(this->filename)) {
84 //suppressing stupid 64/32 errors here (Warning in <TFile::ReadStreamerInfo>: /scratch/buchmann/MC_Spring11_PU_PF/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root: not a TStreamerInfo object)
85 Int_t currlevel=gErrorIgnoreLevel;
86 gErrorIgnoreLevel=5000;
87 this->tfile = new TFile(m_filename.c_str());
88 gErrorIgnoreLevel=currlevel;
89 this->events=(TTree*)(this->tfile)->Get(treename.c_str());
90 if(Verbosity>0) dout << "The " << write_mc_or_data(is_data) << " file " << this->filename << " has been added successfully to the list of samples. " << endl;
91 }
92 else {
93 this->is_active=false;
94 }
95
96 }
97
98 void sample::closeFile()
99 {
100 if(doesROOTFileExist(this->filename)) {
101 (this->tfile)->Close();
102 }
103 else {
104 dout << "SAMPLE " << this->samplename << " cannot be closed as the underlying file (" << this->filename << ") does not exist!" << endl;
105 this->is_active=false;
106 }
107 }
108
109 bool doesROOTFileExist(string filename)
110 {
111 //suppressing stupid 64/32 errors here (Warning in <TFile::ReadStreamerInfo>: /scratch/buchmann/MC_Spring11_PU_PF/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root: not a TStreamerInfo object)
112 Int_t currlevel=gErrorIgnoreLevel;
113 gErrorIgnoreLevel=5000;
114 TFile *f = new TFile(filename.c_str());
115
116 if (f->IsZombie()) {
117 dout << "Error opening file" << filename << endl;
118 return 0;
119 }
120 f->Close();
121 gErrorIgnoreLevel=currlevel;
122 return 1;
123 }
124 //********************************************************
125
126 TCut essentialcut("mll>0");
127 // This will reweight all the events: use "weight" to correct MC for pileUP, "1.0" otherwise
128 TCut cutWeight("weight");
129
130 void setessentialcut(TCut ess) {
131 essentialcut=ess;
132 }
133
134 class samplecollection
135 {
136 public:
137 string name;
138 vector<sample> collection;
139 int nsamples;
140 int ndatasamples;
141 int nmcsamples;
142 int ngroups;
143
144 samplecollection(string m_name);
145 void AddSample(string m_filename,string m_samplename,long m_Nentries, float m_xs,bool is_data, bool m_is_signal, int groupindex, Color_t color);
146
147 //vector quantities
148 TH1F* Draw(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal);
149 TH1F* Draw(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal);
150
151 //minx,maxs quantities
152 TH1F* Draw(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal);
153 TH1F* Draw(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal);
154
155 TH2F* Draw(string m_histoname,string m_var, vector<float> binningx, vector<float> binningy, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal);
156
157 vector<float> get_optimal_binsize(string variable, TCut cut,int nbins, float low, float high);
158
159 THStack DrawStack(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal);
160 THStack DrawStack(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal);
161 vector<int> FindSample(string what);
162 void ListSamples();
163 bool do_sample(int thissample, vector<int> &selected_samples);
164 Color_t GetColor(string filename);
165 Color_t GetColor(int sampleindex);
166
167 TLegend* allbglegend(string title,float x, float y);
168 TLegend* allbglegend(string title, TH1F *data, float x, float y);
169
170 void PickUpFromThisFile(int isamp, string cut, vector<string> &output, vector<string> &pickupfile);
171 void PickUpEvents(string cut);
172 string find_units(string&);
173
174 void RemoveLastSample();
175 };
176
177 samplecollection::samplecollection(string m_name)
178 {
179 this->name=m_name;
180 if(Verbosity>0) dout << "Initiated sample collection " << this->name << " " << endl;
181 }
182
183 void samplecollection::ListSamples()
184 {
185 dout << "---------------------------------------------------------------------------------------------------" << endl;
186 dout << "Listing all " << this->nsamples << " sample(s) of the sample collection called " << this->name << " : " << endl;
187 if(this->ndatasamples>0) {
188 dout << "Data sample(s): " << endl;
189 for (int isamp=0;isamp<this->collection.size();isamp++)
190 {
191 if((this->collection)[isamp].is_data) dout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
192 }
193 }
194 if(this->nmcsamples>0) {
195 dout << "MC sample(s): " << endl;
196 for (int isamp=0;isamp<this->collection.size();isamp++)
197 {
198 if(!(this->collection)[isamp].is_data) dout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
199 }
200 }
201 dout << "---------------------------------------------------------------------------------------------------" << endl;
202 }
203
204 void samplecollection::AddSample(string m_filename,string m_samplename,long m_Nentries, float m_xs,bool is_data, bool m_is_signal, int groupindex, Color_t newcolor)
205 {
206 sample NewSample(m_filename,m_samplename,m_Nentries,m_xs,is_data,m_is_signal,groupindex,newcolor);
207 (this->collection).push_back(NewSample);
208 if((this->collection).size()==1) {
209 this->nmcsamples=0;
210 this->ndatasamples=0;
211 this->nsamples=0;
212 this->ngroups=0;
213 }
214 if(groupindex>this->ngroups) this->ngroups=groupindex;
215 if(is_data) this->ndatasamples++;
216 else this->nmcsamples++;
217 this->nsamples=this->collection.size();
218 }
219
220 bool samplecollection::do_sample(int thissample, vector<int> &selected_samples)
221 {
222 bool drawit=false;
223 for(int isel=0;isel<selected_samples.size();isel++)
224 {
225 if(selected_samples[isel]==thissample) {
226 drawit=true;
227 break;
228 }
229 }
230 return drawit;
231 }
232
233 TH1F* samplecollection::Draw(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal=false)
234 {
235 vector<int> emptyvector;
236 TH1F *histo = this->Draw(m_histoname,m_var, m_nbins, m_minx, m_maxx, m_xlabel, m_ylabel, Cut, m_is_data, luminosity,emptyvector,drawsignal);
237 return histo;
238 }
239
240 TH1F* samplecollection::Draw(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal=false)
241 {
242 vector<int> emptyvector;
243 TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity,emptyvector,drawsignal);
244 return histo;
245 }
246
247 TH1F* samplecollection::Draw(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal=false) {
248 if(Verbosity>0) dout << endl << endl;
249 if(Verbosity>0) dout << "-------------------------------------------------------------------------------------" << endl;
250 if(Verbosity>0) dout << "histoname : " << m_histoname << " , m_var = " << m_var << ", m_xlabel=" << m_xlabel << " m_ylabel=" << m_ylabel << " m_is_data = " << m_is_data << " lumi = " << luminosity << " onlyindex size: " << onlyindex.size() << endl;
251 if(HUSH==0) dout << "Drawing histo called " << m_histoname << "... " << endl;
252 bool do_only_selected_samples=false;
253 if(onlyindex.size()>0&&onlyindex[0]!=-1) {
254 if(Verbosity>0) {dout << "Requested to only draw sample corresponding to the following sample(s) : " << endl;}
255 for(int is=0;is<onlyindex.size();is++) {
256 if(Verbosity>0) dout << " - " << (this->collection)[onlyindex[is]].filename << " (sample index: " << onlyindex[is] << ")" << endl;
257 }
258 do_only_selected_samples=true;
259 }
260 if(onlyindex.size()==1&&onlyindex[0]==-1) {
261 do_only_selected_samples=true; // this is the special case when we request a non-existing sample - we shouldn't draw anything then.
262 onlyindex.clear();
263 }
264 stringstream h_histoname;
265 h_histoname<<"h_"<<m_histoname;
266 float binningarray[binning.size()+1];
267 for(int i=0;i<binning.size();i++) {binningarray[i]=binning[i];}
268 TH1F *histo = new TH1F(m_histoname.c_str(),"",binning.size()-1,binningarray);
269 histo->Sumw2();
270
271 stringstream drawthis;
272 drawthis<<m_var<<">>tempdrawhisto";
273
274 for (unsigned int isample=0;isample<(this->collection).size();isample++) {
275 bool use_this_sample=false;
276 if(!(this->collection)[isample].is_active) continue; // skip inactive samples right away
277 if(((this->collection)[isample].is_data==m_is_data)&&(!do_only_selected_samples)) {
278 if(drawsignal==false&&(this->collection)[isample].is_signal==true) continue;
279 use_this_sample=true;
280 }
281 if(do_only_selected_samples&&this->do_sample(isample,onlyindex)) {
282 use_this_sample=true;
283 }
284 TH1F *tempdrawhisto = new TH1F("tempdrawhisto","tempdrawhisto",binning.size()-1,binningarray);
285 tempdrawhisto->Sumw2();
286 if(use_this_sample) {
287 if(Verbosity>0) dout << "[samplecollection::Draw] : Added contribution from sample " << (this->collection)[isample].samplename << endl;
288 (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*cutWeight);//this weight is based on PU etc. not XS
289 if(addoverunderflowbins) {
290 //now also adding the overflow & underflow bins:
291 tempdrawhisto->SetBinContent(tempdrawhisto->GetNbinsX(),tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()+1)+tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()));
292 tempdrawhisto->SetBinError(tempdrawhisto->GetNbinsX(),TMath::Sqrt(tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX())));
293 tempdrawhisto->SetBinContent(1,tempdrawhisto->GetBinContent(0)+tempdrawhisto->GetBinContent(1));
294 tempdrawhisto->SetBinError(1,TMath::Sqrt(tempdrawhisto->GetBinContent(1)));
295 }
296
297 if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));//weight applied here is XS & N(entries)
298 histo->Add(tempdrawhisto);
299 }
300 tempdrawhisto->Delete();
301 }//end of loop over isample
302 if(Verbosity>0) dout << "Histo has been filled and now contains " << histo->Integral() << " points (integral)" << endl;
303 histo->GetXaxis()->SetTitle(m_xlabel.c_str());
304 // Try to add bin width information: look for units in m_xlabel
305 string units = find_units(m_xlabel);
306 if ( units.length()>0 ) {
307 stringstream ylabel;
308 ylabel << m_ylabel << " / " << histo->GetBinWidth(1) << " " << units;
309 histo->GetYaxis()->SetTitle( ylabel.str().c_str() );
310 } else {
311 histo->GetYaxis()->SetTitle(m_ylabel.c_str());
312 }
313 histo->GetXaxis()->CenterTitle();
314 histo->GetYaxis()->CenterTitle();
315 if(do_only_selected_samples) histo->SetLineColor(this->GetColor(onlyindex[0]));
316 return histo;
317 }
318
319 TH2F* samplecollection::Draw(string m_histoname,string m_var, vector<float> binningx, vector<float> binningy, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal=false) {
320 if(Verbosity>0) dout << endl << endl;
321 if(Verbosity>0) dout << "-------------------------------------------------------------------------------------" << endl;
322 if(Verbosity>0) dout << "histoname : " << m_histoname << " , m_var = " << m_var << ", m_xlabel=" << m_xlabel << " m_ylabel=" << m_ylabel << " m_is_data = " << m_is_data << " lumi = " << luminosity << " onlyindex size: " << onlyindex.size() << endl;
323 if(HUSH==0) dout << "Drawing histo called " << m_histoname << "... " << endl;
324 bool do_only_selected_samples=false;
325 if(onlyindex.size()>0&&onlyindex[0]!=-1) {
326 if(Verbosity>0) {dout << "Requested to only draw sample corresponding to the following sample(s) : " << endl;}
327 for(int is=0;is<onlyindex.size();is++) {
328 if(Verbosity>0) dout << " - " << (this->collection)[onlyindex[is]].filename << " (sample index: " << onlyindex[is] << ")" << endl;
329 }
330 do_only_selected_samples=true;
331 }
332 if(onlyindex.size()==1&&onlyindex[0]==-1) {
333 do_only_selected_samples=true; // this is the special case when we request a non-existing sample - we shouldn't draw anything then.
334 onlyindex.clear();
335 }
336 stringstream h_histoname;
337 h_histoname<<"h_"<<m_histoname;
338 float binningarray[binningx.size()+1];
339 float binningyarray[binningy.size()+1];
340 for(int i=0;i<binningx.size();i++) {binningarray[i]=binningx[i];binningyarray[i]=binningy[i];}
341 TH2F *histo = new TH2F(m_histoname.c_str(),"",binningx.size()-1,binningarray,binningy.size(),binningyarray);
342 histo->Sumw2();
343
344 stringstream drawthis;
345 drawthis<<m_var<<">>tempdrawhisto";
346
347 for (unsigned int isample=0;isample<(this->collection).size();isample++) {
348 bool use_this_sample=false;
349 if(!(this->collection)[isample].is_active) continue; // skip inactive samples right away
350 if(((this->collection)[isample].is_data==m_is_data)&&(!do_only_selected_samples)) {
351 if(drawsignal==false&&(this->collection)[isample].is_signal==true) continue;
352 use_this_sample=true;
353 }
354 if(do_only_selected_samples&&this->do_sample(isample,onlyindex)) {
355 use_this_sample=true;
356 }
357 TH2F *tempdrawhisto = new TH2F("tempdrawhisto","tempdrawhisto",binningx.size()-1,binningarray,binningy.size()-1,binningyarray);
358 tempdrawhisto->Sumw2();
359 if(use_this_sample) {
360 if(Verbosity>0) dout << "[samplecollection::Draw] : Added contribution from sample " << (this->collection)[isample].samplename << endl;
361 (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*cutWeight);//this weight is based on PU etc. not XS
362 if(addoverunderflowbins) {
363 //now also adding the overflow & underflow bins:
364 tempdrawhisto->SetBinContent(tempdrawhisto->GetNbinsX(),tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()+1)+tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()));
365 tempdrawhisto->SetBinError(tempdrawhisto->GetNbinsX(),TMath::Sqrt(tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX())));
366 tempdrawhisto->SetBinContent(1,tempdrawhisto->GetBinContent(0)+tempdrawhisto->GetBinContent(1));
367 tempdrawhisto->SetBinError(1,TMath::Sqrt(tempdrawhisto->GetBinContent(1)));
368 }
369
370 if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));//weight applied here is XS & N(entries)
371 histo->Add(tempdrawhisto);
372 }
373 tempdrawhisto->Delete();
374 }//end of loop over isample
375 if(Verbosity>0) dout << "Histo has been filled and now contains " << histo->Integral() << " points (integral)" << endl;
376 histo->GetXaxis()->SetTitle(m_xlabel.c_str());
377 histo->GetYaxis()->SetTitle(m_ylabel.c_str());
378 histo->GetXaxis()->CenterTitle();
379 histo->GetYaxis()->CenterTitle();
380 if(do_only_selected_samples) histo->SetLineColor(this->GetColor(onlyindex[0]));
381 return histo;
382 }
383
384
385 TH1F* samplecollection::Draw(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, vector<int> onlyindex, bool drawsignal=false)
386 {
387 vector<float> binning;
388 for(int i=0;i<=m_nbins;i++)
389 {
390 binning.push_back(((float)(m_maxx-m_minx)/((float)m_nbins))*i+m_minx);
391 }
392
393 TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity, onlyindex,drawsignal);
394 return histo;
395 }
396
397 THStack samplecollection::DrawStack(string m_histoname,string m_var, vector<float> binning, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool drawsignal=false) {
398 stringstream h_histoname;
399 h_histoname<<"h_"<<m_histoname;
400 THStack thestack("thestack",m_histoname.c_str());
401
402 stringstream drawthis;
403 drawthis<<m_var<<">>tempdrawhisto";
404
405 TH1F *histogroups[this->ngroups+1];
406 int bookedhistos[this->ngroups+1];
407 for(int ih=0;ih<=this->ngroups;ih++) bookedhistos[ih]=0;
408
409 float binningarray[binning.size()+1];
410 for(int i=0;i<binning.size();i++) {binningarray[i]=binning[i];}
411
412 for (unsigned int isample=0;isample<(this->collection).size();isample++) {
413 if(!drawsignal&&(this->collection)[isample].is_signal) continue;
414 if((this->collection)[isample].is_active&&(this->collection)[isample].is_data==m_is_data) {//fills mc if we want mc, else fills data.
415 TH1F *tempdrawhisto = new TH1F("tempdrawhisto","",binning.size()-1,binningarray);
416 (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*cutWeight);
417
418 if(addoverunderflowbins) {
419 //now also adding the overflow & underflow bins:
420 tempdrawhisto->SetBinContent(tempdrawhisto->GetNbinsX(),tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()+1)+tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()));
421 // tempdrawhisto->SetBinError(tempdrawhisto->GetNbinsX(),TMath::Sqrt(tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()))); // errors not necessary for stack!
422 tempdrawhisto->SetBinContent(1,tempdrawhisto->GetBinContent(0)+tempdrawhisto->GetBinContent(1));
423 // tempdrawhisto->SetBinError(1,TMath::Sqrt(tempdrawhisto->GetBinContent(1))); // errors not necessary for stack!
424 }
425
426 if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));
427 tempdrawhisto->SetFillColor(this->GetColor(isample));
428 if(bookedhistos[(this->collection)[isample].groupindex]==0) {
429 tempdrawhisto->SetName(GetNumericHistoName().c_str());
430 histogroups[(this->collection)[isample].groupindex]=(TH1F*)tempdrawhisto->Clone();
431 bookedhistos[(this->collection)[isample].groupindex]=1;
432 }
433 else
434 {
435 histogroups[(this->collection)[isample].groupindex]->Add(tempdrawhisto);
436 bookedhistos[(this->collection)[isample].groupindex]++;
437 }
438 tempdrawhisto->Delete();
439 }
440 }
441 vector<int> ordered_indices;
442 vector<float> ordered_integrals;
443
444 for(int ig=0;ig<=this->ngroups;ig++) {
445 if(bookedhistos[ig]>0) {
446 ordered_indices.push_back(ig);
447 ordered_integrals.push_back(histogroups[ig]->Integral());
448 }
449 }
450
451 /*
452 bubbleSort(ordered_integrals,ordered_indices);
453 // for(int index=ordered_indices.size()-1;index>=0;index--) {
454 for(int index=0;index<ordered_indices.size();index++) {
455 thestack.Add(histogroups[ordered_indices[index]]);
456 }
457 */
458 for(int index=0;index<ordered_indices.size();index++) {
459 thestack.Add(histogroups[ordered_indices[index]]);
460 }
461 return thestack;
462 }
463
464 THStack samplecollection::DrawStack(string m_histoname,string m_var, int m_nbins, float m_minx, float m_maxx, string m_xlabel, string m_ylabel, TCut Cut, bool m_is_data, float luminosity, bool
465 drawsignal=false)
466 {
467 vector<float> binning;
468 for(int i=0;i<=m_nbins;i++)
469 {
470 binning.push_back(((float)(m_maxx-m_minx)/((float)m_nbins))*i+m_minx);
471 }
472 return this->DrawStack(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity, drawsignal);
473 }
474
475
476 vector<int> samplecollection::FindSample(string what)
477 {
478 vector<int> hitcollection;
479 for(unsigned int isam=0;isam<(this->collection).size();isam++)
480 {
481 if(((this->collection)[isam].filename).find(what)!=string::npos) {
482 hitcollection.push_back(isam);
483 } else {
484 }
485 }
486 if(hitcollection.size()==0) {
487 hitcollection.push_back(-1);
488 write_warning(__FUNCTION__,"Couldn't find sample "+string(what)+" using sample collection \""+string(this->name)+"\"");
489 }
490 return hitcollection;
491 }
492
493 Color_t samplecollection::GetColor(string filename) {
494 vector<int> corresponding_samples = this->FindSample(filename);
495 if(corresponding_samples.size()>0) return this->GetColor(corresponding_samples[0]);
496 return kBlack;
497 }
498
499 Color_t samplecollection::GetColor(int sampleindex) {
500 return this->collection[sampleindex].samplecolor;
501 }
502
503 TLegend* samplecollection::allbglegend(string title, TH1F *data,float posx=0.65, float posy=0.60) {
504 // TLegend *leg = new TLegend(0.65,0.60,0.89,0.77);
505 TLegend *leg = new TLegend(posx,posy,0.89,0.89);
506 if(title!="") leg->SetHeader(title.c_str());
507 if(data->GetName()!="nothing") leg->AddEntry(data,"Data","p");
508 leg->SetFillColor(kWhite);
509 leg->SetBorderSize(0);
510 leg->SetLineColor(kWhite);
511
512 TH1F *fakehistos[(this->collection).size()];
513 bool donealready[(this->collection).size()];
514 for(int i=0;i<(this->collection).size();i++) donealready[i]=false;
515 for(int isample=0;isample<(this->collection).size();isample++) {
516 if((this->collection)[isample].is_data||(this->collection)[isample].is_signal) continue;
517
518 if(!donealready[(this->collection)[isample].groupindex]) {
519 fakehistos[(this->collection)[isample].groupindex] = new TH1F(GetNumericHistoName().c_str(),"",1,0,1);
520 fakehistos[(this->collection)[isample].groupindex]->Fill(1);
521 fakehistos[(this->collection)[isample].groupindex]->SetFillColor(this->GetColor(isample));
522 leg->AddEntry(fakehistos[(this->collection)[isample].groupindex],((this->collection)[isample].samplename).c_str(),"f");
523 donealready[(this->collection)[isample].groupindex]=true;
524 }
525 }
526 DrawPrelim();
527 leg->SetTextFont(42);
528 return leg;
529 }
530
531 TLegend* samplecollection::allbglegend(string title="",float posx=0.65, float posy=0.60) {
532 Int_t currlevel=gErrorIgnoreLevel;
533 gErrorIgnoreLevel=5000;
534 TH1F *blub = new TH1F("nothing","nothing",1,0,1);
535 gErrorIgnoreLevel=currlevel;//we know this possibly replaces a previous histo, but we don't care since it's fake anyway.
536 return this->allbglegend(title,blub,posx,posy);
537 }
538
539
540 void set_treename(string treen="events") {
541 dout << "Treename has been set to " << treen << endl;
542 if(treen=="PFevents"||treen=="events") treename=treen;
543 else write_error("sample::set_treename","Setting the treename failed as you provided an invalid tree name.");
544 }
545
546 vector<float> samplecollection::get_optimal_binsize(string variable, TCut cut,int nbins, float low, float high) {
547 TH1F *histo = this->Draw("histo",variable,5000,low,high, "finding_optimal_binsize", "events", cut,0,1000);
548 float eventsperbin=histo->Integral()/nbins;
549 vector<float> binning;
550 binning.push_back(high);
551 float runningsum=0;
552 for(int i=histo->GetNbinsX();i>0;i--) {
553 runningsum+=histo->GetBinContent(i);
554 if(TMath::Abs(runningsum-eventsperbin)<0.05*eventsperbin||runningsum>eventsperbin) {
555 binning.push_back(histo->GetBinLowEdge(i));
556 runningsum=0;
557 }
558 }
559 if(runningsum>0) binning.push_back(low);
560 for(int i=0;i<int(binning.size()/2);i++) {
561 float temp=binning[i];
562 binning[i]=binning[binning.size()-1-i];
563 binning[binning.size()-1-i]=temp;
564 }
565
566 return binning;
567
568 }
569
570 void samplecollection::PickUpFromThisFile(int isamp, string cut, vector<string> &output, vector<string> &pickupfile) {
571 int lumi,eventNum,runNum;
572 float jzb[30];
573 (this->collection)[isamp].events->SetBranchAddress("jzb",&jzb);
574 (this->collection)[isamp].events->SetBranchAddress("lumi",&runNum);
575 (this->collection)[isamp].events->SetBranchAddress("eventNum",&eventNum);
576 (this->collection)[isamp].events->SetBranchAddress("lumi",&lumi);
577 (this->collection)[isamp].events->SetBranchAddress("runNum",&runNum);
578
579
580 TTreeFormula *select = new TTreeFormula("select", cut.c_str()&&essentialcut, (this->collection)[isamp].events);
581 int npickedup=0;
582 for (Int_t entry = 0 ; entry < (this->collection)[isamp].events->GetEntries() ; entry++) {
583 (this->collection)[isamp].events->LoadTree(entry);
584 if (select->EvalInstance()) {
585 (this->collection)[isamp].events->GetEntry(entry);
586 dout << runNum << ":" << lumi << ":" << eventNum << endl;
587 npickedup++;
588 }
589 }
590 dout << "Printed out a total of " << npickedup << " events" << endl;
591 }
592
593
594
595 void samplecollection::PickUpEvents(string cut) {
596 vector<string> output;
597 vector<string> pickupfile;
598 for (int isamp=0;isamp<this->collection.size();isamp++)
599 {
600 if((this->collection)[isamp].is_data) {
601 //we have a data sample !
602 this->PickUpFromThisFile(isamp,cut,output,pickupfile);
603 }
604 }
605 //do something with output and of course the pickup file!
606 }
607
608 //________________________________________________________________________________________
609 // Find units from histogram x-label (looks for '[...]')
610 string samplecollection::find_units(string& xlabel) {
611
612 string units;
613 size_t pos1 = xlabel.find("[");
614 if ( pos1 != string::npos ) {
615 size_t pos2 = xlabel.find("]");
616 units = xlabel.substr(pos1+1,pos2-pos1-1);
617 }// else {
618 // size_t pos1 = xlabel.find("(");
619 // if ( pos1 != string::npos ) {
620 // size_t pos2 = xlabel.find(")");
621 // units = xlabel.substr(pos1+1,pos2-pos1-1);
622 // }
623 // }
624 return units;
625
626 }
627
628 void samplecollection::RemoveLastSample() {
629 dout << "Removing last sample: " << ((this->collection)[(this->collection).size()-1]).filename << std::endl;
630 if(((this->collection)[(this->collection).size()-1]).is_data) this->ndatasamples--;
631 else this->nmcsamples--;
632 this->nsamples--;
633 ((this->collection)[(this->collection).size()-1]).closeFile();
634 (this->collection).pop_back();
635 }
636
637 void switch_overunderflow(bool newpos=false) {
638 addoverunderflowbins=newpos;
639 }
640
641