\(K_s\)粒子

在这一章我们一起来重建\(K_s\)，那在开始之前，你需要在PDG上仔细地研究一下这个粒子。它的质量大概是 \(0.497GeV\)，而寿命 \(\tau_{K_s}\approx8.954\times10^{−11}s\)比一般粒子要长，这是它区别于其他粒子的一个非常显著的特征。再看它的衰变模式，我们发现\(K_s\rightarrow \pi^+\pi^-\)的分支比大概是\(69.2\%\)，因此我们可以在 \(\pi^+\pi^-\)衰变道来寻找 \(K_s\)粒子。

1-Ntuple Maker

由于CMS探测器可以记录带电强子的track，我们是通过对track的寻找和匹配来重建 \(\pi^+\pi^-\)的。

1.1-track数据集

对于track的寻找，我们用的是"packedPFCandidates"数据集。可以用下面的方式查看这个数据集的tag。

$ edmDumpEventContent f8e93985-e14c-4a8a-b28b-f8cceb3c878e.root
...
vector<pat::PackedCandidate>          "packedPFCandidates"        ""                "PAT" 

参照朱峰学长的Multilep learning note中的教程，我们可以把.h文件和.cc文件进行以下修改。

//MuMuEEPat.h
private:
//define toke begin
edm::EDGetTokenT<edm::View<pat::Muon> > gtpatmuonToken_;//获取pat::Muon的token
edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;//获取MagneticField的token
edm::EDGetTokenT<edm::View<pat::PackedCandidate> > trackToken_;//获取pat::PackedCandidate的token
//define toke end

//MuMuEEPat.cc
MuMuEEPat::MuMuEEPat(const edm::ParameterSet& iConfig)
:magneticFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>())
{
  //extract toke
  trackToken_ = consumes<edm::View<pat::PackedCandidate> >(edm::InputTag("packedPFCandidates"));
  //extract toke
}

类似的，这里的trackToken_是一个成员变量，它保存了一个Token。这个Token表示模块需要pat::PackedCandidate类型的对象，数据来源是一个名为“packedPFCandidates”的数据集（由edm::InputTag指定）。这个数据集名称是在数据处理链中定义的，用于标识数据源。

1.2-for循环起始

下面就是最关键的对两个track的循环部分。

//MuMuEEPat.cc
void MuMuEEPat::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
    edm::Handle < edm::View < pat::Muon > > thePATMuonHandle;//这里thePATMuonHandle的作用在for循环中会解释
  iEvent.getByToken(gtpatmuonToken_, thePATMuonHandle );
  edm::Handle < edm::View < pat::PackedCandidate > > theTrackHandle;//定义一个智能指针 theTrackHandle
  iEvent.getByToken(trackToken_,theTrackHandle );//将之前的token提取出来并存储到 theTrackHandle

  if (theTrackHandle->size()>=2 ) {
    //loop for Track1
    //定义指向edm::View < pat::PackedCandidate >中pat::PackedCandidate对象常量迭代器，开始循环工作
    for (edm::View < pat::PackedCandidate >::const_iterator iTrack1 = theTrackHandle->begin(); iTrack1 != theTrackHandle->end(); ++iTrack1) {
      //start对muon循环，防止选中的track是muon的track
      bool matchmuon1 = false;
      if (thePATMuonHandle->size()>0){
        for (edm::View<pat::Muon>::const_iterator iMuonP = thePATMuonHandle->begin(); iMuonP != thePATMuonHandle->end(); ++iMuonP) {
          TrackRef muTrack = iMuonP->track();
          if (muTrack.isNull()) {continue;}
          if (reco::deltaR(iTrack1->momentum(), iMuonP->momentum()) < 0.01){
            matchmuon1 = true;
            break;
          }
        }
      }
      if (matchmuon1) {continue;}
      //end对muon循环，防止选中的track是muon的track

      //start检测选中的track是否有效，检测电荷是否是+/-1
      if (!iTrack1->hasTrackDetails() || iTrack1->charge() == 0){continue;}
      if (!(iTrack1->charge()==1 || iTrack1->charge()==-1)){continue;}
      if (!iTrack1->trackHighPurity()){continue;}
      const reco::Track* trk1 = iTrack1->bestTrack();
      if (!trk1) {continue;}
      //end检测选中的track是否有效，检测电荷是否是+/-1

      // loop for Track2
      for (edm::View < pat::PackedCandidate >::const_iterator iTrack2 = iTrack1 + 1; iTrack2 != theTrackHandle->end(); ++iTrack2) {
        //start对muon循环，防止选中的track是muon的track
        bool matchmuon2 = false;
        if (thePATMuonHandle->size()>0){
          for (edm::View<pat::Muon>::const_iterator iMuonP = thePATMuonHandle->begin(); iMuonP != thePATMuonHandle->end(); ++iMuonP) {
            TrackRef muTrack = iMuonP->track();
            if (muTrack.isNull()) {continue;}
            if (reco::deltaR(iTrack2->momentum(), iMuonP->momentum()) < 0.01){
              matchmuon2 = true;
              break;
            }
          }
        }
        if (matchmuon2) {continue;}
        //end对muon循环，防止选中的track是muon的track

        //start检测选中的track是否有效，检测电荷是否是+/-1
        if (!iTrack2->hasTrackDetails() || iTrack2->charge() == 0){continue;}
        if (!(iTrack2->charge()==1 || iTrack2->charge()==-1)){continue;}
        if (!iTrack2->trackHighPurity()){continue;}
        const reco::Track* trk2 = iTrack2->bestTrack();
        if (!trk2) {continue;}
        //end检测选中的track是否有效，检测电荷是否是+/-1

        if ( (iTrack1->charge()+iTrack2->charge() ) == 0 )
        {
          //next part in "two tracks fit"
        }
      }//loop for Track2
    }//loop for Track1
  }//if two tracks
  .......
  .......
}//analyze

