PageRenderTime 82ms CodeModel.GetById 10ms app.highlight 66ms RepoModel.GetById 1ms app.codeStats 0ms

/ElectroWeakAnalysis/ZMuMu/bin/zChi2Fit.cpp

https://github.com/aivanov-cern/cmssw
C++ | 602 lines | 507 code | 74 blank | 21 comment | 39 complexity | 4c92b2ea0cb0c9e638b1eaf636cdedaf 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/HistoChiSquare.h"
 11#include "PhysicsTools/Utilities/interface/HistoPoissonLikelihoodRatio.h"
 12#include "PhysicsTools/Utilities/interface/RootMinuit.h"
 13#include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
 14#include "PhysicsTools/Utilities/interface/FunctClone.h"
 15#include "PhysicsTools/Utilities/interface/rootPlot.h"
 16#include "PhysicsTools/Utilities/interface/Expression.h"
 17#include "PhysicsTools/Utilities/interface/HistoPdf.h"
 18#include "FWCore/FWLite/interface/AutoLibraryLoader.h"
 19#include "TROOT.h"
 20#include "TSystem.h"
 21#include "TH1.h"
 22#include "TFile.h"
 23#include <fstream>
 24#include <iostream>
 25#include <algorithm> 
 26#include <exception>
 27#include <iterator>
 28#include <string>
 29#include <vector>
 30using namespace std;
 31
 32// A helper function to simplify the main part.
 33template<class T>
 34ostream& operator<<(ostream& os, const vector<T>& v) {
 35  copy(v.begin(), v.end(), ostream_iterator<T>(cout, " ")); 
 36  return os;
 37}
 38
 39// A function that get histogram and sets contents to 0 
 40// if entries are too small
 41TH1 * getHisto(TFile * file, const char * name, unsigned int rebin) {
 42  TObject * h = file->Get(name);
 43  if(h == 0)
 44    throw edm::Exception(edm::errors::Configuration) 
 45      << "Can't find object " << name << "\n";
 46  TH1 * histo = dynamic_cast<TH1*>(h);
 47  if(histo == 0)
 48    throw edm::Exception(edm::errors::Configuration) 
 49      << "Object " << name << " is of type " << h->ClassName() << ", not TH1\n";
 50  histo->Rebin(rebin);  
 51  for(int i = 1; i <= histo->GetNbinsX(); ++i) {
 52    if(histo->GetBinContent(i) < 0.1) {
 53      histo->SetBinContent(i, 0.0);
 54      histo->SetBinError(i, 0.0);
 55    }
 56  }
 57  return histo;
 58}
 59
 60struct sig_tag;
 61struct bkg_tag;
 62
 63typedef funct::FunctExpression Expr;
 64typedef fit::HistoChiSquare<funct::FunctExpression> ExprChi2;
 65typedef fit::HistoPoissonLikelihoodRatio<funct::FunctExpression> ExprPLR;
 66
 67double fMin, fMax;
 68 unsigned int rebinMuMuNoIso ,rebinMuMu =1 , rebinMuMu1HLT , rebinMuMu2HLT , rebinMuTk , rebinMuSa;
 69 // assume that the bin size is 1 GeV!!! 
 70string ext, region;
 71bool nonIsoTemplateFromMC;
 72
 73template<typename T>
 74struct PlotPrefix { };
 75
 76template<>
 77struct PlotPrefix<ExprChi2> {
 78  static string str() { return "chi2"; }
 79};
 80 
 81template<>
 82struct PlotPrefix<ExprPLR> {
 83  static string str() { return "plr"; }
 84};
 85 
 86template<typename T>
 87int main_t(const vector<string> & v_file){
 88  typedef fit::MultiHistoChiSquare<T, T, T, T, T> ChiSquared;
 89  fit::RootMinuitCommands<ChiSquared> commands("zChi2Fit.txt");
 90  
 91  cout << "minuit command file completed" << endl;
 92 
 93  funct::Constant rebinMuMuNoIsoConst(rebinMuMuNoIso), rebinMuMuConst(rebinMuMu), 
 94    rebinMuMu1HLTConst(rebinMuMu1HLT), rebinMuMu2HLTConst(rebinMuMu2HLT), 
 95    rebinMuTkConst(rebinMuTk), rebinMuSaConst(rebinMuSa);
 96  
 97  for(vector<string>::const_iterator it = v_file.begin(); it != v_file.end(); ++it) {
 98    TFile * root_file = new TFile(it->c_str(), "read");
 99    
100    // default when region==all
101    //    TH1 * histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlots/zMass",rebinMuMuNoIso);
102    TH1 * histoZMuMuNoIso = getHisto(root_file, "oneNonIsolatedZToMuMuPlots/zMass",rebinMuMuNoIso);
103    TH1 * histoZMuMu = getHisto(root_file, "goodZToMuMuPlots/zMass",rebinMuMu);
104    TH1 * histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlots/zMass", rebinMuMu1HLT);
105    TH1 * histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlots/zMass", rebinMuMu2HLT);
106    TH1 * histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlots/zMass", rebinMuTk);
107    TH1 * histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlots/zMass", rebinMuSa);
108    TH1 * histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogram/zMass", rebinMuSa);
109   
110    TH1 * histoZMuMuNoIsoTemplateFromMC= histoZMuMu;
111    if (nonIsoTemplateFromMC) {
112      //      histoZMuMuNoIsoTemplateFromMC = getHisto(root_file, "nonIsolatedZToMuMuPlotsMC/zMass",rebinMuMu);
113      histoZMuMuNoIsoTemplateFromMC = getHisto(root_file, "oneNonIsolatedZToMuMuPlotsMC/zMass",rebinMuMu);
114  }    
115    if (region=="barrel"){
116      histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsBarrel/zMass",rebinMuMuNoIso);
117      histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsBarrel/zMass",rebinMuMu);
118      histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsBarrel/zMass", rebinMuMu1HLT);
119      histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsBarrel/zMass", rebinMuMu2HLT);
120      histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsBarrel/zMass", rebinMuTk);
121      histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsBarrel/zMass", rebinMuSa);
122      histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramBarrel/zMass", rebinMuSa);
123    }
124    
125    if (region=="endcap"){
126      histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsEndCap/zMass",rebinMuMuNoIso);
127      histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsEndCap/zMass",rebinMuMu);
128      histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsEndCap/zMass", rebinMuMu1HLT);
129      histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsEndCap/zMass", rebinMuMu2HLT);
130      histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsEndCap/zMass", rebinMuTk);
131      histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsEndCap/zMass", rebinMuSa);
132      histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramEndCap/zMass", rebinMuSa);
133    }
134    
135    if (region=="barrend"){
136      histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsBarrEnd/zMass",rebinMuMuNoIso);
137      histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsBarrEnd/zMass",rebinMuMu);
138      histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsBarrEnd/zMass", rebinMuMu1HLT);
139      histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsBarrEnd/zMass", rebinMuMu2HLT);
140      histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsBarrEnd/zMass", rebinMuTk);
141      histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsBarrEnd/zMass", rebinMuSa);
142      histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramBarrEnd/zMass", rebinMuSa);
143    }
144    
145    if (region!="endcap" && region!="barrel" && region!="barrend" && region!="all"  ){
146      cout<< "not a valid region selected"<< endl;
147      cout << "possible choises are: all, barrel, endcap, barrend " << endl;
148      return 0;
149    }
150    
151    cout << ">>> histogram loaded\n";
152    string f_string = *it + "_" + PlotPrefix<T>::str() + "_";
153    replace(f_string.begin(), f_string.end(), '.', '_');
154    replace(f_string.begin(), f_string.end(), '/', '_');
155    string plot_string = f_string + "." + ext;
156    cout << ">>> Input files loaded\n";
157    
158    const char * kYieldZMuMu = "YieldZMuMu";
159    const char * kEfficiencyTk = "EfficiencyTk";
160    const char * kEfficiencySa = "EfficiencySa";
161    const char * kEfficiencyIso = "EfficiencyIso";
162    const char * kEfficiencyHLT = "EfficiencyHLT"; 
163    const char * kYieldBkgZMuTk = "YieldBkgZMuTk"; 
164    const char * kYieldBkgZMuSa = "YieldBkgZMuSa"; 
165    const char * kYieldBkgZMuMuNotIso = "YieldBkgZMuMuNotIso";
166    const char * kAlpha = "Alpha";
167    const char * kBeta = "Beta";
168    const char * kLambda = "Lambda";
169    const char * kA0 = "A0"; 
170    const char * kA1 = "A1"; 
171    const char * kA2 = "A2"; 
172    const char * kB0 = "B0"; 
173    const char * kB1 = "B1"; 
174    const char * kB2 = "B2"; 
175    const char * kC0 = "C0"; 
176    const char * kC1 = "C1"; 
177    const char * kC2 = "C2"; 
178    
179    funct::Parameter yieldZMuMu(kYieldZMuMu, commands.par(kYieldZMuMu));
180    funct::Parameter effTk(kEfficiencyTk, commands.par(kEfficiencyTk)); 
181    funct::Parameter effSa(kEfficiencySa, commands.par(kEfficiencySa)); 
182    funct::Parameter effIso(kEfficiencyIso, commands.par(kEfficiencyIso)); 
183    funct::Parameter effHLT(kEfficiencyHLT, commands.par(kEfficiencyHLT)); 
184    funct::Parameter yieldBkgZMuTk(kYieldBkgZMuTk, commands.par(kYieldBkgZMuTk));
185    funct::Parameter yieldBkgZMuSa(kYieldBkgZMuSa, commands.par(kYieldBkgZMuSa));
186    funct::Parameter yieldBkgZMuMuNotIso(kYieldBkgZMuMuNotIso, commands.par(kYieldBkgZMuMuNotIso));
187    funct::Parameter alpha(kAlpha, commands.par(kAlpha));
188    funct::Parameter beta(kBeta, commands.par(kBeta));
189    funct::Parameter lambda(kLambda, commands.par(kLambda));
190    funct::Parameter a0(kA0, commands.par(kA0));
191    funct::Parameter a1(kA1, commands.par(kA1));
192    funct::Parameter a2(kA2, commands.par(kA2));
193    funct::Parameter b0(kB0, commands.par(kB0));
194    funct::Parameter b1(kB1, commands.par(kB1));
195    funct::Parameter b2(kB2, commands.par(kB2));
196    funct::Parameter c0(kC0, commands.par(kC0));
197    funct::Parameter c1(kC1, commands.par(kC1));
198    funct::Parameter c2(kC2, commands.par(kC2));
199    funct::Constant cFMin(fMin), cFMax(fMax);
200    
201    // count ZMuMu Yield
202    double nZMuMu = 0, nZMuMu1HLT = 0, nZMuMu2HLT = 0;
203    {
204      unsigned int nBins = histoZMuMu->GetNbinsX();
205      double xMin = histoZMuMu->GetXaxis()->GetXmin();
206      double xMax = histoZMuMu->GetXaxis()->GetXmax();
207      double deltaX =(xMax - xMin) / nBins;
208      for(unsigned int i = 0; i < nBins; ++i) { 
209	double x = xMin + (i +.5) * deltaX;
210	if(x > fMin && x < fMax){
211	  nZMuMu += histoZMuMu->GetBinContent(i+1);
212	  nZMuMu1HLT += histoZMuMu1HLT->GetBinContent(i+1);
213	  nZMuMu2HLT += histoZMuMu2HLT->GetBinContent(i+1);
214	}
215      }
216    }
217    // aggiungi 1HLT 2HLT
218    cout << ">>> count of ZMuMu yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu << endl;
219    cout << ">>> count of ZMuMu (1HLT) yield in the range [" << fMin << ", " << fMax << "]: "  << nZMuMu1HLT << endl;
220    cout << ">>> count of ZMuMu (2HLT) yield in the range [" << fMin << ", " << fMax << "]: "  << nZMuMu2HLT << endl;
221    funct::RootHistoPdf zPdfMuMu(*histoZMuMu, fMin, fMax);
222      //assign ZMuMu as pdf
223    funct::RootHistoPdf zPdfMuMuNonIso = zPdfMuMu;
224    if (nonIsoTemplateFromMC) {
225      funct::RootHistoPdf zPdfMuMuNoIsoFromMC(*histoZMuMuNoIsoTemplateFromMC, fMin, fMax);
226      zPdfMuMuNonIso = zPdfMuMuNoIsoFromMC;
227    }    
228    
229    funct::RootHistoPdf zPdfMuTk = zPdfMuMu;
230    funct::RootHistoPdf zPdfMuMu1HLT = zPdfMuMu;
231    funct::RootHistoPdf zPdfMuMu2HLT = zPdfMuMu;
232    funct::RootHistoPdf zPdfMuSa(*histoZMuSaFromMuMu, fMin, fMax);
233    zPdfMuMuNonIso.rebin(rebinMuMuNoIso/rebinMuMu);
234    zPdfMuTk.rebin(rebinMuTk/rebinMuMu);
235    zPdfMuMu1HLT.rebin(rebinMuMu1HLT/rebinMuMu);
236    zPdfMuMu2HLT.rebin(rebinMuMu2HLT/rebinMuMu);
237    
238    funct::Numerical<2> _2;
239    funct::Numerical<1> _1;
240    
241    //Efficiency term
242    Expr zMuMuEff1HLTTerm = _2 * (effTk ^ _2) *  (effSa ^ _2) * (effIso ^ _2) * effHLT * (_1 - effHLT); 
243    Expr zMuMuEff2HLTTerm = (effTk ^ _2) *  (effSa ^ _2) * (effIso ^ _2) * (effHLT ^ _2) ; 
244    //    Expr zMuMuNoIsoEffTerm = (effTk ^ _2) * (effSa ^ _2) * (_1 - (effIso ^ _2)) * (_1 - ((_1 - effHLT)^_2));
245    // change to both hlt and one not iso
246    Expr zMuMuNoIsoEffTerm = _2 * (effTk ^ _2) * (effSa ^ _2) * effIso * (_1 - effIso) * (effHLT^_2);
247    Expr zMuTkEffTerm = _2 * (effTk ^ _2) * effSa * (_1 - effSa) * (effIso ^ _2) * effHLT;
248    Expr zMuSaEffTerm = _2 * (effSa ^ _2) * effTk * (_1 - effTk) * (effIso ^ _2) * effHLT;
249    
250    Expr zMuMu1HLT = rebinMuMu1HLTConst * zMuMuEff1HLTTerm * yieldZMuMu;
251    Expr zMuMu2HLT = rebinMuMu2HLTConst * zMuMuEff2HLTTerm * yieldZMuMu;
252    
253    Expr zMuTkBkg = yieldBkgZMuTk * funct::Exponential(lambda)* funct::Polynomial<2>(a0, a1, a2);
254    Expr zMuTkBkgScaled = rebinMuTkConst * zMuTkBkg;
255    Expr zMuTk = rebinMuTkConst * (zMuTkEffTerm * yieldZMuMu * zPdfMuTk + zMuTkBkg);
256    
257    Expr zMuMuNoIsoBkg = yieldBkgZMuMuNotIso * funct::Exponential(alpha)* funct::Polynomial<2>(b0, b1, b2);
258    Expr zMuMuNoIsoBkgScaled = rebinMuMuNoIsoConst * zMuMuNoIsoBkg;
259    Expr zMuMuNoIso = rebinMuMuNoIsoConst * (zMuMuNoIsoEffTerm * yieldZMuMu * zPdfMuMuNonIso + zMuMuNoIsoBkg);
260    
261
262    Expr zMuSaBkg = yieldBkgZMuSa * funct::Exponential(beta)* funct::Polynomial<2>(c0, c1, c2);
263    Expr zMuSaBkgScaled = rebinMuSaConst * zMuSaBkg;
264    Expr zMuSa = rebinMuSaConst * (zMuSaEffTerm * yieldZMuMu * zPdfMuSa  + zMuSaBkg );
265    
266    TH1D histoZCount1HLT("histoZCount1HLT", "", 1, fMin, fMax);
267    histoZCount1HLT.Fill(100, nZMuMu1HLT);
268    TH1D histoZCount2HLT("histoZCount2HLT", "", 1, fMin, fMax);
269    histoZCount2HLT.Fill(100, nZMuMu2HLT);
270    
271    ChiSquared chi2(zMuMu1HLT, & histoZCount1HLT,
272		    zMuMu2HLT, & histoZCount2HLT,
273		    zMuTk, histoZMuTk, 
274		    zMuSa, histoZMuSa, 
275		    zMuMuNoIso,histoZMuMuNoIso,
276		    fMin, fMax);
277    cout << "N. bins: " << chi2.numberOfBins() << endl;
278    
279    fit::RootMinuit<ChiSquared> minuit(chi2, true);
280    commands.add(minuit, yieldZMuMu);
281    commands.add(minuit, effTk);
282    commands.add(minuit, effSa);
283    commands.add(minuit, effIso);
284    commands.add(minuit, effHLT);
285    commands.add(minuit, yieldBkgZMuTk);
286    commands.add(minuit, yieldBkgZMuSa);
287    commands.add(minuit, yieldBkgZMuMuNotIso);
288    commands.add(minuit, lambda);
289    commands.add(minuit, alpha);
290    commands.add(minuit, beta);
291    commands.add(minuit, a0);
292    commands.add(minuit, a1);
293    commands.add(minuit, a2);
294    commands.add(minuit, b0);
295    commands.add(minuit, b1);
296    commands.add(minuit, b2);
297    commands.add(minuit, c0);
298    commands.add(minuit, c1);
299    commands.add(minuit, c2);
300    commands.run(minuit);
301    const unsigned int nPar = 20;//WARNIG: this must be updated manually for now
302    ROOT::Math::SMatrix<double, nPar, nPar, ROOT::Math::MatRepSym<double, nPar> > err;
303    minuit.getErrorMatrix(err);
304    
305    std::cout << "error matrix:" << std::endl;
306    for(unsigned int i = 0; i < nPar; ++i) {
307      for(unsigned int j = 0; j < nPar; ++j) {
308	  std::cout << err(i, j) << "\t";
309      }
310      std::cout << std::endl;
311    } 
312    minuit.printFitResults();
313    ofstream myfile;
314    myfile.open ("fitResult.txt", ios::out | ios::app);
315    myfile<<"\n";
316    double Y =  minuit.getParameterError("YieldZMuMu");
317    double dY = minuit.getParameterError("YieldZMuMu", Y);
318    double tk_eff =  minuit.getParameterError("EfficiencyTk");
319    double dtk_eff = minuit.getParameterError("EfficiencyTk", tk_eff);
320    double sa_eff =  minuit.getParameterError("EfficiencySa");
321    double dsa_eff = minuit.getParameterError("EfficiencySa", sa_eff);
322    double iso_eff =  minuit.getParameterError("EfficiencyIso");
323    double diso_eff = minuit.getParameterError("EfficiencyIso", iso_eff);
324    double hlt_eff =  minuit.getParameterError("EfficiencyHLT");
325    double dhlt_eff = minuit.getParameterError("EfficiencyHLT",hlt_eff);
326    myfile<< Y <<" "<< dY <<" "<< tk_eff <<" "<< dtk_eff <<" "<< sa_eff << " " << dsa_eff << " " << iso_eff <<" " << diso_eff<< " " << hlt_eff << " " << dhlt_eff << " " <<chi2()/(chi2.numberOfBins()- minuit.numberOfFreeParameters());
327    
328    myfile.close();
329
330    //Plot
331    double s;
332    s = 0;
333    for(int i = 1; i <= histoZMuMuNoIso->GetNbinsX(); ++i)
334      s += histoZMuMuNoIso->GetBinContent(i);
335    histoZMuMuNoIso->SetEntries(s);
336    s = 0;
337    for(int i = 1; i <= histoZMuMu->GetNbinsX(); ++i)
338      s += histoZMuMu->GetBinContent(i);
339    histoZMuMu->SetEntries(s);
340    s = 0;
341    for(int i = 1; i <= histoZMuMu1HLT->GetNbinsX(); ++i)
342      s += histoZMuMu1HLT->GetBinContent(i);
343    histoZMuMu1HLT->SetEntries(s);
344    s = 0;
345    for(int i = 1; i <= histoZMuMu2HLT->GetNbinsX(); ++i)
346      s += histoZMuMu2HLT->GetBinContent(i);
347    histoZMuMu2HLT->SetEntries(s);
348    s = 0;
349    for(int i = 1; i <= histoZMuTk->GetNbinsX(); ++i)
350      s += histoZMuTk->GetBinContent(i);
351    histoZMuTk->SetEntries(s);
352    s = 0;
353    for(int i = 1; i <= histoZMuSa->GetNbinsX(); ++i)
354      s += histoZMuSa->GetBinContent(i);
355    histoZMuSa->SetEntries(s);
356    
357    string ZMuMu1HLTPlot = "ZMuMu1HLTFit_" + plot_string;
358    root::plot<Expr>(ZMuMu1HLTPlot.c_str(), *histoZMuMu1HLT, zMuMu1HLT, fMin, fMax, 
359		     effTk, effSa, effIso, effHLT, yieldZMuMu, 
360		     kOrange-2, 2, kSolid, 100, 
361		     "Z -> #mu #mu mass", "#mu #mu invariant mass (GeV/c^{2})", 
362		     "Events");
363    
364    string ZMuMu2HLTPlot = "ZMuMu2HLTFit_" + plot_string;
365    root::plot<Expr>(ZMuMu2HLTPlot.c_str(), *histoZMuMu2HLT, zMuMu2HLT, fMin, fMax, 
366		     effTk, effSa, effIso, effHLT, yieldZMuMu, 
367		     kOrange-2, 2, kSolid, 100, 
368		     "Z -> #mu #mu mass", "#mu #mu invariant mass (GeV/c^{2})", 
369		     "Events");
370    
371    
372    string ZMuMuNoIsoPlot = "ZMuMuNoIsoFit_X_" + plot_string;
373    root::plot<Expr>(ZMuMuNoIsoPlot.c_str(), *histoZMuMuNoIso, zMuMuNoIso, fMin, fMax, 
374		     effTk, effSa, effIso, effHLT, yieldZMuMu,
375		     yieldBkgZMuMuNotIso, alpha, b0, b1, b2,
376		     kWhite, 2, kSolid, 100, 
377		     "Z -> #mu #mu Not Iso mass", "#mu #mu invariant mass (GeV/c^{2})", 
378		     "Events");	
379    ZMuMuNoIsoPlot = "ZMuMuNoIsoFit_" + plot_string;
380    TF1 funZMuMuNoIso = root::tf1_t<sig_tag, Expr>("ZMuMuNoIsoFunction", zMuMuNoIso, fMin, fMax, 
381					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
382					      yieldBkgZMuMuNotIso, alpha, b0, b1, b2);
383    funZMuMuNoIso.SetLineColor(kOrange+8);
384    funZMuMuNoIso.SetLineWidth(3);
385    //funZMuMuNoIso.SetLineStyle(kDashed);
386    
387    //funZMuMuNoIso.SetFillColor(kOrange-2);
388    //funZMuMuNoIso.SetFillStyle(3325);
389
390    funZMuMuNoIso.SetNpx(10000);
391    TF1 funZMuMuNoIsoBkg = root::tf1_t<bkg_tag, Expr>("ZMuMuNoIsoBack", zMuMuNoIsoBkgScaled, fMin, fMax, 
392						 yieldBkgZMuMuNotIso, alpha, b0, b1, b2);
393    funZMuMuNoIsoBkg.SetLineColor(kViolet+3);
394    funZMuMuNoIsoBkg.SetLineWidth(2);
395    funZMuMuNoIsoBkg.SetLineStyle(kSolid);
396    funZMuMuNoIsoBkg.SetFillColor(kViolet-5);
397    funZMuMuNoIsoBkg.SetFillStyle(3357);
398
399    funZMuMuNoIsoBkg.SetNpx(10000);
400    histoZMuMuNoIso->SetTitle("Z -> #mu #mu Not Iso mass");
401    histoZMuMuNoIso->SetXTitle("#mu +  #mu invariant mass (GeV/c^{2})");
402    histoZMuMuNoIso->SetYTitle("Events");
403    TCanvas *canvas = new TCanvas("canvas");
404    histoZMuMuNoIso->Draw("e");
405    funZMuMuNoIsoBkg.Draw("same");
406    funZMuMuNoIso.Draw("same");
407    canvas->SaveAs(ZMuMuNoIsoPlot.c_str());
408    canvas->SetLogy();
409    string logZMuMuNoIsoPlot = "log_" + ZMuMuNoIsoPlot;
410    canvas->SaveAs(logZMuMuNoIsoPlot.c_str());
411    
412    double IntSigMMNotIso = ((double) rebinMuMu/  (double) rebinMuMuNoIso) * funZMuMuNoIso.Integral(fMin, fMax);
413    double IntSigMMNotIsoBkg = ((double) rebinMuMu/  (double) rebinMuMuNoIso) * funZMuMuNoIsoBkg.Integral(fMin, fMax);
414    cout << "*********  ZMuMuNoIsoPlot signal yield from the fit ==> " <<  IntSigMMNotIso << endl;
415    cout << "*********  ZMuMuNoIsoPlot background yield from the fit ==> " << IntSigMMNotIsoBkg << endl;
416
417
418
419    string ZMuTkPlot = "ZMuTkFit_X_" + plot_string;
420    root::plot<Expr>(ZMuTkPlot.c_str(), *histoZMuTk, zMuTk, fMin, fMax,
421		     effTk, effSa, effIso, effHLT, yieldZMuMu,
422		     yieldBkgZMuTk, lambda, a0, a1, a2,
423		     kOrange+3, 2, kSolid, 100,
424		     "Z -> #mu + (unmatched) track mass", "#mu #mu invariant mass (GeV/c^{2})",
425		     "Events");
426    ZMuTkPlot = "ZMuTkFit_" + plot_string;
427    TF1 funZMuTk = root::tf1_t<sig_tag, Expr>("ZMuTkFunction", zMuTk, fMin, fMax, 
428					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
429					      yieldBkgZMuTk, lambda, a0, a1, a2);
430    funZMuTk.SetLineColor(kOrange+8);
431    funZMuTk.SetLineWidth(3);
432    funZMuTk.SetLineStyle(kSolid);
433    //  funZMuTk.SetFillColor(kOrange-2);
434    //funZMuTk.SetFillStyle(3325);
435    funZMuTk.SetNpx(10000);
436    TF1 funZMuTkBkg = root::tf1_t<bkg_tag, Expr>("ZMuTkBack", zMuTkBkgScaled, fMin, fMax, 
437						 yieldBkgZMuTk, lambda, a0, a1, a2);
438    funZMuTkBkg.SetLineColor(kViolet+3);
439    funZMuTkBkg.SetLineWidth(2);
440    funZMuTkBkg.SetLineStyle(kSolid);
441    funZMuTkBkg.SetFillColor(kViolet-5);
442    funZMuTkBkg.SetFillStyle(3357);
443    funZMuTkBkg.SetNpx(10000);
444    histoZMuTk->SetTitle("Z -> #mu + (unmatched) track mass");
445    histoZMuTk->SetXTitle("#mu + (unmatched) track invariant mass (GeV/c^{2})");
446    histoZMuTk->SetYTitle("Events");
447    TCanvas *canvas_ = new TCanvas("canvas_");
448    histoZMuTk->Draw("e");
449    funZMuTkBkg.Draw("same");
450    funZMuTk.Draw("same");
451    canvas_->SaveAs(ZMuTkPlot.c_str());
452    canvas_->SetLogy();
453    string logZMuTkPlot = "log_" + ZMuTkPlot;
454    canvas_->SaveAs(logZMuTkPlot.c_str());
455
456
457    double IntSigMT = ((double) rebinMuMu/  (double) rebinMuTk) * funZMuTk.Integral(fMin, fMax);
458    double IntSigMTBkg = ((double) rebinMuMu/  (double) rebinMuTk) * funZMuTkBkg.Integral(fMin, fMax);
459    cout << "*********  ZMuMuTkPlot signal yield from the fit ==> " <<  IntSigMT << endl;
460    cout << "*********  ZMuMuTkPlot background yield from the fit ==> " << IntSigMTBkg << endl;
461
462
463
464    string ZMuSaPlot = "ZMuSaFit_X_" + plot_string;
465    root::plot<Expr>(ZMuSaPlot.c_str(), *histoZMuSa, zMuSa, fMin, fMax, 
466		     effSa, effTk, effIso, yieldZMuMu, 
467		     yieldBkgZMuSa, beta, c0, c1, c2 ,
468		     kOrange+3, 2, kSolid, 100, 
469		     "Z -> #mu + (unmatched) standalone mass", 
470		     "#mu + (unmatched) standalone invariant mass (GeV/c^{2})", 
471		     "Events");
472
473
474
475    ZMuSaPlot = "ZMuSaFit_" + plot_string;
476    TF1 funZMuSa = root::tf1_t<sig_tag, Expr>("ZMuSaFunction", zMuSa, fMin, fMax, 
477					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
478					      yieldBkgZMuSa, beta, c0, c1, c2);
479    funZMuSa.SetLineColor(kOrange+8);
480    funZMuSa.SetLineWidth(3);
481    funZMuSa.SetLineStyle(kSolid);
482    // funZMuSa.SetFillColor(kOrange-2);
483    // funZMuSa.SetFillStyle(3325);
484    funZMuSa.SetNpx(10000);
485    TF1 funZMuSaBkg = root::tf1_t<bkg_tag, Expr>("ZMuSaBack", zMuSaBkgScaled, fMin, fMax, 
486						 yieldBkgZMuSa, beta, c0, c1, c2);
487    funZMuSaBkg.SetLineColor(kViolet+3);
488    funZMuSaBkg.SetLineWidth(2);
489    funZMuSaBkg.SetLineStyle(kSolid);
490    funZMuSaBkg.SetFillColor(kViolet-5);
491    funZMuSaBkg.SetFillStyle(3357);
492    funZMuSaBkg.SetNpx(10000);
493    histoZMuSa->SetTitle("Z -> #mu + (unmatched) standalone mass");
494    histoZMuSa->SetXTitle("#mu + (unmatched) standalone invariant mass (GeV/c^{2})");
495    histoZMuSa->SetYTitle("Events");
496    TCanvas *canvas__ = new TCanvas("canvas__");
497    histoZMuSa->Draw("e");
498    funZMuSaBkg.Draw("same");
499    funZMuSa.Draw("same");
500    canvas__->SaveAs(ZMuSaPlot.c_str());
501    canvas__->SetLogy();
502    string logZMuSaPlot = "log_" + ZMuSaPlot;
503    canvas__->SaveAs(logZMuSaPlot.c_str());
504
505    double IntSigMS = ((double) rebinMuMu/  (double) rebinMuSa) * funZMuSa.Integral(fMin, fMax);
506    double IntSigMSBkg = ((double) rebinMuMu/  (double) rebinMuSa) * funZMuSaBkg.Integral(fMin, fMax);
507    cout << "*********  ZMuMuSaPlot signal yield from the fit ==> " <<  IntSigMS << endl;
508    cout << "*********  ZMuMuSaPlot background yield from the fit ==> " << IntSigMSBkg << endl;
509
510
511  }
512  return 0;
513}
514
515#include <boost/program_options.hpp>
516using namespace boost;
517namespace po = boost::program_options;
518
519int main(int ac, char *av[]) {
520  po::options_description desc("Allowed options");
521  desc.add_options()
522    ("help,h", "produce help message")
523    ("input-file,i", po::value<vector<string> >(), "input file")
524    ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range")
525    ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range")
526    ("rebins,R", po::value<vector<unsigned int> >(), "rebins values: rebinMuMu2HLT , rebinMuMu1HLT , rebinMuMuNoIso , rebinMuSa, rebinMuTk")
527    ("chi2,c", "perform chi-squared fit")
528    ("plr,p", "perform Poisson likelihood-ratio fit")
529    ("nonIsoTemplateFromMC,I", po::value<bool>(&nonIsoTemplateFromMC)->default_value(false) , "take the template for nonIso sample from MC")
530    ("plot-format,f", po::value<string>(&ext)->default_value("eps"), "output plot format")
531    ("detectorRegion,r",po::value<string> (&region)->default_value("all"), "detector region in which muons are detected" );   
532  po::positional_options_description p;
533  p.add("input-file", -1);
534  p.add("rebins", -1);
535
536  
537  po::variables_map vm;
538  po::store(po::command_line_parser(ac, av).
539	    options(desc).positional(p).run(), vm);
540  po::notify(vm);
541  
542  if (vm.count("help")) {
543    cout << "Usage: options_description [options]\n";
544    cout << desc;
545    return 0;
546  }
547  
548  if (!vm.count("input-file")) {
549    return 1;
550  }
551  cout << "Input files are: " 
552       << vm["input-file"].as< vector<string> >() << "\n";
553  vector<string> v_file = vm["input-file"].as< vector<string> >();
554
555
556  if (vm.count("rebins") ) {
557    vector<unsigned int> v_rebin = vm["rebins"].as< vector<unsigned int> >();
558    if (v_rebin.size()!=5){
559      cerr << " please provide 5 numbers in the given order:  rebinMuMu2HLT , rebinMuMu1HLT , rebinMuMuNoIso, rebinMuSa, rebinMuTk \n";
560      return 1;
561    }
562 rebinMuMuNoIso = v_rebin[2], rebinMuMu1HLT = v_rebin[1], rebinMuMu2HLT = v_rebin[0], rebinMuTk = v_rebin[4], rebinMuSa = v_rebin[3];
563  }
564
565
566
567
568  bool chi2Fit = vm.count("chi2"), plrFit = vm.count("plr");
569  
570  if(!(chi2Fit||plrFit))
571    cerr << "Warning: no fit performed. Please, specify either -c or -p options or both" << endl;
572 
573
574  gROOT->SetStyle("Plain");
575
576  int ret = 0;
577  try {
578    if(plrFit) {
579      std::cout << "==================================== " << std::endl;
580      std::cout << "=== Poisson Likelihood Ratio fit === " << std::endl;
581      std::cout << "==================================== " << std::endl;
582      int ret2 = main_t<ExprPLR>(v_file);
583      if(ret2 != 0) ret = 1;
584    }
585    if(chi2Fit) {
586      std::cout << "================= " << std::endl;
587      std::cout << "=== Chi-2 fit === " << std::endl;
588      std::cout << "================= " << std::endl;
589      int ret1 = main_t<ExprChi2>(v_file);
590      if(ret1 != 0) ret = 1;
591    }
592  }
593  catch(std::exception& e) {
594    cerr << "error: " << e.what() << "\n";
595    return 1;
596  }
597  catch(...) {
598    cerr << "Exception of unknown type!\n";
599  }
600  return ret;
601}
602