/ElectroWeakAnalysis/ZMuMu/bin/csa08NewZFit_EtaPtbinned.cpp
C++ | 379 lines | 310 code | 42 blank | 27 comment | 24 complexity | 3d8601e0bc9c0e7edf11317156e3991d MD5 | raw file
1#include "FWCore/Utilities/interface/EDMException.h" 2#include "PhysicsTools/Utilities/interface/Parameter.h" 3#include "PhysicsTools/Utilities/interface/Gaussian.h" 4#include "PhysicsTools/Utilities/interface/Numerical.h" 5#include "PhysicsTools/Utilities/interface/Exponential.h" 6#include "PhysicsTools/Utilities/interface/Polynomial.h" 7#include "PhysicsTools/Utilities/interface/Constant.h" 8#include "PhysicsTools/Utilities/interface/Operations.h" 9#include "PhysicsTools/Utilities/interface/MultiHistoChiSquare.h" 10#include "PhysicsTools/Utilities/interface/RootMinuit.h" 11#include "PhysicsTools/Utilities/interface/RootMinuitCommands.h" 12#include "PhysicsTools/Utilities/interface/FunctClone.h" 13#include "PhysicsTools/Utilities/interface/rootPlot.h" 14#include "PhysicsTools/Utilities/interface/Expression.h" 15#include "PhysicsTools/Utilities/interface/HistoPdf.h" 16#include "TROOT.h" 17#include "TH1.h" 18#include "TFile.h" 19#include <boost/program_options.hpp> 20using namespace boost; 21namespace po = boost::program_options; 22 23#include <iostream> 24#include <algorithm> 25#include <exception> 26#include <iterator> 27#include <string> 28#include <vector> 29using namespace std; 30 31// A helper function to simplify the main part. 32template<class T> 33ostream& operator<<(ostream& os, const vector<T>& v) { 34 copy(v.begin(), v.end(), ostream_iterator<T>(cout, " ")); 35 return os; 36} 37 38//A function that sets istogram contents to 0 39//if they are too small 40void fix(TH1* histo) { 41 for(int i = 1; i <= histo->GetNbinsX(); ++i) { 42 if(histo->GetBinContent(i) < 0.1) { 43 histo->SetBinContent(i, 0.0); 44 histo->SetBinError(i, 0.0); 45 } 46 } 47} 48 49typedef funct::FunctExpression Expr; 50typedef fit::MultiHistoChiSquare<Expr, Expr, Expr, Expr> ChiSquared; 51 52struct sig_tag; 53struct bkg_tag; 54 55int main(int ac, char *av[]) { 56 gROOT->SetStyle("Plain"); 57 try { 58 59 double fMin, fMax; 60 string ext, variable, muCharge; 61 int binNumber; 62 po::options_description desc("Allowed options"); 63 desc.add_options() 64 ("help,h", "produce help message") 65 ("input-file,i", po::value< vector<string> >(), "input file") 66 ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range") 67 ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range") 68 ("eta_or_pt,v", po::value<string>(&variable)->default_value("eta"), "variable to study (eta or pt)") 69 ("charge,q", po::value<string>(&muCharge)->default_value("minus"),"muon charge to study (minus or plus)") 70 ("binNum,b", po::value<int>(&binNumber)->default_value(0), "cynematic bin to fit") 71 ("plot-format,p", po::value<string>(&ext)->default_value("ps"), 72 "output plot format") 73 ; 74 75 po::positional_options_description p; 76 p.add("input-file", -1); 77 78 po::variables_map vm; 79 po::store(po::command_line_parser(ac, av). 80 options(desc).positional(p).run(), vm); 81 po::notify(vm); 82 83 if (vm.count("help")) { 84 cout << "Usage: options_description [options]\n"; 85 cout << desc; 86 return 0; 87 } 88 89 fit::RootMinuitCommands<ChiSquared> commands("csa08NewZFit.txt"); 90 91 const unsigned int rebinMuMu = 1, rebinMuTk = 2, rebinMuSa = 1; 92 // assume that the bin size is 1 GeV!!! 93 funct::Constant rebinMuMuConst(rebinMuMu), rebinMuTkConst(rebinMuTk), rebinMuSaConst(rebinMuSa); 94 95 if (vm.count("input-file")) { 96 cout << "Input files are: " 97 << vm["input-file"].as< vector<string> >() << "\n"; 98 vector<string> v_file = vm["input-file"].as< vector<string> >(); 99 for(vector<string>::const_iterator it = v_file.begin(); 100 it != v_file.end(); ++it) { 101 TFile * root_file = new TFile(it->c_str(),"read"); 102 cout <<"start " << endl; 103 // variable and charge definition at moment in manual way 104 // string variable = "eta"; 105 // string muCharge = "minus"; 106 //////////////////////////////////////// 107 108 stringstream sslabelZMuMu2HLT; 109 sslabelZMuMu2HLT << "zMuMu_efficiencyAnalyzer/" << variable << "Intervals" << "/zmumu2HLT" << muCharge << "_" << variable << "Range" << binNumber; 110 stringstream sslabelZMuMu1HLT; 111 sslabelZMuMu1HLT << "zMuMu_efficiencyAnalyzer/" << variable << "Intervals" << "/zmumu1HLT" << muCharge << "_" << variable << "Range" << binNumber; 112 stringstream sslabelZMuTk; 113 sslabelZMuTk << "zMuMu_efficiencyAnalyzer/" << variable << "Intervals" << "/zmutrack" << muCharge << "_" << variable << "Range" << binNumber; 114 stringstream sslabelZMuSa; 115 sslabelZMuSa << "zMuMu_efficiencyAnalyzer/" << variable << "Intervals" << "/zmusta" << muCharge << "_" << variable << "Range" << binNumber; 116 117 cout << "histoZMuMu2HLT: " << sslabelZMuMu2HLT.str() << endl; 118 TH1D * histoZMuMu2HLT = (TH1D*) root_file->Get(sslabelZMuMu2HLT.str().c_str()); 119 histoZMuMu2HLT->Rebin(rebinMuMu); 120 fix(histoZMuMu2HLT); 121 cout << "histoZMuMu1HLT: " << sslabelZMuMu1HLT.str() << endl; 122 TH1D * histoZMuMu1HLT = (TH1D*) root_file->Get(sslabelZMuMu1HLT.str().c_str()); 123 histoZMuMu1HLT->Rebin(rebinMuMu); 124 fix(histoZMuMu1HLT); 125 cout << "histoZMuTk: " << sslabelZMuTk.str() << endl; 126 TH1D * histoZMuTk = (TH1D*) root_file->Get(sslabelZMuTk.str().c_str()); 127 // histoZMuTk->Rebin(rebinMuTk); 128 fix(histoZMuTk); 129 cout << "histoZMuSa: " << sslabelZMuSa.str() << endl; 130 TH1D * histoZMuSa = (TH1D*) root_file->Get(sslabelZMuSa.str().c_str()); 131 // histoZMuSa->Rebin(rebinMuSa); 132 fix(histoZMuSa); 133 cout << ">>> histogram loaded\n"; 134 135 string f_string = *it; 136 replace(f_string.begin(), f_string.end(), '.', '_'); 137 replace(f_string.begin(), f_string.end(), '/', '_'); 138 string plot_string = f_string + "." + ext; 139 cout << ">>> Input files loaded\n" << f_string << endl; 140 141 const char * kYieldZMuMu = "YieldZMuMu"; 142 const char * kEfficiencyHLT = "EfficiencyHLT"; 143 const char * kEfficiencyTk = "EfficiencyTk"; 144 const char * kEfficiencySa = "EfficiencySa"; 145 const char * kYieldBkgZMuTk = "YieldBkgZMuTk"; 146 const char * kBeta = "Beta"; 147 const char * kLambda = "Lambda"; 148 // const char * kA0 = "A0"; 149 // const char * kA1 = "A1"; 150 // const char * kA2 = "A2"; 151 152 funct::Parameter yieldZMuMu(kYieldZMuMu, commands.par(kYieldZMuMu)); 153 funct::Parameter effHLT(kEfficiencyHLT, commands.par(kEfficiencyHLT)); 154 funct::Parameter effTk(kEfficiencyTk, commands.par(kEfficiencyTk)); 155 funct::Parameter effSa(kEfficiencySa, commands.par(kEfficiencySa)); 156 funct::Parameter yieldBkgZMuTk(kYieldBkgZMuTk, commands.par(kYieldBkgZMuTk)); 157 funct::Parameter beta(kBeta, commands.par(kBeta)); 158 funct::Parameter lambda(kLambda, commands.par(kLambda)); 159 // funct::Parameter a0(kA0, commands.par(kA0)); 160 // funct::Parameter a1(kA1, commands.par(kA1)); 161 // funct::Parameter a2(kA2, commands.par(kA2)); 162 funct::Constant cFMin(fMin), cFMax(fMax); 163 164 // add zMuMu2HLT and zMuMu1HLT to build pdf 165 TH1D *histoZMuMu = (TH1D *) histoZMuMu2HLT->Clone(); 166 histoZMuMu->Sumw2(); 167 histoZMuMu->Add(histoZMuMu2HLT,histoZMuMu1HLT); 168 169 // count ZMuMu Yield 170 double nZMuMu = 0; 171 { 172 unsigned int nBins = histoZMuMu->GetNbinsX(); 173 double xMin = histoZMuMu->GetXaxis()->GetXmin(); 174 double xMax = histoZMuMu->GetXaxis()->GetXmax(); 175 double deltaX =(xMax - xMin) / nBins; 176 for(unsigned int i = 0; i < nBins; ++i) { 177 double x = xMin + (i +.5) * deltaX; 178 if(x > fMin && x < fMax) 179 nZMuMu += histoZMuMu->GetBinContent(i+1); 180 } 181 } 182 183 // count ZMuMu2HLT Yield 184 double nZMuMu2HLT = 0; 185 { 186 unsigned int nBins = histoZMuMu2HLT->GetNbinsX(); 187 double xMin = histoZMuMu2HLT->GetXaxis()->GetXmin(); 188 double xMax = histoZMuMu2HLT->GetXaxis()->GetXmax(); 189 double deltaX =(xMax - xMin) / nBins; 190 for(unsigned int i = 0; i < nBins; ++i) { 191 double x = xMin + (i +.5) * deltaX; 192 if(x > fMin && x < fMax) 193 nZMuMu2HLT += histoZMuMu2HLT->GetBinContent(i+1); 194 } 195 } 196 197 // count ZMuMu1HLT Yield 198 double nZMuMu1HLT = 0; 199 { 200 unsigned int nBins = histoZMuMu1HLT->GetNbinsX(); 201 double xMin = histoZMuMu1HLT->GetXaxis()->GetXmin(); 202 double xMax = histoZMuMu1HLT->GetXaxis()->GetXmax(); 203 double deltaX =(xMax - xMin) / nBins; 204 for(unsigned int i = 0; i < nBins; ++i) { 205 double x = xMin + (i +.5) * deltaX; 206 if(x > fMin && x < fMax) 207 nZMuMu1HLT += histoZMuMu1HLT->GetBinContent(i+1); 208 } 209 } 210 211 // count ZMuSa Yield (too low statistis so we just check the number assuming 0 background) 212 double nZMuSa = 0; 213 { 214 unsigned int nBins = histoZMuSa->GetNbinsX(); 215 double xMin = histoZMuSa->GetXaxis()->GetXmin(); 216 double xMax = histoZMuSa->GetXaxis()->GetXmax(); 217 double deltaX =(xMax - xMin) / nBins; 218 for(unsigned int i = 0; i < nBins; ++i) { 219 double x = xMin + (i +.5) * deltaX; 220 if(x > fMin && x < fMax) 221 nZMuSa += histoZMuSa->GetBinContent(i+1); 222 } 223 } 224 225 cout << ">>> count of ZMuMu yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu << endl; 226 cout << ">>> count of ZMuMu2HLT yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu2HLT << endl; 227 cout << ">>> count of ZMuMu1HLT yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu1HLT << endl; 228 cout << ">>> count of ZMuSa yield in the range [" << fMin << ", " << fMax << "]: " << nZMuSa << endl; 229 230 231 funct::RootHistoPdf zPdfMuTk(*histoZMuMu, fMin, fMax); 232 zPdfMuTk.rebin(rebinMuTk); 233 234 funct::Numerical<2> _2; 235 funct::Numerical<1> _1; 236 237 Expr zMuMu2HLTEffTerm = effTk * effSa * effHLT; 238 Expr zMuMu1HLTEffTerm = effTk * effSa * (_1 - effHLT); 239 Expr zMuTkEffTerm = effTk * (_1 - effSa); 240 Expr zMuSaEffTerm = effSa * (_1 - effTk); 241 242 Expr zMuMu2HLT = rebinMuMuConst * zMuMu2HLTEffTerm * yieldZMuMu; 243 Expr zMuMu1HLT = rebinMuMuConst * zMuMu1HLTEffTerm * yieldZMuMu; 244 245 Expr zMuTkBkg = yieldBkgZMuTk * funct::Exponential(lambda); 246 //* funct::Polynomial<2>(a0, a1, a2); 247 Expr zMuTkBkgScaled = rebinMuTkConst * zMuTkBkg; 248 Expr zMuTk = rebinMuTkConst * (zMuTkEffTerm * yieldZMuMu * zPdfMuTk + zMuTkBkg); 249 Expr zMuSa = rebinMuSaConst * zMuSaEffTerm * yieldZMuMu; 250 251 TH1D histoZMM2HLTCount("histoZMM2HLTCount", "", 1, fMin, fMax); 252 histoZMM2HLTCount.Fill(100, nZMuMu2HLT); 253 TH1D histoZMM1HLTCount("histoZMM1HLTCount", "", 1, fMin, fMax); 254 histoZMM1HLTCount.Fill(100, nZMuMu1HLT); 255 TH1D histoZMSCount("histoZMSCount", "", 1, fMin, fMax); 256 histoZMSCount.Fill(100, nZMuSa); 257 258 ChiSquared chi2(zMuMu2HLT, & histoZMM2HLTCount, 259 zMuMu1HLT, & histoZMM1HLTCount, 260 zMuTk, histoZMuTk, 261 zMuSa, & histoZMSCount, 262 fMin, fMax); 263 cout << "N. deg. of freedom: " << chi2.degreesOfFreedom() << endl; 264 fit::RootMinuit<ChiSquared> minuit(chi2, true); 265 commands.add(minuit, yieldZMuMu); 266 commands.add(minuit, effHLT); 267 commands.add(minuit, effTk); 268 commands.add(minuit, effSa); 269 commands.add(minuit, yieldBkgZMuTk); 270 commands.add(minuit, lambda); 271 commands.add(minuit, beta); 272 // commands.add(minuit, a0); 273 // commands.add(minuit, a1); 274 // commands.add(minuit, a2); 275 commands.run(minuit); 276 const unsigned int nPar = 7;//WARNING: this must be updated manually for now 277 ROOT::Math::SMatrix<double, nPar, nPar, ROOT::Math::MatRepSym<double, nPar> > err; 278 minuit.getErrorMatrix(err); 279 std::cout << "error matrix:" << std::endl; 280 for(unsigned int i = 0; i < nPar; ++i) { 281 for(unsigned int j = 0; j < nPar; ++j) { 282 std::cout << err(i, j) << "\t"; 283 } 284 std::cout << std::endl; 285 } 286 minuit.printFitResults(); 287 288 289 double s; 290 s = 0; 291 for(int i = 1; i <= histoZMuMu2HLT->GetNbinsX(); ++i) 292 s += histoZMuMu2HLT->GetBinContent(i); 293 histoZMuMu2HLT->SetEntries(s); 294 s = 0; 295 for(int i = 1; i <= histoZMuMu1HLT->GetNbinsX(); ++i) 296 s += histoZMuMu1HLT->GetBinContent(i); 297 histoZMuMu1HLT->SetEntries(s); 298 s = 0; 299 for(int i = 1; i <= histoZMuTk->GetNbinsX(); ++i) 300 s += histoZMuTk->GetBinContent(i); 301 histoZMuTk->SetEntries(s); 302 s = 0; 303 for(int i = 1; i <= histoZMuSa->GetNbinsX(); ++i) 304 s += histoZMuSa->GetBinContent(i); 305 histoZMuSa->SetEntries(s); 306 stringstream mybin; 307 mybin << muCharge << "_" << variable << binNumber << "_"; 308 string ZMuMu2HLTPlot = "ZMuMu2HLTFit_muon" + mybin.str() + plot_string; 309 root::plot<Expr>(ZMuMu2HLTPlot.c_str(), *histoZMuMu2HLT, zMuMu2HLT, fMin, fMax, 310 effHLT, effTk, effSa, yieldZMuMu, 311 kRed, 2, kDashed, 100, 312 "Z -> #mu #mu mass (2HLT)", "#mu #mu invariant mass (GeV/c^{2})", 313 "Events"); 314 315 string ZMuMu1HLTPlot = "ZMuMu1HLTFit_muon" + mybin.str() + plot_string; 316 root::plot<Expr>(ZMuMu1HLTPlot.c_str(), *histoZMuMu1HLT, zMuMu1HLT, fMin, fMax, 317 effHLT, effTk, effSa, yieldZMuMu, 318 kRed, 2, kDashed, 100, 319 "Z -> #mu #mu mass (1HLT)", "#mu #mu invariant mass (GeV/c^{2})", 320 "Events"); 321 322 string ZMuTkPlot = "ZMuTkFit_muon" + mybin.str() + plot_string; 323 root::plot<Expr>(ZMuTkPlot.c_str(), *histoZMuTk, zMuTk, fMin, fMax, 324 effHLT, effTk, effSa, yieldZMuMu, 325 yieldBkgZMuTk, lambda, 326 //a0, a1, a2, 327 kRed, 2, kDashed, 100, 328 "Z -> #mu + (unmatched) track mass", "#mu #mu invariant mass (GeV/c^{2})", 329 "Events"); 330 // string ZMuTkPlot = "ZMuTkFit_muon" + muCharge + variable + binNumber + plot_string; 331 TF1 funZMuTk = root::tf1_t<sig_tag, Expr>("ZMuTkFunction", zMuTk, fMin, fMax, 332 effHLT, effTk, effSa, yieldZMuMu, 333 yieldBkgZMuTk, lambda); 334 funZMuTk.SetLineColor(kRed); 335 funZMuTk.SetLineWidth(2); 336 funZMuTk.SetLineStyle(kDashed); 337 funZMuTk.SetNpx(10000); 338 TF1 funZMuTkBkg = root::tf1_t<bkg_tag, Expr>("ZMuTkBack", zMuTkBkgScaled, fMin, fMax, 339 yieldBkgZMuTk, lambda); 340 funZMuTkBkg.SetLineColor(kGreen); 341 funZMuTkBkg.SetLineWidth(2); 342 funZMuTkBkg.SetLineStyle(kDashed); 343 funZMuTkBkg.SetNpx(10000); 344 histoZMuTk->SetTitle("Z -> #mu + (unmatched) track mass"); 345 histoZMuTk->SetXTitle("#mu + (unmatched) track invariant mass (GeV/c^{2})"); 346 histoZMuTk->SetYTitle("Events"); 347 TCanvas *canvas = new TCanvas("canvas"); 348 histoZMuTk->Draw("e"); 349 funZMuTkBkg.Draw("same"); 350 funZMuTk.Draw("same"); 351 canvas->SaveAs(ZMuTkPlot.c_str()); 352 canvas->SetLogy(); 353 string logZMuTkPlot = "log_" + ZMuTkPlot; 354 canvas->SaveAs(logZMuTkPlot.c_str()); 355 string ZMuSaPlot = "ZMuSaFit_muon" + mybin.str() + plot_string; 356 root::plot<Expr>(ZMuSaPlot.c_str(), *histoZMuSa, zMuSa, fMin, fMax, 357 effHLT, effSa, effTk, yieldZMuMu, 358 kRed, 2, kDashed, 10000, 359 "Z -> #mu + (unmatched) standalone mass", 360 "#mu + (unmatched) standalone invariant mass (GeV/c^{2})", 361 "Events"); 362 363 } 364 365 } 366 } 367 catch(exception& e) { 368 cerr << "error: " << e.what() << "\n"; 369 return 1; 370 } 371 catch(...) { 372 cerr << "Exception of unknown type!\n"; 373 } 374 375 return 0; 376} 377 378 379