1.3-two tracks fit

//MuMuEEPat.cc
void MuMuEEPat::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  const double myPimass=0.13957039;
  const double myPimasserr = myPimass * 1e-6;
  const MagneticField &bFieldHandle = iSetup.getData(magneticFieldToken_);//在获取track信息时需要磁场信息
  ......
  if ( (iTrack1->charge()+iTrack2->charge() ) == 0 )
  {
    // Get two tracks information
    TransientTrack trackTT1(*iTrack1->bestTrack(), &(bFieldHandle));
    TransientTrack trackTT2(*iTrack2->bestTrack(), &(bFieldHandle));
    //粒子加工厂，结合track、mass、sigma就可以构造出我们想要的粒子
    KinematicParticleFactoryFromTransientTrack PiPiFactory;
    // pdg mass
    ParticleMass pion_mass = myPimass;
    float pion_sigma = myPimasserr;
    ParticleMass proton_mass = myProtonmass;
    float proton_sigma = myProtonmasserr;
    // initial chi2 and ndf before kinematic fits.
    float chi = 0.;
    float ndf = 0.;
    //构建一个RefCountedKinematicParticle类型的vector，把我们想要拟合的粒子都放进去
    vector < RefCountedKinematicParticle > PiPiParticles;
    PiPiParticles.push_back(PiPiFactory.particle(trackTT1, pion_mass, chi, ndf, pion_sigma));
    PiPiParticles.push_back(PiPiFactory.particle(trackTT2, pion_mass, chi, ndf, pion_sigma));
    //start two tracks vertex fit
    KinematicParticleVertexFitter fitter;
    RefCountedKinematicTree PiPiVertexFitTree;
    PiPiVertexFitTree = fitter.fit(PiPiParticles);
    //顶点拟合有效
    if (PiPiVertexFitTree->isValid())
    {
      //对这个tree提取数据，这个tree的结构是Ks->pi+ pi-
      PiPiVertexFitTree->movePointerToTheTop();
      RefCountedKinematicParticle PiPi_vFit_noMC = PiPiVertexFitTree->currentParticle();
      RefCountedKinematicVertex PiPi_vFit_vertex_noMC = PiPiVertexFitTree->currentDecayVertex();
      KinematicParameters mypipipara = PiPi_vFit_noMC->currentState().kinematicParameters();

      if (PiPi_vFit_vertex_noMC->vertexIsValid())
      {
        //输出我们想要的变量
      }
    }//顶点拟合有效
  }//if two track charge is 0
  ......
}//analyze

