PageRenderTime 43ms CodeModel.GetById 6ms app.highlight 32ms RepoModel.GetById 1ms app.codeStats 0ms

/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>
 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