C++ TRandom3::BreitWigner方法代码示例

void TSVDUnfoldExample()
{
   gROOT->Reset();
   gROOT->SetStyle("Plain");
   gStyle->SetOptStat(0);

   TRandom3 R;

   const Double_t cutdummy= -99999.0;

   // --- Data/MC toy generation -----------------------------------

   // The MC input
   Int_t nbins = 40;
   TH1D *xini = new TH1D("xini", "MC truth", nbins, -10.0, 10.0);
   TH1D *bini = new TH1D("bini", "MC reco", nbins, -10.0, 10.0);
   TH2D *Adet = new TH2D("Adet", "detector response", nbins, -10.0, 10.0, nbins, -10.0, 10.0);

   // Data
   TH1D *data = new TH1D("data", "data", nbins, -10.0, 10.0);
   // Data "truth" distribution to test the unfolding
   TH1D *datatrue = new TH1D("datatrue", "data truth", nbins, -10.0, 10.0);
   // Statistical covariance matrix
   TH2D *statcov = new TH2D("statcov", "covariance matrix", nbins, -10.0, 10.0, nbins, -10.0, 10.0);

   // Fill the MC using a Breit-Wigner, mean 0.3 and width 2.5.
   for (Int_t i= 0; i<100000; i++) {
      Double_t xt = R.BreitWigner(0.3, 2.5);
      xini->Fill(xt);
      Double_t x = Reconstruct( xt, R );
      if (x != cutdummy) {
         Adet->Fill(x, xt);
         bini->Fill(x);
      }
   }

   // Fill the "data" with a Gaussian, mean 0 and width 2.
   for (Int_t i=0; i<10000; i++) {
      Double_t xt = R.Gaus(0.0, 2.0);
      datatrue->Fill(xt);
      Double_t x = Reconstruct( xt, R );
      if (x != cutdummy)
      data->Fill(x);
   }

   cout << "Created toy distributions and errors for: " << endl;
   cout << "... \"true MC\"   and \"reconstructed (smeared) MC\"" << endl;
   cout << "... \"true data\" and \"reconstructed (smeared) data\"" << endl;
   cout << "... the \"detector response matrix\"" << endl;

   // Fill the data covariance matrix
   for (int i=1; i<=data->GetNbinsX(); i++) {
       statcov->SetBinContent(i,i,data->GetBinError(i)*data->GetBinError(i));
   }

   // --- Here starts the actual unfolding -------------------------

   // Create TSVDUnfold object and initialise
   TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );

   // It is possible to normalise unfolded spectrum to unit area
   tsvdunf->SetNormalize( kFALSE ); // no normalisation here

   // Perform the unfolding with regularisation parameter kreg = 13
   // - the larger kreg, the finer grained the unfolding, but the more fluctuations occur
   // - the smaller kreg, the stronger is the regularisation and the bias
   TH1D* unfres = tsvdunf->Unfold( 13 );

   // Get the distribution of the d to cross check the regularization
   // - choose kreg to be the point where |d_i| stop being statistically significantly >>1
   TH1D* ddist = tsvdunf->GetD();

   // Get the distribution of the singular values
   TH1D* svdist = tsvdunf->GetSV();

   // Compute the error matrix for the unfolded spectrum using toy MC
   // using the measured covariance matrix as input to generate the toys
   // 100 toys should usually be enough
   // The same method can be used for different covariance matrices separately.
   TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );

   // Now compute the error matrix on the unfolded distribution originating
   // from the finite detector matrix statistics
   TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );

   // Sum up the two (they are uncorrelated)
   ustatcov->Add( uadetcov );

   //Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics.
   TH2D* utaucov = tsvdunf->GetXtau();
   utaucov->Add( uadetcov );

   //Get the computed inverse of the covariance matrix
   TH2D* uinvcov = tsvdunf->GetXinv();


   // --- Only plotting stuff below ------------------------------

   for (int i=1; i<=unfres->GetNbinsX(); i++) {
      unfres->SetBinError(i, TMath::Sqrt(utaucov->GetBinContent(i,i)));
//.........这里部分代码省略.........