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