///////////////////////////////////////////////////////////////////////// // // 'VALIDATION AND MC STUDIES' RooFit tutorial macro #801 // // A Toy Monte Carlo study that perform cycles of // event generation and fittting // // ///////////////////////////////////////////////////////////////////////// #ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooConstVar.h" #include "RooChebychev.h" #include "RooAddPdf.h" #include "RooMCStudy.h" #include "RooPlot.h" #include "TCanvas.h" #include "TAxis.h" #include "TH2.h" #include "RooFitResult.h" #include "TStyle.h" #include "TDirectory.h" #include #include "RooAddition.h" #include "TH1F.h" using namespace RooFit ; using namespace std; void GetEta4Histos(vector &myvect); void rf801_mcstudy_MyModel_eta4_PythiaB() { // C r e a t e m o d e l // ----------------------- // Declare observable x RooRealVar x("ptrel","p_{T}^{rel} [GeV]",0.,4.) ; x.setBins(4); vector hhs; GetEta4Histos(hhs); // Create PDFs and their paramaters //all filtered MB 0 1 2 3, with PythiaB 0 1 5 6 RooDataHist * dhH = new RooDataHist ("dhH","dhH ",x,Import( *hhs.at(0) )); RooDataHist * dhG = new RooDataHist ("dhC","dhG ",x,Import( *hhs.at(1) )); RooDataHist * dhC = new RooDataHist ("dhC","dhC ",x,Import( *hhs.at(5) )); RooDataHist * dhB = new RooDataHist ("dhB","dhB ",x,Import( *hhs.at(6) )); RooHistPdf histpdfH ("histpdfH" , "histpdfH" ,x,*dhH,0) ; RooHistPdf histpdfG ("histpdfG" , "histpdfG" ,x,*dhG,0) ; RooHistPdf histpdfC ("histpdfC" , "histpdfC" ,x,*dhC,0) ; RooHistPdf histpdfB ("histpdfB" , "histpdfB" ,x,*dhB,0) ; // Sum the composite signal and background (from "data" fit log) double nmax = 2467 ; //for "eta4" eta bin double nResB = 1257.73 ; // with chad from MB frac: 1257.73; with h=0, only C in chad: 1388.1; double nResChad = 586.036 ; // 586.036; 315.306; double nResG = 623.493 ; // 623.493; 763.785; double hfracMC = 0.223368; // 0.223368; 0.0; RooRealVar hadfrac("hadfrac","fraction of hadrons in C+Hadrons",hfracMC); hadfrac.setConstant(kTRUE);/// later change this? for systematics RooAddPdf chad ("chad","PDF for C electrons and Hadrons",RooArgList(histpdfH,histpdfC),hadfrac); RooRealVar nbsig("nbsig","number of B #rightarrow e " , nResB , 0.0 , nmax); // the max should be the total integral of data RooRealVar nchad("nchad","number of C #rightarrow e + hadron ", nResChad, 0.0 , nmax); RooRealVar nconv("nconv","number of #gamma conversions" , nResG , 0.0 , nmax); RooAddition ntot("ntot","total sum of components", RooArgSet(nbsig,nchad,nconv)); // RooAddPdf bgfull ("bgfull","BG electrons full pdf conversions, hadrons and C",RooArgList(histpdfG,chad),nconv,nchad); RooAddPdf model("fullModel","full model B + CHad + Gamma",RooArgList(histpdfB,chad,histpdfG),RooArgList(nbsig,nchad,nconv)); // C r e a t e m a n a g e r // --------------------------- // Instantiate RooMCStudy manager on model with x as observable and given choice of fit options // // The Silence() option kills all messages below the PROGRESS level, leaving only a single message // per sample executed, and any error message that occur during fitting // // The Extended() option has two effects: // 1) The extended ML term is included in the likelihood and // 2) A poisson fluctuation is introduced on the number of generated events // // The FitOptions() given here are passed to the fitting stage of each toy experiment. // If Save() is specified, the fit result of each experiment is saved by the manager // // A Binned() option is added in this example to bin the data between generation and fitting // to speed up the study at the expemse of some precision RooMCStudy* mcstudy = new RooMCStudy(model,x,Binned(kTRUE),Silence(),Extended(), FitOptions(Save(kTRUE),PrintEvalErrors(0),Minos(kTRUE))) ; // G e n e r a t e a n d f i t e v e n t s // --------------------------------------------- int nPE=1500; // Generate and fit 1000 samples of Poisson(nExpected) events mcstudy->generateAndFit(nPE/*,nmax,kTRUE*/) ;// hopefully saved //no poisson?? // E x p l o r e r e s u l t s o f s t u d y // ------------------------------------------------ // Make plots of the distributions of mean, the error on mean and the pull of mean RooPlot* frame1 = mcstudy->plotParam(nbsig,Bins(30)) ; RooPlot* frame2 = mcstudy->plotError(nbsig,Bins(30)) ; RooPlot* frame3 = mcstudy->plotPull(nbsig,Bins(30),FrameRange(-8,8),FitGauss(kTRUE)); //nbsig,-8,8,40); ; // Plot distribution of minimized likelihood RooPlot* frame4 = mcstudy->plotNLL(Bins(30)) ; // // Make some histograms from the parameter dataset // TH1* hh_cor_a0_s1f = mcstudy->fitParDataSet().createHistogram("hh",a1,YVar(sig1frac)) ; // TH1* hh_cor_a0_a1 = mcstudy->fitParDataSet().createHistogram("hh",a0,YVar(a1)) ; // // Access some of the saved fit results from individual toys // TH2* corrHist000 = mcstudy->fitResult(0)->correlationHist("c000") ; // TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ; // TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ; // Draw all plots on a canvas gStyle->SetPalette(1) ; gStyle->SetOptStat(0) ; TCanvas* c = new TCanvas("nbsig","rf801_mcstudy",600,600) ; c->Divide(2,2) ; c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ; c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ; c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3->Draw() ; c->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.4) ; frame4->Draw() ; // c->cd(5) ; gPad->SetLeftMargin(0.15) ; hh_cor_a0_s1f->GetYaxis()->SetTitleOffset(1.4) ; hh_cor_a0_s1f->Draw("box") ; // c->cd(6) ; gPad->SetLeftMargin(0.15) ; hh_cor_a0_a1->GetYaxis()->SetTitleOffset(1.4) ; hh_cor_a0_a1->Draw("box") ; // c->cd(7) ; gPad->SetLeftMargin(0.15) ; corrHist000->GetYaxis()->SetTitleOffset(1.4) ; corrHist000->Draw("colz") ; // c->cd(8) ; gPad->SetLeftMargin(0.15) ; corrHist127->GetYaxis()->SetTitleOffset(1.4) ; corrHist127->Draw("colz") ; // c->cd(9) ; gPad->SetLeftMargin(0.15) ; corrHist953->GetYaxis()->SetTitleOffset(1.4) ; corrHist953->Draw("colz") ; RooPlot* frame1ch = mcstudy->plotParam(nchad,Bins(30)) ; RooPlot* frame2ch = mcstudy->plotError(nchad,Bins(30)) ; RooPlot* frame3ch = mcstudy->plotPull(nchad,FrameBins(35),FrameRange(-8,8),FitGauss(kTRUE)) ; // Plot distribution of minimized likelihood //RooPlot* frame4ch = mcstudy->plotNLL(Bins(50)) ; // // Make some histograms from the parameter dataset //TH1* hh_cor_nbsig_nchad = mcstudy->fitParDataSet().createHistogram("hh1",nbsig,YVar(nchad)) ; //TH1* hh_cor_nbsig_nconv = mcstudy->fitParDataSet().createHistogram("hh2",a0,YVar(a1)) ; // // Access some of the saved fit results from individual toys TH2* corrHist000 = mcstudy->fitResult(0)->correlationHist("c000") ; // TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ; // TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ; // Draw all plots on a canvas TCanvas* cchad = new TCanvas("cchad","rf801_mcstudy",600,600) ; cchad->Divide(2,2) ; cchad->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1ch->Draw() ; cchad->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2ch->Draw() ; cchad->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3ch->Draw() ; cchad->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.4) ; corrHist000->Draw("text") ; /////////////////////// RooPlot* frame1cv = mcstudy->plotParam(nconv,Bins(30)) ; RooPlot* frame2cv = mcstudy->plotError(nconv,Bins(30)) ; RooPlot* frame3cv = mcstudy->plotPull(nconv,Bins(30),FrameRange(-8,8),FitGauss(kTRUE)) ; // RooPlot* frame4cv = mcstudy->plotParam(ntot,Bins(50)); // Draw all plots on a canvas // // Access some of the saved fit results from individual toys // TH2* corrHist5871 = mcstudy->fitResult(0)->correlationHist("c5871") ; // TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ; // TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ; TH1F* hh_ntot = new TH1F("hh_ntot","sum of components",nmax/5,0.,2*nmax) ; for (int j=0; jfitResult(j); RooArgSet * ral=(RooArgSet*)fr->floatParsFinal(); double ntot0=ral->getRealValue("nbsig")+ral->getRealValue("nchad")+ral->getRealValue("nconv"); // cout<<"ntot0= "<Fill(ntot0); } TCanvas* cconv = new TCanvas("cconv","mcstudy conversions",600,600) ; cconv->Divide(2,2) ; cconv->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1cv->Draw() ; cconv->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2cv->Draw() ; cconv->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3cv->Draw() ; // cconv->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.4) ; corrHist5871->Draw("text") ; cconv->cd(4); gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.4) ; hh_ntot->Draw("p") ; // Make RooMCStudy object available on command line after // macro finishes gDirectory->Add(mcstudy) ; gDirectory->Add(hh_ntot) ; } void GetEta4Histos(vector & myvect){ TFile infile("FMBtemplates_wPythiaB_4trackEtaBins_wTaus_new.root","READ"); TH1D * h_B = (TH1D*)infile.Get("histBeta4"); TH1D * h_C = (TH1D*)infile.Get("histCeta4"); TH1D * h_BP = (TH1D*)infile.Get("histBPeta4"); TH1D * h_CP = (TH1D*)infile.Get("histCPeta4"); TH1D * h_G = (TH1D*)infile.Get("histGeta4"); TH1D * h_H = (TH1D*)infile.Get("histHeta4"); //TH2D * h2D_T=(TH2D*)infile.Get("hall"); double bins[5]={ 0., 1., 2., 3., 4. }; // // string hnameB="hB_ET10_12"; // TH1F* histB =(TH1F*) h2D_B->ProjectionX(hnameB.c_str(),10,12,"e"); string hnameB="histBR__eta4"; TH1F* histBR=h_B->Rebin(4,hnameB.c_str(),bins); int overflow=histBR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histBR "<SetBinContent(4,histBR->GetBinContent(4)+overflow); histBR->SetBinContent(5,0.0); } // string hnameC="hC_ET10_12"; // TH1F* histC =(TH1F*) h2D_C->ProjectionX(hnameC.c_str(),10,12,"e"); string hnameC="histCR__eta4"; TH1F* histCR=h_C->Rebin(4,hnameC.c_str(),bins); overflow=histCR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histCR "<SetBinContent(4,histCR->GetBinContent(4)+overflow); histCR->SetBinContent(5,0.0); } // string hnameG="hG_ET10_12"; // TH1F* histG =(TH1F*) h2D_G->ProjectionX(hnameG.c_str(),10,12,"e"); string hnameG="histGR__eta4"; TH1F* histGR=h_G->Rebin(4,hnameG.c_str(),bins); overflow=histGR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histGR "<SetBinContent(4,histGR->GetBinContent(4)+overflow); histGR->SetBinContent(5,0.0); } // string hnameH="hH_ET10_12"; // TH1F* histH =(TH1F*) h2D_H->ProjectionX(hnameH.c_str(),10,12,"e"); string hnameH="histHR__eta4"; TH1F* histHR=h_H->Rebin(4,hnameH.c_str(),bins); overflow=histHR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histHR "<SetBinContent(4,histHR->GetBinContent(4)+overflow); histHR->SetBinContent(5,0.0); } //pythiaB string hnameCP="histCPR__eta4"; TH1F* histCPR=h_CP->Rebin(4,hnameCP.c_str(),bins); overflow=histCPR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histCPR "<SetBinContent(4,histCPR->GetBinContent(4)+overflow); histCPR->SetBinContent(5,0.0); } string hnameBP="histBPR__eta4"; TH1F* histBPR=h_BP->Rebin(4,hnameBP.c_str(),bins); overflow=histBPR->GetBinContent(5); if(overflow>0){ //overflow cout<<"Overflow in histBPR "<SetBinContent(4,histBPR->GetBinContent(4)+overflow); histBPR->SetBinContent(5,0.0); } // /// add them up TH1F* histTot=histBR->Clone("hTot"); histTot->Add(histCR); histTot->Add(histGR); histTot->Add(histHR); histTot->SetTitle("all electrons in filtered MB"); TH1F* histTotPB=histHR->Clone("hTotPB"); histTotPB->Add(histGR); //add C and B from PythiaB to h_TP double scaleC=histCR.Integral()/histCPR.Integral(); histTotPB->Add(histCPR,scaleC); double scaleB=histBR.Integral()/histBPR.Integral(); histTotPB->Add(histBPR,scaleB); histTotPB->SetTitle("all electrons in filtered MB with PythiaB B and C"); cout<<"histTot.Integral()= "<Integral()<