/ElectroWeakAnalysis/ZMuMu/bin/zMuMuRooFit.cpp

https://github.com/aivanov-cern/cmssw · C++ · 548 lines · 402 code · 78 blank · 68 comment · 29 complexity · 6ccd34cc425f9d171d9e9e9498d99584 MD5 · raw file

  1. #include "RooRealVar.h"
  2. #include "RooRealConstant.h"
  3. #include "RooGaussian.h"
  4. #include "RooExponential.h"
  5. #include "RooAddPdf.h"
  6. #include "RooAddition.h"
  7. #include "RooDataSet.h"
  8. #include "RooDataHist.h"
  9. #include "RooHistPdf.h"
  10. #include "RooChebychev.h"
  11. #include "RooExponential.h"
  12. #include "RooProdPdf.h"
  13. #include "RooChi2Var.h"
  14. #include "RooGlobalFunc.h"
  15. #include "RooPlot.h"
  16. #include "RooMinuit.h"
  17. #include "RooFitResult.h"
  18. #include "RooFormulaVar.h"
  19. #include "RooGenericPdf.h"
  20. #include "RooExtendPdf.h"
  21. #include "TCanvas.h"
  22. #include "TROOT.h"
  23. #include "TFile.h"
  24. #include "TH1.h"
  25. #include "RooNLLVar.h"
  26. #include "RooRandom.h"
  27. #include "TRandom3.h"
  28. #include <iostream>
  29. #include <fstream>
  30. #include <boost/program_options.hpp>
  31. using namespace boost;
  32. namespace po = boost::program_options;
  33. #include <vector>
  34. using namespace std;
  35. using namespace RooFit;
  36. // A function that get histogram, restricting it on fit range and sets contents to 0.0 if entries are too small
  37. TH1F * getHisto(TFile * file, const char * name, double fMin, double fMax, unsigned int rebin) {
  38. TObject * h = file->Get(name);
  39. if(h == 0)
  40. cout << "Can't find object " << name << "\n";
  41. TH1F * histo = dynamic_cast<TH1F*>(h);
  42. if(histo == 0)
  43. cout << "Object " << name << " is of type " << h->ClassName() << ", not TH1\n";
  44. TH1F * new_histo = new TH1F(name, name, (int) (fMax-fMin), fMin, fMax);
  45. int bin_num=0;
  46. for (int i = (int)fMin; i <= (int)fMax; ++i ) {
  47. bin_num= (i - (int)fMin + 1);
  48. new_histo->SetBinContent( bin_num, histo->GetBinContent(i) );
  49. }
  50. delete histo;
  51. new_histo->Sumw2();
  52. new_histo->Rebin(rebin);
  53. for(int i = 1; i <= new_histo->GetNbinsX(); ++i) {
  54. if(new_histo->GetBinContent(i) == 0.00) {
  55. cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;
  56. }
  57. if(new_histo->GetBinContent(i) < 0.1) {
  58. new_histo->SetBinContent(i, 0.0);
  59. new_histo->SetBinError(i, 0.0);
  60. cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;
  61. }
  62. }
  63. return new_histo;
  64. }
  65. // a function that create fromm a model pdf a toy RooDataHist
  66. RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model, RooRealVar *var, double ScaleLumi, int range, int rebin, int seed ) {
  67. double genEvents = model->expectedEvents(*var);
  68. TRandom3 *rndm = new TRandom3();
  69. rndm->SetSeed(seed);
  70. double nEvt = rndm->PoissonD( genEvents) ;
  71. int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  72. RooDataSet * data = model->generate(*var , intEvt );
  73. cout<< " expected events for " << name << " = "<< genEvents << endl;
  74. cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  75. // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  76. //cout<< " cast of nEvt for" << name << " == " << intEvt<< endl;
  77. RooAbsData *binned_data = data->binnedClone();
  78. TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin ) );
  79. for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
  80. toy_hist->SetBinError( i, sqrt( toy_hist->GetBinContent(i)) );
  81. if(toy_hist->GetBinContent(i) == 0.00) {
  82. cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;
  83. }
  84. if(toy_hist->GetBinContent(i) < 0.1) {
  85. toy_hist->SetBinContent(i, 0.0);
  86. toy_hist->SetBinError(i, 0.0);
  87. cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;
  88. }
  89. }
  90. RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  91. return toy_rooHist;
  92. }
  93. // a function to create the pdf used for the fit, need the histo model, should be zmm except for zmusta case.....
  94. RooHistPdf * createHistPdf( const char * name, TH1F *model, RooRealVar *var, int rebin){
  95. TH1F * model_clone = new TH1F(*model);
  96. model_clone->Sumw2();
  97. model_clone-> Rebin( rebin );
  98. RooDataHist * model_dataHist = new RooDataHist(name, name , RooArgList(*var), model_clone );
  99. RooHistPdf * HistPdf = new RooHistPdf(name, name, *var, *model_dataHist, 0);
  100. delete model_clone;
  101. return HistPdf;
  102. }
  103. void fit( RooAbsReal & chi2, int numberOfBins, const char * outFileNameWithFitResult ){
  104. TFile * out_root_file = new TFile(outFileNameWithFitResult , "recreate");
  105. RooMinuit m_tot(chi2) ;
  106. m_tot.migrad();
  107. // m_tot.hesse();
  108. RooFitResult* r_chi2 = m_tot.save() ;
  109. cout << "==> Chi2 Fit results " << endl ;
  110. r_chi2->Print("v") ;
  111. // r_chi2->floatParsFinal().Print("v") ;
  112. int NumberOfFreeParameters = r_chi2->floatParsFinal().getSize() ;
  113. for (int i =0; i< NumberOfFreeParameters; ++i){
  114. r_chi2->floatParsFinal()[i].Print();
  115. }
  116. cout<<"chi2:" <<chi2.getVal() << ", numberOfBins: " << numberOfBins << ", NumberOfFreeParameters: " << NumberOfFreeParameters << endl;
  117. cout<<"Normalized Chi2 = " << chi2.getVal()/ (numberOfBins - NumberOfFreeParameters)<<endl;
  118. r_chi2->Write( ) ;
  119. delete out_root_file;
  120. }
  121. int main(int argc, char** argv){
  122. gROOT->SetStyle("Plain");
  123. double fMin , fMax, lumi, scaleLumi =1;
  124. int seed;
  125. Bool_t toy = kFALSE;
  126. Bool_t fitFromData = kFALSE;
  127. string infile, outfile;
  128. int rebinZMuMu =1, rebinZMuSa = 1, rebinZMuTk = 1, rebinZMuMuNoIso=1, rebinZMuMuHlt = 1;
  129. po::options_description desc("Allowed options");
  130. desc.add_options()
  131. ("help,h", "produce help message")
  132. ("toy,t", "toy enabled")
  133. ("seed,s", po::value<int>(&seed)->default_value(34567), "seed value for toy")
  134. ("luminosity,l", po::value<double>(&lumi)->default_value(45.), "luminosity value for toy ")
  135. ("fit,f", "fit from data enabled")
  136. ("rebin,r", po::value< vector<int> > (), "rebin value: r_mutrk r mumuNotIso r_musa r _muhlt")
  137. ("input-file,i", po::value< string> (&infile), "input file")
  138. ("output-file,o", po::value< string> (&outfile), "output file with fit results")
  139. ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range")
  140. ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range");
  141. po::positional_options_description p;
  142. p.add("rebin", -1);
  143. po::variables_map vm;
  144. po::store(po::command_line_parser(argc, argv).
  145. options(desc).positional(p).run(), vm);
  146. po::notify(vm);
  147. if (vm.count("help")) {
  148. cout << "Usage: options_description [options]\n";
  149. cout << desc;
  150. return 0;
  151. }
  152. if (vm.count("toy")) {
  153. cout << "toy enabled with seed " << seed << "\n";
  154. toy = kTRUE;
  155. //RooRandom::randomGenerator()->SetSeed(seed) ;
  156. // lumi should be intented as pb-1 and passed from outside
  157. scaleLumi= lumi / 45.0 ; // 45 is the current lumi correspondent to the histogram.....
  158. }
  159. if (vm.count("fit")) {
  160. cout << "fit from data enabled \n";
  161. fitFromData = kTRUE;
  162. }
  163. if (!vm.count("toy") && !vm.count("fit")) {
  164. cerr << "Choose one beetween fitting form data or with a toy MC \n";
  165. return 1;
  166. }
  167. if ( toy == fitFromData ){
  168. cerr << "Choose if fit from data or with a toy MC \n";
  169. return 1;
  170. }
  171. if (vm.count("rebin") ) {
  172. vector<int> v_rebin = vm["rebin"].as< vector<int> >();
  173. if (v_rebin.size()!=4){
  174. cerr << " please provide 4 numbers in the given order:r_mutrk r mumuNotIso r_musa r _muhlt \n";
  175. return 1;
  176. }
  177. rebinZMuTk = v_rebin[0]; rebinZMuMuNoIso=v_rebin[1]; rebinZMuSa = v_rebin[2];rebinZMuMuHlt = v_rebin[3];
  178. }
  179. RooRealVar * mass = new RooRealVar("mass", "mass (GeV/c^{2})", fMin, fMax );
  180. TFile * root_file = new TFile(infile.c_str(), "read");
  181. int range = (int)(fMax - fMin);
  182. int numberOfBins = range/rebinZMuSa + range/rebinZMuTk + range/rebinZMuMuNoIso + 2* range/rebinZMuMuHlt ;
  183. // zmm histograms used for pdf
  184. TH1F * zmm = getHisto(root_file, "goodZToMuMuPlots/zMass",fMin, fMax, rebinZMuMu );
  185. // zmsta used for pdf
  186. TH1F *zmsta = getHisto(root_file, "zmumuSaMassHistogram/zMass",fMin, fMax, 1); // histogramms to fit.....
  187. TH1F *zmmsta = getHisto(root_file, "goodZToMuMuOneStandAloneMuonPlots/zMass",fMin, fMax, rebinZMuSa/rebinZMuMu);
  188. TH1F *zmt = getHisto(root_file, "goodZToMuMuOneTrackPlots/zMass", fMin, fMax, rebinZMuTk / rebinZMuMu);
  189. TH1F *zmmNotIso = getHisto(root_file,"nonIsolatedZToMuMuPlots/zMass", fMin, fMax, rebinZMuMuNoIso / rebinZMuMu);
  190. TH1F *zmm1hlt = getHisto(root_file, "goodZToMuMu1HLTPlots/zMass", fMin, fMax, rebinZMuMuHlt / rebinZMuMu);
  191. TH1F *zmm2hlt = getHisto(root_file, "goodZToMuMu2HLTPlots/zMass", fMin, fMax, rebinZMuMuHlt / rebinZMuMu);
  192. // creating a pdf for Zmt
  193. RooHistPdf * ZmtPdf = createHistPdf( "ZmtPdf", zmm, mass, rebinZMuTk/ rebinZMuMu );
  194. // creating a pdf for Zmm not iso
  195. RooHistPdf * ZmmNoIsoPdf = createHistPdf( "ZmmNoIsoPdf", zmm, mass, rebinZMuMuNoIso/ rebinZMuMu );
  196. // creating a pdf for Zms from zmsta!!!
  197. RooHistPdf * ZmsPdf = createHistPdf( "ZmsPdf", zmsta, mass, rebinZMuSa/ rebinZMuMu );
  198. // creating a pdf for Zmmhlt
  199. RooHistPdf * ZmmHltPdf = createHistPdf( "ZmmHltPdf", zmm, mass, rebinZMuMuHlt/ rebinZMuMu );
  200. // creating the variable with random init values
  201. RooRealVar Yield("Yield","Yield", 15000.,345.,3567890.) ;
  202. RooRealVar nbkg_mutrk("nbkg_mutrk","background _mutrk fraction",500,0.,100000.) ;
  203. RooRealVar nbkg_mumuNotIso("nbkg_mumuNotIso","background fraction",500,0.,100000.) ;
  204. RooRealVar nbkg_musa("nbkg_musa","background fraction",50,0.,100000.) ;
  205. RooRealVar eff_tk("eff_tk","signal _mutrk fraction",.99,0.8,1.0) ;
  206. RooRealVar eff_sa("eff_sa","eff musta",0.99,0.8,1.0) ;
  207. RooRealVar eff_iso("eff_iso","eff mumuNotIso",.99,0.8,1.0) ;
  208. RooRealVar eff_hlt("eff_hlt","eff 1hlt",0.99, 0.8,1.0) ;
  209. RooRealVar alpha ("alpha","coefficient alpha", -0.01 , -1000, 1000.) ;
  210. RooRealVar a0 ("a0","coefficient 0", 1,-1000.,1000.) ;
  211. RooRealVar a1 ("a1","coefficient 1", -0.001,-1000,1000.) ;
  212. RooRealVar a2 ("a2","coefficient 2", 0.0,-1000.,1000.) ;
  213. RooRealVar beta ("beta","coefficient beta", -0.01,-1000 , 1000. ) ;
  214. RooRealVar b0 ("b0","coefficient 0", 1,-1000.,1000.) ;
  215. RooRealVar b1 ("b1","coefficient 1", -0.001,-1000,1000.) ;
  216. RooRealVar b2("b2","coefficient 2", 0.0,-1000.,1000.) ;
  217. RooRealVar gamma ("gamma","coefficient gamma", -0.01,-1000 , 1000. ) ;
  218. RooRealVar c0 ("c0","coefficient 0", 1,-1000.,1000.) ;
  219. RooRealVar c1 ("c1","coefficient 1", -0.001,-1000,1000.) ;
  220. // fit parameters setted from datacard
  221. filebuf fb;
  222. fb.open ("zMuMuRooFit.txt",ios::in);
  223. istream is(&fb);
  224. char line[1000];
  225. Yield.readFromStream(is.getline(line, 1000), kFALSE);
  226. nbkg_mutrk.readFromStream(is.getline(line,1000), kFALSE);
  227. nbkg_mumuNotIso.readFromStream(is.getline(line,1000), kFALSE);
  228. nbkg_musa.readFromStream(is.getline(line,1000), kFALSE);
  229. eff_tk.readFromStream(is.getline(line,1000), kFALSE);
  230. eff_sa.readFromStream(is.getline(line,1000), kFALSE);
  231. eff_iso.readFromStream(is.getline(line,1000), kFALSE);
  232. eff_hlt.readFromStream(is.getline(line,1000), kFALSE);
  233. alpha.readFromStream(is.getline(line,1000), kFALSE);
  234. a0.readFromStream(is.getline(line,1000), kFALSE);
  235. a1.readFromStream(is.getline(line,1000), kFALSE );
  236. a2.readFromStream(is.getline(line,1000), kFALSE);
  237. beta.readFromStream(is.getline(line,1000), kFALSE);
  238. b0.readFromStream(is.getline(line,1000), kFALSE);
  239. b1.readFromStream(is.getline(line,1000), kFALSE);
  240. b2.readFromStream(is.getline(line,1000), kFALSE);
  241. gamma.readFromStream(is.getline(line,1000), kFALSE);
  242. c0.readFromStream(is.getline(line,1000), kFALSE);
  243. c1.readFromStream(is.getline(line,1000), kFALSE);
  244. fb.close();
  245. // scaling to lumi if toy is enabled...
  246. if (vm.count("toy")) {
  247. Yield.setVal(scaleLumi * (Yield.getVal()));
  248. nbkg_mutrk.setVal(scaleLumi * (nbkg_mutrk.getVal()));
  249. nbkg_mumuNotIso.setVal(scaleLumi * (nbkg_mumuNotIso.getVal()));
  250. }
  251. //efficiency term
  252. //zMuMuEff1HLTTerm = _2 * (effTk ^ _2) * (effSa ^ _2) * (effIso ^ _2) * effHLT * (_1 - effHLT);
  253. RooFormulaVar zMuMu1HLTEffTerm("zMuMu1HLTEffTerm", "Yield * (2.* (eff_tk)^2 * (eff_sa)^2 * (eff_iso)^2 * eff_hlt *(1.- eff_hlt))",
  254. RooArgList(eff_tk, eff_sa,eff_iso, eff_hlt, Yield)); //zMuMuEff2HLTTerm = (effTk ^ _2) * (effSa ^ _2) * (effIso ^ _2) * (effHLT ^ _2) ;
  255. RooFormulaVar zMuMu2HLTEffTerm("zMuMu2HLTEffTerm", "Yield * ((eff_tk)^2 * (eff_sa)^2 * (eff_iso)^2 * (eff_hlt)^2)",
  256. RooArgList(eff_tk, eff_sa,eff_iso, eff_hlt, Yield));
  257. //zMuMuNoIsoEffTerm = (effTk ^ _2) * (effSa ^ _2) * (_1 - (effIso ^ _2)) * (_1 - ((_1 - effHLT)^_2));
  258. RooFormulaVar zMuMuNoIsoEffTerm("zMuMuNoIsoEffTerm", "Yield * ((eff_tk)^2 * (eff_sa)^2 * (1.- (eff_iso)^2) * (1.- ((1.-eff_hlt)^2)))",
  259. RooArgList(eff_tk, eff_sa,eff_iso, eff_hlt, Yield));
  260. //zMuTkEffTerm = _2 * (effTk ^ _2) * effSa * (_1 - effSa) * (effIso ^ _2) * effHLT;
  261. RooFormulaVar zMuTkEffTerm("zMuTkEffTerm", "Yield * (2. *(eff_tk)^2 * eff_sa * (1.-eff_sa)* (eff_iso)^2 * eff_hlt)",
  262. RooArgList(eff_tk, eff_sa,eff_iso, eff_hlt, Yield));
  263. //zMuSaEffTerm = _2 * (effSa ^ _2) * effTk * (_1 - effTk) * (effIso ^ _2) * effHLT;
  264. RooFormulaVar zMuSaEffTerm("zMuSaEffTerm", "Yield * (2. *(eff_sa)^2 * eff_tk * (1.-eff_tk)* (eff_iso)^2 * eff_hlt)",
  265. RooArgList(eff_tk, eff_sa,eff_iso, eff_hlt, Yield));
  266. // creating model for the fit
  267. // z mu track
  268. RooGenericPdf *bkg_mutrk = new RooGenericPdf("bkg_mutrk","zmt bkg_model", "exp(alpha*mass) * ( a0 + a1 * mass + a2 * mass^2 )", RooArgSet( *mass, alpha, a0, a1,a2) );
  269. // RooFormulaVar fracSigMutrk("fracSigMutrk", "@0 / (@0 + @1)", RooArgList(zMuTkEffTerm, nbkg_mutrk ));
  270. RooAddPdf * model_mutrk= new RooAddPdf("model_mutrk","model_mutrk",RooArgList(*ZmtPdf,*bkg_mutrk),RooArgList(zMuTkEffTerm , nbkg_mutrk)) ;
  271. // z mu mu not Iso
  272. // creating background pdf for zmu mu not Iso
  273. RooGenericPdf *bkg_mumuNotIso = new RooGenericPdf("bkg_mumuNotIso","zmumuNotIso bkg_model", "exp(beta * mass) * (b0 + b1 * mass + b2 * mass^2)", RooArgSet( *mass, beta, b0, b1,b2) );
  274. // RooFormulaVar fracSigMuMuNoIso("fracSigMuMuNoIso", "@0 / (@0 + @1)", RooArgList(zMuMuNoIsoEffTerm, nbkg_mumuNotIso ));
  275. RooAddPdf * model_mumuNotIso= new RooAddPdf("model_mumuNotIso","model_mumuNotIso",RooArgList(*ZmmNoIsoPdf,*bkg_mumuNotIso), RooArgList(zMuMuNoIsoEffTerm, nbkg_mumuNotIso));
  276. // z mu sta
  277. // RooGenericPdf model_musta("model_musta", " ZmsPdf * zMuSaEffTerm ", RooArgSet( *ZmsPdf, zMuSaEffTerm)) ;
  278. RooGenericPdf *bkg_musa = new RooGenericPdf("bkg_musa","zmusa bkg_model", "exp(gamma * mass) * (c0 + c1 * mass )", RooArgSet( *mass, gamma, c0, c1) );
  279. // RooAddPdf * eZmsSig= new RooAddPdf("eZmsSig","eZmsSig",RooArgList(*,*bkg_mumuNotIso), RooArgList(zMuMuNoIsoEffTerm, nbkg_mumuNotIso));
  280. RooAddPdf *eZmsSig = new RooAddPdf("eZmsSig","eZmsSig",RooArgList(*ZmsPdf,*bkg_musa),RooArgList(zMuSaEffTerm , nbkg_musa)) ;
  281. //RooExtendPdf * eZmsSig= new RooExtendPdf("eZmsSig","extended signal p.d.f for zms ",*ZmsPdf, zMuSaEffTerm ) ;
  282. // z mu mu HLT
  283. // count ZMuMu Yield
  284. double nZMuMu = 0.;
  285. double nZMuMu1HLT = 0.;
  286. double nZMuMu2HLT = 0.;
  287. unsigned int nBins = zmm->GetNbinsX();
  288. double xMin = zmm->GetXaxis()->GetXmin();
  289. double xMax = zmm->GetXaxis()->GetXmax();
  290. double deltaX =(xMax - xMin) / nBins;
  291. for(unsigned int i = 0; i < nBins; ++i) {
  292. double x = xMin + (i +.5) * deltaX;
  293. if(x > fMin && x < fMax){
  294. nZMuMu += zmm->GetBinContent(i+1);
  295. nZMuMu1HLT += zmm1hlt->GetBinContent(i+1);
  296. nZMuMu2HLT += zmm2hlt->GetBinContent(i+1);
  297. }
  298. }
  299. cout << ">>> count of ZMuMu yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu << endl;
  300. cout << ">>> count of ZMuMu (1HLT) yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu1HLT << endl;
  301. cout << ">>> count of ZMuMu (2HLT) yield in the range [" << fMin << ", " << fMax << "]: " << nZMuMu2HLT << endl;
  302. // we set eff_hlt
  303. //eff_hlt.setVal( 1. / (1. + (nZMuMu1HLT/ (2 * nZMuMu2HLT))) ) ;
  304. // creating the pdf for z mu mu 1hlt
  305. RooExtendPdf * eZmm1hltSig= new RooExtendPdf("eZmm1hltSig","extended signal p.d.f for zmm 1hlt",*ZmmHltPdf, zMuMu1HLTEffTerm ) ;
  306. // creating the pdf for z mu mu 2hlt
  307. RooExtendPdf * eZmm2hltSig= new RooExtendPdf("eZmm2hltSig","extended signal p.d.f for zmm 2hlt",*ZmmHltPdf, zMuMu2HLTEffTerm ) ;
  308. // getting the data if fit otherwise constructed the data for model if toy....
  309. RooDataHist *zmtMass, * zmmNotIsoMass, *zmsMass, *zmm1hltMass, *zmm2hltMass;
  310. if (toy){
  311. zmtMass = genHistFromModelPdf("zmtMass", model_mutrk, mass, scaleLumi, range,rebinZMuTk, seed);
  312. zmmNotIsoMass = genHistFromModelPdf("zmmNotIsoMass", model_mumuNotIso, mass, scaleLumi, range,rebinZMuMuNoIso, seed);
  313. zmsMass = genHistFromModelPdf("zmsMass", eZmsSig, mass, scaleLumi, range,rebinZMuSa , seed);
  314. zmm1hltMass = genHistFromModelPdf("zmm1hltMass", eZmm1hltSig, mass, scaleLumi, range,rebinZMuMuHlt, seed );
  315. zmm2hltMass = genHistFromModelPdf("zmm2hltMass", eZmm2hltSig, mass, scaleLumi, range,rebinZMuMuHlt, seed);
  316. } else {// if fit from data....
  317. zmtMass = new RooDataHist("zmtMass", "good z mu track" , RooArgList(*mass), zmt );
  318. zmmNotIsoMass = new RooDataHist("ZmmNotIso", "good z mu mu not isolated" , RooArgList(*mass), zmmNotIso );
  319. zmsMass = new RooDataHist("zmsMass", "good z mu sta mass" , RooArgList(*mass), zmmsta );
  320. zmm1hltMass = new RooDataHist("zmm1hltMass", "good Z mu mu 1hlt" , RooArgList(*mass), zmm1hlt );
  321. zmm2hltMass = new RooDataHist("zmm2hltMass", "good Z mu mu 2hlt" , RooArgList(*mass), zmm2hlt );
  322. }
  323. // creting the chi2s
  324. RooChi2Var *chi2_mutrk = new RooChi2Var("chi2_mutrk","chi2_mutrk",*model_mutrk,*zmtMass, Extended(kTRUE),DataError(RooAbsData::SumW2) ) ;
  325. RooChi2Var *chi2_mumuNotIso = new RooChi2Var("chi2_mumuNotIso","chi2_mumuNotIso",*model_mumuNotIso,*zmmNotIsoMass,Extended(kTRUE), DataError(RooAbsData::SumW2)) ;
  326. RooChi2Var *chi2_musta = new RooChi2Var("chi2_musta","chi2_musta",*eZmsSig, *zmsMass, Extended(kTRUE) ,DataError(RooAbsData::SumW2) ) ;
  327. // uncomment this line if you want to use logLik for mu sta
  328. // RooNLLVar *chi2_musta = new RooNLLVar("chi2_musta","chi2_musta",*eZmsSig, *zmsMass, Extended(kTRUE), DataError(RooAbsData::SumW2) ) ;
  329. RooChi2Var *chi2_mu1hlt = new RooChi2Var("chi2_mu1hlt","chi2_mu1hlt", *eZmm1hltSig, *zmm1hltMass, Extended(kTRUE) , DataError(RooAbsData::SumW2)) ;
  330. RooChi2Var *chi2_mu2hlt = new RooChi2Var("chi2_mu2hlt","chi2_mu2hlt", *eZmm2hltSig, *zmm2hltMass, Extended(kTRUE), DataError(RooAbsData::SumW2) ) ;
  331. // adding the chi2
  332. RooAddition totChi2("totChi2","chi2_mutrk + chi2_mumuNotIso + chi2_musta + chi2_mu1hlt + chi2_mu2hlt " , RooArgSet(*chi2_mutrk ,*chi2_mumuNotIso , *chi2_musta , *chi2_mu1hlt , *chi2_mu2hlt )) ;
  333. // printing out the model integral befor fitting
  334. double N_zMuMu1hlt = eZmm1hltSig->expectedEvents(*mass);
  335. double N_bkgTk = bkg_mutrk->expectedEvents(*mass);
  336. double N_bkgIso = bkg_mumuNotIso->expectedEvents(*mass);
  337. double N_zMuTk = zMuTkEffTerm.getVal();
  338. double e_hlt = eff_hlt.getVal();
  339. double e_tk = eff_tk.getVal();
  340. double e_sa = eff_sa.getVal();
  341. double e_iso = eff_iso.getVal();
  342. double Y_hlt = N_zMuMu1hlt / ((2.* (e_tk * e_tk) * (e_sa * e_sa) * (e_iso* e_iso) * e_hlt *(1.- e_hlt))) ;
  343. double Y_mutk = N_zMuTk / ((2.* (e_tk * e_tk) * e_sa *(1.- e_sa) * (e_iso* e_iso) * e_hlt )) ;
  344. cout << "Yield prediction from mumu1hlt integral befor fitting: " << Y_hlt << endl
  345. ;
  346. cout << "Yield prediction from mutk integral befor fitting: " << Y_mutk << endl;
  347. cout << "Bkg for mutk prediction from mutk integral after fitting: " << N_bkgTk << endl;
  348. cout << "Bkg for mumuNotIso prediction from mumuNotIso integral after fitting: " << N_bkgIso << endl;
  349. fit( totChi2, numberOfBins, outfile.c_str());
  350. N_zMuMu1hlt = eZmm1hltSig->expectedEvents(*mass);
  351. N_zMuTk = zMuTkEffTerm.getVal();
  352. double N_zMuMuNoIso = zMuMuNoIsoEffTerm.getVal();
  353. double N_Tk = model_mutrk->expectedEvents(*mass);
  354. double N_Iso = model_mumuNotIso->expectedEvents(*mass);
  355. e_hlt = eff_hlt.getVal();
  356. e_tk = eff_tk.getVal();
  357. e_sa = eff_sa.getVal();
  358. e_iso = eff_iso.getVal();
  359. Y_hlt = N_zMuMu1hlt / ((2.* (e_tk * e_tk) * (e_sa * e_sa) * (e_iso* e_iso) * e_hlt *(1.- e_hlt))) ;
  360. Y_mutk = N_zMuTk / ((2.* (e_tk * e_tk) * e_sa *(1.- e_sa) * (e_iso* e_iso) * e_hlt )) ;
  361. cout << "Yield prediction from mumu1hlt integral after fitting: " << Y_hlt << endl;
  362. //cout << "Yield prediction from mutk integral after fitting: " << Y_mutk << endl;
  363. cout << "N + B prediction from mutk integral after fitting: " << N_Tk << endl;
  364. cout << "zMuTkEffTerm " << N_zMuTk << endl;
  365. cout << "N + B prediction from mumuNotIso integral after fitting: " << N_Iso << endl;
  366. cout << "zMuMuNoIsoEffTerm " << N_zMuMuNoIso << endl;
  367. cout << "chi2_mutrk:" << chi2_mutrk->getVal()<< endl;
  368. cout << "chi2_mumuNotIso:" << chi2_mumuNotIso->getVal() << endl;
  369. cout << "chi2_musta:" << chi2_musta->getVal() << endl;
  370. cout << "chi2_mumu1hlt:" << chi2_mu1hlt->getVal()<< endl;
  371. cout << "chi2_mumu2hlt:" << chi2_mu2hlt->getVal()<< endl;
  372. //plotting
  373. RooPlot * massFrame_mutrk = mass->frame() ;
  374. RooPlot * massFrame_mumuNotIso = mass->frame() ;
  375. RooPlot * massFrame_musta = mass->frame() ;
  376. RooPlot * massFrame_mumu1hlt = mass->frame() ;
  377. RooPlot * massFrame_mumu2hlt = mass->frame() ;
  378. TCanvas * canv1 = new TCanvas("canvas");
  379. TCanvas * canv2 = new TCanvas("new_canvas");
  380. canv1->Divide(2,3);
  381. canv1->cd(1);
  382. canv1->SetLogy(kTRUE);
  383. zmtMass->plotOn(massFrame_mutrk, LineColor(kBlue)) ;
  384. model_mutrk->plotOn(massFrame_mutrk,LineColor(kRed)) ;
  385. model_mutrk->plotOn(massFrame_mutrk,Components(*bkg_mutrk),LineColor(kGreen)) ;
  386. massFrame_mutrk->SetTitle("Z -> #mu track");
  387. // massFrame_mutrk->GetYaxis()->SetLogScale();
  388. massFrame_mutrk->Draw();
  389. gPad->SetLogy(1);
  390. canv2->cd();
  391. canv2->SetLogy(kTRUE);
  392. massFrame_mutrk->Draw();
  393. canv2->SaveAs("LogZMuTk.eps");
  394. canv2->SetLogy(kFALSE);
  395. canv2->SaveAs("LinZMuTk.eps");
  396. canv1->cd(2);
  397. zmmNotIsoMass->plotOn(massFrame_mumuNotIso, LineColor(kBlue)) ;
  398. model_mumuNotIso->plotOn(massFrame_mumuNotIso,LineColor(kRed)) ;
  399. model_mumuNotIso->plotOn(massFrame_mumuNotIso,Components(*bkg_mumuNotIso), LineColor(kGreen)) ;
  400. massFrame_mumuNotIso->SetTitle("Z -> #mu #mu not isolated");
  401. massFrame_mumuNotIso->Draw();
  402. gPad->SetLogy(1);
  403. canv2->cd();
  404. canv2->Clear();
  405. canv2->SetLogy(kTRUE);
  406. massFrame_mumuNotIso->Draw();
  407. canv2->SaveAs("LogZMuMuNotIso.eps");
  408. canv2->SetLogy(kFALSE);
  409. canv2->SaveAs("LinZMuMuNotIso.eps");
  410. canv1->cd(3);
  411. zmsMass->plotOn(massFrame_musta, LineColor(kBlue)) ;
  412. eZmsSig->plotOn(massFrame_musta,Components(*bkg_musa), LineColor(kGreen)) ;
  413. eZmsSig->plotOn(massFrame_musta,LineColor(kRed)) ;
  414. massFrame_musta->SetTitle("Z -> #mu sta");
  415. massFrame_musta->Draw();
  416. canv2->cd();
  417. canv2->Clear();
  418. canv2->SetLogy(kTRUE);
  419. massFrame_musta->Draw();
  420. canv2->SaveAs("LogZMuSa.eps");
  421. canv2->SetLogy(kFALSE);
  422. canv2->SaveAs("LinZMuSa.eps");
  423. canv1->cd(4);
  424. zmm1hltMass->plotOn(massFrame_mumu1hlt, LineColor(kBlue)) ;
  425. eZmm1hltSig->plotOn(massFrame_mumu1hlt,LineColor(kRed)) ;
  426. massFrame_mumu1hlt->SetTitle("Z -> #mu #mu 1hlt");
  427. massFrame_mumu1hlt->Draw();
  428. canv2->cd();
  429. canv2->Clear();
  430. canv2->SetLogy(kTRUE);
  431. massFrame_mumu1hlt->Draw();
  432. canv2->SaveAs("LogZMuMu1Hlt.eps");
  433. canv2->SetLogy(kFALSE);
  434. canv2->SaveAs("LinZMuMu1Hlt.eps");
  435. canv1->cd(5);
  436. zmm2hltMass->plotOn(massFrame_mumu2hlt, LineColor(kBlue)) ;
  437. eZmm2hltSig->plotOn(massFrame_mumu2hlt,LineColor(kRed)) ;
  438. massFrame_mumu2hlt->SetTitle("Z -> #mu #mu 2hlt");
  439. massFrame_mumu2hlt->Draw();
  440. canv1->SaveAs("mass.eps");
  441. canv2->cd();
  442. canv2->Clear();
  443. canv2->SetLogy(kTRUE);
  444. massFrame_mumu2hlt->Draw();
  445. canv2->SaveAs("LogZMuMu2Hlt.eps");
  446. canv2->SetLogy(kFALSE);
  447. canv2->SaveAs("LinZMuMu2Hlt.eps");
  448. /* how to read the fit result in root
  449. TH1D h_Yield("h_Yield", "h_Yield", 100, 10000, 30000)
  450. for (int i =0: i < 100; i++){
  451. RooFitResult* r = gDirectory->Get(Form("toy_totChi2;%d)",i)
  452. //r->floatParsFinal().Print("s");
  453. // without s return a list, can we get the number?
  454. RooFitResult* r = gDirectory->Get("toy_totChi2;1")
  455. // chi2
  456. r->minNll();
  457. //distamce form chi2.....
  458. //r->edm();
  459. // yield
  460. r->floatParsFinal()[0]->Print();
  461. //RooAbsReal * l = r->floatParsFinal()->first()
  462. RooAbsReal * y = r->floatParsFinal()->find("Yield");
  463. h_Yield->Fill(y->getVal());
  464. }
  465. */
  466. delete root_file;
  467. //delete out_root_file;
  468. return 0;
  469. }