在输出变量之前，我们需要仔细思考自己需要用到哪些变量。\(K_s\)粒子的一个非常重要的特征是寿命比较长，我们要抓住这个特点筛选掉干扰我们寻找 \(K_s\)粒子的本底，所以 \(c\tau\)是一个必须要输出的变量，关于 \(c\tau\)的定义在Multilep learning note中的10-mass constraint结尾可以找到。另外，我们是通过两个track来重建\(K_s\)的，因此必须保证两个track是从同一个顶点出来的，而对应的物理量是vertex ChiSquared Probability，这个值越接近1就表示两个track从同一个顶点出来的可能性越大，如果没有Vertex Probability，我们重建出来的 \(\pi^+\pi^-\)质量绝大部分都是本底（垃圾）。

试着把下面这些变量在tree里面都存储下来。

//MuMuEEPat.cc
if (PiPi_vFit_vertex_noMC->vertexIsValid())
{
  float mypipionlyctau=GetcTau(PiPi_vFit_vertex_noMC,PiPi_vFit_noMC,theBeamSpotV);
  float mypipionlyctauerr=GetcTauErr(PiPi_vFit_vertex_noMC,PiPi_vFit_noMC,theBeamSpotV);

  pipionlyctau->push_back(mypipionlyctau);
  pipionlyctauerr->push_back(mypipionlyctauerr);
  //我们必须保证two tracks是从同一个顶点出来的，因此对Vertex Probability的cut是必须的
  double chi2 = PiPi_vFit_vertex_noMC->chiSquared();
  double ndof = PiPi_vFit_vertex_noMC->degreesOfFreedom();
  if (chi2 >=0.0 && ndof >= 0.0)
  {
    pipionlyVtxCL->push_back( ChiSquaredProbability((double)(PiPi_vFit_vertex_noMC->chiSquared()),(double)(PiPi_vFit_vertex_noMC->degreesOfFreedom())));
  }
  else
  {
    pipionlyVtxCL->push_back(-9);
  }

  pipionlyMass->push_back(PiPi_vFit_noMC->currentState().mass());
  if(  PiPi_vFit_noMC->currentState().kinematicParametersError().matrix()(6,6)>0)
  {
    pipionlyMassErr->push_back( sqrt(PiPi_vFit_noMC->currentState().kinematicParametersError().matrix()(6,6)) );
  } 
  else
  {
    pipionlyMassErr->push_back( -9);
  }
  pipionlyPx->push_back( mypipipara.momentum().x() );
  pipionlyPy->push_back( mypipipara.momentum().y() );
  pipionlyPz->push_back( mypipipara.momentum().z() );
  TLorentzVector pipionlyp4;
  pipionlyp4.SetXYZM(mypipipara.momentum().x(), mypipipara.momentum().y(), mypipipara.momentum().z(), PiPi_vFit_noMC->currentState().mass());
  pipionlyPt->push_back( pipionlyp4.Pt());
  pion1Px->push_back( iTrack1->momentum().x() );
  pion1Py->push_back( iTrack1->momentum().y() );
  pion1Pz->push_back( iTrack1->momentum().z() );
  TLorentzVector pion1p4;
  pion1p4.SetXYZM(iTrack1->momentum().x(), iTrack1->momentum().y(), iTrack1->momentum().z(), iTrack1->mass());
  pion1Pt->push_back( pion1p4.Pt());
  pion2Px->push_back( iTrack2->momentum().x() );
  pion2Py->push_back( iTrack2->momentum().y() );
  pion2Pz->push_back( iTrack2->momentum().z() );
  TLorentzVector pion2p4;
  pion2p4.SetXYZM(iTrack2->momentum().x(), iTrack2->momentum().y(), iTrack2->momentum().z(), iTrack2->mass());
  pion2Pt->push_back( pion2p4.Pt());
  ++npipionly;
}//if (PiPi_vFit_vertex_noMC->vertexIsValid())

请记住在完成.h和.cc文件的编写之后需要编译。

scramv1 b ProjectRename 
scramv1 b clean #当你进入到一个新的工作目录的时候你需要上面两行去清空一些缓存
scramv1 b           #当你修改了你的代码的时候用这一句就可以编译了

2-产生multilep.root文件

2.1-cmsRun本地产生

