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