/ElectroWeakAnalysis/ZMuMu/bin/zChi2Fit.cpp
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> (®ion)->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