PageRenderTime 109ms CodeModel.GetById 17ms app.highlight 85ms RepoModel.GetById 1ms app.codeStats 0ms

/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>
 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()