ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/SampleClass.C
Revision: 1.5
Committed: Mon May 14 15:14:53 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.4: +1 -1 lines
Log Message:
Added names to histograms in stack so they are more easily identifiable

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