为了检测你编写好的.h和.cc文件是否正确，我们需要先在本地运行处理一定事例。在runMultiLepPAT_dataRun3_miniAOD.py中传入所需要的原始root文件cms Data Aggregation System ，并且保证你的CMSSW Release版本和Global Tag与你传入的原始root文件是一致的，运行

cmsRun runMultiLepPAT_dataRun3_miniAOD.py

如果你运行之后出现了crash的报错，那可能是你的程序的内部逻辑出现了问题，而不是语法错误，这在编译的时候是查不出来的。

如果你的cmsRun可以成功运行，你会发现本地处理event的速度非常慢，三天三夜大概只能处理五万左右的事例，这是因为一个event中包括的track数量是非常庞大的。显然本地处理原始root行不通，因此你需要交到服务器上，让服务器帮助你处理。

2.2-submit crab jobs

这在Multilep learning note中的pre 1-CRAB jobs submit已经讲的非常详细了。

在我们提交crab jobs之后，可能需要等待一两天才能跑完所有的数据，但是你不能等它完全finish才开始下一步，一旦有multilep.root产生就立刻开始下一步myntuple分析，这样可以检查你的Ntuple Maker是否仍然存在其他错误，以节省时间，提高效率。

3-myntuple分析root文件

3.1-本地运行

一个ntuple的job需要myntuple.C.rootmap，myntuple.C，myntuple.h和一个运行的Runjobs.C文件。详细的介绍可以参见Multilep learning note中的5-myntuple分析root文件。这里仅列出我的myntuple.C作为参考。

//myntuple.C
void myntuple::Loop()
{
  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntries();
  Long64_t nbytes = 0, nb = 0;

  //fstream foutFinalEvents;
  //foutFinalEvents.open((outputname + "KsEvents.txt").Data(), std::ios::out);
  TFile* myhbk = new TFile ("myhbk.root","recreate");
  TH1F* pipimass = new TH1F("pipimass","pipimass",100,0.4,0.6);

  for (Long64_t jentry=0; jentry<nentries;jentry++)
  {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    for (unsigned int myi = 0; myi < pipionlyMass->size(); myi++) 
    { 
      if ( (*pipionlyctau)[myi] > 5 && (*pipionlyVtxCL)[myi] > 0.05 && (*pion1Pt)[myi] > 1.5 && (*pion2Pt)[myi] > 1.5) 
      {
        pipimass->Fill((*pipionlyMass)[myi]);

        /*foutFinalEvents <<std::fixed
                << (*pipionlyMass)[myi] << " "
                << (*pipionlyctau)[myi] << " "
                << (*pipionlyVtxCL)[myi] << " "
                << (*pion1Pt)[myi] << " "
                << (*pion2Pt)[myi] << " "
                << std::endl;*/
      }
    }
    //这句话其实非常重要，它可以告诉你处理了多少事例数
    if (jentry%1000 == 0) std::cout << "I am running " << jentry << "th entries out of " << nentries << " total entries" << std::endl;
  }
//foutFinalEvents.close();
myhbk->Write();
}

请记住写好myntuple.C之后，需要进行编译，否则你运行结束产生的图就和上次一模一样。

$ root -l 
.L myntuple.C++

编译完成之后运行

$ root -l -b -q Runjobs.C #-l 是不显示root启动界面，-b 以批处理模式运行，不显示图形界面， -q 处理完命令行宏文件后退出。这个命令通常在需要跑很多文件的时候配合脚本使用

刚开始你产生的候选事例还比较少，看到的图可能长成下面这个样子。这个时候其实是看不出来到底是否存在信号的，需要用更多的multilep.root。

但是在本地myntuple处理数据的过程中，你会发现跑的很慢，大概一万个事例要跑半个小时，效率非常低，因为一个事例的track数目非常庞大。因此需要交condor作业用服务器来帮助你跑。

3.2-submit condor jobs

详细的介绍可以参见Multilep learning note中的pre2-Condor Jobs submit。

4-Where is the \(K_s\)

当你拥有了足够多的事例之后，你可能会看到下面左侧这样的图，但是你会感觉左侧图的bin分的太细了，有点看不清楚是不是有信号，这里有一个小trick

$ pipimass->Rebin(2)

将bin的数量砍成原来的一半，可以得到右侧的图。

我们发现在0.5GeV左右的地方好像是存在 \(K_s\)，但是还不能完全确定。

