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

  2. #include "PhysicsTools/Utilities/interface/Parameter.h"
    3. 

  3. #include "PhysicsTools/Utilities/interface/Gaussian.h"
    4. 

  4. #include "PhysicsTools/Utilities/interface/Numerical.h"
    5. 

  5. #include "PhysicsTools/Utilities/interface/Exponential.h"
    6. 

  6. #include "PhysicsTools/Utilities/interface/Polynomial.h"
    7. 

  7. #include "PhysicsTools/Utilities/interface/Constant.h"
    8. 

  8. #include "PhysicsTools/Utilities/interface/Operations.h"
    9. 

  9. #include "PhysicsTools/Utilities/interface/MultiHistoChiSquare.h"
    10. 

  10. #include "PhysicsTools/Utilities/interface/HistoChiSquare.h"
    11. 

  11. #include "PhysicsTools/Utilities/interface/HistoPoissonLikelihoodRatio.h"
    12. 

  12. #include "PhysicsTools/Utilities/interface/RootMinuit.h"
    13. 

  13. #include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
    14. 

  14. #include "PhysicsTools/Utilities/interface/FunctClone.h"
    15. 

  15. #include "PhysicsTools/Utilities/interface/rootPlot.h"
    16. 

  16. #include "PhysicsTools/Utilities/interface/Expression.h"
    17. 

  17. #include "PhysicsTools/Utilities/interface/HistoPdf.h"
    18. 

  18. #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
    19. 

  19. #include "TROOT.h"
    20. 

  20. #include "TSystem.h"
    21. 

  21. #include "TH1.h"
    22. 

  22. #include "TFile.h"
    23. 

  23. #include <fstream>
    24. 

  24. #include <iostream>
    25. 

  25. #include <algorithm> 
    26. 

  26. #include <exception>
    27. 

  27. #include <iterator>
    28. 

  28. #include <string>
    29. 

  29. #include <vector>
    30. 

  30. using namespace std;
    31. 

  31. 

  32. // A helper function to simplify the main part.
    33. 

  33. template<class T>
    34. 

  34. ostream& operator<<(ostream& os, const vector<T>& v) {
    35. 

  35.   copy(v.begin(), v.end(), ostream_iterator<T>(cout, " ")); 
    36. 

  36.   return os;
    37. 

  37. }
    38. 

  38. 

  39. // A function that get histogram and sets contents to 0 
    40. 

  40. // if entries are too small
    41. 

  41. TH1 * getHisto(TFile * file, const char * name, unsigned int rebin) {
    42. 

  42.   TObject * h = file->Get(name);
    43. 

  43.   if(h == 0)
    44. 

  44.     throw edm::Exception(edm::errors::Configuration) 
    45. 

  45.       << "Can't find object " << name << "\n";
    46. 

  46.   TH1 * histo = dynamic_cast<TH1*>(h);
    47. 

  47.   if(histo == 0)
    48. 

  48.     throw edm::Exception(edm::errors::Configuration) 
    49. 

  49.       << "Object " << name << " is of type " << h->ClassName() << ", not TH1\n";
    50. 

  50.   histo->Rebin(rebin);  
    51. 

  51.   for(int i = 1; i <= histo->GetNbinsX(); ++i) {
    52. 

  52.     if(histo->GetBinContent(i) < 0.1) {
    53. 

  53.       histo->SetBinContent(i, 0.0);
    54. 

  54.       histo->SetBinError(i, 0.0);
    55. 

  55.     }
    56. 

  56.   }
    57. 

  57.   return histo;
    58. 

  58. }
    59. 

  59. 

  60. struct sig_tag;
    61. 

  61. struct bkg_tag;
    62. 

  62. 

  63. typedef funct::FunctExpression Expr;
    64. 

  64. typedef fit::HistoChiSquare<funct::FunctExpression> ExprChi2;
    65. 

  65. typedef fit::HistoPoissonLikelihoodRatio<funct::FunctExpression> ExprPLR;
    66. 

  66. 

  67. double fMin, fMax;
    68. 

  68.  unsigned int rebinMuMuNoIso ,rebinMuMu =1 , rebinMuMu1HLT , rebinMuMu2HLT , rebinMuTk , rebinMuSa;
    69. 

  69.  // assume that the bin size is 1 GeV!!! 
    70. 

  70. string ext, region;
    71. 

  71. bool nonIsoTemplateFromMC;
    72. 

  72. 

  73. template<typename T>
    74. 

  74. struct PlotPrefix { };
    75. 

  75. 

  76. template<>
    77. 

  77. struct PlotPrefix<ExprChi2> {
    78. 

  78.   static string str() { return "chi2"; }
    79. 

  79. };
    80. 

  80.  
    81. 

  81. template<>
    82. 

  82. struct PlotPrefix<ExprPLR> {
    83. 

  83.   static string str() { return "plr"; }
    84. 

  84. };
    85. 

  85.  
    86. 

  86. template<typename T>
    87. 

  87. int main_t(const vector<string> & v_file){
    88. 

  88.   typedef fit::MultiHistoChiSquare<T, T, T, T, T> ChiSquared;
    89. 

  89.   fit::RootMinuitCommands<ChiSquared> commands("zChi2Fit.txt");
    90. 

  90.   
    91. 

  91.   cout << "minuit command file completed" << endl;
    92. 

  92.  
    93. 

  93.   funct::Constant rebinMuMuNoIsoConst(rebinMuMuNoIso), rebinMuMuConst(rebinMuMu), 
    94. 

  94.     rebinMuMu1HLTConst(rebinMuMu1HLT), rebinMuMu2HLTConst(rebinMuMu2HLT), 
    95. 

  95.     rebinMuTkConst(rebinMuTk), rebinMuSaConst(rebinMuSa);
    96. 

  96.   
    97. 

  97.   for(vector<string>::const_iterator it = v_file.begin(); it != v_file.end(); ++it) {
    98. 

  98.     TFile * root_file = new TFile(it->c_str(), "read");
    99. 

  99.     
    100. 

  100.     // default when region==all
    101. 

  101.     //    TH1 * histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlots/zMass",rebinMuMuNoIso);
    102. 

  102.     TH1 * histoZMuMuNoIso = getHisto(root_file, "oneNonIsolatedZToMuMuPlots/zMass",rebinMuMuNoIso);
    103. 

  103.     TH1 * histoZMuMu = getHisto(root_file, "goodZToMuMuPlots/zMass",rebinMuMu);
    104. 

  104.     TH1 * histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlots/zMass", rebinMuMu1HLT);
    105. 

  105.     TH1 * histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlots/zMass", rebinMuMu2HLT);
    106. 

  106.     TH1 * histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlots/zMass", rebinMuTk);
    107. 

  107.     TH1 * histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlots/zMass", rebinMuSa);
    108. 

  108.     TH1 * histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogram/zMass", rebinMuSa);
    109. 

  109.    
    110. 

  110.     TH1 * histoZMuMuNoIsoTemplateFromMC= histoZMuMu;
    111. 

  111.     if (nonIsoTemplateFromMC) {
    112. 

  112.       //      histoZMuMuNoIsoTemplateFromMC = getHisto(root_file, "nonIsolatedZToMuMuPlotsMC/zMass",rebinMuMu);
    113. 

  113.       histoZMuMuNoIsoTemplateFromMC = getHisto(root_file, "oneNonIsolatedZToMuMuPlotsMC/zMass",rebinMuMu);
    114. 

  114.   }    
    115. 

  115.     if (region=="barrel"){
    116. 

  116.       histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsBarrel/zMass",rebinMuMuNoIso);
    117. 

  117.       histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsBarrel/zMass",rebinMuMu);
    118. 

  118.       histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsBarrel/zMass", rebinMuMu1HLT);
    119. 

  119.       histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsBarrel/zMass", rebinMuMu2HLT);
    120. 

  120.       histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsBarrel/zMass", rebinMuTk);
    121. 

  121.       histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsBarrel/zMass", rebinMuSa);
    122. 

  122.       histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramBarrel/zMass", rebinMuSa);
    123. 

  123.     }
    124. 

  124.     
    125. 

  125.     if (region=="endcap"){
    126. 

  126.       histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsEndCap/zMass",rebinMuMuNoIso);
    127. 

  127.       histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsEndCap/zMass",rebinMuMu);
    128. 

  128.       histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsEndCap/zMass", rebinMuMu1HLT);
    129. 

  129.       histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsEndCap/zMass", rebinMuMu2HLT);
    130. 

  130.       histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsEndCap/zMass", rebinMuTk);
    131. 

  131.       histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsEndCap/zMass", rebinMuSa);
    132. 

  132.       histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramEndCap/zMass", rebinMuSa);
    133. 

  133.     }
    134. 

  134.     
    135. 

  135.     if (region=="barrend"){
    136. 

  136.       histoZMuMuNoIso = getHisto(root_file, "nonIsolatedZToMuMuPlotsBarrEnd/zMass",rebinMuMuNoIso);
    137. 

  137.       histoZMuMu = getHisto(root_file, "goodZToMuMuPlotsBarrEnd/zMass",rebinMuMu);
    138. 

  138.       histoZMuMu1HLT = getHisto(root_file, "goodZToMuMu1HLTPlotsBarrEnd/zMass", rebinMuMu1HLT);
    139. 

  139.       histoZMuMu2HLT = getHisto(root_file, "goodZToMuMu2HLTPlotsBarrEnd/zMass", rebinMuMu2HLT);
    140. 

  140.       histoZMuTk = getHisto(root_file, "goodZToMuMuOneTrackPlotsBarrEnd/zMass", rebinMuTk);
    141. 

  141.       histoZMuSa = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlotsBarrEnd/zMass", rebinMuSa);
    142. 

  142.       histoZMuSaFromMuMu = getHisto(root_file, "zmumuSaMassHistogramBarrEnd/zMass", rebinMuSa);
    143. 

  143.     }
    144. 

  144.     
    145. 

  145.     if (region!="endcap" && region!="barrel" && region!="barrend" && region!="all"  ){
    146. 

  146.       cout<< "not a valid region selected"<< endl;
    147. 

  147.       cout << "possible choises are: all, barrel, endcap, barrend " << endl;
    148. 

  148.       return 0;
    149. 

  149.     }
    150. 

  150.     
    151. 

  151.     cout << ">>> histogram loaded\n";
    152. 

  152.     string f_string = *it + "_" + PlotPrefix<T>::str() + "_";
    153. 

  153.     replace(f_string.begin(), f_string.end(), '.', '_');
    154. 

  154.     replace(f_string.begin(), f_string.end(), '/', '_');
    155. 

  155.     string plot_string = f_string + "." + ext;
    156. 

  156.     cout << ">>> Input files loaded\n";
    157. 

  157.     
    158. 

  158.     const char * kYieldZMuMu = "YieldZMuMu";
    159. 

  159.     const char * kEfficiencyTk = "EfficiencyTk";
    160. 

  160.     const char * kEfficiencySa = "EfficiencySa";
    161. 

  161.     const char * kEfficiencyIso = "EfficiencyIso";
    162. 

  162.     const char * kEfficiencyHLT = "EfficiencyHLT"; 
    163. 

  163.     const char * kYieldBkgZMuTk = "YieldBkgZMuTk"; 
    164. 

  164.     const char * kYieldBkgZMuSa = "YieldBkgZMuSa"; 
    165. 

  165.     const char * kYieldBkgZMuMuNotIso = "YieldBkgZMuMuNotIso";
    166. 

  166.     const char * kAlpha = "Alpha";
    167. 

  167.     const char * kBeta = "Beta";
    168. 

  168.     const char * kLambda = "Lambda";
    169. 

  169.     const char * kA0 = "A0"; 
    170. 

  170.     const char * kA1 = "A1"; 
    171. 

  171.     const char * kA2 = "A2"; 
    172. 

  172.     const char * kB0 = "B0"; 
    173. 

  173.     const char * kB1 = "B1"; 
    174. 

  174.     const char * kB2 = "B2"; 
    175. 

  175.     const char * kC0 = "C0"; 
    176. 

  176.     const char * kC1 = "C1"; 
    177. 

  177.     const char * kC2 = "C2"; 
    178. 

  178.     
    179. 

  179.     funct::Parameter yieldZMuMu(kYieldZMuMu, commands.par(kYieldZMuMu));
    180. 

  180.     funct::Parameter effTk(kEfficiencyTk, commands.par(kEfficiencyTk)); 
    181. 

  181.     funct::Parameter effSa(kEfficiencySa, commands.par(kEfficiencySa)); 
    182. 

  182.     funct::Parameter effIso(kEfficiencyIso, commands.par(kEfficiencyIso)); 
    183. 

  183.     funct::Parameter effHLT(kEfficiencyHLT, commands.par(kEfficiencyHLT)); 
    184. 

  184.     funct::Parameter yieldBkgZMuTk(kYieldBkgZMuTk, commands.par(kYieldBkgZMuTk));
    185. 

  185.     funct::Parameter yieldBkgZMuSa(kYieldBkgZMuSa, commands.par(kYieldBkgZMuSa));
    186. 

  186.     funct::Parameter yieldBkgZMuMuNotIso(kYieldBkgZMuMuNotIso, commands.par(kYieldBkgZMuMuNotIso));
    187. 

  187.     funct::Parameter alpha(kAlpha, commands.par(kAlpha));
    188. 

  188.     funct::Parameter beta(kBeta, commands.par(kBeta));
    189. 

  189.     funct::Parameter lambda(kLambda, commands.par(kLambda));
    190. 

  190.     funct::Parameter a0(kA0, commands.par(kA0));
    191. 

  191.     funct::Parameter a1(kA1, commands.par(kA1));
    192. 

  192.     funct::Parameter a2(kA2, commands.par(kA2));
    193. 

  193.     funct::Parameter b0(kB0, commands.par(kB0));
    194. 

  194.     funct::Parameter b1(kB1, commands.par(kB1));
    195. 

  195.     funct::Parameter b2(kB2, commands.par(kB2));
    196. 

  196.     funct::Parameter c0(kC0, commands.par(kC0));
    197. 

  197.     funct::Parameter c1(kC1, commands.par(kC1));
    198. 

  198.     funct::Parameter c2(kC2, commands.par(kC2));
    199. 

  199.     funct::Constant cFMin(fMin), cFMax(fMax);
    200. 

  200.     
    201. 

  201.     // count ZMuMu Yield
    202. 

  202.     double nZMuMu = 0, nZMuMu1HLT = 0, nZMuMu2HLT = 0;
    203. 

  203.     {
    204. 

  204.       unsigned int nBins = histoZMuMu->GetNbinsX();
    205. 

  205.       double xMin = histoZMuMu->GetXaxis()->GetXmin();
    206. 

  206.       double xMax = histoZMuMu->GetXaxis()->GetXmax();
    207. 

  207.       double deltaX =(xMax - xMin) / nBins;
    208. 

  208.       for(unsigned int i = 0; i < nBins; ++i) { 
    209. 

  209. 	double x = xMin + (i +.5) * deltaX;
    210. 

  210. 	if(x > fMin && x < fMax){
    211. 

  211. 	  nZMuMu += histoZMuMu->GetBinContent(i+1);
    212. 

  212. 	  nZMuMu1HLT += histoZMuMu1HLT->GetBinContent(i+1);
    213. 

  213. 	  nZMuMu2HLT += histoZMuMu2HLT->GetBinContent(i+1);
    214. 

  214. 	}
    215. 

  215.       }
    216. 

  216.     }
    217. 

  217.     // aggiungi 1HLT 2HLT
    218. 

  218.     cout << ">>> count of ZMuMu yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu << endl;
    219. 

  219.     cout << ">>> count of ZMuMu (1HLT) yield in the range [" << fMin << ", " << fMax << "]: "  << nZMuMu1HLT << endl;
    220. 

  220.     cout << ">>> count of ZMuMu (2HLT) yield in the range [" << fMin << ", " << fMax << "]: "  << nZMuMu2HLT << endl;
    221. 

  221.     funct::RootHistoPdf zPdfMuMu(*histoZMuMu, fMin, fMax);
    222. 

  222.       //assign ZMuMu as pdf
    223. 

  223.     funct::RootHistoPdf zPdfMuMuNonIso = zPdfMuMu;
    224. 

  224.     if (nonIsoTemplateFromMC) {
    225. 

  225.       funct::RootHistoPdf zPdfMuMuNoIsoFromMC(*histoZMuMuNoIsoTemplateFromMC, fMin, fMax);
    226. 

  226.       zPdfMuMuNonIso = zPdfMuMuNoIsoFromMC;
    227. 

  227.     }    
    228. 

  228.     
    229. 

  229.     funct::RootHistoPdf zPdfMuTk = zPdfMuMu;
    230. 

  230.     funct::RootHistoPdf zPdfMuMu1HLT = zPdfMuMu;
    231. 

  231.     funct::RootHistoPdf zPdfMuMu2HLT = zPdfMuMu;
    232. 

  232.     funct::RootHistoPdf zPdfMuSa(*histoZMuSaFromMuMu, fMin, fMax);
    233. 

  233.     zPdfMuMuNonIso.rebin(rebinMuMuNoIso/rebinMuMu);
    234. 

  234.     zPdfMuTk.rebin(rebinMuTk/rebinMuMu);
    235. 

  235.     zPdfMuMu1HLT.rebin(rebinMuMu1HLT/rebinMuMu);
    236. 

  236.     zPdfMuMu2HLT.rebin(rebinMuMu2HLT/rebinMuMu);
    237. 

  237.     
    238. 

  238.     funct::Numerical<2> _2;
    239. 

  239.     funct::Numerical<1> _1;
    240. 

  240.     
    241. 

  241.     //Efficiency term
    242. 

  242.     Expr zMuMuEff1HLTTerm = _2 * (effTk ^ _2) *  (effSa ^ _2) * (effIso ^ _2) * effHLT * (_1 - effHLT); 
    243. 

  243.     Expr zMuMuEff2HLTTerm = (effTk ^ _2) *  (effSa ^ _2) * (effIso ^ _2) * (effHLT ^ _2) ; 
    244. 

  244.     //    Expr zMuMuNoIsoEffTerm = (effTk ^ _2) * (effSa ^ _2) * (_1 - (effIso ^ _2)) * (_1 - ((_1 - effHLT)^_2));
    245. 

  245.     // change to both hlt and one not iso
    246. 

  246.     Expr zMuMuNoIsoEffTerm = _2 * (effTk ^ _2) * (effSa ^ _2) * effIso * (_1 - effIso) * (effHLT^_2);
    247. 

  247.     Expr zMuTkEffTerm = _2 * (effTk ^ _2) * effSa * (_1 - effSa) * (effIso ^ _2) * effHLT;
    248. 

  248.     Expr zMuSaEffTerm = _2 * (effSa ^ _2) * effTk * (_1 - effTk) * (effIso ^ _2) * effHLT;
    249. 

  249.     
    250. 

  250.     Expr zMuMu1HLT = rebinMuMu1HLTConst * zMuMuEff1HLTTerm * yieldZMuMu;
    251. 

  251.     Expr zMuMu2HLT = rebinMuMu2HLTConst * zMuMuEff2HLTTerm * yieldZMuMu;
    252. 

  252.     
    253. 

  253.     Expr zMuTkBkg = yieldBkgZMuTk * funct::Exponential(lambda)* funct::Polynomial<2>(a0, a1, a2);
    254. 

  254.     Expr zMuTkBkgScaled = rebinMuTkConst * zMuTkBkg;
    255. 

  255.     Expr zMuTk = rebinMuTkConst * (zMuTkEffTerm * yieldZMuMu * zPdfMuTk + zMuTkBkg);
    256. 

  256.     
    257. 

  257.     Expr zMuMuNoIsoBkg = yieldBkgZMuMuNotIso * funct::Exponential(alpha)* funct::Polynomial<2>(b0, b1, b2);
    258. 

  258.     Expr zMuMuNoIsoBkgScaled = rebinMuMuNoIsoConst * zMuMuNoIsoBkg;
    259. 

  259.     Expr zMuMuNoIso = rebinMuMuNoIsoConst * (zMuMuNoIsoEffTerm * yieldZMuMu * zPdfMuMuNonIso + zMuMuNoIsoBkg);
    260. 

  260.     
    261. 

  261. 

  262.     Expr zMuSaBkg = yieldBkgZMuSa * funct::Exponential(beta)* funct::Polynomial<2>(c0, c1, c2);
    263. 

  263.     Expr zMuSaBkgScaled = rebinMuSaConst * zMuSaBkg;
    264. 

  264.     Expr zMuSa = rebinMuSaConst * (zMuSaEffTerm * yieldZMuMu * zPdfMuSa  + zMuSaBkg );
    265. 

  265.     
    266. 

  266.     TH1D histoZCount1HLT("histoZCount1HLT", "", 1, fMin, fMax);
    267. 

  267.     histoZCount1HLT.Fill(100, nZMuMu1HLT);
    268. 

  268.     TH1D histoZCount2HLT("histoZCount2HLT", "", 1, fMin, fMax);
    269. 

  269.     histoZCount2HLT.Fill(100, nZMuMu2HLT);
    270. 

  270.     
    271. 

  271.     ChiSquared chi2(zMuMu1HLT, & histoZCount1HLT,
    272. 

  272. 		    zMuMu2HLT, & histoZCount2HLT,
    273. 

  273. 		    zMuTk, histoZMuTk, 
    274. 

  274. 		    zMuSa, histoZMuSa, 
    275. 

  275. 		    zMuMuNoIso,histoZMuMuNoIso,
    276. 

  276. 		    fMin, fMax);
    277. 

  277.     cout << "N. bins: " << chi2.numberOfBins() << endl;
    278. 

  278.     
    279. 

  279.     fit::RootMinuit<ChiSquared> minuit(chi2, true);
    280. 

  280.     commands.add(minuit, yieldZMuMu);
    281. 

  281.     commands.add(minuit, effTk);
    282. 

  282.     commands.add(minuit, effSa);
    283. 

  283.     commands.add(minuit, effIso);
    284. 

  284.     commands.add(minuit, effHLT);
    285. 

  285.     commands.add(minuit, yieldBkgZMuTk);
    286. 

  286.     commands.add(minuit, yieldBkgZMuSa);
    287. 

  287.     commands.add(minuit, yieldBkgZMuMuNotIso);
    288. 

  288.     commands.add(minuit, lambda);
    289. 

  289.     commands.add(minuit, alpha);
    290. 

  290.     commands.add(minuit, beta);
    291. 

  291.     commands.add(minuit, a0);
    292. 

  292.     commands.add(minuit, a1);
    293. 

  293.     commands.add(minuit, a2);
    294. 

  294.     commands.add(minuit, b0);
    295. 

  295.     commands.add(minuit, b1);
    296. 

  296.     commands.add(minuit, b2);
    297. 

  297.     commands.add(minuit, c0);
    298. 

  298.     commands.add(minuit, c1);
    299. 

  299.     commands.add(minuit, c2);
    300. 

  300.     commands.run(minuit);
    301. 

  301.     const unsigned int nPar = 20;//WARNIG: this must be updated manually for now
    302. 

  302.     ROOT::Math::SMatrix<double, nPar, nPar, ROOT::Math::MatRepSym<double, nPar> > err;
    303. 

  303.     minuit.getErrorMatrix(err);
    304. 

  304.     
    305. 

  305.     std::cout << "error matrix:" << std::endl;
    306. 

  306.     for(unsigned int i = 0; i < nPar; ++i) {
    307. 

  307.       for(unsigned int j = 0; j < nPar; ++j) {
    308. 

  308. 	  std::cout << err(i, j) << "\t";
    309. 

  309.       }
    310. 

  310.       std::cout << std::endl;
    311. 

  311.     } 
    312. 

  312.     minuit.printFitResults();
    313. 

  313.     ofstream myfile;
    314. 

  314.     myfile.open ("fitResult.txt", ios::out | ios::app);
    315. 

  315.     myfile<<"\n";
    316. 

  316.     double Y =  minuit.getParameterError("YieldZMuMu");
    317. 

  317.     double dY = minuit.getParameterError("YieldZMuMu", Y);
    318. 

  318.     double tk_eff =  minuit.getParameterError("EfficiencyTk");
    319. 

  319.     double dtk_eff = minuit.getParameterError("EfficiencyTk", tk_eff);
    320. 

  320.     double sa_eff =  minuit.getParameterError("EfficiencySa");
    321. 

  321.     double dsa_eff = minuit.getParameterError("EfficiencySa", sa_eff);
    322. 

  322.     double iso_eff =  minuit.getParameterError("EfficiencyIso");
    323. 

  323.     double diso_eff = minuit.getParameterError("EfficiencyIso", iso_eff);
    324. 

  324.     double hlt_eff =  minuit.getParameterError("EfficiencyHLT");
    325. 

  325.     double dhlt_eff = minuit.getParameterError("EfficiencyHLT",hlt_eff);
    326. 

  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. 

  327.     
    328. 

  328.     myfile.close();
    329. 

  329. 

  330.     //Plot
    331. 

  331.     double s;
    332. 

  332.     s = 0;
    333. 

  333.     for(int i = 1; i <= histoZMuMuNoIso->GetNbinsX(); ++i)
    334. 

  334.       s += histoZMuMuNoIso->GetBinContent(i);
    335. 

  335.     histoZMuMuNoIso->SetEntries(s);
    336. 

  336.     s = 0;
    337. 

  337.     for(int i = 1; i <= histoZMuMu->GetNbinsX(); ++i)
    338. 

  338.       s += histoZMuMu->GetBinContent(i);
    339. 

  339.     histoZMuMu->SetEntries(s);
    340. 

  340.     s = 0;
    341. 

  341.     for(int i = 1; i <= histoZMuMu1HLT->GetNbinsX(); ++i)
    342. 

  342.       s += histoZMuMu1HLT->GetBinContent(i);
    343. 

  343.     histoZMuMu1HLT->SetEntries(s);
    344. 

  344.     s = 0;
    345. 

  345.     for(int i = 1; i <= histoZMuMu2HLT->GetNbinsX(); ++i)
    346. 

  346.       s += histoZMuMu2HLT->GetBinContent(i);
    347. 

  347.     histoZMuMu2HLT->SetEntries(s);
    348. 

  348.     s = 0;
    349. 

  349.     for(int i = 1; i <= histoZMuTk->GetNbinsX(); ++i)
    350. 

  350.       s += histoZMuTk->GetBinContent(i);
    351. 

  351.     histoZMuTk->SetEntries(s);
    352. 

  352.     s = 0;
    353. 

  353.     for(int i = 1; i <= histoZMuSa->GetNbinsX(); ++i)
    354. 

  354.       s += histoZMuSa->GetBinContent(i);
    355. 

  355.     histoZMuSa->SetEntries(s);
    356. 

  356.     
    357. 

  357.     string ZMuMu1HLTPlot = "ZMuMu1HLTFit_" + plot_string;
    358. 

  358.     root::plot<Expr>(ZMuMu1HLTPlot.c_str(), *histoZMuMu1HLT, zMuMu1HLT, fMin, fMax, 
    359. 

  359. 		     effTk, effSa, effIso, effHLT, yieldZMuMu, 
    360. 

  360. 		     kOrange-2, 2, kSolid, 100, 
    361. 

  361. 		     "Z -> #mu #mu mass", "#mu #mu invariant mass (GeV/c^{2})", 
    362. 

  362. 		     "Events");
    363. 

  363.     
    364. 

  364.     string ZMuMu2HLTPlot = "ZMuMu2HLTFit_" + plot_string;
    365. 

  365.     root::plot<Expr>(ZMuMu2HLTPlot.c_str(), *histoZMuMu2HLT, zMuMu2HLT, fMin, fMax, 
    366. 

  366. 		     effTk, effSa, effIso, effHLT, yieldZMuMu, 
    367. 

  367. 		     kOrange-2, 2, kSolid, 100, 
    368. 

  368. 		     "Z -> #mu #mu mass", "#mu #mu invariant mass (GeV/c^{2})", 
    369. 

  369. 		     "Events");
    370. 

  370.     
    371. 

  371.     
    372. 

  372.     string ZMuMuNoIsoPlot = "ZMuMuNoIsoFit_X_" + plot_string;
    373. 

  373.     root::plot<Expr>(ZMuMuNoIsoPlot.c_str(), *histoZMuMuNoIso, zMuMuNoIso, fMin, fMax, 
    374. 

  374. 		     effTk, effSa, effIso, effHLT, yieldZMuMu,
    375. 

  375. 		     yieldBkgZMuMuNotIso, alpha, b0, b1, b2,
    376. 

  376. 		     kWhite, 2, kSolid, 100, 
    377. 

  377. 		     "Z -> #mu #mu Not Iso mass", "#mu #mu invariant mass (GeV/c^{2})", 
    378. 

  378. 		     "Events");	
    379. 

  379.     ZMuMuNoIsoPlot = "ZMuMuNoIsoFit_" + plot_string;
    380. 

  380.     TF1 funZMuMuNoIso = root::tf1_t<sig_tag, Expr>("ZMuMuNoIsoFunction", zMuMuNoIso, fMin, fMax, 
    381. 

  381. 					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
    382. 

  382. 					      yieldBkgZMuMuNotIso, alpha, b0, b1, b2);
    383. 

  383.     funZMuMuNoIso.SetLineColor(kOrange+8);
    384. 

  384.     funZMuMuNoIso.SetLineWidth(3);
    385. 

  385.     //funZMuMuNoIso.SetLineStyle(kDashed);
    386. 

  386.     
    387. 

  387.     //funZMuMuNoIso.SetFillColor(kOrange-2);
    388. 

  388.     //funZMuMuNoIso.SetFillStyle(3325);
    389. 

  389. 

  390.     funZMuMuNoIso.SetNpx(10000);
    391. 

  391.     TF1 funZMuMuNoIsoBkg = root::tf1_t<bkg_tag, Expr>("ZMuMuNoIsoBack", zMuMuNoIsoBkgScaled, fMin, fMax, 
    392. 

  392. 						 yieldBkgZMuMuNotIso, alpha, b0, b1, b2);
    393. 

  393.     funZMuMuNoIsoBkg.SetLineColor(kViolet+3);
    394. 

  394.     funZMuMuNoIsoBkg.SetLineWidth(2);
    395. 

  395.     funZMuMuNoIsoBkg.SetLineStyle(kSolid);
    396. 

  396.     funZMuMuNoIsoBkg.SetFillColor(kViolet-5);
    397. 

  397.     funZMuMuNoIsoBkg.SetFillStyle(3357);
    398. 

  398. 

  399.     funZMuMuNoIsoBkg.SetNpx(10000);
    400. 

  400.     histoZMuMuNoIso->SetTitle("Z -> #mu #mu Not Iso mass");
    401. 

  401.     histoZMuMuNoIso->SetXTitle("#mu +  #mu invariant mass (GeV/c^{2})");
    402. 

  402.     histoZMuMuNoIso->SetYTitle("Events");
    403. 

  403.     TCanvas *canvas = new TCanvas("canvas");
    404. 

  404.     histoZMuMuNoIso->Draw("e");
    405. 

  405.     funZMuMuNoIsoBkg.Draw("same");
    406. 

  406.     funZMuMuNoIso.Draw("same");
    407. 

  407.     canvas->SaveAs(ZMuMuNoIsoPlot.c_str());
    408. 

  408.     canvas->SetLogy();
    409. 

  409.     string logZMuMuNoIsoPlot = "log_" + ZMuMuNoIsoPlot;
    410. 

  410.     canvas->SaveAs(logZMuMuNoIsoPlot.c_str());
    411. 

  411.     
    412. 

  412.     double IntSigMMNotIso = ((double) rebinMuMu/  (double) rebinMuMuNoIso) * funZMuMuNoIso.Integral(fMin, fMax);
    413. 

  413.     double IntSigMMNotIsoBkg = ((double) rebinMuMu/  (double) rebinMuMuNoIso) * funZMuMuNoIsoBkg.Integral(fMin, fMax);
    414. 

  414.     cout << "*********  ZMuMuNoIsoPlot signal yield from the fit ==> " <<  IntSigMMNotIso << endl;
    415. 

  415.     cout << "*********  ZMuMuNoIsoPlot background yield from the fit ==> " << IntSigMMNotIsoBkg << endl;
    416. 

  416. 

  417. 

  418. 

  419.     string ZMuTkPlot = "ZMuTkFit_X_" + plot_string;
    420. 

  420.     root::plot<Expr>(ZMuTkPlot.c_str(), *histoZMuTk, zMuTk, fMin, fMax,
    421. 

  421. 		     effTk, effSa, effIso, effHLT, yieldZMuMu,
    422. 

  422. 		     yieldBkgZMuTk, lambda, a0, a1, a2,
    423. 

  423. 		     kOrange+3, 2, kSolid, 100,
    424. 

  424. 		     "Z -> #mu + (unmatched) track mass", "#mu #mu invariant mass (GeV/c^{2})",
    425. 

  425. 		     "Events");
    426. 

  426.     ZMuTkPlot = "ZMuTkFit_" + plot_string;
    427. 

  427.     TF1 funZMuTk = root::tf1_t<sig_tag, Expr>("ZMuTkFunction", zMuTk, fMin, fMax, 
    428. 

  428. 					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
    429. 

  429. 					      yieldBkgZMuTk, lambda, a0, a1, a2);
    430. 

  430.     funZMuTk.SetLineColor(kOrange+8);
    431. 

  431.     funZMuTk.SetLineWidth(3);
    432. 

  432.     funZMuTk.SetLineStyle(kSolid);
    433. 

  433.     //  funZMuTk.SetFillColor(kOrange-2);
    434. 

  434.     //funZMuTk.SetFillStyle(3325);
    435. 

  435.     funZMuTk.SetNpx(10000);
    436. 

  436.     TF1 funZMuTkBkg = root::tf1_t<bkg_tag, Expr>("ZMuTkBack", zMuTkBkgScaled, fMin, fMax, 
    437. 

  437. 						 yieldBkgZMuTk, lambda, a0, a1, a2);
    438. 

  438.     funZMuTkBkg.SetLineColor(kViolet+3);
    439. 

  439.     funZMuTkBkg.SetLineWidth(2);
    440. 

  440.     funZMuTkBkg.SetLineStyle(kSolid);
    441. 

  441.     funZMuTkBkg.SetFillColor(kViolet-5);
    442. 

  442.     funZMuTkBkg.SetFillStyle(3357);
    443. 

  443.     funZMuTkBkg.SetNpx(10000);
    444. 

  444.     histoZMuTk->SetTitle("Z -> #mu + (unmatched) track mass");
    445. 

  445.     histoZMuTk->SetXTitle("#mu + (unmatched) track invariant mass (GeV/c^{2})");
    446. 

  446.     histoZMuTk->SetYTitle("Events");
    447. 

  447.     TCanvas *canvas_ = new TCanvas("canvas_");
    448. 

  448.     histoZMuTk->Draw("e");
    449. 

  449.     funZMuTkBkg.Draw("same");
    450. 

  450.     funZMuTk.Draw("same");
    451. 

  451.     canvas_->SaveAs(ZMuTkPlot.c_str());
    452. 

  452.     canvas_->SetLogy();
    453. 

  453.     string logZMuTkPlot = "log_" + ZMuTkPlot;
    454. 

  454.     canvas_->SaveAs(logZMuTkPlot.c_str());
    455. 

  455. 

  456. 

  457.     double IntSigMT = ((double) rebinMuMu/  (double) rebinMuTk) * funZMuTk.Integral(fMin, fMax);
    458. 

  458.     double IntSigMTBkg = ((double) rebinMuMu/  (double) rebinMuTk) * funZMuTkBkg.Integral(fMin, fMax);
    459. 

  459.     cout << "*********  ZMuMuTkPlot signal yield from the fit ==> " <<  IntSigMT << endl;
    460. 

  460.     cout << "*********  ZMuMuTkPlot background yield from the fit ==> " << IntSigMTBkg << endl;
    461. 

  461. 

  462. 

  463. 

  464.     string ZMuSaPlot = "ZMuSaFit_X_" + plot_string;
    465. 

  465.     root::plot<Expr>(ZMuSaPlot.c_str(), *histoZMuSa, zMuSa, fMin, fMax, 
    466. 

  466. 		     effSa, effTk, effIso, yieldZMuMu, 
    467. 

  467. 		     yieldBkgZMuSa, beta, c0, c1, c2 ,
    468. 

  468. 		     kOrange+3, 2, kSolid, 100, 
    469. 

  469. 		     "Z -> #mu + (unmatched) standalone mass", 
    470. 

  470. 		     "#mu + (unmatched) standalone invariant mass (GeV/c^{2})", 
    471. 

  471. 		     "Events");
    472. 

  472. 

  473. 

  474. 

  475.     ZMuSaPlot = "ZMuSaFit_" + plot_string;
    476. 

  476.     TF1 funZMuSa = root::tf1_t<sig_tag, Expr>("ZMuSaFunction", zMuSa, fMin, fMax, 
    477. 

  477. 					      effTk, effSa, effIso, effHLT, yieldZMuMu, 
    478. 

  478. 					      yieldBkgZMuSa, beta, c0, c1, c2);
    479. 

  479.     funZMuSa.SetLineColor(kOrange+8);
    480. 

  480.     funZMuSa.SetLineWidth(3);
    481. 

  481.     funZMuSa.SetLineStyle(kSolid);
    482. 

  482.     // funZMuSa.SetFillColor(kOrange-2);
    483. 

  483.     // funZMuSa.SetFillStyle(3325);
    484. 

  484.     funZMuSa.SetNpx(10000);
    485. 

  485.     TF1 funZMuSaBkg = root::tf1_t<bkg_tag, Expr>("ZMuSaBack", zMuSaBkgScaled, fMin, fMax, 
    486. 

  486. 						 yieldBkgZMuSa, beta, c0, c1, c2);
    487. 

  487.     funZMuSaBkg.SetLineColor(kViolet+3);
    488. 

  488.     funZMuSaBkg.SetLineWidth(2);
    489. 

  489.     funZMuSaBkg.SetLineStyle(kSolid);
    490. 

  490.     funZMuSaBkg.SetFillColor(kViolet-5);
    491. 

  491.     funZMuSaBkg.SetFillStyle(3357);
    492. 

  492.     funZMuSaBkg.SetNpx(10000);
    493. 

  493.     histoZMuSa->SetTitle("Z -> #mu + (unmatched) standalone mass");
    494. 

  494.     histoZMuSa->SetXTitle("#mu + (unmatched) standalone invariant mass (GeV/c^{2})");
    495. 

  495.     histoZMuSa->SetYTitle("Events");
    496. 

  496.     TCanvas *canvas__ = new TCanvas("canvas__");
    497. 

  497.     histoZMuSa->Draw("e");
    498. 

  498.     funZMuSaBkg.Draw("same");
    499. 

  499.     funZMuSa.Draw("same");
    500. 

  500.     canvas__->SaveAs(ZMuSaPlot.c_str());
    501. 

  501.     canvas__->SetLogy();
    502. 

  502.     string logZMuSaPlot = "log_" + ZMuSaPlot;
    503. 

  503.     canvas__->SaveAs(logZMuSaPlot.c_str());
    504. 

  504. 

  505.     double IntSigMS = ((double) rebinMuMu/  (double) rebinMuSa) * funZMuSa.Integral(fMin, fMax);
    506. 

  506.     double IntSigMSBkg = ((double) rebinMuMu/  (double) rebinMuSa) * funZMuSaBkg.Integral(fMin, fMax);
    507. 

  507.     cout << "*********  ZMuMuSaPlot signal yield from the fit ==> " <<  IntSigMS << endl;
    508. 

  508.     cout << "*********  ZMuMuSaPlot background yield from the fit ==> " << IntSigMSBkg << endl;
    509. 

  509. 

  510. 

  511.   }
    512. 

  512.   return 0;
    513. 

  513. }
    514. 

  514. 

  515. #include <boost/program_options.hpp>
    516. 

  516. using namespace boost;
    517. 

  517. namespace po = boost::program_options;
    518. 

  518. 

  519. int main(int ac, char *av[]) {
    520. 

  520.   po::options_description desc("Allowed options");
    521. 

  521.   desc.add_options()
    522. 

  522.     ("help,h", "produce help message")
    523. 

  523.     ("input-file,i", po::value<vector<string> >(), "input file")
    524. 

  524.     ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range")
    525. 

  525.     ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range")
    526. 

  526.     ("rebins,R", po::value<vector<unsigned int> >(), "rebins values: rebinMuMu2HLT , rebinMuMu1HLT , rebinMuMuNoIso , rebinMuSa, rebinMuTk")
    527. 

  527.     ("chi2,c", "perform chi-squared fit")
    528. 

  528.     ("plr,p", "perform Poisson likelihood-ratio fit")
    529. 

  529.     ("nonIsoTemplateFromMC,I", po::value<bool>(&nonIsoTemplateFromMC)->default_value(false) , "take the template for nonIso sample from MC")
    530. 

  530.     ("plot-format,f", po::value<string>(&ext)->default_value("eps"), "output plot format")
    531. 

  531.     ("detectorRegion,r",po::value<string> (&region)->default_value("all"), "detector region in which muons are detected" );   
    532. 

  532.   po::positional_options_description p;
    533. 

  533.   p.add("input-file", -1);
    534. 

  534.   p.add("rebins", -1);
    535. 

  535. 

  536.   
    537. 

  537.   po::variables_map vm;
    538. 

  538.   po::store(po::command_line_parser(ac, av).
    539. 

  539. 	    options(desc).positional(p).run(), vm);
    540. 

  540.   po::notify(vm);
    541. 

  541.   
    542. 

  542.   if (vm.count("help")) {
    543. 

  543.     cout << "Usage: options_description [options]\n";
    544. 

  544.     cout << desc;
    545. 

  545.     return 0;
    546. 

  546.   }
    547. 

  547.   
    548. 

  548.   if (!vm.count("input-file")) {
    549. 

  549.     return 1;
    550. 

  550.   }
    551. 

  551.   cout << "Input files are: " 
    552. 

  552.        << vm["input-file"].as< vector<string> >() << "\n";
    553. 

  553.   vector<string> v_file = vm["input-file"].as< vector<string> >();
    554. 

  554. 

  555. 

  556.   if (vm.count("rebins") ) {
    557. 

  557.     vector<unsigned int> v_rebin = vm["rebins"].as< vector<unsigned int> >();
    558. 

  558.     if (v_rebin.size()!=5){
    559. 

  559.       cerr << " please provide 5 numbers in the given order:  rebinMuMu2HLT , rebinMuMu1HLT , rebinMuMuNoIso, rebinMuSa, rebinMuTk \n";
    560. 

  560.       return 1;
    561. 

  561.     }
    562. 

  562.  rebinMuMuNoIso = v_rebin[2], rebinMuMu1HLT = v_rebin[1], rebinMuMu2HLT = v_rebin[0], rebinMuTk = v_rebin[4], rebinMuSa = v_rebin[3];
    563. 

  563.   }
    564. 

  564. 

  565. 

  566. 

  567. 

  568.   bool chi2Fit = vm.count("chi2"), plrFit = vm.count("plr");
    569. 

  569.   
    570. 

  570.   if(!(chi2Fit||plrFit))
    571. 

  571.     cerr << "Warning: no fit performed. Please, specify either -c or -p options or both" << endl;
    572. 

  572.  
    573. 

  573. 

  574.   gROOT->SetStyle("Plain");
    575. 

  575. 

  576.   int ret = 0;
    577. 

  577.   try {
    578. 

  578.     if(plrFit) {
    579. 

  579.       std::cout << "==================================== " << std::endl;
    580. 

  580.       std::cout << "=== Poisson Likelihood Ratio fit === " << std::endl;
    581. 

  581.       std::cout << "==================================== " << std::endl;
    582. 

  582.       int ret2 = main_t<ExprPLR>(v_file);
    583. 

  583.       if(ret2 != 0) ret = 1;
    584. 

  584.     }
    585. 

  585.     if(chi2Fit) {
    586. 

  586.       std::cout << "================= " << std::endl;
    587. 

  587.       std::cout << "=== Chi-2 fit === " << std::endl;
    588. 

  588.       std::cout << "================= " << std::endl;
    589. 

  589.       int ret1 = main_t<ExprChi2>(v_file);
    590. 

  590.       if(ret1 != 0) ret = 1;
    591. 

  591.     }
    592. 

  592.   }
    593. 

  593.   catch(std::exception& e) {
    594. 

  594.     cerr << "error: " << e.what() << "\n";
    595. 

  595.     return 1;
    596. 

  596.   }
    597. 

  597.   catch(...) {
    598. 

  598.     cerr << "Exception of unknown type!\n";
    599. 

  599.   }
    600. 

  600.   return ret;
    601. 

  601. }
    602.