| 1 |
#include <iostream>
|
| 2 |
#include <fstream>
|
| 3 |
#include <sstream>
|
| 4 |
|
| 5 |
|
| 6 |
#include <TROOT.h>
|
| 7 |
#include <TStyle.h>
|
| 8 |
#include <TTree.h>
|
| 9 |
#include <TFile.h>
|
| 10 |
#include <THStack.h>
|
| 11 |
#include <TCut.h>
|
| 12 |
#include <TF1.h>
|
| 13 |
#include <TCanvas.h>
|
| 14 |
#include <TLegend.h>
|
| 15 |
#include <TText.h>
|
| 16 |
#include <TLatex.h>
|
| 17 |
#include <TSystem.h>
|
| 18 |
#include <TMath.h>
|
| 19 |
#include <TPad.h>
|
| 20 |
#include <TLine.h>
|
| 21 |
#include <TColor.h>
|
| 22 |
|
| 23 |
using namespace std;
|
| 24 |
|
| 25 |
float get_error_R(float nemout, float nemin, float nemouterr, float neminerr);
|
| 26 |
float give_n_eemm_in_correct_error(float nouteemm, float noutem, float ninem, float ninemerror, float noutemerror, float nouteemmerror);
|
| 27 |
float get_new_stat_error(TH1 *totalhisto,float a, float b, float c, float d);
|
| 28 |
|
| 29 |
void compute_new_results(TH1* hdataeemmout, TH1* hdataeemmin, TH1* hdataemout, TH1* hdataemin, TH1* hdataeeout, TH1* hdataeein, TH1* hdatammout, TH1* hdatammin, TH1* hemout, TH1* hemin, TH1* heemmout, TH1* heemmin, TH1* heeout, TH1* heein, TH1* hmmout, TH1* hmmin, TH1* httbareemmout, TH1* httbareemmin, float A, float B, float C, float D, float E, float F, float results[2][3][3][2], float errors[2][3][3][2], float *correctresult, int nsel)
|
| 30 |
{
|
| 31 |
|
| 32 |
/*
|
| 33 |
cout << "hdataeemmout :" << hdataeemmout->GetName() << endl;
|
| 34 |
cout << "hdataeemmin :" << hdataeemmin->GetName() << endl;
|
| 35 |
cout << "hdataemout :" << hdataemout->GetName() << endl;
|
| 36 |
cout << "hdataemin :" << hdataemin->GetName() << endl;
|
| 37 |
cout << "hdataeeout :" << hdataeeout->GetName() << endl;
|
| 38 |
cout << "hdataeein :" << hdataeein->GetName() << endl;
|
| 39 |
cout << "hdatammout :" << hdatammout->GetName() << endl;
|
| 40 |
cout << "hdatammin :" << hdatammin->GetName() << endl;
|
| 41 |
cout << "hemout :" << hemout->GetName() << endl;
|
| 42 |
cout << "hemin :" << hemin->GetName() << endl;
|
| 43 |
cout << "heemmout :" << heemmout->GetName() << endl;
|
| 44 |
cout << "heemmin :" << heemmin->GetName() << endl;
|
| 45 |
cout << "heeout :" << heeout->GetName() << endl;
|
| 46 |
cout << "heein :" << heein->GetName() << endl;
|
| 47 |
cout << "hmmout :" << hmmout->GetName() << endl;
|
| 48 |
cout << "hmmin :" << hmmin->GetName() << endl;
|
| 49 |
cout << "httbareemmout :" << httbareemmout->GetName() << endl;
|
| 50 |
*/
|
| 51 |
|
| 52 |
|
| 53 |
bool debug_compute_results=true; //if you want to debug the method, this will inform you about a lot of the steps.
|
| 54 |
cout << "working on compute_results ";
|
| 55 |
if (nsel==0) cout << "with all samples" << endl;
|
| 56 |
if (nsel==1) {cout << "with all nonZ samples" << endl;}
|
| 57 |
|
| 58 |
//MC
|
| 59 |
float em_inside_MC = hemin->Integral();
|
| 60 |
float em_outside_MC = hemout->Integral();
|
| 61 |
float ee_inside_MC = heein->Integral();
|
| 62 |
float ee_outside_MC = heeout->Integral();
|
| 63 |
float mm_inside_MC = hmmin->Integral();
|
| 64 |
float mm_outside_MC = hmmout->Integral();
|
| 65 |
float eemm_outside_MC = heemmout->Integral();
|
| 66 |
|
| 67 |
float em_inside_MC_error = get_new_stat_error(hemin,E,F,0,0);
|
| 68 |
float em_outside_MC_error = get_new_stat_error(hemout,A,B,C,D);
|
| 69 |
float eemm_outside_MC_error = get_new_stat_error(heemmout,A,B,C,D);
|
| 70 |
|
| 71 |
float ee_inside_MC_error = get_new_stat_error(heein,E,F,0.0,0.0);
|
| 72 |
float ee_outside_MC_error = get_new_stat_error(heeout,A,B,C,D);
|
| 73 |
|
| 74 |
float mm_inside_MC_error = get_new_stat_error(hmmin,E,F,0.0,0.0);
|
| 75 |
float mm_outside_MC_error = get_new_stat_error(hmmout,A,B,C,D);
|
| 76 |
|
| 77 |
//data
|
| 78 |
float em_inside_data = hdataemin->Integral();
|
| 79 |
float em_outside_data = hdataemout->Integral();
|
| 80 |
float ee_inside_data = hdataeein->Integral();
|
| 81 |
float ee_outside_data = hdataeeout->Integral();
|
| 82 |
float mm_inside_data = hdatammin->Integral();
|
| 83 |
float mm_outside_data = hdatammout->Integral();
|
| 84 |
float eemm_outside_data = hdataeemmout->Integral();
|
| 85 |
|
| 86 |
if (em_outside_data==0)
|
| 87 |
{
|
| 88 |
cout << "The estimate for data cannot be computed as there are no data em events in the sidebands (!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!)" << endl;
|
| 89 |
}
|
| 90 |
if (em_outside_MC==0)
|
| 91 |
{
|
| 92 |
cout << "The estimate for MC cannot be computed as there are no MC em events in the sidebands (!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!)" << endl;
|
| 93 |
}
|
| 94 |
if ((em_outside_data==0)||(em_outside_MC==0)) return;
|
| 95 |
|
| 96 |
float eemm_inside_data = eemm_outside_data * (em_inside_data / em_outside_data);
|
| 97 |
|
| 98 |
float em_inside_data_error = TMath::Sqrt(em_inside_data);
|
| 99 |
float em_outside_data_error = TMath::Sqrt(em_outside_data);
|
| 100 |
//float eemm_inside_data_error = TMath::Sqrt(eemm_inside_data);//the error is NOT calculated in this way (it's the method's job to calculate it!) so this is commented out to illustrate it!
|
| 101 |
float eemm_outside_data_error = TMath::Sqrt(eemm_outside_data);
|
| 102 |
|
| 103 |
float ee_inside_data_error = TMath::Sqrt(ee_inside_data);
|
| 104 |
float ee_outside_data_error = TMath::Sqrt(ee_outside_data);
|
| 105 |
|
| 106 |
float mm_inside_data_error = TMath::Sqrt(mm_inside_data);
|
| 107 |
float mm_outside_data_error = TMath::Sqrt(mm_outside_data);
|
| 108 |
|
| 109 |
if(debug_compute_results){
|
| 110 |
cout << endl << endl;
|
| 111 |
cout << "BEFORE PREDICTING ANYTHING, THESE ARE OUR NUMBERS:" << endl;
|
| 112 |
cout << "type \t inside \t outside" << endl;
|
| 113 |
cout << "em/MC \t" << em_inside_MC << "+/-" << em_inside_MC_error << "\t" << em_outside_MC << "+/-" << em_outside_MC_error << endl;
|
| 114 |
cout << "em/data \t" << em_inside_data << "+/-" << em_inside_data_error << "\t" << em_outside_data << "+/-" << em_outside_data_error << endl;
|
| 115 |
cout << "ee/MC \t" << ee_inside_MC << "+/-" << ee_inside_MC_error << "\t" << ee_outside_MC << "+/-" << ee_outside_MC_error << endl;
|
| 116 |
cout << "ee/data" << "\t" << ee_inside_data << "+/-" << ee_inside_data_error << "\t" << ee_outside_data << "+/-" << ee_outside_data_error << endl;
|
| 117 |
cout << "mm/MC \t" << mm_inside_MC << "+/-" << mm_inside_MC_error << "\t" << mm_outside_MC << "+/-" << mm_outside_MC_error << endl;
|
| 118 |
cout << "mm/data \t" << mm_inside_data << "+/-" << mm_inside_data_error << "\t" << mm_outside_data << "+/-" << mm_outside_data_error << endl;
|
| 119 |
cout << endl << endl;
|
| 120 |
}
|
| 121 |
|
| 122 |
//predictions
|
| 123 |
|
| 124 |
float eemm_inside_MC_pred = eemm_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 125 |
float ee_inside_MC_pred = ee_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 126 |
float mm_inside_MC_pred = mm_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 127 |
|
| 128 |
float eemm_inside_data_pred = eemm_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 129 |
float ee_inside_data_pred = ee_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 130 |
float mm_inside_data_pred = mm_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 131 |
|
| 132 |
|
| 133 |
|
| 134 |
/*
|
| 135 |
TCanvas *test_distributions = new TCanvas("test_distributions","distribution verification");
|
| 136 |
test_distributions->cd();
|
| 137 |
hem->Draw();
|
| 138 |
cout << endl;
|
| 139 |
*/
|
| 140 |
/*
|
| 141 |
cout << "mc abs args: " << endl;
|
| 142 |
cout << eemm_outside_MC << endl;
|
| 143 |
cout << em_outside_MC << endl;
|
| 144 |
cout << em_inside_MC << endl;
|
| 145 |
cout << em_inside_MC_error << endl;
|
| 146 |
cout << em_outside_MC_error << endl;
|
| 147 |
cout << eemm_outside_MC_error << endl;
|
| 148 |
*/
|
| 149 |
|
| 150 |
|
| 151 |
float R_data = em_inside_data / em_outside_data;
|
| 152 |
float R_data_error = get_error_R(em_outside_data,em_inside_data, em_outside_data_error,em_inside_data_error);
|
| 153 |
float R_MC = em_inside_MC / em_outside_MC;
|
| 154 |
float R_MC_error = get_error_R(em_outside_MC,em_inside_MC, em_outside_MC_error,em_inside_MC_error);
|
| 155 |
|
| 156 |
float eemm_inside_MC_pred_error = give_n_eemm_in_correct_error(eemm_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error,em_outside_MC_error,eemm_outside_MC_error);
|
| 157 |
float eemm_inside_data_pred_error = give_n_eemm_in_correct_error(eemm_outside_data, em_outside_data, em_inside_data, em_inside_data_error, em_outside_data_error,eemm_outside_data_error);
|
| 158 |
if(debug_compute_results){
|
| 159 |
cout << "RESULTS " << endl;
|
| 160 |
cout << "TOTALS : " << endl;
|
| 161 |
cout << "Data: \t" << eemm_inside_data_pred << "+/- " << eemm_inside_data_pred_error << "\t R = " << R_data << " +/- " << R_data_error;
|
| 162 |
cout << endl;
|
| 163 |
cout << "MC \t " << eemm_inside_MC_pred << " +/- " << eemm_inside_MC_pred_error << "\t R = " << R_MC << " +/- " << R_MC_error << endl;
|
| 164 |
cout << endl;
|
| 165 |
}
|
| 166 |
|
| 167 |
float ee_inside_MC_pred_error=give_n_eemm_in_correct_error(ee_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error, em_outside_MC_error, ee_outside_MC_error);
|
| 168 |
float ee_inside_data_pred_error=give_n_eemm_in_correct_error(ee_outside_data, em_outside_data, em_inside_data, em_inside_data_error,em_outside_data_error,ee_outside_data_error);
|
| 169 |
if(debug_compute_results){
|
| 170 |
cout << "FOR ELECTRONS: " << endl;
|
| 171 |
cout << "Data \t " << ee_inside_data_pred << " +/- " << ee_inside_data_pred_error << "\t R = " << R_data << endl;
|
| 172 |
cout << "MC \t " << ee_inside_MC_pred << " +/- " << ee_inside_MC_pred_error << "\t R = " << R_MC << endl;
|
| 173 |
}
|
| 174 |
float mm_inside_MC_pred_error=give_n_eemm_in_correct_error(mm_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error,em_outside_MC_error,mm_outside_MC_error);
|
| 175 |
float mm_inside_data_pred_error=give_n_eemm_in_correct_error(mm_outside_data, em_outside_data, em_inside_data, em_inside_data_error,em_outside_data_error,mm_outside_data_error);
|
| 176 |
|
| 177 |
if(debug_compute_results){
|
| 178 |
cout << "FOR MUONS: " << endl;
|
| 179 |
cout << "Data \t " << mm_inside_data_pred << " +/- " << ee_inside_data_pred_error << "\t R = " << em_inside_data / em_outside_data << endl;
|
| 180 |
cout << "MC \t " << mm_inside_MC_pred << " +/- " << mm_inside_MC_pred_error << "\t R = " << em_inside_MC / em_outside_MC << endl;
|
| 181 |
|
| 182 |
cout << endl << endl;
|
| 183 |
cout << "the correct result would be " << httbareemmin->Integral() << endl;
|
| 184 |
}
|
| 185 |
//first: MC
|
| 186 |
results[nsel][0][0][0]=em_outside_MC; errors[nsel][0][0][0]=em_outside_MC_error;
|
| 187 |
results[nsel][0][1][0]=em_inside_MC; errors[nsel][0][1][0]=em_inside_MC_error;
|
| 188 |
results[nsel][0][2][0]=R_MC; errors[nsel][0][2][0]=R_MC_error;
|
| 189 |
|
| 190 |
results[nsel][1][0][0]=ee_outside_MC; errors[nsel][1][0][0]=ee_outside_MC_error;
|
| 191 |
results[nsel][1][1][0]=ee_inside_MC; errors[nsel][1][1][0]=ee_inside_MC_error;
|
| 192 |
results[nsel][1][2][0]=ee_inside_MC_pred; errors[nsel][1][2][0]=ee_inside_MC_pred_error;
|
| 193 |
|
| 194 |
results[nsel][2][0][0]=mm_outside_MC; errors[nsel][2][0][0]=mm_outside_MC_error;
|
| 195 |
results[nsel][2][1][0]=mm_inside_MC; errors[nsel][2][1][0]=mm_inside_MC_error;
|
| 196 |
results[nsel][2][2][0]=mm_inside_MC_pred; errors[nsel][2][2][0]=mm_inside_MC_pred_error;
|
| 197 |
|
| 198 |
//second: data
|
| 199 |
results[nsel][0][0][1]=em_outside_data; errors[nsel][0][0][1]=em_outside_data_error;
|
| 200 |
results[nsel][0][1][1]=em_inside_data; errors[nsel][0][1][1]=em_inside_data_error;
|
| 201 |
results[nsel][0][2][1]=R_data; errors[nsel][0][2][1]=R_data_error;
|
| 202 |
|
| 203 |
results[nsel][1][0][1]=ee_outside_data; errors[nsel][1][0][1]=ee_outside_data_error;
|
| 204 |
results[nsel][1][1][1]=ee_inside_data; errors[nsel][1][1][1]=ee_inside_data_error;
|
| 205 |
results[nsel][1][2][1]=ee_inside_data_pred; errors[nsel][1][2][1]=ee_inside_data_pred_error;
|
| 206 |
|
| 207 |
results[nsel][2][0][1]=mm_outside_data; errors[nsel][2][0][1]=mm_outside_data_error;
|
| 208 |
results[nsel][2][1][1]=mm_inside_data; errors[nsel][2][1][1]=mm_inside_data_error;
|
| 209 |
results[nsel][2][2][1]=mm_inside_data_pred; errors[nsel][2][2][1]=mm_inside_data_pred_error;
|
| 210 |
|
| 211 |
*correctresult=httbareemmin->Integral();
|
| 212 |
|
| 213 |
}
|
| 214 |
|
| 215 |
|
| 216 |
//**********************************************************************
|
| 217 |
|
| 218 |
float get_new_stat_error(TH1 *totalhisto,float a, float b, float c, float d)
|
| 219 |
{
|
| 220 |
float err_sum=0;
|
| 221 |
for (int i=totalhisto->FindBin(a);i<totalhisto->FindBin(b);i++)
|
| 222 |
{
|
| 223 |
float currerr = totalhisto->GetBinError(i);
|
| 224 |
err_sum += currerr * currerr;
|
| 225 |
}
|
| 226 |
|
| 227 |
if (TMath::Abs(c-d)>0.01)//if c&d are different (or sufficiently different -- they are floats!) we calculate the integral in between the two as well and add it!
|
| 228 |
{
|
| 229 |
for (int i=totalhisto->FindBin(c);i<totalhisto->FindBin(d);i++)
|
| 230 |
{
|
| 231 |
float currerr = totalhisto->GetBinError(i);
|
| 232 |
err_sum += currerr * currerr;
|
| 233 |
}
|
| 234 |
}
|
| 235 |
float err_result = TMath::Sqrt(err_sum);
|
| 236 |
return err_result;
|
| 237 |
}
|
| 238 |
|
| 239 |
|
| 240 |
|
| 241 |
void compute_wrong_results(TH1* hdataeemm, TH1* hdataem, TH1* hdataee, TH1* hdatamm, TH1* hem, TH1* heemm, TH1* hee, TH1* hmm, TH1* httbareemm, float A, float B, float C, float D, float E, float F, float results[2][3][3], float errors[2][3][3], float *correctresult, int nsel)
|
| 242 |
{
|
| 243 |
|
| 244 |
bool debug_compute_results=true; //if you want to debug the method, this will inform you about a lot of the steps.
|
| 245 |
cout << "working on compute_results ";
|
| 246 |
if (nsel==0) cout << "with all samples" << endl;
|
| 247 |
if (nsel==1) cout << "with all nonZ samples" << endl;
|
| 248 |
|
| 249 |
if(debug_compute_results)
|
| 250 |
{
|
| 251 |
cout << hdataeemm->Integral() << endl;
|
| 252 |
cout << hdataeemm->Integral() << endl;
|
| 253 |
cout << hdataem->Integral() << endl;
|
| 254 |
cout << hdataee->Integral() << endl;
|
| 255 |
cout << hdatamm->Integral() << endl;
|
| 256 |
cout << "hem" << hem->Integral() << endl;
|
| 257 |
cout << "heemm" << heemm->Integral() << endl;
|
| 258 |
cout << "hee" << hee->Integral() << endl;
|
| 259 |
cout << "hmm" << hmm->Integral() << endl;
|
| 260 |
cout << httbareemm->Integral() << endl;
|
| 261 |
|
| 262 |
|
| 263 |
}
|
| 264 |
/* //MC
|
| 265 |
float em_inside_MC = (hem->Integral(hem->FindBin(E),hem->FindBin(F)));
|
| 266 |
float em_outside_MC = (hem->Integral(hem->FindBin(A),hem->FindBin(B))+hem->Integral(hem->FindBin(C),hem->FindBin(D)));
|
| 267 |
float ee_inside_MC = (hee->Integral(hee->FindBin(E),hee->FindBin(F)));
|
| 268 |
float ee_outside_MC = (hee->Integral(hee->FindBin(A),hee->FindBin(B))+hee->Integral(hee->FindBin(C),hee->FindBin(D)));
|
| 269 |
float mm_inside_MC = (hmm->Integral(hmm->FindBin(E),hmm->FindBin(F)));
|
| 270 |
float mm_outside_MC = (hmm->Integral(hmm->FindBin(A),hmm->FindBin(B))+hmm->Integral(hmm->FindBin(C),hmm->FindBin(D)));
|
| 271 |
float eemm_outside_MC = heemm->Integral(heemm->FindBin(A),heemm->FindBin(B)) + heemm->Integral(heemm->FindBin(C),heemm->FindBin(D));
|
| 272 |
|
| 273 |
float em_inside_MC_error = get_new_stat_error(hem,E,F,0,0);
|
| 274 |
float em_outside_MC_error = get_new_stat_error(hem,A,B,C,D);
|
| 275 |
float eemm_outside_MC_error = get_new_stat_error(heemm,A,B,C,D);
|
| 276 |
|
| 277 |
float ee_inside_MC_error = get_new_stat_error(hee,E,F,0.0,0.0);
|
| 278 |
float ee_outside_MC_error = get_new_stat_error(hee,A,B,C,D);
|
| 279 |
|
| 280 |
float mm_inside_MC_error = get_new_stat_error(hmm,E,F,0.0,0.0);
|
| 281 |
float mm_outside_MC_error = get_new_stat_error(hmm,A,B,C,D);
|
| 282 |
|
| 283 |
//data
|
| 284 |
float em_inside_data = hdataem->Integral(hdataem->FindBin(E),hdataem->FindBin(F));
|
| 285 |
float em_outside_data = hdataem->Integral(hdataem->FindBin(A),hdataem->FindBin(B)) + hdataem->Integral(hdataem->FindBin(C),hdataem->FindBin(D));
|
| 286 |
float ee_inside_data = hdataee->Integral(hdataee->FindBin(E),hdataee->FindBin(F));
|
| 287 |
float ee_outside_data = hdataee->Integral(hdataee->FindBin(A),hdataee->FindBin(B)) + hdataee->Integral(hdataee->FindBin(C),hdataee->FindBin(D));
|
| 288 |
float mm_inside_data = hdatamm->Integral(hdatamm->FindBin(E),hdatamm->FindBin(F));
|
| 289 |
float mm_outside_data = hdatamm->Integral(hdatamm->FindBin(A),hdatamm->FindBin(B)) + hdatamm->Integral(hdatamm->FindBin(C),hdatamm->FindBin(D));
|
| 290 |
float eemm_outside_data = hdataeemm->Integral(hdataeemm->FindBin(A),hdataeemm->FindBin(B)) + hdataeemm->Integral(hdataeemm->FindBin(C),heemm->FindBin(D));
|
| 291 |
float eemm_inside_data = eemm_outside_data * (em_inside_data / em_outside_data);
|
| 292 |
|
| 293 |
float em_inside_data_error = TMath::Sqrt(em_inside_data);
|
| 294 |
float em_outside_data_error = TMath::Sqrt(em_outside_data);
|
| 295 |
//float eemm_inside_data_error = TMath::Sqrt(eemm_inside_data);//the error is NOT calculated in this way (it's the method's job to calculate it!) so this is commented out to illustrate it!
|
| 296 |
float eemm_outside_data_error = TMath::Sqrt(eemm_outside_data);
|
| 297 |
|
| 298 |
float ee_inside_data_error = TMath::Sqrt(ee_inside_data);
|
| 299 |
float ee_outside_data_error = TMath::Sqrt(ee_outside_data);
|
| 300 |
|
| 301 |
float mm_inside_data_error = TMath::Sqrt(mm_inside_data);
|
| 302 |
float mm_outside_data_error = TMath::Sqrt(mm_outside_data);
|
| 303 |
|
| 304 |
//predictions
|
| 305 |
|
| 306 |
float eemm_inside_MC_pred = eemm_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 307 |
float ee_inside_MC_pred = ee_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 308 |
float mm_inside_MC_pred = mm_outside_MC * (em_inside_MC/em_outside_MC); //this is our prediction!
|
| 309 |
|
| 310 |
float eemm_inside_data_pred = eemm_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 311 |
float ee_inside_data_pred = ee_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 312 |
float mm_inside_data_pred = mm_outside_data * (em_inside_data/em_outside_data); //this is our prediction!
|
| 313 |
|
| 314 |
|
| 315 |
if(debug_compute_results){
|
| 316 |
cout << endl << endl;
|
| 317 |
cout << "BEFORE PREDICTING ANYTHING, THESE ARE OUR NUMBERS:" << endl;
|
| 318 |
cout << "em inside (MC) " << em_inside_MC << "+/-" << em_inside_MC_error << endl;
|
| 319 |
cout << "em outside (MC) " << em_outside_MC << "+/-" << em_outside_MC_error << endl;
|
| 320 |
cout << "em inside (data) " << em_inside_data << "+/-" << em_inside_data_error << endl;
|
| 321 |
cout << "em outside (data) " << em_outside_data << "+/-" << em_outside_data_error << endl;
|
| 322 |
cout << endl;
|
| 323 |
cout << "ee inside (MC) " << ee_inside_MC << "+/-" << ee_inside_MC_error << endl;
|
| 324 |
cout << "ee outside (MC) " << ee_outside_MC << "+/-" << ee_outside_MC_error << endl;
|
| 325 |
cout << "ee inside (data) " << ee_inside_data << "+/-" << ee_inside_data_error << endl;
|
| 326 |
cout << "ee outside (data) " << ee_outside_data << "+/-" << ee_outside_data_error << endl;
|
| 327 |
cout << endl;
|
| 328 |
cout << "mm inside (MC) " << mm_inside_MC << "+/-" << mm_inside_MC_error << endl;
|
| 329 |
cout << "mm outside (MC) " << mm_outside_MC << "+/-" << mm_outside_MC_error << endl;
|
| 330 |
cout << "mm inside (data) " << mm_inside_data << "+/-" << mm_inside_data_error << endl;
|
| 331 |
cout << "mm outside (data) " << mm_outside_data << "+/-" << mm_outside_data_error << endl;
|
| 332 |
cout << endl << endl << "MC: " << endl;
|
| 333 |
}
|
| 334 |
|
| 335 |
|
| 336 |
float R_data = em_inside_data / em_outside_data;
|
| 337 |
float R_data_error = get_error_R(em_outside_data,em_inside_data, em_outside_data_error,em_inside_data_error);
|
| 338 |
float R_MC = em_inside_MC / em_outside_MC;
|
| 339 |
float R_MC_error = get_error_R(em_outside_MC,em_inside_MC, em_outside_MC_error,em_inside_MC_error);
|
| 340 |
|
| 341 |
float eemm_inside_MC_pred_error = give_n_eemm_in_correct_error(eemm_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error,em_outside_MC_error,eemm_outside_MC_error);
|
| 342 |
float eemm_inside_data_pred_error = give_n_eemm_in_correct_error(eemm_outside_data, em_outside_data, em_inside_data, em_inside_data_error, em_outside_data_error,eemm_outside_data_error);
|
| 343 |
if(debug_compute_results){
|
| 344 |
cout << "RESULTS " << endl;
|
| 345 |
cout << "TOTALS : " << endl;
|
| 346 |
cout << "Data: \t" << eemm_inside_data_pred << "+/- " << eemm_inside_data_pred_error << "\t R = " << R_data << " +/- " << R_data_error;
|
| 347 |
cout << endl;
|
| 348 |
cout << "MC \t " << eemm_inside_MC_pred << " +/- " << eemm_inside_MC_pred_error << "\t R = " << R_MC << " +/- " << R_MC_error << endl;
|
| 349 |
cout << endl << endl;
|
| 350 |
}
|
| 351 |
|
| 352 |
float ee_inside_MC_pred_error=give_n_eemm_in_correct_error(ee_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error, em_outside_MC_error, ee_outside_MC_error);
|
| 353 |
float ee_inside_data_pred_error=give_n_eemm_in_correct_error(ee_outside_data, em_outside_data, em_inside_data, em_inside_data_error,em_outside_data_error,ee_outside_data_error);
|
| 354 |
if(debug_compute_results){
|
| 355 |
cout << "FOR ELECTRONS: " << endl;
|
| 356 |
cout << "Data \t " << ee_inside_data_pred << " +/- " << ee_inside_data_pred_error << "\t R = " << R_data << endl;
|
| 357 |
cout << "MC \t " << ee_inside_MC_pred << " +/- " << ee_inside_MC_pred_error << "\t R = " << R_MC << endl;
|
| 358 |
}
|
| 359 |
float mm_inside_MC_pred_error=give_n_eemm_in_correct_error(mm_outside_MC, em_outside_MC, em_inside_MC, em_inside_MC_error,em_outside_MC_error,mm_outside_MC_error);
|
| 360 |
float ee_inside_data_pred_error=give_n_eemm_in_correct_error(mm_outside_data, em_outside_data, em_inside_data, em_inside_data_error,em_outside_data_error,mm_outside_data_error);
|
| 361 |
|
| 362 |
if(debug_compute_results){
|
| 363 |
cout << "FOR MUONS: " << endl;
|
| 364 |
cout << "Data \t " << mm_inside_data_pred << " +/- " << ee_inside_data_pred_error << "\t R = " << em_inside_data / em_outside_data << endl;
|
| 365 |
cout << "MC \t " << mm_inside_MC_pred << " +/- " << mm_inside_MC_pred_error << "\t R = " << em_inside_MC / em_outside_MC << endl;
|
| 366 |
|
| 367 |
cout << endl << endl;
|
| 368 |
cout << "the correct result would be " << httbareemm->Integral(httbareemm->FindBin(E),httbareemm->FindBin(F)) << endl;
|
| 369 |
}
|
| 370 |
|
| 371 |
results[nsel][0][0]=em_outside_MC; errors[nsel][0][0]=em_outside_MC_error;
|
| 372 |
results[nsel][0][1]=em_inside_MC; errors[nsel][0][1]=em_inside_MC_error;
|
| 373 |
results[nsel][0][2]=R_MC; errors[nsel][0][2]=R_MC_error;
|
| 374 |
|
| 375 |
results[nsel][1][0]=ee_outside_MC; errors[nsel][1][0]=ee_outside_MC_error;
|
| 376 |
results[nsel][1][1]=ee_inside_MC; errors[nsel][1][1]=ee_inside_MC_error;
|
| 377 |
results[nsel][1][2]=ee_inside_MC_pred; errors[nsel][1][2]=ee_inside_MC_pred_error;
|
| 378 |
|
| 379 |
results[nsel][2][0]=mm_outside_MC; errors[nsel][2][0]=mm_outside_MC_error;
|
| 380 |
results[nsel][2][1]=mm_inside_MC; errors[nsel][2][1]=mm_inside_MC_error;
|
| 381 |
results[nsel][2][2]=mm_inside_MC_pred; errors[nsel][2][2]=mm_inside_MC_pred_error;
|
| 382 |
|
| 383 |
*correctresult=httbareemm->Integral(httbareemm->FindBin(E),httbareemm->FindBin(F));
|
| 384 |
*/
|
| 385 |
}
|
| 386 |
|
| 387 |
|
| 388 |
void do_validation_step1(TH1F *eehisto,TH1F *mmhisto,TH1F *eemmhisto,TH1F *emhisto, TLine *oneline)
|
| 389 |
{
|
| 390 |
//validation step 1: plot the em, ee, mm, and eemm distributions to show that em:eemm ~ 1 and the shapes agree.
|
| 391 |
//cout << "validation step one" << endl;
|
| 392 |
|
| 393 |
//first of all: set up the canvas with two pads
|
| 394 |
TCanvas *methodval1 = new TCanvas("methodval1","Method Validation Canvas 1",600,600);
|
| 395 |
|
| 396 |
methodval1->Range(0,0,1,1);
|
| 397 |
methodval1->SetBorderSize(0);
|
| 398 |
methodval1->SetFrameFillColor(0);
|
| 399 |
TPad *method_validation_canvas1_pad1 = new TPad("method_validation_canvas1_pad1","method_validation_canvas1_pad1",0,0.33,1,1);
|
| 400 |
method_validation_canvas1_pad1->Draw();
|
| 401 |
method_validation_canvas1_pad1->cd();
|
| 402 |
method_validation_canvas1_pad1->Range(0,0,1,1);
|
| 403 |
method_validation_canvas1_pad1->SetFillColor(kWhite);
|
| 404 |
method_validation_canvas1_pad1->SetBorderSize(0);
|
| 405 |
method_validation_canvas1_pad1 ->SetFrameFillColor(0);
|
| 406 |
method_validation_canvas1_pad1 ->Modified();
|
| 407 |
methodval1->cd();
|
| 408 |
|
| 409 |
TPad *method_validation_canvas1_pad2 = new TPad("method_validation_canvas1_pad2", "method_validation_canvas1_pad2",0,0.0,1,0.33);
|
| 410 |
method_validation_canvas1_pad2->Draw();
|
| 411 |
method_validation_canvas1_pad2->cd();
|
| 412 |
method_validation_canvas1_pad2->Range(0,0,1,1);
|
| 413 |
method_validation_canvas1_pad2->SetFillColor(kWhite);
|
| 414 |
method_validation_canvas1_pad2->SetBorderSize(0);
|
| 415 |
method_validation_canvas1_pad2->SetFrameFillColor(0);
|
| 416 |
method_validation_canvas1_pad2->Modified();
|
| 417 |
methodval1->cd();
|
| 418 |
methodval1->Modified();
|
| 419 |
methodval1->cd();
|
| 420 |
|
| 421 |
//and now let's get plotting!
|
| 422 |
|
| 423 |
method_validation_canvas1_pad1->cd();
|
| 424 |
|
| 425 |
method_validation_canvas1_pad1->SetLogy();
|
| 426 |
|
| 427 |
eehisto->SetLineColor(kBlue);//BnoZ, ee
|
| 428 |
mmhisto->SetLineColor(kGreen);//BnoZ, mm
|
| 429 |
eemmhisto->SetLineColor(kBlack);//BnoZ, eemm
|
| 430 |
emhisto->SetLineColor(kRed);//BnoZ, em
|
| 431 |
|
| 432 |
emhisto->SetStats(0);
|
| 433 |
emhisto->GetXaxis()->SetTitle("m_{ll} (GeV)");
|
| 434 |
emhisto->GetXaxis()->CenterTitle();
|
| 435 |
emhisto->GetYaxis()->SetTitle("events");
|
| 436 |
emhisto->GetYaxis()->CenterTitle();
|
| 437 |
|
| 438 |
emhisto->Draw("histo");//BnoZ, em
|
| 439 |
eehisto->Draw("histo,same");//BnoZ, ee
|
| 440 |
mmhisto->Draw("histo,same");//BnoZ, mm
|
| 441 |
eemmhisto->Draw("histo,same");//BnoZ, eemm
|
| 442 |
|
| 443 |
TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
|
| 444 |
leg->AddEntry(emhisto,"e#mu, BnoZ (MC)","L");
|
| 445 |
leg->AddEntry(eemmhisto,"ee or#mu#mu, BnoZ (MC)","L");
|
| 446 |
leg->AddEntry(eehisto,"ee, BnoZ (MC)","L");
|
| 447 |
leg->AddEntry(mmhisto,"#mu#mu, BnoZ (MC)","L");
|
| 448 |
leg->SetFillColor(kWhite);
|
| 449 |
leg->SetLineColor(kWhite);
|
| 450 |
leg->SetBorderSize(0);
|
| 451 |
leg->Draw();
|
| 452 |
|
| 453 |
method_validation_canvas1_pad2->cd();
|
| 454 |
|
| 455 |
TH1F *ratio_one;
|
| 456 |
ratio_one = create_ratio_plot(emhisto,eemmhisto);
|
| 457 |
ratio_one->GetYaxis()->SetRangeUser(0,4);
|
| 458 |
ratio_one->SetFillColor(8);
|
| 459 |
ratio_one->GetYaxis()->SetTitle("ratio ee#mu#mu/e#mu");
|
| 460 |
ratio_one->GetYaxis()->CenterTitle();
|
| 461 |
ratio_one->Draw("e5");
|
| 462 |
ratio_one->SetMarkerSize(0);
|
| 463 |
ratio_one->GetXaxis()->SetTitle("");
|
| 464 |
oneline->Draw();
|
| 465 |
|
| 466 |
methodval1->SaveAs("em_method_output/method_validation_canvas1.png");
|
| 467 |
methodval1->SaveAs("em_method_output/method_validation_canvas1.C");
|
| 468 |
|
| 469 |
method_validation_canvas1_pad1->SaveAs("em_method_output/method_validation_canvas1_pad1.png");
|
| 470 |
method_validation_canvas1_pad2->SaveAs("em_method_output/method_validation_canvas1_pad2.png");
|
| 471 |
method_validation_canvas1_pad1->SaveAs("em_method_output/method_validation_canvas1_pad1.C");
|
| 472 |
method_validation_canvas1_pad2->SaveAs("em_method_output/method_validation_canvas1_pad2.C");
|
| 473 |
|
| 474 |
}//end of do_validation_step2
|
| 475 |
|
| 476 |
void do_validation_step2(TH1F *eehisto,TH1F *mmhisto,TH1F *emhisto,TLine *oneline, float A, float B, float C, float D, float E, float F)
|
| 477 |
{
|
| 478 |
//validation step 2: plot ee BnoZ vs. em BnoZ on one pad (with a ratio pad) and mm BnoZ vs. em BnoZ on a second (incl. a ratio pad)
|
| 479 |
TCanvas *methodval2 = new TCanvas("methodval2","Method Validation Canvas 2",1200,600);
|
| 480 |
methodval2->Range(0,0,1,1);
|
| 481 |
methodval2->SetBorderSize(0);
|
| 482 |
methodval2->SetFrameFillColor(0);
|
| 483 |
TPad *method_validation_canvas2_pad1 = new TPad("method_validation_canvas2_pad1","method_validation_canvas2_pad1",0,0.33,0.5,1);
|
| 484 |
method_validation_canvas2_pad1->Draw();
|
| 485 |
method_validation_canvas2_pad1->cd();
|
| 486 |
method_validation_canvas2_pad1->Range(0,0,1,1);
|
| 487 |
method_validation_canvas2_pad1->SetFillColor(kWhite);
|
| 488 |
method_validation_canvas2_pad1->SetBorderSize(0);
|
| 489 |
method_validation_canvas2_pad1 ->SetFrameFillColor(0);
|
| 490 |
method_validation_canvas2_pad1 ->Modified();
|
| 491 |
methodval2->cd();
|
| 492 |
|
| 493 |
TPad *method_validation_canvas2_pad2 = new TPad("method_validation_canvas2_pad2", "method_validation_canvas2_pad2",0,0.0,0.5,0.33);
|
| 494 |
method_validation_canvas2_pad2->Draw();
|
| 495 |
method_validation_canvas2_pad2->cd();
|
| 496 |
method_validation_canvas2_pad2->Range(0,0,1,1);
|
| 497 |
method_validation_canvas2_pad2->SetFillColor(kWhite);
|
| 498 |
method_validation_canvas2_pad2->SetBorderSize(0);
|
| 499 |
method_validation_canvas2_pad2->SetFrameFillColor(0);
|
| 500 |
method_validation_canvas2_pad2->Modified();
|
| 501 |
methodval2->cd();
|
| 502 |
methodval2->Modified();
|
| 503 |
methodval2->cd();
|
| 504 |
|
| 505 |
TPad *method_validation_canvas2_pad3 = new TPad("method_validation_canvas2_pad3","method_validation_canvas2_pad3",0.5,0.33,1,1);
|
| 506 |
method_validation_canvas2_pad3->Draw();
|
| 507 |
method_validation_canvas2_pad3->cd();
|
| 508 |
method_validation_canvas2_pad3->Range(0,0,1,1);
|
| 509 |
method_validation_canvas2_pad3->SetFillColor(kWhite);
|
| 510 |
method_validation_canvas2_pad3->SetBorderSize(0);
|
| 511 |
method_validation_canvas2_pad3 ->SetFrameFillColor(0);
|
| 512 |
method_validation_canvas2_pad3 ->Modified();
|
| 513 |
methodval2->cd();
|
| 514 |
|
| 515 |
TPad *method_validation_canvas2_pad4 = new TPad("method_validation_canvas2_pad4", "method_validation_canvas2_pad4",0.5,0.0,1,0.33);
|
| 516 |
method_validation_canvas2_pad4->Draw();
|
| 517 |
method_validation_canvas2_pad4->cd();
|
| 518 |
method_validation_canvas2_pad4->Range(0,0,1,1);
|
| 519 |
method_validation_canvas2_pad4->SetFillColor(kWhite);
|
| 520 |
method_validation_canvas2_pad4->SetBorderSize(0);
|
| 521 |
method_validation_canvas2_pad4->SetFrameFillColor(0);
|
| 522 |
method_validation_canvas2_pad4->Modified();
|
| 523 |
methodval2->cd();
|
| 524 |
methodval2->Modified();
|
| 525 |
methodval2->cd();
|
| 526 |
|
| 527 |
method_validation_canvas2_pad1->cd();
|
| 528 |
TH1F *eebnoznormalized = (TH1F*) eehisto->Clone("eebnoznormalized");
|
| 529 |
TH1F *embnoznormalized = (TH1F*) emhisto->Clone("embnoznormalized");
|
| 530 |
TH1F *mmbnoznormalized = (TH1F*) mmhisto->Clone("mmbnoznormalized");
|
| 531 |
|
| 532 |
// eebnoznormalized->Scale(1/(eebnoznormalized->Integral()));
|
| 533 |
// embnoznormalized->Scale(0.5);
|
| 534 |
// mmbnoznormalized->Scale(1/(mmbnoznormalized->Integral()));
|
| 535 |
|
| 536 |
eebnoznormalized->SetLineColor(kBlue);
|
| 537 |
mmbnoznormalized->SetLineColor(kBlue);
|
| 538 |
embnoznormalized->SetLineColor(kBlack);
|
| 539 |
embnoznormalized->SetLineStyle(3);
|
| 540 |
|
| 541 |
method_validation_canvas2_pad1->cd();
|
| 542 |
eebnoznormalized->GetYaxis()->SetRangeUser(0.8*(eebnoznormalized->GetMinimum()),1.1*(eebnoznormalized->GetMaximum()));
|
| 543 |
embnoznormalized->GetYaxis()->SetRangeUser(0.8*(embnoznormalized->GetMinimum()),1.1*(embnoznormalized->GetMaximum()));
|
| 544 |
mmbnoznormalized->GetYaxis()->SetRangeUser(0.8*(mmbnoznormalized->GetMinimum()),1.1*(mmbnoznormalized->GetMaximum()));
|
| 545 |
|
| 546 |
eebnoznormalized->Draw("histo");
|
| 547 |
embnoznormalized->Draw("histo,same");
|
| 548 |
|
| 549 |
|
| 550 |
//**************************************************************************************** first set
|
| 551 |
|
| 552 |
float dia1_min=0;
|
| 553 |
if(embnoznormalized->GetMinimum()<eebnoznormalized->GetMinimum()){
|
| 554 |
dia1_min=embnoznormalized->GetMinimum();
|
| 555 |
}else{
|
| 556 |
dia1_min=eebnoznormalized->GetMinimum();
|
| 557 |
}
|
| 558 |
method_validation_canvas2_pad1->cd()->SetLogy(1);//log - kill this!
|
| 559 |
float dia1_max=0;
|
| 560 |
if(embnoznormalized->GetMaximum()>eebnoznormalized->GetMaximum()){
|
| 561 |
dia1_max=embnoznormalized->GetMaximum();
|
| 562 |
}else{
|
| 563 |
dia1_max=eebnoznormalized->GetMaximum();
|
| 564 |
}
|
| 565 |
|
| 566 |
TBox *box1 = new TBox(A,dia1_min,B,dia1_max);
|
| 567 |
box1->SetFillColor(kAzure-9);
|
| 568 |
box1->Draw();
|
| 569 |
|
| 570 |
TBox *box2 = new TBox(C,dia1_min,D,dia1_max);
|
| 571 |
box2->SetFillColor(kAzure-9);
|
| 572 |
box2->Draw();
|
| 573 |
|
| 574 |
TBox *box3 = new TBox(E,dia1_min,F,dia1_max);
|
| 575 |
box3->SetFillColor(kGreen-9);
|
| 576 |
box3->Draw();
|
| 577 |
|
| 578 |
embnoznormalized->Draw("histo,same");
|
| 579 |
eebnoznormalized->Draw("histo,same");
|
| 580 |
method_validation_canvas2_pad1->RedrawAxis();
|
| 581 |
|
| 582 |
TText *atext_sr = new TText((E+F)/2,(dia1_max-dia1_min)*0.05+dia1_min,"SR");
|
| 583 |
atext_sr->SetTextColor(kGreen+2);
|
| 584 |
atext_sr->SetTextAlign(21);
|
| 585 |
atext_sr->Draw();
|
| 586 |
|
| 587 |
TText *atext_sb1 = new TText((A+B)/2,(dia1_max-dia1_min)*0.05+dia1_min,"SB");
|
| 588 |
atext_sb1->SetTextColor(kBlue-4);
|
| 589 |
atext_sb1->SetTextAlign(21);
|
| 590 |
atext_sb1->Draw();
|
| 591 |
|
| 592 |
TText *atext_sb2 = new TText((C+D)/2,(dia1_max-dia1_min)*0.05+dia1_min,"SB");
|
| 593 |
atext_sb2->SetTextColor(kBlue-4);
|
| 594 |
atext_sb2->SetTextAlign(21);
|
| 595 |
atext_sb2->Draw();
|
| 596 |
|
| 597 |
TLegend *leg2 = new TLegend(0.7,0.7,0.89,0.89);
|
| 598 |
leg2->AddEntry(eebnoznormalized,"e^{+}e^{-} (MC)","L");
|
| 599 |
leg2->AddEntry(embnoznormalized,"0.5 x e^{#pm}#mu^{#bar{+}} (MC)","L");
|
| 600 |
leg2->SetFillColor(kAzure-9);
|
| 601 |
leg2->SetLineColor(kAzure-9);
|
| 602 |
leg2->SetBorderSize(0);
|
| 603 |
leg2->Draw();
|
| 604 |
|
| 605 |
method_validation_canvas2_pad2->cd();
|
| 606 |
TH1F *ratio_two;
|
| 607 |
ratio_two = create_ratio_plot(eebnoznormalized,embnoznormalized);
|
| 608 |
ratio_two->GetYaxis()->SetRangeUser(0,4);
|
| 609 |
ratio_two->GetYaxis()->SetTitle("ratio ee/e#mu");
|
| 610 |
ratio_two->GetYaxis()->CenterTitle();
|
| 611 |
ratio_two->SetFillColor(8);
|
| 612 |
ratio_two->SetMarkerSize(0);
|
| 613 |
|
| 614 |
ratio_two->Draw("e5");
|
| 615 |
oneline->Draw();
|
| 616 |
|
| 617 |
//**************************************************************************************** second set
|
| 618 |
method_validation_canvas2_pad3->cd();
|
| 619 |
mmbnoznormalized->Draw("histo");
|
| 620 |
embnoznormalized->Draw("histo,same");
|
| 621 |
|
| 622 |
|
| 623 |
float dia2_min=0;
|
| 624 |
if(embnoznormalized->GetMinimum()<eebnoznormalized->GetMinimum()){
|
| 625 |
dia2_min=embnoznormalized->GetMinimum();
|
| 626 |
}else{
|
| 627 |
dia2_min=mmbnoznormalized->GetMinimum();
|
| 628 |
}
|
| 629 |
|
| 630 |
float dia2_max=0;
|
| 631 |
if(embnoznormalized->GetMaximum()>eebnoznormalized->GetMaximum()){
|
| 632 |
dia2_max=embnoznormalized->GetMaximum();
|
| 633 |
}else{
|
| 634 |
dia2_max=mmbnoznormalized->GetMaximum();
|
| 635 |
}
|
| 636 |
|
| 637 |
method_validation_canvas2_pad3->cd()->SetLogy(1); //log kill this!
|
| 638 |
|
| 639 |
TBox *box4 = new TBox(A,dia2_min,B,dia2_max);
|
| 640 |
box4->SetFillColor(kAzure-9);
|
| 641 |
box4->Draw();
|
| 642 |
|
| 643 |
TBox *box5 = new TBox(C,dia2_min,D,dia2_max);
|
| 644 |
box5->SetFillColor(kAzure-9);
|
| 645 |
box5->Draw();
|
| 646 |
|
| 647 |
TBox *box6 = new TBox(E,dia2_min,F,dia2_max);
|
| 648 |
box6->SetFillColor(kGreen-9);
|
| 649 |
box6->Draw();
|
| 650 |
|
| 651 |
embnoznormalized->Draw("histo,same");
|
| 652 |
mmbnoznormalized->Draw("histo,same");
|
| 653 |
method_validation_canvas2_pad3->RedrawAxis();
|
| 654 |
|
| 655 |
TText *text_sr = new TText((E+F)/2,(dia2_max-dia2_min)*0.05+dia2_min,"SR");
|
| 656 |
text_sr->SetTextColor(kGreen+2);
|
| 657 |
text_sr->SetTextAlign(21);
|
| 658 |
text_sr->Draw();
|
| 659 |
|
| 660 |
TText *text_sb1 = new TText((B+A)/2,(dia2_max-dia2_min)*0.05+dia2_min,"SB");
|
| 661 |
text_sb1->SetTextColor(kBlue-4);
|
| 662 |
text_sb1->SetTextAlign(21);
|
| 663 |
text_sb1->Draw();
|
| 664 |
|
| 665 |
TText *text_sb2 = new TText(151,(dia2_max-dia2_min)*0.05+dia2_min,"SB");
|
| 666 |
text_sb2->SetTextColor(kBlue-4);
|
| 667 |
text_sb2->SetTextAlign(21);
|
| 668 |
text_sb2->Draw();
|
| 669 |
|
| 670 |
TLegend *leg3 = new TLegend(0.7,0.7,0.89,0.89);
|
| 671 |
leg3->AddEntry(mmbnoznormalized,"#mu^{+}#mu^{-} (MC)","L");
|
| 672 |
leg3->AddEntry(embnoznormalized,"0.5 x e^{#pm}#mu^{#bar{+}} (MC)","L");
|
| 673 |
leg3->SetFillColor(kAzure-9);
|
| 674 |
leg3->SetLineColor(kAzure-9);
|
| 675 |
leg3->SetBorderSize(0);
|
| 676 |
leg3->Draw();
|
| 677 |
|
| 678 |
method_validation_canvas2_pad4->cd();
|
| 679 |
TH1F *ratio_three;
|
| 680 |
ratio_three = create_ratio_plot(mmbnoznormalized,embnoznormalized);
|
| 681 |
ratio_three->GetYaxis()->SetRangeUser(0,4);
|
| 682 |
ratio_three->GetYaxis()->SetTitle("ratio #mu#mu/e#mu");
|
| 683 |
ratio_three->GetYaxis()->CenterTitle();
|
| 684 |
ratio_three->SetFillColor(8);
|
| 685 |
ratio_three->SetMarkerSize(0);
|
| 686 |
|
| 687 |
ratio_three->Draw("e5");
|
| 688 |
oneline->Draw();
|
| 689 |
methodval2->SaveAs("em_method_output/method_validation_canvas2.png");
|
| 690 |
methodval2->SaveAs("em_method_output/method_validation_canvas2.C");
|
| 691 |
method_validation_canvas2_pad1->SaveAs("em_method_output/method_validation_canvas2_pad1.png");
|
| 692 |
method_validation_canvas2_pad2->SaveAs("em_method_output/method_validation_canvas2_pad2.png");
|
| 693 |
method_validation_canvas2_pad3->SaveAs("em_method_output/method_validation_canvas2_pad3.png");
|
| 694 |
method_validation_canvas2_pad4->SaveAs("em_method_output/method_validation_canvas2_pad4.png");
|
| 695 |
|
| 696 |
method_validation_canvas2_pad1->SaveAs("em_method_output/method_validation_canvas2_pad1.eps");
|
| 697 |
method_validation_canvas2_pad2->SaveAs("em_method_output/method_validation_canvas2_pad2.eps");
|
| 698 |
method_validation_canvas2_pad3->SaveAs("em_method_output/method_validation_canvas2_pad3.eps");
|
| 699 |
method_validation_canvas2_pad4->SaveAs("em_method_output/method_validation_canvas2_pad4.eps");
|
| 700 |
|
| 701 |
method_validation_canvas2_pad1->SaveAs("em_method_output/method_validation_canvas2_pad1.C");
|
| 702 |
method_validation_canvas2_pad2->SaveAs("em_method_output/method_validation_canvas2_pad2.C");
|
| 703 |
method_validation_canvas2_pad3->SaveAs("em_method_output/method_validation_canvas2_pad3.C");
|
| 704 |
method_validation_canvas2_pad4->SaveAs("em_method_output/method_validation_canvas2_pad4.C");
|
| 705 |
|
| 706 |
}
|
| 707 |
|
| 708 |
void do_validation_step3(TH1F *ee_sb,TH1F *ee_sr,TH1F *mm_sb,TH1F *mm_sr,TLine *oneline)
|
| 709 |
{
|
| 710 |
TCanvas *methodval3 = new TCanvas("methodval3","Method Validation Canvas 3",1200,600);
|
| 711 |
methodval3->Range(0,0,1,1);
|
| 712 |
methodval3->SetBorderSize(0);
|
| 713 |
methodval3->SetFrameFillColor(0);
|
| 714 |
TPad *method_validation_canvas3_pad1 = new TPad("method_validation_canvas3_pad1","method_validation_canvas3_pad1",0,0.33,0.5,1);
|
| 715 |
method_validation_canvas3_pad1->Draw();
|
| 716 |
method_validation_canvas3_pad1->cd();
|
| 717 |
method_validation_canvas3_pad1->Range(0,0,1,1);
|
| 718 |
method_validation_canvas3_pad1->SetFillColor(kWhite);
|
| 719 |
method_validation_canvas3_pad1->SetBorderSize(0);
|
| 720 |
method_validation_canvas3_pad1 ->SetFrameFillColor(0);
|
| 721 |
method_validation_canvas3_pad1 ->Modified();
|
| 722 |
methodval3->cd();
|
| 723 |
|
| 724 |
TPad *method_validation_canvas3_pad2 = new TPad("method_validation_canvas3_pad2", "method_validation_canvas3_pad2",0,0.0,0.5,0.33);
|
| 725 |
method_validation_canvas3_pad2->Draw();
|
| 726 |
method_validation_canvas3_pad2->cd();
|
| 727 |
method_validation_canvas3_pad2->Range(0,0,1,1);
|
| 728 |
method_validation_canvas3_pad2->SetFillColor(kWhite);
|
| 729 |
method_validation_canvas3_pad2->SetBorderSize(0);
|
| 730 |
method_validation_canvas3_pad2->SetFrameFillColor(0);
|
| 731 |
method_validation_canvas3_pad2->Modified();
|
| 732 |
methodval3->cd();
|
| 733 |
methodval3->Modified();
|
| 734 |
methodval3->cd();
|
| 735 |
|
| 736 |
TPad *method_validation_canvas3_pad3 = new TPad("method_validation_canvas3_pad3","method_validation_canvas3_pad3",0.5,0.33,1,1);
|
| 737 |
method_validation_canvas3_pad3->Draw();
|
| 738 |
method_validation_canvas3_pad3->cd();
|
| 739 |
method_validation_canvas3_pad3->Range(0,0,1,1);
|
| 740 |
method_validation_canvas3_pad3->SetFillColor(kWhite);
|
| 741 |
method_validation_canvas3_pad3->SetBorderSize(0);
|
| 742 |
method_validation_canvas3_pad3 ->SetFrameFillColor(0);
|
| 743 |
method_validation_canvas3_pad3 ->Modified();
|
| 744 |
methodval3->cd();
|
| 745 |
|
| 746 |
TPad *method_validation_canvas3_pad4 = new TPad("method_validation_canvas3_pad4", "method_validation_canvas3_pad4",0.5,0.0,1,0.33);
|
| 747 |
method_validation_canvas3_pad4->Draw();
|
| 748 |
method_validation_canvas3_pad4->cd();
|
| 749 |
method_validation_canvas3_pad4->Range(0,0,1,1);
|
| 750 |
method_validation_canvas3_pad4->SetFillColor(kWhite);
|
| 751 |
method_validation_canvas3_pad4->SetBorderSize(0);
|
| 752 |
method_validation_canvas3_pad4->SetFrameFillColor(0);
|
| 753 |
method_validation_canvas3_pad4->Modified();
|
| 754 |
methodval3->cd();
|
| 755 |
methodval3->Modified();
|
| 756 |
methodval3->cd();
|
| 757 |
|
| 758 |
method_validation_canvas3_pad1->cd();
|
| 759 |
TH1F *eesb = (TH1F*) ee_sb->Clone("eesb");
|
| 760 |
TH1F *eesr = (TH1F*) ee_sr->Clone("eesr");
|
| 761 |
|
| 762 |
TH1F *mmsb = (TH1F*) mm_sb->Clone("mmsb");
|
| 763 |
TH1F *mmsr = (TH1F*) mm_sr->Clone("mmsr");
|
| 764 |
|
| 765 |
eesb->Scale(1/(eesb->Integral()));
|
| 766 |
eesr->Scale(1/(eesr->Integral()));
|
| 767 |
|
| 768 |
mmsb->Scale(1/(mmsb->Integral()));
|
| 769 |
mmsr->Scale(1/(mmsr->Integral()));
|
| 770 |
|
| 771 |
eesb->SetLineColor(kBlue);
|
| 772 |
eesr->SetLineColor(kBlack);
|
| 773 |
eesr->SetLineStyle(3);
|
| 774 |
|
| 775 |
mmsb->SetLineColor(kBlue);
|
| 776 |
mmsr->SetLineColor(kBlack);
|
| 777 |
mmsr->SetLineStyle(3);
|
| 778 |
|
| 779 |
|
| 780 |
method_validation_canvas3_pad1->cd();
|
| 781 |
if (eesr->GetMaximum()>eesb->GetMaximum()) {
|
| 782 |
eesr->Draw("histo");
|
| 783 |
eesb->Draw("histo,same");
|
| 784 |
}
|
| 785 |
else
|
| 786 |
{
|
| 787 |
eesb->Draw("histo");
|
| 788 |
eesr->Draw("histo,same");
|
| 789 |
}
|
| 790 |
TLegend *leg3a = new TLegend(0.6,0.6,0.89,0.89);
|
| 791 |
leg3a->AddEntry(eesb,"ee, Sidebands (MC)","L");
|
| 792 |
leg3a->AddEntry(eesr,"ee, Signal Region (MC)","L");
|
| 793 |
leg3a->SetFillColor(kWhite);
|
| 794 |
leg3a->SetLineColor(kWhite);
|
| 795 |
leg3a->SetBorderSize(0);
|
| 796 |
leg3a->Draw();
|
| 797 |
|
| 798 |
method_validation_canvas3_pad2->cd();
|
| 799 |
TH1F *ratio_3a;
|
| 800 |
ratio_3a = create_ratio_plot(eesb,eesr);
|
| 801 |
ratio_3a->GetYaxis()->SetRangeUser(0,4);
|
| 802 |
ratio_3a->GetYaxis()->SetTitle("ratio ee SB/SR");
|
| 803 |
ratio_3a->GetYaxis()->CenterTitle();
|
| 804 |
ratio_3a->SetFillColor(8);
|
| 805 |
ratio_3a->SetMarkerSize(0);
|
| 806 |
ratio_3a->Draw("e5");
|
| 807 |
|
| 808 |
TLine *twoline = new TLine(51,1,200,1);
|
| 809 |
twoline->SetLineColor(kBlue);
|
| 810 |
twoline->SetLineStyle(2);
|
| 811 |
twoline->Draw();
|
| 812 |
//oneline->Draw();
|
| 813 |
|
| 814 |
|
| 815 |
method_validation_canvas3_pad3->cd();
|
| 816 |
if (mmsr->GetMaximum()>mmsb->GetMaximum())
|
| 817 |
{
|
| 818 |
mmsr->Draw("histo");
|
| 819 |
mmsb->Draw("histo,same");
|
| 820 |
}
|
| 821 |
else
|
| 822 |
{
|
| 823 |
mmsb->Draw("histo");
|
| 824 |
mmsr->Draw("histo,same");
|
| 825 |
}
|
| 826 |
TLegend *leg3b = new TLegend(0.6,0.6,0.89,0.89);
|
| 827 |
leg3b->AddEntry(mmsb,"#mu#mu, Sidebands (MC)","L");
|
| 828 |
leg3b->AddEntry(mmsr,"#mu#mu, Signal Region (MC)","L");
|
| 829 |
leg3b->SetFillColor(kWhite);
|
| 830 |
leg3b->SetLineColor(kWhite);
|
| 831 |
leg3b->SetBorderSize(0);
|
| 832 |
leg3b->Draw();
|
| 833 |
|
| 834 |
method_validation_canvas3_pad4->cd();
|
| 835 |
TH1F *ratio_3b;
|
| 836 |
ratio_3b = create_ratio_plot(mmsb,mmsr);
|
| 837 |
ratio_3b->GetYaxis()->SetRangeUser(0,4);
|
| 838 |
ratio_3b->GetYaxis()->SetTitle("ratio #mu#mu SB/SR");
|
| 839 |
ratio_3b->GetYaxis()->CenterTitle();
|
| 840 |
ratio_3b->SetFillColor(8);
|
| 841 |
ratio_3b->SetMarkerSize(0);
|
| 842 |
ratio_3b->Draw("e5");
|
| 843 |
// oneline->Draw();
|
| 844 |
twoline->Draw();
|
| 845 |
|
| 846 |
methodval3->SaveAs("em_method_output/method_validation_canvas3.png");
|
| 847 |
methodval3->SaveAs("em_method_output/method_validation_canvas3.C");
|
| 848 |
|
| 849 |
method_validation_canvas3_pad1->SaveAs("em_method_output/method_validation_canvas3_pad1.png");
|
| 850 |
method_validation_canvas3_pad2->SaveAs("em_method_output/method_validation_canvas3_pad2.png");
|
| 851 |
method_validation_canvas3_pad3->SaveAs("em_method_output/method_validation_canvas3_pad3.png");
|
| 852 |
method_validation_canvas3_pad4->SaveAs("em_method_output/method_validation_canvas3_pad4.png");
|
| 853 |
|
| 854 |
method_validation_canvas3_pad1->SaveAs("em_method_output/method_validation_canvas3_pad1.C");
|
| 855 |
method_validation_canvas3_pad2->SaveAs("em_method_output/method_validation_canvas3_pad2.C");
|
| 856 |
method_validation_canvas3_pad3->SaveAs("em_method_output/method_validation_canvas3_pad3.C");
|
| 857 |
method_validation_canvas3_pad4->SaveAs("em_method_output/method_validation_canvas3_pad4.C");
|
| 858 |
|
| 859 |
}
|
| 860 |
|
| 861 |
string get_selection( int i)
|
| 862 |
{
|
| 863 |
string selection="";
|
| 864 |
if (i==0) selection="(id1==id2)&&(id1==0)";//ee
|
| 865 |
if (i==1) selection="(id1==id2)&&(id1==1)";//mm
|
| 866 |
if (i==2) selection="(id1==id2)";//eemm
|
| 867 |
if (i==3) selection="(id1!=id2)";//em
|
| 868 |
return selection;
|
| 869 |
}
|
| 870 |
|
| 871 |
string get_area(int i, float A, float B, float C, float D, float E, float F)
|
| 872 |
{
|
| 873 |
|
| 874 |
if (i==0) return "";//all
|
| 875 |
stringstream selection;
|
| 876 |
if (i==1)//outside (i.e. sidebands)
|
| 877 |
{
|
| 878 |
selection << "((mll>" << A << "&&mll<" << B << ")||(mll>" << C << "&&mll<" << D << "))";
|
| 879 |
return selection.str();
|
| 880 |
}
|
| 881 |
if (i==2) //inside
|
| 882 |
{
|
| 883 |
selection << "(mll>" << E << "&&mll<" << F << ")";
|
| 884 |
return selection.str();
|
| 885 |
}
|
| 886 |
return "";
|
| 887 |
}
|
| 888 |
|
| 889 |
string get_jzb(int i, float peak_pos)
|
| 890 |
{
|
| 891 |
if (i==0) return "";
|
| 892 |
if (i==1)
|
| 893 |
{
|
| 894 |
char selection[1000];
|
| 895 |
sprintf(selection,"jzb[1]>%f",peak_pos);
|
| 896 |
return selection;
|
| 897 |
}
|
| 898 |
return "";
|
| 899 |
}
|
| 900 |
|
| 901 |
|
| 902 |
//******************************
|
| 903 |
|
| 904 |
string get_selection_name( int i)
|
| 905 |
{
|
| 906 |
string selection="";
|
| 907 |
if (i==0) return "ee";
|
| 908 |
if (i==1) return "mm";
|
| 909 |
if (i==2) return "eemm";
|
| 910 |
if (i==3) return "em";
|
| 911 |
return "get_selection_name_error";
|
| 912 |
}
|
| 913 |
|
| 914 |
string get_area_name(int i)
|
| 915 |
{
|
| 916 |
|
| 917 |
if (i==0) return "all";
|
| 918 |
if (i==1) return "sidebands";
|
| 919 |
if (i==2) return "signalregion";
|
| 920 |
return "get_area_name_error";
|
| 921 |
}
|
| 922 |
|
| 923 |
string get_jzb_name(int i)
|
| 924 |
{
|
| 925 |
if (i==0) return "no_jzb_condition";
|
| 926 |
if (i==1) return "with_jzb_condition";
|
| 927 |
return "get_jzb_name_error";
|
| 928 |
}
|
| 929 |
|
| 930 |
string get_var_name(int i)
|
| 931 |
{
|
| 932 |
if (i==0) return "met4";
|
| 933 |
if (i==1) return "mll";
|
| 934 |
return "get_var_name_error";
|
| 935 |
}
|
| 936 |
|
| 937 |
string get_sample_name(int i)
|
| 938 |
{
|
| 939 |
if (i==0) return "data";
|
| 940 |
if (i==1) return "zjets";
|
| 941 |
if (i==2) return "ttbar";
|
| 942 |
if (i==3) return "astar";
|
| 943 |
if (i==4) return "znunu";
|
| 944 |
if (i==5) return "wjets";
|
| 945 |
if (i==6) return "vvjets";
|
| 946 |
if (i==7) return "singletop";
|
| 947 |
if (i==8) return "qcd1";
|
| 948 |
if (i==9) return "qcd2";
|
| 949 |
if (i==10) return "qcd3";
|
| 950 |
if (i==11) return "qcd4";
|
| 951 |
if (i==12) return "lm4";
|
| 952 |
return "get_sample_name_error";
|
| 953 |
}
|
| 954 |
|
| 955 |
//*************************************
|
| 956 |
float get_weight(int isample)
|
| 957 |
{
|
| 958 |
if (isample==0) return 1.0;
|
| 959 |
|
| 960 |
float lumi=34.3; // normalize MC to lumi
|
| 961 |
|
| 962 |
float ZJetsCrossSection = 3048.0; //NNLO
|
| 963 |
float TTbarCrossSection = 165.0; // approx. NNLO
|
| 964 |
float WJetsCrossSection = 3.131e4; //NNLO
|
| 965 |
float ZnunuCrossSection = 4.5e+3; //(LO);
|
| 966 |
float SingleTopCrossSection = 64.6; // NLO;
|
| 967 |
float VVJetsCrossSection = 4.8; // LO;
|
| 968 |
float LM4CrossSection = 2.537; // k*LO
|
| 969 |
float QCD100to250CrossSection = 7e+6 ; // LO
|
| 970 |
float QCD250to500CrossSection = 1.71e+5 ; // LO
|
| 971 |
float QCD500to1000CrossSection = 5.2e+3 ; // LO
|
| 972 |
float QCD1000toInfCrossSection = 8.3e+1 ; // LO
|
| 973 |
|
| 974 |
float totEventsZjets = 1.07407e+06;
|
| 975 |
float totEventsTTbar = 1.4389e+06;
|
| 976 |
float totEventsWJets = 9867519.0;
|
| 977 |
float totEventsZnunu = 1530040.0;
|
| 978 |
float totEventsVVJets = 96685.0;
|
| 979 |
float totEventsSingleTop = 528593.0;
|
| 980 |
float totEventsQCD100to250 = 10820757.0;
|
| 981 |
float totEventsQCD250to500 = 4642821.;
|
| 982 |
float totEventsQCD500to1000 = 3152612.;
|
| 983 |
float totEventsQCD1000toInf = 479830.;
|
| 984 |
|
| 985 |
|
| 986 |
float totEventsLM4=10e15;//to make it zero
|
| 987 |
float totEventsAstar=189069.;//from before
|
| 988 |
float AstarCrossSection=310.0;//from before
|
| 989 |
|
| 990 |
|
| 991 |
float lm4_weight=(lumi/totEventsLM4)*LM4CrossSection;
|
| 992 |
float astar_weight=(lumi/totEventsAstar)*AstarCrossSection;
|
| 993 |
float znunu_weight=(lumi/totEventsZnunu)*ZnunuCrossSection;
|
| 994 |
float zjets_weight=(lumi/totEventsZjets)*ZJetsCrossSection;
|
| 995 |
float ttbar_weight=(lumi/totEventsTTbar)*TTbarCrossSection;
|
| 996 |
float wjet_weight=(lumi/totEventsWJets)*WJetsCrossSection;
|
| 997 |
float vvjets_weight=(lumi/totEventsVVJets)*VVJetsCrossSection;
|
| 998 |
float singletop_weight=(lumi/totEventsSingleTop)*SingleTopCrossSection;
|
| 999 |
float qcd1000toinf_weight=(lumi/totEventsQCD1000toInf)*QCD1000toInfCrossSection;
|
| 1000 |
float qcd500to1000_weight=(lumi/totEventsQCD500to1000)*QCD500to1000CrossSection;
|
| 1001 |
float qcd250to500_weight=(lumi/totEventsQCD250to500)*QCD250to500CrossSection;
|
| 1002 |
float qcd100to250_weight=(lumi/totEventsQCD100to250)*QCD100to250CrossSection;
|
| 1003 |
// float data_weight=1.0;
|
| 1004 |
|
| 1005 |
if (isample==1) return zjets_weight;
|
| 1006 |
if (isample==2) return ttbar_weight;
|
| 1007 |
if (isample==3) return astar_weight;
|
| 1008 |
if (isample==4) return znunu_weight;
|
| 1009 |
if (isample==5) return wjet_weight;
|
| 1010 |
if (isample==6) return vvjets_weight;
|
| 1011 |
if (isample==7) return singletop_weight;
|
| 1012 |
if (isample==8) return qcd1000toinf_weight;
|
| 1013 |
if (isample==9) return qcd500to1000_weight;
|
| 1014 |
if (isample==10) return qcd250to500_weight;
|
| 1015 |
if (isample==11) return qcd100to250_weight;
|
| 1016 |
if (isample==12) return lm4_weight;
|
| 1017 |
|
| 1018 |
return 0.0; // if there is no match, we return 0.
|
| 1019 |
}
|
| 1020 |
|
| 1021 |
string get_histo_name(int isample,int ivar,int iarea,int isel,int ijzb)
|
| 1022 |
{
|
| 1023 |
stringstream histoname;
|
| 1024 |
histoname << get_sample_name(isample) << "_" << get_var_name(ivar) << "_" << get_area_name(iarea) << "_" << get_selection_name(isel) << "_" << get_jzb_name(ijzb);
|
| 1025 |
// cout << "get_histo_name has been called with " << isample << ivar << iarea << isel << ijzb << " --> " << histoname.str() << endl;
|
| 1026 |
return histoname.str();
|
| 1027 |
}
|
| 1028 |
|
| 1029 |
string get_sum_sample_name(int isample)
|
| 1030 |
{
|
| 1031 |
if (isample==0) return "all";
|
| 1032 |
if (isample==1) return "BwithZ";
|
| 1033 |
if (isample==2) return "BnoZ";
|
| 1034 |
if (isample==3) return "Data";
|
| 1035 |
return "get_sum_sample_name_error";
|
| 1036 |
}
|
| 1037 |
|
| 1038 |
string get_sum_histo_name(int isample,int ivar,int iarea,int isel,int ijzb)
|
| 1039 |
{
|
| 1040 |
stringstream histoname;
|
| 1041 |
histoname << "SUM_" << get_sum_sample_name(isample) << "_" << get_var_name(ivar) << "_" << get_area_name(iarea) << "_" << get_selection_name(isel) << "_" << get_jzb_name(ijzb);
|
| 1042 |
// cout << "get_histo_name has been called with " << isample << ivar << iarea << isel << ijzb << " --> " << histoname.str() << endl;
|
| 1043 |
return histoname.str();
|
| 1044 |
}
|
| 1045 |
|
| 1046 |
|
| 1047 |
float get_stat_error(TH1F *hzjethisto,TH1F *httbarhisto, float a, float b, float c, float d, TH1F *totalhisto)//since only one diagram is now processed, the method was renamed "get_new_stat_error" (top)
|
| 1048 |
{
|
| 1049 |
float err_sum=0;
|
| 1050 |
for (int i=totalhisto->FindBin(a);i<totalhisto->FindBin(b);i++)
|
| 1051 |
{
|
| 1052 |
float currerr = totalhisto->GetBinError(i);
|
| 1053 |
err_sum += currerr * currerr;
|
| 1054 |
}
|
| 1055 |
|
| 1056 |
if (c!=d)
|
| 1057 |
{
|
| 1058 |
for (int i=totalhisto->FindBin(c);i<totalhisto->FindBin(d);i++)
|
| 1059 |
{
|
| 1060 |
float currerr = totalhisto->GetBinError(i);
|
| 1061 |
err_sum += currerr * currerr;
|
| 1062 |
}
|
| 1063 |
}
|
| 1064 |
float err_result = TMath::Sqrt(err_sum);
|
| 1065 |
/*
|
| 1066 |
float sum = hzjethisto->Integral(hzjethisto->FindBin(a),hzjethisto->FindBin(b))*zjets_weight*zjets_weight;
|
| 1067 |
sum += httbarhisto->Integral(httbarhisto->FindBin(a),httbarhisto->FindBin(b))*ttbar_weight*ttbar_weight;
|
| 1068 |
if (c!=d)
|
| 1069 |
{
|
| 1070 |
sum += hzjethisto->Integral(hzjethisto->FindBin(c),hzjethisto->FindBin(d))*zjets_weight*zjets_weight;
|
| 1071 |
sum += httbarhisto->Integral(httbarhisto->FindBin(c),httbarhisto->FindBin(d))*ttbar_weight*ttbar_weight;
|
| 1072 |
}
|
| 1073 |
float result = TMath::Sqrt(sum);
|
| 1074 |
*/
|
| 1075 |
return err_result;
|
| 1076 |
}
|
| 1077 |
|
| 1078 |
|
| 1079 |
|
| 1080 |
|
| 1081 |
float get_error_R(float nemout, float nemin, float nemouterr, float neminerr)
|
| 1082 |
{
|
| 1083 |
float sum =0;
|
| 1084 |
sum += (neminerr*neminerr)/(nemout*nemout);
|
| 1085 |
sum += ((nemin)/(nemout*nemout))*((nemin)/(nemout*nemout))*(neminerr*neminerr);
|
| 1086 |
float result = TMath::Sqrt(sum);
|
| 1087 |
return result;
|
| 1088 |
}
|
| 1089 |
|
| 1090 |
float give_n_eemm_in_correct_error(float nouteemm, float noutem, float ninem, float ninemerror, float noutemerror, float nouteemmerror)
|
| 1091 |
{
|
| 1092 |
if (noutem==0) {cout << "WATCH OUT ERROR IS INFINITE SINCE NOUTEM IS ZERO " << endl; return -9999;}
|
| 1093 |
// float r = ninem/noutem;
|
| 1094 |
// float error = TMath::Sqrt(nouteemm * r * (nouteemm/noutem + r + r*nouteemm/noutem));
|
| 1095 |
|
| 1096 |
float sum = ((nouteemm*nouteemm)/(noutem*noutem))*(ninemerror*ninemerror);
|
| 1097 |
sum += ((ninem*ninem*nouteemm*nouteemm)/(noutem*noutem*noutem*noutem))*(noutemerror*noutemerror);
|
| 1098 |
sum += ((ninem*ninem)/(noutem*noutem))*(nouteemmerror*nouteemmerror);
|
| 1099 |
float error = TMath::Sqrt(sum);
|
| 1100 |
return error;
|
| 1101 |
}
|
| 1102 |
/*
|
| 1103 |
TH1F *create_ratio_plot(TH1F *one, TH1F *two)
|
| 1104 |
{
|
| 1105 |
// returns (one)/(two)
|
| 1106 |
int nbin_one = one->GetNbinsX();
|
| 1107 |
int nbin_two = two->GetNbinsX();
|
| 1108 |
if (nbin_one!=nbin_two)
|
| 1109 |
{
|
| 1110 |
cerr << "ERROR WHEN TRYING TO GENERATE RATIO PLOT -- BINNING IS DIFFERENT!" << endl;
|
| 1111 |
return 0;
|
| 1112 |
}
|
| 1113 |
|
| 1114 |
const int nbin=nbin_one;
|
| 1115 |
float bincenter[nbin+1];
|
| 1116 |
|
| 1117 |
for (int hist_i=1;hist_i<(nbin_one+1);hist_i++)
|
| 1118 |
{
|
| 1119 |
bincenter[hist_i-1]=one->GetBinLowEdge(hist_i);
|
| 1120 |
}
|
| 1121 |
|
| 1122 |
|
| 1123 |
bincenter[nbin_one]=one->GetBinLowEdge(nbin_one)+one->GetBinWidth(nbin_one);
|
| 1124 |
|
| 1125 |
TH1F *ratio_plot = new TH1F("ratio_plot","",nbin_one,bincenter);
|
| 1126 |
|
| 1127 |
float ratio = 0;
|
| 1128 |
|
| 1129 |
for (int i=1;i<(nbin_one+1);i++)
|
| 1130 |
{
|
| 1131 |
if(two->GetBinContent(i)!=0)
|
| 1132 |
{
|
| 1133 |
float val1 = (float)one->GetBinContent(i);
|
| 1134 |
float val2 = (float)two->GetBinContent(i);
|
| 1135 |
float e1 = (float)one->GetBinError(i);
|
| 1136 |
float e2 = (float)two->GetBinError(i);
|
| 1137 |
float ratio = (val1/val2);
|
| 1138 |
float error = TMath::Sqrt((e1*e1)/(val2*val2)+(e2*e2*val1*val1)/(val2*val2*val2*val2));
|
| 1139 |
ratio_plot->SetBinContent(i,ratio);
|
| 1140 |
ratio_plot->SetBinError(i,error);
|
| 1141 |
}
|
| 1142 |
|
| 1143 |
}
|
| 1144 |
ratio_plot->SetFillColor(kBlue);
|
| 1145 |
ratio_plot->SetStats(0);
|
| 1146 |
ratio_plot->GetYaxis()->SetTitle("ratio");
|
| 1147 |
ratio_plot->GetXaxis()->SetTitle(one->GetXaxis()->GetTitle());
|
| 1148 |
ratio_plot->GetYaxis()->SetRangeUser(0,4);
|
| 1149 |
return ratio_plot;
|
| 1150 |
}
|
| 1151 |
*/
|
| 1152 |
void compute_Gaussian(TH1F *hist, TF1 *g, Double_t &err)
|
| 1153 |
{
|
| 1154 |
bool debug=false; //if you encounter an error ("could not fit function correctly") you can activate this to see the last fit that was attempted.
|
| 1155 |
if (hist->GetEntries()<10)
|
| 1156 |
{
|
| 1157 |
cout << "WATCH OUT YOU PROVIDED AN EMPTY HISTOGRAM!" << endl;
|
| 1158 |
float a,b,c;
|
| 1159 |
a=-999;
|
| 1160 |
b=-999;
|
| 1161 |
c=-999;
|
| 1162 |
g->SetParameters(a,b,c);
|
| 1163 |
}
|
| 1164 |
|
| 1165 |
|
| 1166 |
float top_edge=hist->GetBinCenter(hist->GetNbinsX())+hist->GetBinWidth(1)/2;
|
| 1167 |
float low_edge=hist->GetBinLowEdge(1);
|
| 1168 |
float inverse_bin_width = (top_edge-low_edge)/((float)hist->GetNbinsX());
|
| 1169 |
int steps=100;
|
| 1170 |
float stepwidth=top_edge/((float)steps); //we start at zero because the gaussian is suppose to have mean 0.
|
| 1171 |
bool peak_out_of_range=1; //indicates whether the peak of the fitted Gaussian is within the range used to fit the Gaussian
|
| 1172 |
bool sufficient_events=0; //indicates whether we have used at least 65% of the events for the fit.
|
| 1173 |
int counter=0;
|
| 1174 |
float a = -999, b = -999, c = -999; //parameters of the fitted Gaussian
|
| 1175 |
while (peak_out_of_range || !sufficient_events )
|
| 1176 |
{
|
| 1177 |
hist->Fit("gaus","Q","same",low_edge,counter*stepwidth);
|
| 1178 |
TF1 *fit = hist->GetFunction("gaus");
|
| 1179 |
if (fit->GetParameter(1)<((counter-2)*stepwidth))
|
| 1180 |
{
|
| 1181 |
a=fit->GetParameter(0);
|
| 1182 |
b=fit->GetParameter(1);
|
| 1183 |
c=fit->GetParameter(2);
|
| 1184 |
peak_out_of_range=0;
|
| 1185 |
if (((hist->Integral(1,hist->FindBin(counter*stepwidth)))/(hist->Integral(1,hist->GetNbinsX())))>0.9)
|
| 1186 |
{
|
| 1187 |
sufficient_events=1;
|
| 1188 |
}
|
| 1189 |
}
|
| 1190 |
if (counter>=steps)
|
| 1191 |
{
|
| 1192 |
cout << "could not fit function correctly (the mean was always outside the sub-range used to fit; please check what happened..." << endl;
|
| 1193 |
cerr << "ERROR REPORTED BY compute_Gaussian : The mean was always outside the sub-range used to fit it. Please check what happened." << endl;
|
| 1194 |
peak_out_of_range=0;
|
| 1195 |
sufficient_events=1;
|
| 1196 |
}
|
| 1197 |
if(debug&&(!peak_out_of_range)&&sufficient_events)
|
| 1198 |
{
|
| 1199 |
TCanvas *compute_Gaussian_canvas = new TCanvas("compute_Gaussian_canvas","compute_Gaussian_canvas");
|
| 1200 |
hist->Draw();
|
| 1201 |
hist->GetFunction("gaus")->Draw("same");
|
| 1202 |
compute_Gaussian_canvas->SaveAs("compute_Gaussian_debug_canvas.eps");
|
| 1203 |
}
|
| 1204 |
|
| 1205 |
counter++;
|
| 1206 |
}
|
| 1207 |
g->SetParameters(a,b,c);
|
| 1208 |
err = c/sqrt((hist->Integral(1,hist->FindBin((counter-1)*stepwidth)))*inverse_bin_width);
|
| 1209 |
|
| 1210 |
|
| 1211 |
// this corresponds to sigma/sqrt(N)
|
| 1212 |
}
|
| 1213 |
|
| 1214 |
float get_peak(TH1F *histo, Double_t *err)
|
| 1215 |
{
|
| 1216 |
float peak=-999;
|
| 1217 |
if (histo->GetEntries()>9)//only fit if there are enough points!
|
| 1218 |
{
|
| 1219 |
float low=histo->GetBinLowEdge(1);
|
| 1220 |
float hi=histo->GetBinCenter(histo->GetNbinsX())+histo->GetBinWidth(1)/2;
|
| 1221 |
TF1 *f = new TF1("f","gaus",low,hi); // MARCO HACK
|
| 1222 |
compute_Gaussian(histo,f,*err); //MARCO HACK
|
| 1223 |
peak=f->GetParameter(1);//this is the peak //MARCO HACK
|
| 1224 |
}
|
| 1225 |
else
|
| 1226 |
{
|
| 1227 |
cout << "WATCH OUT YOU PROVIDED A HISTOGRAM WITH VERY LOW STATISTICS (<10 EVENTS) ! [" << histo->GetName() << "]" << endl;
|
| 1228 |
}
|
| 1229 |
return peak;
|
| 1230 |
}
|
| 1231 |
|
| 1232 |
|
| 1233 |
float give_n_eemm_in_error(float nouteemm, float noutem, float ninem)
|
| 1234 |
{
|
| 1235 |
if (noutem==0) {cout << "WATCH OUT ERROR IS INFINITE SINCE NOUTEM IS ZERO " << endl; return -9999;}
|
| 1236 |
float r = ninem/noutem;
|
| 1237 |
float error = TMath::Sqrt(nouteemm * r * (nouteemm/noutem + r + r*nouteemm/noutem));
|
| 1238 |
return error;
|
| 1239 |
}
|
| 1240 |
|
| 1241 |
|
