/ElectroWeakAnalysis/ZMuMu/bin/csa08NewZFit_EtaPtbinned.cpp

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