C++ TH1D::GetRMSError方法代码示例

void GaussianProfile::getTruncatedMeanRMS(TH1* hist, float& mean, float& mean_error, float& rms, float& rms_error) {
  int nBins = hist->GetNbinsX();
  double xMin = hist->GetXaxis()->GetXmin();
  double xMax = hist->GetXaxis()->GetXmax();
  //double binWidth = (xMax - xMin) / (double) nBins; //WARNING: this works only if bins are of the same size
  double integral = hist->Integral();

  int maxBin = 0;
  TF1* gaussian = new TF1("gaussian", "gaus");
  fitProjection(hist, gaussian, 1.5, "RQN");
  //maxBin = (int) ceil((gaussian->GetParameter(1) - xMin) / binWidth);
  maxBin = hist->FindBin(gaussian->GetParameter(1));
  delete gaussian;

  TH1D* newHisto = new TH1D("newHisto", "", nBins, xMin, xMax);
  newHisto->SetBinContent(maxBin, hist->GetBinContent(maxBin));
  newHisto->SetBinError(maxBin, hist->GetBinError(maxBin));
  int iBin = maxBin;
  int delta_iBin = 1;
  int sign  = 1;

  while (newHisto->Integral() < 0.99 * integral) {
    iBin += sign * delta_iBin;

    newHisto->SetBinContent(iBin, hist->GetBinContent(iBin));
    newHisto->SetBinError(iBin, hist->GetBinError(iBin));

    delta_iBin += 1;
    sign *= -1;
  }

  rms = newHisto->GetRMS();
  rms_error = newHisto->GetRMSError();

  while (newHisto->Integral() < 0.99 * integral) {
    iBin += sign * delta_iBin;

    newHisto->SetBinContent(iBin, hist->GetBinContent(iBin));
    newHisto->SetBinError(iBin, hist->GetBinError(iBin));

    delta_iBin += 1;
    sign *= -1;
  }

  mean = newHisto->GetMean();
  mean_error = newHisto->GetMeanError();

  delete newHisto;
}

