LightFlavorRatios/mass_histograms/generate_sWeights.C

0001 void generate_sWeights()
0002 {
0003   TFile* Ks_file = TFile::Open("Kshort_3runs.root");
0004   TFile* lambda_file = TFile::Open("Lambda_3runs.root");
0005
0006   TTree* Ks_tree = (TTree*)Ks_file->Get("DecayTree");
0007   TTree* lambda_tree = (TTree*)lambda_file->Get("DecayTree");
0008
0009   TFile* Ks_SWfile = new TFile("Kshort_SW_3runs.root","RECREATE");
0010   TTree* Ks_SWtree = Ks_tree->CloneTree(-1);
0011
0012   TFile* lambda_SWfile = new TFile("lambda_SW_3runs.root","RECREATE");
0013   TTree* lambda_SWtree = lambda_tree->CopyTree("Lambda0_mass>=1.1&&Lambda0_mass<1.13");
0014
0015   RooRealVar m_ks("K_S0_mass","K_S0_mass",0.45,0.55);
0016   RooRealVar m_lambda("Lambda0_mass","Lambda0_mass",1.1,1.13);
0017
0018   RooDataSet Ks_ds("K_S0","K_S0",RooArgSet(m_ks),RooFit::Import(*Ks_tree));
0019   RooDataSet lambda_ds("Lambda0","Lambda0",RooArgSet(m_lambda),RooFit::Import(*lambda_tree));
0020
0021   // Kshort signal and background models
0022   RooRealVar ks_mean("ks_mean","mean",0.49,0.48,0.5);
0023   RooRealVar ks_width("ks_width","width",0.01,0.0001,0.2);
0024   RooGaussian ks_signal("ks_signal","signal",m_ks,ks_mean,ks_width);
0025
0026   RooRealVar ks_lin_coef("ks_lin_coef","lin_coef",0.,0.,1.);
0027   RooRealVar ks_quad_coef("ks_quad_coef","quad_coef",0.,0.,1.);
0028   RooRealVar ks_cubic_coef("ks_cubic_coef","cubic_coef",0.,0.,1.);
0029   RooRealVar ks_quartic_coef("ks_quartic_coef","quartic_coef",0.,0.,1.);
0030   RooChebychev ks_bkg("ks_bkg","background",m_ks,RooArgList(ks_lin_coef,ks_quad_coef,ks_cubic_coef,ks_quartic_coef));
0031
0032   RooRealVar ks_nsignal("ks_nsignal","nSignal",0.,2*Ks_ds.numEntries());
0033   RooRealVar ks_nbkg("ks_nbkg","nBackground",0.,2*Ks_ds.numEntries());
0034   RooAddPdf ks_model("ks_model","combined signal and background",RooArgList(ks_signal,ks_bkg),RooArgList(ks_nsignal,ks_nbkg));
0035
0036   // Lambda signal and background models
0037   RooRealVar lambda_mean("lambda_mean","mean",1.113,1.11,1.12);
0038   RooRealVar lambda_width("lambda_width","width",0.01,0.0001,0.2);
0039   RooGaussian lambda_signal("lambda_signal","signal",m_lambda,lambda_mean,lambda_width);
0040
0041   RooRealVar lambda_lin_coef("lambda_lin_coef","lin_coef",0.,0.,1.);
0042   RooRealVar lambda_quad_coef("lambda_quad_coef","quad_coef",0.,0.,1.);
0043   RooRealVar lambda_cubic_coef("lambda_cubic_coef","cubic_coef",0.,0.,1.);
0044   RooRealVar lambda_quartic_coef("lambda_quartic_coef","quartic_coef",0.,0.,1.);
0045   RooChebychev lambda_bkg("lambda_bkg","background",m_lambda,RooArgList(lambda_lin_coef,lambda_quad_coef,lambda_cubic_coef));
0046
0047   RooRealVar lambda_nsignal("lambda_nsignal","nSignal",0.,2*Ks_ds.numEntries());
0048   RooRealVar lambda_nbkg("lambda_nbkg","nBackground",0.,2*Ks_ds.numEntries());
0049   RooAddPdf lambda_model("lambda_model","combined signal and background",RooArgList(lambda_signal,lambda_bkg),RooArgList(lambda_nsignal,lambda_nbkg));
0050
0051   ks_model.fitTo(Ks_ds);
0052   lambda_model.fitTo(lambda_ds);
0053
0054   RooStats::SPlot* Ks_sData = new RooStats::SPlot("Ks_sData","Ks_sPlot",Ks_ds,&ks_model,RooArgList(ks_nsignal,ks_nbkg));
0055   RooStats::SPlot* lambda_sData = new RooStats::SPlot("lambda_sData","lambda_sPlot",lambda_ds,&lambda_model,RooArgList(lambda_nsignal,lambda_nbkg));
0056
0057   double ks_SW;
0058   double lambda_SW;
0059
0060   TBranch* Ks_SWbranch = Ks_SWtree->Branch("sWeight",&ks_SW,"sWeight/D");
0061   TBranch* lambda_SWbranch = lambda_SWtree->Branch("sWeight",&lambda_SW,"sWeight/D");
0062
0063   std::cout << "Check SWeights:" << std::endl;
0064
0065   std::cout << "Yield of signal is "  << ks_nsignal.getVal() << ", " << lambda_nsignal.getVal()
0066               << ".  From sWeights it is " << Ks_sData->GetYieldFromSWeight("ks_nsignal") << ", " << lambda_sData->GetYieldFromSWeight("lambda_nsignal") << std::endl;
0067
0068   for(int i=0; i<Ks_ds.numEntries(); i++)
0069   {
0070     const RooArgSet* row = Ks_ds.get(i);
0071     ks_SW = (double)row->getRealValue("ks_nsignal_sw");
0072     Ks_SWbranch->Fill();
0073   }
0074
0075   std::cout << "Kshort done" << std::endl;
0076
0077   for(int i=0; i<lambda_ds.numEntries(); i++)
0078   {
0079     const RooArgSet* row = lambda_ds.get(i);
0080     lambda_SW = (double)row->getRealValue("lambda_nsignal_sw");
0081     lambda_SWbranch->Fill();
0082   }
0083
0084   Ks_SWfile->Write();
0085   Ks_SWfile->Close();
0086
0087   lambda_SWfile->Write();
0088   lambda_SWfile->Close();
0089 }