那到底怎么才能看到非常确定的 \(K_s\)呢？

注意到我们在myntuple.C里设置的筛选条件是

if ( (*pipionlyctau)[myi] > 5 && (*pipionlyVtxCL)[myi] > 0.05 && (*pion1Pt)[myi] > 1.5 && (*pion2Pt)[myi] > 1.5) 

可以把条件设的稍微宽松一些

if ( (*pipionlyctau)[myi] > 2 && (*pipionlyVtxCL)[myi] > 0.05 && (*pion1Pt)[myi] > 1.2 && (*pion2Pt)[myi] > 1.2) 

再加之交更多的crab作业，产生更多的multilep.root，我们最终可以得到下面的图

多么pretty的峰！这下我们心里终于有底了，这确确实实就是我们要找的 \(K_s\)！

5-Fit

当然，找到信号可以说是已经成功了，但是这还没完，下面我们要对信号进行拟合，看看它的质量、宽度、显著度。

易老师：做事就要把它做完整才算成功。

5.1-直接用直方图拟合

拟合的教程在这里可以找到。由于 \(K_s\)的宽度很窄，所以我们采用高斯函数拟合信号，一阶切比雪夫多项式拟合本底。下面仅列出我的fitKs.C作为参考。

#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooPlot.h"
using namespace RooFit;

void fitKs() {
  TFile *f = new TFile("outKs30Oct.root", "READ");
  if (!f || f->IsZombie()) {
      std::cout << "Error: cannot open file!" << std::endl;
      return;
  }
  TH1F *h = (TH1F*) f->Get("pipimass");
  if (!h) {
      std::cout << "Error: histogram not found!" << std::endl;
      return;
  }
  h->SetMinimum(0);
  double xmin = h->GetXaxis()->GetXmin();
  double xmax = h->GetXaxis()->GetXmax();
  TF1 *fitFunc = new TF1("fitFunc", "gaus(0)+cheb1(3)", xmin, xmax);
  double mean  = h->GetMean();
  double sigma = h->GetRMS();
  double peak  = h->GetMaximum();
  fitFunc->SetParameters(peak, mean, sigma, 0.0, 0.0);

  h->Fit(fitFunc, "RQM");
  TCanvas *c1 = new TCanvas("c1", "Gaussian + Chebyshev Fit", 800, 600);
  h->Draw();
  fitFunc->SetLineColor(kRed);
  fitFunc->Draw("same");

  TF1 *gausPart = new TF1("gausPart", "gaus", xmin, xmax);
  gausPart->SetParameters(fitFunc->GetParameter(0),  // [0] amplitude
                          fitFunc->GetParameter(1),  // [1] mean
                          fitFunc->GetParameter(2)); // [2] sigma
  gausPart->SetLineColor(kBlue);
  gausPart->SetLineStyle(2); 
  gausPart->Draw("same");

  TF1 *bkgPart = new TF1("bkgPart", "cheb1", xmin, xmax);
  bkgPart->SetParameters(fitFunc->GetParameter(3),  // [3] constant
                         fitFunc->GetParameter(4)); // [4] linear
  bkgPart->SetLineColor(kGreen+2);
  bkgPart->SetLineStyle(3);
  bkgPart->Draw("same");

  std::cout << "Fit Results:" << std::endl;
  std::cout << "Mean  = " << fitFunc->GetParameter(1) << " ± " << fitFunc->GetParError(1) << std::endl;
  std::cout << "Sigma = " << fitFunc->GetParameter(2) << " ± " << fitFunc->GetParError(2) << std::endl;
}

拟合结果是

Fit Results:
Mean  = 0.497576 ± 0.00027899
Sigma = 0.00618197 ± 0.000371872

看到Mean值，已经没有任何人敢质疑你的 \(K_s\)了。而Sigma和我们看到的信号“宽度”是有关的，这里的“宽度”大概是 \(15MeV\)，但是我们在PDG上却发现 \(\tau_{K_s}\approx0.8954\times10^{−10}s\)，对应的实际宽度 \(\Gamma_{K_s}\approx7.35\times10^{−12}MeV\)，这与我们看到的“宽度”相比少了将近12个数量级，因此我们看到的“宽度”其实不是信号真正的宽度，完全是由探测器的resolution(分辨率) 决定的。

