void convolution_hist () { // macro to create a MC dataset from an exponential function convoluted // with a smeared distribution for the tau constant // steps // 1) load the histogram containg the possible tau values // 2) build a RooHistPdf from a RooHistDataSet based on the histogram // 3) create an exponential pdf // 4) convolve the exp pdf and the smearing function (use RooFFTConvPdf) // 5) create a MC sample based on the convolved Pdf // gSystem->Load("libRooFit.so"); using namespace RooFit; TCanvas *c1 = new TCanvas("c1","",1024,768); c1->cd(); // 1) load histogram TFile* f= new TFile("tests_140224-fei4-250um-fl3.8e15-extract_traps-V400-taue-in-ns.root","READ"); f->cd(); TH1D* hTauE = (TH1D*)f->Get("hTauE"); if ( !hTauE) { std::cout << "hTauE not found.\Exiting...\n"; exit(1); } // determining tau range int nbins = hTauE->GetNbinsX(); double xmin = hTauE->GetBinLowEdge(1); double xmax = hTauE->GetBinLowEdge(nbins+1); double mean = hTauE->GetMean(); std::cout << "mean = " << mean << "\n"; std::cout << "xmin = " << xmin << "\n"; std::cout << "xmax = " << xmax << "\n"; // the tau parameter RooRealVar tau("tau","#tau [ns]",mean,xmin,xmax); // the data set RooDataHist* tauEDH = new RooDataHist("tauEDH","tauEDH",tau,Import(*hTauE)); // 2) build the RooHistPdf RooHistPdf tauHPDF("tauHPDF","tauHPDF",tau,*tauEDH,0); // check the result RooPlot* frame1 = tau.frame(Title("Tau histogram pdf")); tauEDH->plotOn(frame1); tauHPDF.plotOn(frame1); frame1->Draw(); // 3) The exponential pdf RooRealVar t("t","t [ns]",0,fabs(xmax)*5); t.Print(); tau.Print(); RooExponential expPDF("expPDF","expPDF",t,tau); RooPlot* frame2 = t.frame(Title("Trapping time")); expPDF.plotOn(frame2); TCanvas *c2 = new TCanvas("c2","",1024,768); c2->cd(); frame2->Draw(); // 4) convolve the exp pdf and the smearing function (use RooFFTConvPdf) // inspired by rf211_paramconv.C // Convolution in tau parameter trapConv = expPDF(t,tau) (x) tauHPDF(tau) t.setBins(1000,"cache"); tau.setBins(50,"cache"); RooFFTConvPdf trapConv("trapConv","trapConv",tau,expPDF,tauHPDF); //trapConv.Print(); // Configure convolution to construct a 2-D cache in (t,tau) // rather than a 1-d cache in tau that needs to be recalculated // for each value of t trapConv.setCacheObservables(t) ; trapConv.setBufferFraction(1.0) ; // Integrate trapConv over tau; projModel = Int trapConv dtau RooAbsPdf* projtrapConv = trapConv.createProjection(tau) ; // Generate 1000 toy events RooDataHist* d = projtrapConv->generateBinned(t,1000) ; // Fit p.d.f. to toy data projtrapConv->fitTo(*d,Verbose()) ; // Plot data and fitted p.d.f. RooPlot* frame3 = t.frame(Bins(25)) ; d->plotOn(frame3) ; projtrapConv->plotOn(frame3) ; TCanvas *c3 = new TCanvas("c3","",1024,768); c3->cd(); frame3->Draw(); // Make 2d histogram of trapConv(x;mean) TH1* hh = trapConv.createHistogram("hh",t,Binning(50),YVar(tau,Binning(50)),ConditionalObservables(tau)) ; hh->SetTitle("histogram of trapConv(t|tau)") ; hh->SetLineColor(kBlue) ; // Draw frame on canvas TCanvas* c = new TCanvas("rf211_paramconv","rf211_paramconv",800,400) ; c->Divide(2) ; c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3->Draw() ; c->cd(2) ; gPad->SetLeftMargin(0.20) ; hh->GetZaxis()->SetTitleOffset(2.5) ; hh->Draw("surf") ; // closing the file // f->Close(); }