ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/SampleClass.C
Revision: 1.1
Committed: Wed Jun 22 11:07:37 2011 UTC (13 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Log Message:
Initial commit of Plotting tools

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 <TCut.h>
10     #include <THStack.h>
11     #include <TColor.h>
12     #include <TCanvas.h>
13     #include <TError.h>
14     #include <TText.h>
15     #include <TLegend.h>
16     #include <TError.h>
17    
18     #define SampleClassLoaded
19     #ifndef Verbosity
20     #define Verbosity 0
21     #endif
22     #ifndef HUSH
23     #define HUSH 1
24     #endif
25     #ifndef GeneralToolBoxLoaded
26     #include "GeneralToolBox.C"
27     #endif
28     using namespace std;
29    
30     bool doesROOTFileExist(string filename);
31    
32     string treename="PFevents";
33    
34     class sample
35     {
36     public:
37     string filename;
38     string samplename;
39     int Nentries;
40     float xs;//cross section
41     float weight;
42     bool is_data;
43     bool is_signal;
44     bool is_active;
45     int groupindex;
46     Color_t samplecolor;
47    
48     TFile *tfile;
49     TTree *events;
50    
51     sample(string m_filename,string m_samplename,int m_Nentries, float m_xs,bool is_data, bool m_is_signal, int m_groupindex, Color_t color);
52     void closeFile();
53     };
54    
55     string write_mc_or_data(bool is_data)
56     {
57     if(is_data) return "data";
58     return "MC";
59     }
60    
61     sample::sample(string m_filename, string m_samplename, int m_Nentries, float m_xs,bool m_is_data, bool m_is_signal, int m_groupindex, Color_t mycolor)
62     {
63     this->filename=m_filename;
64     this->samplename=m_samplename;
65     this->Nentries=m_Nentries;
66     this->xs=m_xs;
67     this->is_data=m_is_data;
68     this->is_signal=m_is_signal;
69     this->weight=(xs/(float)Nentries);
70     this->groupindex=m_groupindex;
71     this->is_active=true;
72     this->samplecolor=mycolor;
73     if(!doesROOTFileExist(this->filename)) {
74     stringstream message;
75     message << "The " << write_mc_or_data(is_data) << " sample " << this->samplename << " is invalid because the associated file path, " << this->filename << " is incorrect.";
76     write_error("Sample::Sample",message.str());
77     this->is_active=false;
78     }
79    
80     if(doesROOTFileExist(this->filename)) {
81     //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)
82     Int_t currlevel=gErrorIgnoreLevel;
83     gErrorIgnoreLevel=5000;
84     this->tfile = new TFile(m_filename.c_str());
85     gErrorIgnoreLevel=currlevel;
86     this->events=(TTree*)(this->tfile)->Get(treename.c_str());
87     if(Verbosity>0) cout << "The " << write_mc_or_data(is_data) << " file " << this->filename << " has been added successfully to the list of samples. " << endl;
88     }
89     else {
90     this->is_active=false;
91     }
92    
93     }
94    
95     void sample::closeFile()
96     {
97     if(doesROOTFileExist(this->filename)) {
98     (this->tfile)->Close();
99     }
100     else {
101     cout << "SAMPLE " << this->samplename << " cannot be closed as the underlying file (" << this->filename << ") does not exist!" << endl;
102     this->is_active=false;
103     }
104     }
105    
106     bool doesROOTFileExist(string filename)
107     {
108     //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)
109     Int_t currlevel=gErrorIgnoreLevel;
110     gErrorIgnoreLevel=5000;
111     TFile *f = new TFile(filename.c_str());
112    
113     if (f->IsZombie()) {
114     cout << "Error opening file" << filename << endl;
115     return 0;
116     }
117     f->Close();
118     gErrorIgnoreLevel=currlevel;
119     return 1;
120     }
121     //********************************************************
122    
123     TCut essentialcut("mll>0");
124    
125     void setessentialcut(TCut ess) {
126     essentialcut=ess;
127     }
128    
129     class samplecollection
130     {
131     public:
132     string name;
133     vector<sample> collection;
134     int nsamples;
135     int ndatasamples;
136     int nmcsamples;
137     int ngroups;
138    
139     samplecollection(string m_name);
140     void AddSample(string m_filename,string m_samplename,int m_Nentries, float m_xs,bool is_data, bool m_is_signal, int groupindex, Color_t color);
141    
142     //vector quantities
143     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);
144     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);
145    
146     //minx,maxs quantities
147     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);
148     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);
149    
150     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);
151     vector<int> FindSample(string what);
152     void ListSamples();
153     bool do_sample(int thissample, vector<int> &selected_samples);
154     Color_t GetColor(string filename);
155     Color_t GetColor(int sampleindex);
156    
157     TLegend* allbglegend(string title);
158     TLegend* allbglegend(string title, TH1F *data);
159     };
160    
161     samplecollection::samplecollection(string m_name)
162     {
163     this->name=m_name;
164     if(Verbosity>0) cout << "Initiated sample collection " << this->name << " " << endl;
165     }
166    
167     void samplecollection::ListSamples()
168     {
169     cout << "---------------------------------------------------------------------------------------------------" << endl;
170     cout << "Listing all " << this->nsamples << " sample(s) of the sample collection called " << this->name << " : " << endl;
171     if(this->ndatasamples>0) {
172     cout << "Data sample(s): " << endl;
173     for (int isamp=0;isamp<this->collection.size();isamp++)
174     {
175     if((this->collection)[isamp].is_data) cout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
176     }
177     }
178     if(this->nmcsamples>0) {
179     cout << "MC sample(s): " << endl;
180     for (int isamp=0;isamp<this->collection.size();isamp++)
181     {
182     if(!(this->collection)[isamp].is_data) cout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
183     }
184     }
185     cout << "---------------------------------------------------------------------------------------------------" << endl;
186     }
187    
188     void samplecollection::AddSample(string m_filename,string m_samplename,int m_Nentries, float m_xs,bool is_data, bool m_is_signal, int groupindex, Color_t newcolor)
189     {
190     sample NewSample(m_filename,m_samplename,m_Nentries,m_xs,is_data,m_is_signal,groupindex,newcolor);
191     (this->collection).push_back(NewSample);
192     if((this->collection).size()==1) {
193     this->nmcsamples=0;
194     this->ndatasamples=0;
195     this->nsamples=0;
196     this->ngroups=0;
197     }
198     if(groupindex>this->ngroups) this->ngroups=groupindex;
199     if(is_data) this->ndatasamples++;
200     else this->nmcsamples++;
201     this->nsamples=this->collection.size();
202     }
203    
204     bool samplecollection::do_sample(int thissample, vector<int> &selected_samples)
205     {
206     bool drawit=false;
207     for(int isel=0;isel<selected_samples.size();isel++)
208     {
209     if(selected_samples[isel]==thissample) {
210     drawit=true;
211     break;
212     }
213     }
214     return drawit;
215     }
216    
217     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)
218     {
219     vector<int> emptyvector;
220     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);
221     return histo;
222     }
223    
224     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)
225     {
226     vector<int> emptyvector;
227     TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity,emptyvector,drawsignal);
228     return histo;
229     }
230    
231     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) {
232     if(Verbosity>0) cout << endl << endl;
233     if(Verbosity>0) cout << "-------------------------------------------------------------------------------------" << endl;
234     if(Verbosity>0) cout << "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;
235     if(HUSH==0) cout << "Drawing histo called " << m_histoname << "... " << endl;
236     bool do_only_selected_samples=false;
237     if(onlyindex.size()>0&&onlyindex[0]!=-1) {
238     if(Verbosity>0) {cout << "Requested to only draw sample corresponding to the following sample(s) : " << endl;}
239     for(int is=0;is<onlyindex.size();is++) {
240     if(Verbosity>0) cout << " - " << (this->collection)[onlyindex[is]].filename << " (sample index: " << onlyindex[is] << ")" << endl;
241     }
242     do_only_selected_samples=true;
243     }
244     if(onlyindex.size()==1&&onlyindex[0]==-1) {
245     do_only_selected_samples=true; // this is the special case when we request a non-existing sample - we shouldn't draw anything then.
246     onlyindex.clear();
247     }
248     stringstream h_histoname;
249     h_histoname<<"h_"<<m_histoname;
250     float binningarray[binning.size()+1];
251     for(int i=0;i<binning.size();i++) {binningarray[i]=binning[i];}
252     TH1F *histo = new TH1F(m_histoname.c_str(),"",binning.size()-1,binningarray);
253     histo->Sumw2();
254    
255     stringstream drawthis;
256     drawthis<<m_var<<">>tempdrawhisto";
257    
258     for (unsigned int isample=0;isample<(this->collection).size();isample++) {
259     bool use_this_sample=false;
260     if(!(this->collection)[isample].is_active) continue; // skip inactive samples right away
261     if(((this->collection)[isample].is_data==m_is_data)&&(!do_only_selected_samples)) {
262     if(drawsignal==false&&(this->collection)[isample].is_signal==true) continue;
263     use_this_sample=true;
264     }
265     if(do_only_selected_samples&&this->do_sample(isample,onlyindex)) {
266     use_this_sample=true;
267     }
268     TH1F *tempdrawhisto = new TH1F("tempdrawhisto","tempdrawhisto",binning.size()-1,binningarray);
269     tempdrawhisto->Sumw2();
270     if(use_this_sample) {
271     if(Verbosity>0) cout << "[samplecollection::Draw] : Added contribution from sample " << (this->collection)[isample].samplename << endl;
272     (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*"weight");
273     if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));
274     histo->Add(tempdrawhisto);
275     }
276     tempdrawhisto->Delete();
277     }//end of loop over isample
278     if(Verbosity>0) cout << "Histo has been filled and now contains " << histo->Integral() << " points (integral)" << endl;
279     histo->GetXaxis()->SetTitle(m_xlabel.c_str());
280     histo->GetYaxis()->SetTitle(m_ylabel.c_str());
281     histo->GetXaxis()->CenterTitle();
282     histo->GetYaxis()->CenterTitle();
283     if(do_only_selected_samples) histo->SetLineColor(this->GetColor(onlyindex[0]));
284     return histo;
285     }
286    
287    
288     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)
289     {
290     vector<float> binning;
291     for(int i=0;i<=m_nbins;i++)
292     {
293     binning.push_back(((float)(m_maxx-m_minx)/((float)m_nbins))*i+m_minx);
294     }
295    
296     TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity, onlyindex,drawsignal);
297     return histo;
298     }
299     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 drawsignal=false)
300     {
301     stringstream h_histoname;
302     h_histoname<<"h_"<<m_histoname;
303     THStack thestack("thestack",m_histoname.c_str());
304    
305     stringstream drawthis;
306     drawthis<<m_var<<">>tempdrawhisto";
307    
308     TH1F *histogroups[this->ngroups+1];
309     int bookedhistos[this->ngroups+1];
310     for(int ih=0;ih<=this->ngroups;ih++) bookedhistos[ih]=0;
311    
312     for (unsigned int isample=0;isample<(this->collection).size();isample++) {
313     if(!drawsignal&&(this->collection)[isample].is_signal) continue;
314     if((this->collection)[isample].is_active&&(this->collection)[isample].is_data==m_is_data) {//fills mc if we want mc, else fills data.
315     TH1F *tempdrawhisto = new TH1F("tempdrawhisto","",m_nbins,m_minx,m_maxx);
316     (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*"weight");
317     if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));
318     tempdrawhisto->SetFillColor(this->GetColor(isample));
319     if(bookedhistos[(this->collection)[isample].groupindex]==0) {
320     tempdrawhisto->SetName(GetNumericHistoName().c_str());
321     histogroups[(this->collection)[isample].groupindex]=(TH1F*)tempdrawhisto->Clone();
322     bookedhistos[(this->collection)[isample].groupindex]=1;
323     }
324     else
325     {
326     histogroups[(this->collection)[isample].groupindex]->Add(tempdrawhisto);
327     bookedhistos[(this->collection)[isample].groupindex]++;
328     }
329     tempdrawhisto->Delete();
330     }
331     }
332     vector<int> ordered_indices;
333     vector<float> ordered_integrals;
334    
335     for(int ig=0;ig<=this->ngroups;ig++) {
336     if(bookedhistos[ig]>0) {
337     ordered_indices.push_back(ig);
338     ordered_integrals.push_back(histogroups[ig]->Integral());
339     }
340     }
341    
342     /*
343     bubbleSort(ordered_integrals,ordered_indices);
344     // for(int index=ordered_indices.size()-1;index>=0;index--) {
345     for(int index=0;index<ordered_indices.size();index++) {
346     thestack.Add(histogroups[ordered_indices[index]]);
347     }
348     */
349     for(int index=0;index<ordered_indices.size();index++) {
350     thestack.Add(histogroups[ordered_indices[index]]);
351     }
352     return thestack;
353     }
354    
355    
356     vector<int> samplecollection::FindSample(string what)
357     {
358     vector<int> hitcollection;
359     for(unsigned int isam=0;isam<(this->collection).size();isam++)
360     {
361     if(((this->collection)[isam].filename).find(what)!=string::npos) {
362     hitcollection.push_back(isam);
363     }
364     else {
365     }
366     }
367     if(hitcollection.size()==0) hitcollection.push_back(-1);
368     return hitcollection;
369     }
370    
371     Color_t samplecollection::GetColor(string filename) {
372     vector<int> corresponding_samples = this->FindSample(filename);
373     if(corresponding_samples.size()>0) return this->GetColor(corresponding_samples[0]);
374     return kBlack;
375     }
376    
377     Color_t samplecollection::GetColor(int sampleindex) {
378     return this->collection[sampleindex].samplecolor;
379     }
380    
381     TLegend* samplecollection::allbglegend(string title, TH1F *data) {
382     TLegend *leg = new TLegend(0.65,0.60,0.89,0.77);
383     if(title!="") leg->SetHeader(title.c_str());
384     if(data->GetName()!="nothing") leg->AddEntry(data,"Data","p");
385     leg->SetFillColor(kWhite);
386     leg->SetBorderSize(0);
387     leg->SetLineColor(kWhite);
388    
389     TH1F *fakehistos[(this->collection).size()];
390     bool donealready[(this->collection).size()];
391     for(int i=0;i<(this->collection).size();i++) donealready[i]=false;
392     for(int isample=0;isample<(this->collection).size();isample++) {
393     if((this->collection)[isample].is_data||(this->collection)[isample].is_signal) continue;
394    
395     if(!donealready[(this->collection)[isample].groupindex]) {
396     fakehistos[(this->collection)[isample].groupindex] = new TH1F(GetNumericHistoName().c_str(),"",1,0,1);
397     fakehistos[(this->collection)[isample].groupindex]->Fill(1);
398     fakehistos[(this->collection)[isample].groupindex]->SetFillColor(this->GetColor(isample));
399     leg->AddEntry(fakehistos[(this->collection)[isample].groupindex],((this->collection)[isample].samplename).c_str(),"f");
400     donealready[(this->collection)[isample].groupindex]=true;
401     }
402     }
403     TText *writeline1 = write_text(0.77,0.87,"CMS Preliminary 2011");
404     stringstream lumitext;
405     lumitext<<"#sqrt{s}=7, L="<<lumi<<" pb^{-1}";
406     TText *writeline2 = write_text(0.77,0.83,lumitext.str().c_str());
407     writeline1->SetTextSize(0.03);
408     writeline2->SetTextSize(0.03);
409     writeline1->Draw();
410     writeline2->Draw();
411     leg->SetTextFont(42);
412     return leg;
413     }
414    
415     TLegend* samplecollection::allbglegend(string title="") {
416     Int_t currlevel=gErrorIgnoreLevel;
417     gErrorIgnoreLevel=5000;
418     TH1F *blub = new TH1F("nothing","nothing",1,0,1);
419     gErrorIgnoreLevel=currlevel;//we know this possibly replaces a previous histo, but we don't care since it's fake anyway.
420     return this->allbglegend(title,blub);
421     }
422    
423     void set_treename(string treen="events") {
424     cout << "Treename has been set to " << treen << endl;
425     if(treen=="PFevents"||treen=="events") treename=treen;
426     else write_error("sample::set_treename","Setting the treename failed as you provided an invalid tree name.");
427     }
428