#ifndef __CINT__ #endif #include #include #include #include using namespace std; TTree *generate_data(Int_t nevts=10000) { // generate the data sample gROOT->Delete("data"); TTree *tree = new TTree("data","dummy data"); // variables Double_t mass; Int_t dtype; // create tree branches tree->Branch("mass",&mass,"mass/D"); tree->Branch("type",&dtype,"type/I"); // data extraction loop for (Int_t i=0;i!=10000;++i) { Double_t type = gRandom->Uniform(); if (type<.6) { // generate signal A dtype = 1; mass = gRandom->Gaus(.42,.1); } else if (type<.7) { // generate signal B dtype = 2; mass = gRandom->Gaus(.7,.08); } else { // generate background events dtype = 0; Double_t shape = -.1; mass = (-1+TMath::Sqrt(1+gRandom->Uniform()*2*(1+shape/2)*shape))/shape; } tree->Fill(); } // end extractino loop return tree; } void rooTest(Bool_t binned=kTRUE) { // define the fit variables RooRealVar Mass("mass","Inv. mass [GeV/c^{2}]",0,1); // define fit parameters RooRealVar Peak1("Peak1","Peak 1",.4,.3,.5); RooRealVar Width1("Weak1","Width 1",.1,.01,.3); RooRealVar Peak2("Peak2","Peak 2",.8,.65,.82); RooRealVar Width2("Weak2","Width 2",.1,.01,.3); RooRealVar bkg_shape("bkg_shape","shape",0,-1,0); RooRealVar frac_bkg("frac_bkg","Bkg frac",.1,0,.5); RooRealVar frac_sgn("frac_bkg","Bkg frac",.3,0,.5); // define the fit model RooGaussian sgn1("sgn1","gaus(Mass,Peak1,Width1)", Mass,Peak1,Width1); RooGaussian sgn2("sgn1","gaus(Mass,Peak2,Width2)", Mass,Peak2,Width2); RooAddPdf sgn("sgn","Signals",sgn2,sgn1,frac_sgn); RooPolynomial bkg("bkg","Combinatoric Background", Mass,bkg_shape,1); RooAddPdf model_complete("model","Mass model",bkg,sgn,frac_bkg); RooAddPdf model_partial("model","Mass model",bkg,sgn1,frac_bkg); // take data TTree *data = generate_data(); // create the dataset gROOT->Delete("histo_mass"); TH1F *histo_mass = new TH1F("histo_mass","Mass",50,0,1); data->Draw("mass>>histo_mass",0,"goff"); RooDataHist data_histo("data_histo","Binned dataset",Mass,histo_mass); RooDataSet data_unbinned("data_unbinned","Dataset",data,Mass); // fit the data model_partial.printCompactTree(); Mass.setRange("low",0,.55); Mass.setRange("high",.8,1); RooFitResult *res_low; RooFitResult *res_high; RooFitResult *res_low_high; if (binned) { res_low = model_partial.fitTo(data_histo,RooFit::Range("low"),RooFit::Save()); res_high = model_partial.fitTo(data_histo,RooFit::Range("high"),RooFit::Save()); res_low_high = model_partial.fitTo(data_histo,RooFit::Range("low,high"),RooFit::Save()); } else { res_low = model_partial.fitTo(data_unbinned,RooFit::Range("low"),RooFit::Save()); res_high = model_partial.fitTo(data_unbinned,RooFit::Range("high"),RooFit::Save()); res_low_high = model_partial.fitTo(data_unbinned,RooFit::Range("low,high"),RooFit::Save()); } // print fit results cout << "*** low sideband ***" << endl; res_low->Print(); cout << "*** high sideband ***" << endl; res_high->Print(); cout << "*** low+high sidebands ***" << endl; res_low_high->Print(); // show the projections RooPlot *Mass_proj = Mass.frame(); TCanvas *canvas_bin = new TCanvas("c1","Mass projection (binned)",1); if (binned) data_histo.plotOn(Mass_proj); else data_unbinned.plotOn(Mass_proj); model_partial.plotOn(Mass_proj); model_partial.plotOn(Mass_proj,RooFit::Components(sgn1), RooFit::LineStyle(kDashed)); model_partial.plotOn(Mass_proj,RooFit::Components(bkg), RooFit::FillStyle(1001),RooFit::FillColor(kGray)); Mass_proj->Draw(); }