//.........这里部分代码省略.........

	      //if (fabs(y-yt)>5) std::cout << " --- iL=" << iL << " iS=" << iS 
	      //<< " x=" << x << " xt=" << xt 
	      //<< " y=" << y << " yt=" << yt 
	      //<< std::endl;
	      p_posx[ipu][iL][iS]->Fill(x-xt);
	      p_posy[ipu][iL][iS]->Fill(y-yt);
	    }

	  }//loop on layers
	  if (savePoint) p_intercalibSigmaSquare[ipu]->Fill(Etotsq);
	}//loop on entries
    
      }//loop on pt

      //fill first point with linear weighting
      TLatex lat;
      char buf[500];

      for (unsigned iL(0);iL<nLayers;++iL){
	mycFit->cd();
	TH1D *projy = p_deltavsreco_x[ipu][iL]->ProjectionY();
	projy->Draw();
	projy->Fit("gaus","0+Q");
	TF1 *fitx = projy->GetFunction("gaus");
	fitx->SetLineColor(6);
	fitx->Draw("same");
	sprintf(buf,"Layer %d, linear weighting, pu=%d",iL,pu[ipu]);
	lat.DrawLatexNDC(0.1,0.96,buf);
	
	grX[ipu][iL]->SetPoint(0,0.5,fitx->GetParameter(2));
	grX[ipu][iL]->SetPointError(0,0,fitx->GetParError(2));
	grXrms[ipu][iL]->SetPoint(0,0.5,projy->GetRMS());
	grXrms[ipu][iL]->SetPointError(0,0,projy->GetRMSError());
	mycFit->Update();
	if (savePoint) mycFit->Print("PLOTS/fits_x.pdf");
	
	mycFit->cd();
	projy = p_deltavsreco_y[ipu][iL]->ProjectionY();
	projy->Draw();
	projy->Fit("gaus","0+Q");
	TF1 *fity = projy->GetFunction("gaus");
	fitx->SetLineColor(6);
	fitx->Draw("same");
	sprintf(buf,"Layer %d, linear weighting, pu=%d",iL,pu[ipu]);
	lat.DrawLatexNDC(0.1,0.96,buf);
	
	grY[ipu][iL]->SetPoint(0,0.5,fity->GetParameter(2));
	grY[ipu][iL]->SetPointError(0,0,fity->GetParError(2));
	grYrms[ipu][iL]->SetPoint(0,0.5,projy->GetRMS());
	grYrms[ipu][iL]->SetPointError(0,0,projy->GetRMSError());
	mycFit->Update();
	if (savePoint) mycFit->Print("PLOTS/fits_y.pdf");
      }

      //fit vs w0, get w0min
      for (unsigned iL(0);iL<nLayers;++iL){
	for (unsigned iS(0); iS<nScans;++iS){
	  mycFit->cd();
	  p_posx[ipu][iL][iS]->Draw();
	  double w0 = wStart+iS*wStep;
	  //myGaus->SetParameters();
	  p_posx[ipu][iL][iS]->Fit("gaus","0+Q","",-2.,2.);
	  TF1 *fitx = p_posx[ipu][iL][iS]->GetFunction("gaus");
	  fitx->SetLineColor(6);
	  fitx->Draw("same");

void fitTools::getTruncatedMeanAndRMS(TH1D* h1_projection, Float_t& mean, Float_t& mean_err, Float_t& rms, Float_t& rms_err, Double_t percentIntegral_MEAN, Double_t percentIntegral_RMS) {
//TCanvas* getTruncatedMeanAndRMS(TH1D* h1_projection, Float_t& mean, Float_t& mean_err, Float_t& rms, Float_t& rms_err, Double_t percentIntegral_MEAN=0.9, Double_t percentIntegral_RMS=0.68) {

   bool useMode = false;


   if( percentIntegral_MEAN<0. || percentIntegral_MEAN>1. ) {
     std::cout << "WARNING! percentIntegral_MEAN is " << percentIntegral_MEAN << "!! Setting it to 90%." << std::endl;
     percentIntegral_MEAN = 0.9;
   }

   if( percentIntegral_RMS<0. || percentIntegral_RMS>1. ) {
     std::cout << "WARNING! percentIntegral_RMS is " << percentIntegral_RMS << "!! Setting it to 68%." << std::endl;
     percentIntegral_RMS = 0.68;
   }

   Int_t nBins = h1_projection->GetNbinsX();
   Double_t xMin = h1_projection->GetXaxis()->GetXmin();
   Double_t xMax = h1_projection->GetXaxis()->GetXmax();
   Double_t binWidth = (xMax-xMin)/(Double_t)nBins; //WARNING: this works only if bins are of the same size
   Double_t integral = h1_projection->Integral();
//  std::cout << "xmax: " << xMax << "\txMin: " << xMin << std::endl;
  
   //first: find maximum
//  std::cout << "N: " << gaussian->GetParameter(0) << "\tmu: " << gaussian->GetParameter(1) << "\tsigma: " << gaussian->GetParameter(2) << std::endl;
   Int_t maxBin;
   if( useMode ) {
     maxBin = h1_projection->GetMaximumBin();
   } else {
     TF1* gaussian = new TF1("gaussian", "gaus");
     gaussian->SetLineColor(kGreen);
     fitProjection(h1_projection, gaussian, 1.5, "RQN");
     maxBin = (Int_t)ceil((gaussian->GetParameter(1)-xMin)/binWidth);
     delete gaussian;
   }

//  std::cout << "maxBin: " << maxBin << "\tbin center: " << h1_projection->GetXaxis()->GetBinCenter(maxBin) << "\t gauss mu: " << gaussian->GetParameter(1) << std::endl;
   TH1D* newHisto = new TH1D("newHisto", "", nBins, xMin, xMax);
   newHisto->SetBinContent( maxBin, h1_projection->GetBinContent(maxBin) );
   newHisto->SetBinError( maxBin, h1_projection->GetBinError(maxBin) );
   Int_t iBin = maxBin;
   Int_t delta_iBin = 1;
   Int_t sign  = 1;
//  std::cout << "iBin: " << iBin << "\tint: " << newHisto->Integral()/integral << std::endl;

   while( newHisto->Integral() < percentIntegral_RMS*integral ) {

     iBin += sign*delta_iBin; 
     
//  std::cout << "iBin: " << iBin << "\tint: " << newHisto->Integral()/integral << std::endl;
     newHisto->SetBinContent( iBin, h1_projection->GetBinContent(iBin) );
     newHisto->SetBinError( iBin, h1_projection->GetBinError(iBin) );

     delta_iBin += 1;
     sign *= -1;

   }

//  std::cout << "done with rms." << std::endl;
//    TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
//    c1->cd();
//    h1_projection->Draw();
//    newHisto->SetFillColor(kRed);
//    newHisto->DrawClone("HISTO same");

   rms = newHisto->GetRMS();
   rms_err = newHisto->GetRMSError();

//std::cout << "rms: " << rms << std::endl;
   while( newHisto->Integral() < percentIntegral_MEAN*integral ) {
//  std::cout << "iBin: " << iBin << "\tint: " << newHisto->Integral()/integral << std::endl;

     iBin += sign*delta_iBin; 
     
     newHisto->SetBinContent( iBin, h1_projection->GetBinContent(iBin) );
     newHisto->SetBinError( iBin, h1_projection->GetBinError(iBin) );

     delta_iBin += 1;
     sign *= -1;

   }

//    newHisto->SetFillStyle(3004);
//    newHisto->SetFillColor(kBlue);
//    newHisto->DrawClone("HISTO same");

   mean = newHisto->GetMean();
   mean_err = newHisto->GetMeanError();

   delete newHisto;

//    return c1;
}

// ============================================================================
/// finalize the algorithm
// ============================================================================
StatusCode Aida2Root::finalize()
{

  always() << "Get the native ROOT representation of histograms!" << endmsg ;

  {  // loop over all 1D-histograms
    for ( List::const_iterator ipath = m_1Ds.begin() ;
          m_1Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram1D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 1D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
      TH1D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 1D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 1D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;

      info () << " |  Compare       | AIDA/HistoStats |     ROOT/TH1    |      Delta      | "  << endmsg ;
      const std::string format = " | %1$-14.14s | %2$ 15.8g | %3$- 15.8g | %4$= 15.8g | "  ;
      info () << print
        ( Gaudi::Utils::HistoStats::mean        ( aida ) ,
          root->GetMean      ()    , "'mean'"        , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::meanErr     ( aida ) ,
          root->GetMeanError ()    , "'meanErr'"     , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::rms         ( aida ) ,
          root->GetRMS       ()    , "'rms'"         , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::rmsErr      ( aida ) ,
          root->GetRMSError  ()    , "'rmsErr'"      , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::skewness    ( aida ) ,
          root->GetSkewness ()     , "'skewness'"    , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::skewnessErr ( aida ) ,
          root->GetSkewness ( 11 ) , "'skewnessErr'" , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::kurtosis  ( aida ) ,
          root->GetKurtosis ()     , "'kurtosis'"    , format  ) << endmsg ;
      info () << print
        ( Gaudi::Utils::HistoStats::kurtosisErr ( aida ) ,
          root->GetKurtosis ( 11 ) , "'kurtosisErr'" , format  ) << endmsg ;
    }
  }

  { // loop over all 2D-histograms
    for ( List::const_iterator ipath = m_2Ds.begin() ;
          m_2Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram2D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 2D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
      TH2D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 2D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 2D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;
    }
  }

  { // loop over all 3D-histograms
    for ( List::const_iterator ipath = m_3Ds.begin() ;
          m_3Ds.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IHistogram3D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
      if ( sc.isFailure() || 0 == aida )
      { return Error ( "Unable to retrieve 3D-histogram '" + (*ipath) + "'"  ) ; }
      /// convert it to ROOT
	TH3D* root = Gaudi::Utils::Aida2ROOT::aida2root ( aida ) ;
      if ( 0 == root )
      { return Error ( "Unable to convert to ROOT the 3D-histogram '"+(*ipath)+"'") ; }
      /// use the native printout from ROOT
      info() << "The native ROOT printout for 3D-histogram '" << (*ipath) << "':" << endmsg ;
      root->Print() ;
    }
  }


  { // loop over all 1D-profiles
    for ( List::const_iterator ipath = m_1Ps.begin() ;
          m_1Ps.end() != ipath ; ++ipath )
    {
      /// retrieve the historam by full path:
      AIDA::IProfile1D* aida = 0 ;
      StatusCode sc = histoSvc()->retrieveObject( *ipath , aida ) ;
//.........这里部分代码省略.........