ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/SampleClass.C
Revision: 1.3
Committed: Fri Jul 8 06:20:31 2011 UTC (13 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.2: +4 -0 lines
Log Message:
Introduced some very early and preliminary sideband checks; fixed bug in Create_All_Plots (global canvas was passed to a lot of functions and at some point wasn't saved correctly anymore); added response correction function (correcting for |pfmet-pt|/pt deficiency); changed the fitting function (cb_fit_...) to return the result (functions) as a vector with 0 being the down variation, 1 the central and 2 the up variation instead of using global tf1s that could be (falsely) misused; updated to sample with 1068 /pb; introduced very simplistic histo integral function (to be improved); other changes ...

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 buchmann 1.2 bool addoverunderflowbins=false;
32 buchmann 1.1 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 buchmann 1.2 vector<float> get_optimal_binsize(string variable, TCut cut,int nbins, float low, float high);
151    
152 buchmann 1.1 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);
153     vector<int> FindSample(string what);
154     void ListSamples();
155     bool do_sample(int thissample, vector<int> &selected_samples);
156     Color_t GetColor(string filename);
157     Color_t GetColor(int sampleindex);
158    
159     TLegend* allbglegend(string title);
160     TLegend* allbglegend(string title, TH1F *data);
161     };
162    
163     samplecollection::samplecollection(string m_name)
164     {
165     this->name=m_name;
166     if(Verbosity>0) cout << "Initiated sample collection " << this->name << " " << endl;
167     }
168    
169     void samplecollection::ListSamples()
170     {
171     cout << "---------------------------------------------------------------------------------------------------" << endl;
172     cout << "Listing all " << this->nsamples << " sample(s) of the sample collection called " << this->name << " : " << endl;
173     if(this->ndatasamples>0) {
174     cout << "Data sample(s): " << endl;
175     for (int isamp=0;isamp<this->collection.size();isamp++)
176     {
177     if((this->collection)[isamp].is_data) cout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
178     }
179     }
180     if(this->nmcsamples>0) {
181     cout << "MC sample(s): " << endl;
182     for (int isamp=0;isamp<this->collection.size();isamp++)
183     {
184     if(!(this->collection)[isamp].is_data) cout << " - " << (this->collection)[isamp].samplename << " from " << (this->collection)[isamp].filename << endl;
185     }
186     }
187     cout << "---------------------------------------------------------------------------------------------------" << endl;
188     }
189    
190     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)
191     {
192     sample NewSample(m_filename,m_samplename,m_Nentries,m_xs,is_data,m_is_signal,groupindex,newcolor);
193     (this->collection).push_back(NewSample);
194     if((this->collection).size()==1) {
195     this->nmcsamples=0;
196     this->ndatasamples=0;
197     this->nsamples=0;
198     this->ngroups=0;
199     }
200     if(groupindex>this->ngroups) this->ngroups=groupindex;
201     if(is_data) this->ndatasamples++;
202     else this->nmcsamples++;
203     this->nsamples=this->collection.size();
204     }
205    
206     bool samplecollection::do_sample(int thissample, vector<int> &selected_samples)
207     {
208     bool drawit=false;
209     for(int isel=0;isel<selected_samples.size();isel++)
210     {
211     if(selected_samples[isel]==thissample) {
212     drawit=true;
213     break;
214     }
215     }
216     return drawit;
217     }
218    
219     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)
220     {
221     vector<int> emptyvector;
222     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);
223     return histo;
224     }
225    
226     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)
227     {
228     vector<int> emptyvector;
229     TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity,emptyvector,drawsignal);
230     return histo;
231     }
232    
233     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) {
234     if(Verbosity>0) cout << endl << endl;
235     if(Verbosity>0) cout << "-------------------------------------------------------------------------------------" << endl;
236     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;
237     if(HUSH==0) cout << "Drawing histo called " << m_histoname << "... " << endl;
238     bool do_only_selected_samples=false;
239     if(onlyindex.size()>0&&onlyindex[0]!=-1) {
240     if(Verbosity>0) {cout << "Requested to only draw sample corresponding to the following sample(s) : " << endl;}
241     for(int is=0;is<onlyindex.size();is++) {
242     if(Verbosity>0) cout << " - " << (this->collection)[onlyindex[is]].filename << " (sample index: " << onlyindex[is] << ")" << endl;
243     }
244     do_only_selected_samples=true;
245     }
246     if(onlyindex.size()==1&&onlyindex[0]==-1) {
247     do_only_selected_samples=true; // this is the special case when we request a non-existing sample - we shouldn't draw anything then.
248     onlyindex.clear();
249     }
250     stringstream h_histoname;
251     h_histoname<<"h_"<<m_histoname;
252     float binningarray[binning.size()+1];
253     for(int i=0;i<binning.size();i++) {binningarray[i]=binning[i];}
254     TH1F *histo = new TH1F(m_histoname.c_str(),"",binning.size()-1,binningarray);
255     histo->Sumw2();
256    
257     stringstream drawthis;
258     drawthis<<m_var<<">>tempdrawhisto";
259    
260     for (unsigned int isample=0;isample<(this->collection).size();isample++) {
261     bool use_this_sample=false;
262     if(!(this->collection)[isample].is_active) continue; // skip inactive samples right away
263     if(((this->collection)[isample].is_data==m_is_data)&&(!do_only_selected_samples)) {
264     if(drawsignal==false&&(this->collection)[isample].is_signal==true) continue;
265     use_this_sample=true;
266     }
267     if(do_only_selected_samples&&this->do_sample(isample,onlyindex)) {
268     use_this_sample=true;
269     }
270     TH1F *tempdrawhisto = new TH1F("tempdrawhisto","tempdrawhisto",binning.size()-1,binningarray);
271     tempdrawhisto->Sumw2();
272     if(use_this_sample) {
273     if(Verbosity>0) cout << "[samplecollection::Draw] : Added contribution from sample " << (this->collection)[isample].samplename << endl;
274 buchmann 1.2 (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*"weight");//this weight is based on PU etc. not XS
275     if(addoverunderflowbins) {
276     //now also adding the overflow & underflow bins:
277     tempdrawhisto->SetBinContent(tempdrawhisto->GetNbinsX(),tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()+1));
278     tempdrawhisto->SetBinError(tempdrawhisto->GetNbinsX(),TMath::Sqrt(tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX())));
279     tempdrawhisto->SetBinContent(1,tempdrawhisto->GetBinContent(0));
280     tempdrawhisto->SetBinError(1,TMath::Sqrt(tempdrawhisto->GetBinContent(1)));
281     }
282    
283     if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));//weight applied here is XS & N(entries)
284 buchmann 1.1 histo->Add(tempdrawhisto);
285     }
286     tempdrawhisto->Delete();
287     }//end of loop over isample
288     if(Verbosity>0) cout << "Histo has been filled and now contains " << histo->Integral() << " points (integral)" << endl;
289     histo->GetXaxis()->SetTitle(m_xlabel.c_str());
290     histo->GetYaxis()->SetTitle(m_ylabel.c_str());
291     histo->GetXaxis()->CenterTitle();
292     histo->GetYaxis()->CenterTitle();
293     if(do_only_selected_samples) histo->SetLineColor(this->GetColor(onlyindex[0]));
294     return histo;
295     }
296    
297    
298     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)
299     {
300     vector<float> binning;
301     for(int i=0;i<=m_nbins;i++)
302     {
303     binning.push_back(((float)(m_maxx-m_minx)/((float)m_nbins))*i+m_minx);
304     }
305    
306     TH1F *histo = this->Draw(m_histoname,m_var, binning, m_xlabel, m_ylabel, Cut, m_is_data, luminosity, onlyindex,drawsignal);
307     return histo;
308     }
309     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)
310     {
311     stringstream h_histoname;
312     h_histoname<<"h_"<<m_histoname;
313     THStack thestack("thestack",m_histoname.c_str());
314    
315     stringstream drawthis;
316     drawthis<<m_var<<">>tempdrawhisto";
317    
318     TH1F *histogroups[this->ngroups+1];
319     int bookedhistos[this->ngroups+1];
320     for(int ih=0;ih<=this->ngroups;ih++) bookedhistos[ih]=0;
321    
322     for (unsigned int isample=0;isample<(this->collection).size();isample++) {
323     if(!drawsignal&&(this->collection)[isample].is_signal) continue;
324     if((this->collection)[isample].is_active&&(this->collection)[isample].is_data==m_is_data) {//fills mc if we want mc, else fills data.
325     TH1F *tempdrawhisto = new TH1F("tempdrawhisto","",m_nbins,m_minx,m_maxx);
326     (this->collection)[isample].events->Draw(drawthis.str().c_str(),(essentialcut&&Cut)*"weight");
327 buchmann 1.2
328     if(addoverunderflowbins) {
329     //now also adding the overflow & underflow bins:
330     tempdrawhisto->SetBinContent(tempdrawhisto->GetNbinsX(),tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX()+1));
331     tempdrawhisto->SetBinError(tempdrawhisto->GetNbinsX(),TMath::Sqrt(tempdrawhisto->GetBinContent(tempdrawhisto->GetNbinsX())));
332     tempdrawhisto->SetBinContent(1,tempdrawhisto->GetBinContent(0));
333     tempdrawhisto->SetBinError(1,TMath::Sqrt(tempdrawhisto->GetBinContent(1)));
334     }
335    
336 buchmann 1.1 if(!(this->collection)[isample].is_data) tempdrawhisto->Scale(luminosity*((this->collection)[isample].weight));
337     tempdrawhisto->SetFillColor(this->GetColor(isample));
338     if(bookedhistos[(this->collection)[isample].groupindex]==0) {
339     tempdrawhisto->SetName(GetNumericHistoName().c_str());
340     histogroups[(this->collection)[isample].groupindex]=(TH1F*)tempdrawhisto->Clone();
341     bookedhistos[(this->collection)[isample].groupindex]=1;
342     }
343     else
344     {
345     histogroups[(this->collection)[isample].groupindex]->Add(tempdrawhisto);
346     bookedhistos[(this->collection)[isample].groupindex]++;
347     }
348     tempdrawhisto->Delete();
349     }
350     }
351     vector<int> ordered_indices;
352     vector<float> ordered_integrals;
353    
354     for(int ig=0;ig<=this->ngroups;ig++) {
355     if(bookedhistos[ig]>0) {
356     ordered_indices.push_back(ig);
357     ordered_integrals.push_back(histogroups[ig]->Integral());
358     }
359     }
360    
361     /*
362     bubbleSort(ordered_integrals,ordered_indices);
363     // for(int index=ordered_indices.size()-1;index>=0;index--) {
364     for(int index=0;index<ordered_indices.size();index++) {
365     thestack.Add(histogroups[ordered_indices[index]]);
366     }
367     */
368     for(int index=0;index<ordered_indices.size();index++) {
369     thestack.Add(histogroups[ordered_indices[index]]);
370     }
371     return thestack;
372     }
373    
374    
375     vector<int> samplecollection::FindSample(string what)
376     {
377     vector<int> hitcollection;
378     for(unsigned int isam=0;isam<(this->collection).size();isam++)
379     {
380     if(((this->collection)[isam].filename).find(what)!=string::npos) {
381     hitcollection.push_back(isam);
382     }
383     else {
384     }
385     }
386     if(hitcollection.size()==0) hitcollection.push_back(-1);
387     return hitcollection;
388     }
389    
390     Color_t samplecollection::GetColor(string filename) {
391     vector<int> corresponding_samples = this->FindSample(filename);
392     if(corresponding_samples.size()>0) return this->GetColor(corresponding_samples[0]);
393     return kBlack;
394     }
395    
396     Color_t samplecollection::GetColor(int sampleindex) {
397     return this->collection[sampleindex].samplecolor;
398     }
399    
400     TLegend* samplecollection::allbglegend(string title, TH1F *data) {
401     TLegend *leg = new TLegend(0.65,0.60,0.89,0.77);
402     if(title!="") leg->SetHeader(title.c_str());
403     if(data->GetName()!="nothing") leg->AddEntry(data,"Data","p");
404     leg->SetFillColor(kWhite);
405     leg->SetBorderSize(0);
406     leg->SetLineColor(kWhite);
407    
408     TH1F *fakehistos[(this->collection).size()];
409     bool donealready[(this->collection).size()];
410     for(int i=0;i<(this->collection).size();i++) donealready[i]=false;
411     for(int isample=0;isample<(this->collection).size();isample++) {
412     if((this->collection)[isample].is_data||(this->collection)[isample].is_signal) continue;
413    
414     if(!donealready[(this->collection)[isample].groupindex]) {
415     fakehistos[(this->collection)[isample].groupindex] = new TH1F(GetNumericHistoName().c_str(),"",1,0,1);
416     fakehistos[(this->collection)[isample].groupindex]->Fill(1);
417     fakehistos[(this->collection)[isample].groupindex]->SetFillColor(this->GetColor(isample));
418     leg->AddEntry(fakehistos[(this->collection)[isample].groupindex],((this->collection)[isample].samplename).c_str(),"f");
419     donealready[(this->collection)[isample].groupindex]=true;
420     }
421     }
422     TText *writeline1 = write_text(0.77,0.87,"CMS Preliminary 2011");
423     stringstream lumitext;
424     lumitext<<"#sqrt{s}=7, L="<<lumi<<" pb^{-1}";
425     TText *writeline2 = write_text(0.77,0.83,lumitext.str().c_str());
426     writeline1->SetTextSize(0.03);
427     writeline2->SetTextSize(0.03);
428     writeline1->Draw();
429     writeline2->Draw();
430     leg->SetTextFont(42);
431     return leg;
432     }
433    
434     TLegend* samplecollection::allbglegend(string title="") {
435     Int_t currlevel=gErrorIgnoreLevel;
436     gErrorIgnoreLevel=5000;
437     TH1F *blub = new TH1F("nothing","nothing",1,0,1);
438     gErrorIgnoreLevel=currlevel;//we know this possibly replaces a previous histo, but we don't care since it's fake anyway.
439     return this->allbglegend(title,blub);
440     }
441    
442     void set_treename(string treen="events") {
443     cout << "Treename has been set to " << treen << endl;
444     if(treen=="PFevents"||treen=="events") treename=treen;
445     else write_error("sample::set_treename","Setting the treename failed as you provided an invalid tree name.");
446     }
447 buchmann 1.2
448     vector<float> samplecollection::get_optimal_binsize(string variable, TCut cut,int nbins, float low, float high) {
449     TH1F *histo = this->Draw("histo",variable,5000,low,high, "finding_optimal_binsize", "events", cut,0,1000);
450     float eventsperbin=histo->Integral()/nbins;
451     vector<float> binning;
452     binning.push_back(high);
453     float runningsum=0;
454     for(int i=histo->GetNbinsX();i>0;i--) {
455     runningsum+=histo->GetBinContent(i);
456     if(TMath::Abs(runningsum-eventsperbin)<0.05*eventsperbin||runningsum>eventsperbin) {
457     binning.push_back(histo->GetBinLowEdge(i));
458     runningsum=0;
459     }
460     }
461     if(runningsum>0) binning.push_back(low);
462     for(int i=0;i<int(binning.size()/2);i++) {
463     float temp=binning[i];
464     binning[i]=binning[binning.size()-1-i];
465     binning[binning.size()-1-i]=temp;
466     }
467    
468     return binning;
469    
470 buchmann 1.3 }
471    
472     void switch_overunderflow(bool newpos=false) {
473     addoverunderflowbins=newpos;
474 buchmann 1.2 }