/ElectroWeakAnalysis/ZMuMu/bin/zMuMuExpFit.cpp

  1#include "PhysicsTools/Utilities/interface/Parameter.h"
  2#include "PhysicsTools/Utilities/interface/ZLineShape.h"
  3#include "PhysicsTools/Utilities/interface/Gaussian.h"
  4#include "PhysicsTools/Utilities/interface/Exponential.h"
  5#include "PhysicsTools/Utilities/interface/Polynomial.h"
  6#include "PhysicsTools/Utilities/interface/Convolution.h"
  7#include "PhysicsTools/Utilities/interface/Operations.h"
  8#include "PhysicsTools/Utilities/interface/MultiHistoChiSquare.h"
  9#include "PhysicsTools/Utilities/interface/RootMinuit.h"
 10#include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
 11#include "PhysicsTools/Utilities/interface/rootPlot.h"
 12#include "TROOT.h"
 13#include "TH1.h"
 14#include "TFile.h"
 15#include <boost/program_options.hpp>
 16using namespace boost;
 17namespace po = boost::program_options;
 18
 19#include <iostream>
 20#include <algorithm> 
 21#include <exception>
 22#include <iterator>
 23#include <string>
 24#include <vector>
 25using namespace std;
 26
 27// A helper function to simplify the main part.
 28template<class T>
 29ostream& operator<<(ostream& os, const vector<T>& v) {
 30  copy(v.begin(), v.end(), ostream_iterator<T>(cout, " ")); 
 31  return os;
 32}
 33
 34//A function that sets istogram contents to 0 
 35//if they are too small
 36void fix(TH1* histo) {
 37  for(int i = 1; i <= histo->GetNbinsX(); ++i) {
 38    if(histo->GetBinContent(i) < 0.1) {
 39      histo->SetBinContent(i, 0.0);
 40      histo->SetBinError(i, 0.0);
 41    }
 42  }
 43}
 44
 45typedef funct::GaussIntegrator IntegratorConv;
 46
 47int main(int ac, char *av[]) {
 48  gROOT->SetStyle("Plain");
 49  try {
 50    typedef funct::Product<funct::Exponential, 
 51                           funct::Convolution<funct::ZLineShape, funct::Gaussian, IntegratorConv>::type >::type ZPeak;
 52    typedef funct::Product<funct::Parameter, ZPeak>::type ZMuMuSig;
 53    typedef ZMuMuSig ZMuMu;
 54    typedef ZMuMuSig ZMuTkSig;
 55    typedef funct::Product<funct::Parameter, 
 56                           funct::Product<funct::Exponential, funct::Polynomial<2> >::type >::type ZMuTkBkg;
 57    typedef funct::Sum<ZMuTkSig, ZMuTkBkg>::type ZMuTk;
 58    typedef funct::Product<funct::Parameter, funct::Gaussian>::type ZMuSaSig;
 59    typedef funct::Parameter ZMuSaBkg;
 60    typedef funct::Sum<ZMuSaSig, ZMuSaBkg>::type ZMuSa;
 61    
 62    typedef fit::MultiHistoChiSquare<ZMuMu, ZMuTk, ZMuSa> ChiSquared;
 63    fit::RootMinuitCommands<ChiSquared> commands("ElectroWeakAnalysis/ZMuMu/test/zMuMuExpFit.txt");
 64    
 65    double fMin, fMax;
 66    string ext;
 67    po::options_description desc("Allowed options");
 68    desc.add_options()
 69      ("help", "produce help message")
 70      ("include-path,I", po::value< vector<string> >(), 
 71       "include path")
 72      ("input-file", po::value< vector<string> >(), "input file")
 73      ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range")
 74      ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range")
 75      ("output-file,O", po::value<string>(&ext)->default_value(".ps"), 
 76       "output file format")
 77      ;
 78    
 79    po::positional_options_description p;
 80    p.add("input-file", -1);
 81    
 82    po::variables_map vm;
 83    po::store(po::command_line_parser(ac, av).
 84	    options(desc).positional(p).run(), vm);
 85    po::notify(vm);
 86    
 87    if (vm.count("help")) {
 88      cout << "Usage: options_description [options]\n";
 89      cout << desc;
 90      return 0;
 91      }
 92    
 93    if (vm.count("include-path")) {
 94      cout << "Include paths are: " 
 95	   << vm["include-path"].as< vector<string> >() << "\n";
 96    }
 97    
 98    if (vm.count("input-file")) {
 99      cout << "Input files are: " 
100	   << vm["input-file"].as< vector<string> >() << "\n";
101      vector<string> v_file = vm["input-file"].as< vector<string> >();
102      for(vector<string>::const_iterator it = v_file.begin(); 
103	  it != v_file.end(); ++it) {
104	TFile * root_file = new TFile(it->c_str(),"read");
105	TH1D * histoZMuMu = (TH1D*) root_file->Get("zToMM");
106	fix(histoZMuMu);
107	TH1D * histoZMuTk = (TH1D*) root_file->Get("zToMTk");
108	fix(histoZMuTk);
109	TH1D * histoZMuSa = (TH1D*) root_file->Get("zToMS");
110	fix(histoZMuSa);
111	cout << ">>> histogram loaded\n";
112	string f_string = *it;
113	replace(f_string.begin(), f_string.end(), '.', '_');
114	string plot_string = f_string + ext;
115	cout << ">>> Input files loaded\n";
116	
117	const char * kYieldZMuMu = "YieldZMuMu";
118	const char * kYieldZMuTk = "YieldZMuTk";
119	const char * kYieldZMuSa = "YieldZMuSa";
120	const char * kYieldBkgZMuTk = "YieldBkgZMuTk"; 
121	const char * kYieldBkgZMuSa = "YieldBkgZMuSa"; 
122	const char * kLambdaZMuMu = "LambdaZMuMu";
123	const char * kMass = "Mass";
124	const char * kGamma = "Gamma";
125	const char * kPhotonFactorZMuMu = "PhotonFactorZMuMu";
126	const char * kInterferenceFactorZMuMu = "InterferenceFactorZMuMu";
127	//const char * kPhotonFactorZMuTk = "PhotonFactorZMuTk";
128	//const char * kInterferenceFactorZMuTk = "InterferenceFactorZMuTk";
129	const char * kMeanZMuMu = "MeanZMuMu";
130	const char * kSigmaZMuMu = "SigmaZMuMu";
131	const char * kLambda = "Lambda";
132	const char * kA0 = "A0"; 
133	const char * kA1 = "A1"; 
134	const char * kA2 = "A2"; 
135	const char * kSigmaZMuSa = "SigmaZMuSa";
136
137	IntegratorConv integratorConv(1.e-5);
138	
139	funct::Parameter lambdaZMuMu(kLambdaZMuMu, commands.par(kLambdaZMuMu));
140	funct::Parameter mass(kMass, commands.par(kMass));
141	funct::Parameter gamma(kGamma, commands.par(kGamma));
142	funct::Parameter photonFactorZMuMu(kPhotonFactorZMuMu, commands.par(kPhotonFactorZMuMu)); 
143	funct::Parameter interferenceFactorZMuMu(kInterferenceFactorZMuMu, commands.par(kInterferenceFactorZMuMu)); 
144	//funct::Parameter photonFactorZMuTk(kPhotonFactorZMuTk, commands.par(kPhotonFactorZMuTk)); 
145	//funct::Parameter interferenceFactorZMuTk(kInterferenceFactorZMuTk, commands.par(kInterferenceFactorZMuTk)); 
146	funct::Parameter yieldZMuMu(kYieldZMuMu, commands.par(kYieldZMuMu));
147	funct::Parameter yieldZMuTk(kYieldZMuTk, commands.par(kYieldZMuTk)); 
148	funct::Parameter yieldZMuSa(kYieldZMuSa, commands.par(kYieldZMuSa)); 
149	funct::Parameter yieldBkgZMuTk(kYieldBkgZMuTk, commands.par(kYieldBkgZMuTk));
150	funct::Parameter yieldBkgZMuSa(kYieldBkgZMuSa, commands.par(kYieldBkgZMuSa));
151	funct::Parameter meanZMuMu(kMeanZMuMu, commands.par(kMeanZMuMu));
152	funct::Parameter sigmaZMuMu(kSigmaZMuMu, commands.par(kSigmaZMuMu)); 
153	funct::Parameter sigmaZMuSa(kSigmaZMuSa, commands.par(kSigmaZMuSa)); 
154	funct::Parameter lambda(kLambda, commands.par(kLambda));
155	funct::Parameter a0(kA0, commands.par(kA0));
156	funct::Parameter a1(kA1, commands.par(kA1));
157	funct::Parameter a2(kA2, commands.par(kA2));
158	                
159	ZPeak zPeak = funct::Exponential(lambdaZMuMu) * 
160	              funct::conv(funct::ZLineShape(mass, gamma, photonFactorZMuMu, interferenceFactorZMuMu), 
161			 	  funct::Gaussian(meanZMuMu, sigmaZMuMu), 
162				  -3*sigmaZMuMu.value(), 3*sigmaZMuMu.value(), integratorConv);
163	ZMuMu zMuMu = yieldZMuMu * zPeak;
164	ZMuTkBkg zMuTkBkg = 
165	  yieldBkgZMuTk * (funct::Exponential(lambda) * funct::Polynomial<2>(a0, a1, a2));
166	ZMuTk zMuTk = yieldZMuTk * zPeak + zMuTkBkg;
167	ZMuSa zMuSa = yieldZMuSa * funct::Gaussian(mass, sigmaZMuSa) + yieldBkgZMuSa;
168	
169	ChiSquared chi2(zMuMu, histoZMuMu, 
170			zMuTk, histoZMuTk, 
171			zMuSa, histoZMuSa, 
172			fMin, fMax);
173	cout << "N. deg. of freedom: " << chi2.numberOfBins() << endl;
174	fit::RootMinuit<ChiSquared> minuit(chi2, true);
175	commands.add(minuit, yieldZMuMu);
176	commands.add(minuit, yieldZMuTk);
177	commands.add(minuit, yieldZMuSa);
178	commands.add(minuit, yieldBkgZMuTk);
179	commands.add(minuit, yieldBkgZMuSa);
180	commands.add(minuit, lambdaZMuMu);
181	commands.add(minuit, mass);
182	commands.add(minuit, gamma);
183	commands.add(minuit, photonFactorZMuMu);
184	commands.add(minuit, interferenceFactorZMuMu);
185	//commands.add(minuit, photonFactorZMuTk);
186	//commands.add(minuit, interferenceFactorZMuTk);
187	commands.add(minuit, meanZMuMu);
188	commands.add(minuit, sigmaZMuMu);
189	commands.add(minuit, lambda);
190	commands.add(minuit, a0);
191	commands.add(minuit, a1);
192	commands.add(minuit, a2);
193	commands.add(minuit, sigmaZMuSa);
194	commands.run(minuit);
195	const unsigned int nPar = 17;
196	ROOT::Math::SMatrix<double, nPar, nPar, ROOT::Math::MatRepSym<double, nPar> > err;
197	minuit.getErrorMatrix(err);
198	std::cout << "error matrix:" << std::endl;
199	for(unsigned int i = 0; i < nPar; ++i) {
200	  for(unsigned int j = 0; j < nPar; ++j) {
201	    std::cout << err(i, j) << "\t";
202	  }
203	  std::cout << std::endl;
204	} 
205	minuit.printFitResults();
206	string ZMuMuPlot = "ZMuMuFit" + plot_string;
207	root::plot<ZMuMu>(ZMuMuPlot.c_str(), *histoZMuMu, zMuMu, fMin, fMax, 
208			  yieldZMuMu, lambdaZMuMu, mass, gamma, photonFactorZMuMu, interferenceFactorZMuMu, 
209			  meanZMuMu, sigmaZMuMu, 
210			  kRed, 2, kDashed, 10000, 
211			  "Z -> #mu #mu mass with isolation cut", "#mu #mu invariant mass (GeV/c^{2})", 
212			  "Events");
213	string ZMuTkPlot = "ZMuTkFit" + plot_string;
214	TF1 funZMuTk = root::tf1<ZMuTk>("ZMuTkFunction", zMuTk, fMin, fMax, 
215					yieldZMuTk, lambdaZMuMu, mass, gamma, photonFactorZMuMu, interferenceFactorZMuMu, 
216					meanZMuMu, sigmaZMuMu, 
217					yieldBkgZMuTk, lambda, a0, a1, a2);
218	funZMuTk.SetLineColor(kRed);
219	funZMuTk.SetLineWidth(2);
220	funZMuTk.SetLineStyle(kDashed);
221	funZMuTk.SetNpx(10000);
222	TF1 funZMuTkBkg = root::tf1<ZMuTkBkg>("ZMuTkBack", zMuTkBkg, fMin, fMax, 
223					      yieldBkgZMuTk, lambda, a0, a1, a2);
224	funZMuTkBkg.SetLineColor(kGreen);
225	funZMuTkBkg.SetLineWidth(2);
226	funZMuTkBkg.SetLineStyle(kDashed);
227	funZMuTkBkg.SetNpx(10000);
228	histoZMuTk->SetTitle("Z -> #mu + (unmatched) track mass with isolation cut");
229	histoZMuTk->SetXTitle("#mu + (unmatched) track invariant mass (GeV/c^{2})");
230	histoZMuTk->SetYTitle("Events");
231	TCanvas *canvas = new TCanvas("canvas");
232	histoZMuTk->Draw("e");
233	funZMuTk.Draw("same");
234	funZMuTkBkg.Draw("same");
235	canvas->SaveAs(ZMuTkPlot.c_str());
236	canvas->SetLogy();
237	string logZMuTkPlot = "log_" + ZMuTkPlot;
238	canvas->SaveAs(logZMuTkPlot.c_str());
239	string ZMuSaPlot = "ZMuSaFit" + plot_string;
240	root::plot<ZMuSa>(ZMuSaPlot.c_str(), *histoZMuSa, zMuSa, fMin, fMax, 
241			  yieldZMuSa, mass, sigmaZMuSa, yieldBkgZMuSa, 
242			  kRed, 2, kDashed, 10000, 
243			  "Z -> #mu + (unmatched) standalone mass with isolation cut", 
244			  "#mu + (unmatched) standalone invariant mass (GeV/c^{2})", 
245			  "Events");
246      }
247    }
248    
249  }
250  catch(std::exception& e) {
251    cerr << "error: " << e.what() << "\n";
252    return 1;
253  }
254  catch(...) {
255    cerr << "Exception of unknown type!\n";
256  }
257  return 0;
258}
259