5.2-用RooFit拟合

在上面的拟合中峰的位置其实并没有处在正中间，感觉不太美观，峰的右侧长啥样也还不知道，我们想把 \(\pi^+\pi^-\)的质量设置在 \([0.45, 0.55]GeV\)之间，但是直方图根本没有把 \(0.52GeV\)以上的数据信息囊括进来。所以我们在实际处理的过程中，一般都会在myntuple.C把信息输出到txt文件里，设置的范围较宽松，然后再在做拟合的时候进一步限制变量范围，这样可以方便随时修改取用。

除此之外，我们还想知道这个信号的Significance是多少，这个可以通过做RooFit拟合计算出来。

RooFit拟合在这里可以找到案例。下面仅列出我的fitKs.C作为参考，我采用的是Unbinned fit。

#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooPlot.h"
using namespace RooFit;

void fitKs() {
  int STRATEGY = 1; //0, 1, 2
  int NCPU = 4; //num of cpus
  int FFT_BINS = 10000;
  bool MINOS = false;

  double mxMin = 0.45, mxMax = 0.55;
  int mxBins = 200;
  double mxBinWidth = (mxMax - mxMin) / mxBins;

  TString XTitle, YTitle;
  XTitle.Form("m(#pi^{+}#pi^{-}) (GeV)");
  YTitle.Form("Candidates / %.2f MeV", mxBinWidth * 1000);

  ofstream fitresult("fitoutput/fitResult.txt");

  RooRealVar mx("mx", "mx", mxMin, mxMax);
  RooArgSet variables;
  variables.add(mx);

  RooDataSet data_sig = *RooDataSet::read("outKs30Oct.txt",variables, "Q");//Unbinned fit

  RooRealVar mean("mean", "mean", 0.497, 0.45, 0.55);
  RooRealVar sigma("sigma", "sigma", 0.006, 0.002, 0.01);
  RooAbsPdf* gaus = new RooGaussian("gaus", "gaus", mx, mean, sigma);
  RooRealVar nsig("nsig","nsig",0,0,50000000);
  RooExtendPdf* sig=new RooExtendPdf("sig","sig",*gaus,nsig);

  RooRealVar a0("a0","a0",0.005,-1,1);
  RooAbsPdf* shev=new RooChebychev("shev","shev",mx,a0);
  RooRealVar nbkg("nbkg","nbkg",0,0,50000000);
  RooExtendPdf* bkg=new RooExtendPdf("bkg","bkg",*shev,nbkg);

  RooAddPdf null_Pdf("null_Pdf","null_Pdf",RooArgList(*bkg),RooArgList(nbkg));
  RooAddPdf sig_Pdf("sig_Pdf","sig_Pdf",RooArgList(*sig,*bkg),RooArgList(nsig,nbkg));

  //fit
  RooFitResult *null_result= null_Pdf.fitTo(data_sig, Save(kTRUE), Minos(MINOS), Strategy(STRATEGY), NumCPU(NCPU));
  double edm_null = null_result->edm();
  double minNll_null = null_result->minNll(); //Return minimized -log(L) value.
  int status_null = null_result->status();
  int covQual_null = null_result->covQual();
  double significance;
  if(edm_null<0.001 && status_null==0 && covQual_null==3)
  {
    RooFitResult *sig_result= sig_Pdf.fitTo(data_sig, Save(kTRUE), Minos(MINOS), Strategy(STRATEGY), NumCPU(NCPU));
    double edm_sig = sig_result->edm();
    double minNll_sig = sig_result->minNll(); //Return minimized -log(L) value.
    int status_sig = sig_result->status();
    int covQual_sig = sig_result->covQual();
    if(edm_sig<0.001 && status_sig==0 && covQual_sig==3)
    {  
      significance = sqrt(2*fabs(minNll_null-minNll_sig));

      RooPlot *xframe = mx.frame(RooFit::Title("Fit K_{s}"), Range(mxMin, mxMax), Bins(mxBins));
      data_sig.plotOn(xframe, Name("data"));
      sig_Pdf.plotOn(xframe, Name("model"), LineColor(kBlue), LineStyle(1));
      sig_Pdf.plotOn(xframe, Components(*sig),LineStyle(kDashed),LineColor(2), Name("signal"));
      sig_Pdf.plotOn(xframe, Components(*bkg),LineStyle(kDashed),LineColor(6), Name("background"));

      //=======write fit parameters ======//
      RooArgSet params(mean, sigma);
      RooPlot* paramBox = sig_Pdf.paramOn(xframe,Parameters(params),Layout(0.12, 0.46, 0.9),Format("NEU", AutoPrecision(2)));
      TPaveText *pt = new TPaveText(0.12, 0.65, 0.4, 0.75, "NDC");
      pt->SetFillColor(0);
      pt->SetTextColor(kRed);
      pt->SetTextSize(0.045);
      pt->SetTextFont(42);
      pt->AddText(Form("Significance = %.2f #sigma", significance));
      //=======write fit parameters ======//

      xframe->GetXaxis()->SetTitle(XTitle.Data());
      xframe->GetXaxis()->SetLabelColor(0, 0);
      xframe->GetYaxis()->SetTitle(YTitle.Data());
      xframe->GetYaxis()->SetTitleSize(0.05);
      xframe->GetYaxis()->SetLabelSize(0.05);
      xframe->GetYaxis()->SetTitleOffset(0.8);

      TLegend leg1(0.7, 0.7, 0.9, 0.9);
      leg1.AddEntry(xframe->findObject("data"), "Data", "pe");
      leg1.AddEntry(xframe->findObject("model"), "Fit", "l");
      leg1.AddEntry(xframe->findObject("signal"), "Signal", "l");
      leg1.AddEntry(xframe->findObject("background"), "Bkg", "l");

      RooPlot *pullFrame1 = mx.frame(RooFit::Title(" "), Range(mxMin, mxMax), Bins(mxBins));
      RooHist *pull1 = xframe->pullHist("data", "model");
      pullFrame1->addPlotable(pull1, "p");
      pullFrame1->GetXaxis()->SetTitle(XTitle.Data());
      pullFrame1->GetXaxis()->SetTitleSize(0.15);
      pullFrame1->GetXaxis()->SetLabelSize(0.12);
      pullFrame1->GetYaxis()->SetTitle("Pull");
      pullFrame1->GetYaxis()->SetTitleSize(0.15);
      pullFrame1->GetYaxis()->SetLabelSize(0.1);
      pullFrame1->GetYaxis()->SetTitleOffset(0.2);

      TCanvas c1("c1", "c1", 800, 600);
      c1.cd();
      TPad pad11("pad11", "pad11", 0, 0.3, 1, 1.0);
      pad11.SetTopMargin(0.08);
      pad11.SetBottomMargin(0.017);
      pad11.Draw();
      pad11.cd();
      xframe->Draw();
      leg1.Draw();
      pt->Draw("same");//write fit parameters
      c1.cd();
      TPad pad12("pad12", "pad12", 0.0, 0.0, 1, 0.3);
      pad12.SetTopMargin(0.03);
      pad12.SetBottomMargin(0.325);
      pad12.SetGridx();
      pad12.SetGridy(2);
      pad12.Draw();
      pad12.cd();
      pullFrame1->Draw();
      c1.Update();
      c1.SaveAs("fitoutput/figure/fit_sig.pdf");
    }
  }
  fitresult << std::fixed
     << "NumSig = " << nsig.getVal() << " ± "<< nsig.getError() << " "
     << "Mean  = " << mean.getVal() << " ± " << mean.getError() << " "
     << "Sigma = " << sigma.getVal() << " ±" << sigma.getError() << " "
     << "Significance = " << significance << " "
     << std::endl;

  std::cout << "Mean  = " << mean.getVal() << " ± " << mean.getError() << std::endl;
  std::cout << "Sigma = " << sigma.getVal() << " ± " <<sigma.getError() << std::endl;
  std::cout << "NumSig = " << nsig.getVal() << " ± " << nsig.getError() << std::endl;
  std::cout << "Significance = " << significance << std::endl;
  fitresult.close();
} 

拟合结果已经部分地标在图上了，一个professional的图必须包括Title、X-axis title、Y-axis title、legend、pull分布，所有字体、颜色要合适。

恭喜你！到这儿就已经成功地完整地打通了 \(K_s\)这一关！有没有觉得信心大增？在下一章我们会继续挑战重建 \(\Lambda\)粒子！