ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/SampleClass.C
Revision: 1.8
Committed: Wed May 23 16:01:42 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.7: +0 -1 lines
Log Message:
Avoiding redeclaration

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