LightFlavorRatios/swimming_correction/fit_pythia_lifetime.C

0001 void fit_pythia_lifetime()
0002 {
0003   TFile* f = TFile::Open("output/merged_fout.root");
0004
0005   TH1F* truth_lifetime = (TH1F*)f->Get("lifetime_labframe");
0006   std::cout << "h entries " << truth_lifetime->GetEntries() << std::endl;
0007   RooRealVar tau("tau","tau",0.,0.5);
0008
0009   RooDataHist tlh("tlh","truth lifetime hist",tau,RooFit::Import(*truth_lifetime));
0010   float lifetime_guess = 0.26; // [ns]
0011   RooRealVar k("k","k",-1./lifetime_guess,-1000.,-0.01);
0012   RooExponential expmodel("expmodel","exponential fit",tau,k);
0013
0014   RooPlot* frame = tau.frame(RooFit::Title("Truth lifetime"));
0015   tlh.plotOn(frame);
0016
0017   expmodel.fitTo(tlh,RooFit::SumW2Error(true));
0018
0019   std::cout << "fit lifetime: " << -1./k.getVal() << " ns" << std::endl;
0020
0021   TCanvas* ctruth = new TCanvas("ctruth","ctruth",600,600);
0022   expmodel.plotOn(frame);
0023   frame->Draw();
0024
0025
0026   TH1F* truth_smeared = (TH1F*)f->Get("lifetime_smeared");
0027   RooDataHist lsh("lsh","smeared lifetime hist",tau,RooFit::Import(*truth_smeared));
0028
0029   RooPlot* sm_frame = tau.frame(RooFit::Title("Resolution-smeared lifetime"));
0030   lsh.plotOn(sm_frame);
0031
0032   RooRealVar gmean("gmean","gmean",0.02,0.,0.1);
0033   RooRealVar gsigma("gsigma","gsigma",0.02,0.,0.1);
0034   RooRealVar lt("lt","lifetime",0.26,0.,10.);
0035   RooRealVar ltd("ltd","dlifetime",0.26,0.2,0.3);
0036
0037   RooGExpModel gm("gm","gm",tau,gmean,gsigma,lt,RooFit::RooConst(1.),RooFit::RooConst(1.),RooFit::RooConst(1.),true);
0038   RooDecay dc("dc","dc",tau,ltd,gm,RooDecay::DecayType::DoubleSided);
0039
0040   TCanvas* cmod = new TCanvas("cmod","cmod",600,600);
0041   dc.fitTo(lsh,RooFit::SumW2Error(true));
0042   dc.plotOn(sm_frame);
0043   sm_frame->Draw();
0044
0045   TH1F* t_res = (TH1F*)f->Get("lifetime_smear_amount");
0046   RooRealVar dt("dt","dt",-0.05,0.05);
0047   RooDataHist lres("lres","lifetime resolution",dt,RooFit::Import(*t_res));
0048
0049   RooRealVar rmean("rmean","rmean",0.,-0.1,0.1);
0050   RooRealVar rsigma("rmean","rmean",0.,-0.1,0.1);
0051   RooRealVar rk("rk","rk",0.,-100.,100.);
0052   RooGExpModel gmr("gmr","gmr",dt,rmean,rsigma,rk,RooFit::RooConst(1.),RooFit::RooConst(1.),RooFit::RooConst(1.),true);
0053   gmr.fitTo(lres);
0054
0055   TCanvas* cres = new TCanvas("cres","cres",600,600);
0056   RooPlot* res_frame = dt.frame(RooFit::Title("lifetime resolution"));
0057   lres.plotOn(res_frame);
0058   gmr.plotOn(res_frame);
0059   res_frame->Draw();
0060
0061
0062 /*
0063   RooRealVar gmean("gmean","gmean",0.06,0.04,0.08);
0064   RooRealVar gsigma("gsigma","gsigma",0.08,0.04,0.12);
0065   RooRealVar lt("lt","lifetime",0.05,0.,1.);
0066   RooGExpModel gexpmodel("gexpmodel","Exponential-modified Gaussian model",tau,gmean,gsigma,lt,RooFit::RooConst(1.),RooFit::RooConst(1.),RooFit::RooConst(1.),true);
0067
0068   gexpmodel.fitTo(lsh,RooFit::SumW2Error(true));
0069   gexpmodel.plotOn(sm_frame);
0070 */
0071   //sm_frame->Draw();
0072 }