ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/SampleClass.C
Revision: 1.9
Committed: Tue Jun 5 09:18:29 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.8: +6 -3 lines
Log Message:
Writing out weight normalization factor as well

File Contents

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