C++ TH2::GetBinContent方法代码示例

TH2* doSmooth (TH2* hRaw, bool useLog = true) {
    //
    // smooth histogram
    //
    TH2* hSmooth = (TH2*)hRaw->Clone("hSmooth");
    hSmooth->SetTitle("hSmooth");
    if ( useLog ) {
        for ( int ix=1; ix<=hSmooth->GetNbinsX(); ++ix ) {
            for ( int iy=1; iy<=hSmooth->GetNbinsY(); ++iy ) {
                double c = hSmooth->GetBinContent(ix,iy);
                c = c>0. ? log(c) : -100.;
                hSmooth->SetBinContent(ix,iy,c);
            }
        }
    }
    fitLinear(hSmooth,hRaw);

    if ( useLog ) {
        for ( int ix=1; ix<=hSmooth->GetNbinsX(); ++ix ) {
            for ( int iy=1; iy<=hSmooth->GetNbinsY(); ++iy ) {
                double craw = hRaw->GetBinContent(ix,iy);
                double c = hSmooth->GetBinContent(ix,iy);
                c = craw>0. ? exp(c) : 0.;
                hSmooth->SetBinContent(ix,iy,c);
            }
        }
    }
    hSmooth->SetName(hRaw->GetName());

    return hSmooth;
}

void smoothHisto(TH2* hist, int Niter) {

  for(int iter=0; iter<Niter; iter++) {
    if(iter>2) {
      TH2* clone = (TH2*) hist->Clone("clone");
      for(int ix=2; ix<hist->GetNbinsX()-1; ix++) {
	for(int iy=2; iy<hist->GetNbinsY()-1; iy++) {
	  double nn1 = clone->GetBinContent(ix-1,iy);
	  double nn2 = clone->GetBinContent(ix,iy-1);
	  double nn3 = clone->GetBinContent(ix+1,iy);
	  double nn4 = clone->GetBinContent(ix,iy+1);
	  hist->SetBinContent(ix,iy,0.25*(nn1+nn2+nn3+nn4));
	}
      }
      delete clone;
    } else {
      TH2* clone = (TH2*) hist->Clone("clone");
      for(int ix=2; ix<hist->GetNbinsX()-1; ix++) {
	for(int iy=2; iy<hist->GetNbinsY()-1; iy++) {
	  double nn1 = clone->GetBinContent(ix-1,iy-1);
	  double nn2 = clone->GetBinContent(ix-1,iy+1);
	  double nn3 = clone->GetBinContent(ix+1,iy-1);
	  double nn4 = clone->GetBinContent(ix+1,iy+1);
	  if(nn1+nn2+nn3+nn4>2 && hist->GetBinContent(ix,iy))
	    hist->SetBinContent(ix,iy,1.);
	  else
	    hist->SetBinContent(ix,iy,0.25*(nn1+nn2+nn3+nn4));
	}
      }
      delete clone;      
    }
  }
}

double BTagEffService::GetEfficiency(BTagger const &bTagger, double pt, double eta,
  unsigned flavour) const
{
    // Find the appropriate efficiency histogram. Load it if needed
    auto histGroupIt = effHists.find(bTagger);
    TH2 *hist = nullptr;
    
    if (histGroupIt == effHists.end())
    {
        // Try to load histograms for the given b-tagger
        const_cast<BTagEffService *>(this)->LoadEfficiencies(bTagger);
        hist = effHists.at(bTagger).at(flavour).get();
    }
    else
        hist = histGroupIt->second.at(flavour).get();
    
    
    // Make sure the histogram exists
    if (not hist)
    {
        std::ostringstream message;
        message << "BTagEffService[\"" << GetName() << "\"]::GetEfficiency: " <<
          "Failed to find an efficiency histogram for b tagger " << bTagger.GetTextCode() <<
          ", efficiency label \"" << curEffLabel << "\", jet flavour " << flavour << ".";
        throw std::runtime_error(message.str());
    }

    
    // Return the efficiency
    return hist->GetBinContent(hist->FindFixBin(pt, eta));
}

  void getPlotData() {
    m_name = "TbData -";
    TH2 * h = (TH2*) m_file->Get(m_direc.c_str());

    for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
      m_xs.push_back(h->GetXaxis()->GetBinCenter(i));

    for (int i=0; i<h->GetYaxis()->GetNbins(); i++)
      m_ys.push_back(h->GetYaxis()->GetBinCenter(i));

    for (int i=0; i<h->GetXaxis()->GetNbins(); i++) {
      QVector<double> col;
      for (int j=0; j<h->GetYaxis()->GetNbins(); j++)
        col.push_back(h->GetBinContent(h->GetBin(i, j)));

      m_zs.push_back(col);
    }

    m_plot->m_xAxisTitle = std::string(h->GetXaxis()->GetTitle());
    m_plot->m_yAxisTitle = std::string(h->GetYaxis()->GetTitle());
    m_plot->m_zAxisTitle = std::string(h->GetZaxis()->GetTitle());
    m_plot->m_title = std::string(h->GetTitle());

    std::stringstream ssN;
    m_statsTitles.push_back("N: ");
    m_statsValues.push_back(ssN.str());
  }

void SF_TH2F_And_Eff::init(string filename,string effdata,string effmc,string errordata,string errormc)
{
    f = TFile::Open(filename.c_str() ) ;

    if (f == NULL){
        Log(__FUNCTION__,"ERROR","file '" + filename + "' does not exist");
        throw abortException() ;
    }

    TH2 * hDataEff = getHisto(effdata);
    TH2 * hDataErr = NULL;
    if (errordata != "" ) hDataErr = getHisto(errordata);
    TH2 * hMcEff = getHisto(effmc);
    TH2 * hMcErr = NULL;
    if (errormc != "" ) hMcErr = getHisto(errormc);


    for( int aetabin =1; aetabin <= hDataEff->GetNbinsX() ; ++aetabin)
    for( int ptbin =1; ptbin <= hDataEff->GetNbinsY() ; ++ptbin)
    {
        float ptmin = hDataEff->GetYaxis()->GetBinLowEdge(ptbin); 
        float ptmax = hDataEff->GetYaxis()->GetBinLowEdge(ptbin+1); 
        float aetamin = hDataEff->GetXaxis()->GetBinLowEdge(aetabin); 
        float aetamax = hDataEff->GetXaxis()->GetBinLowEdge(aetabin+1); 
        float effData = hDataEff->GetBinContent(aetabin,ptbin);
        float errData = 0.0 ;
        if (hDataErr) err=hDataErr->GetBinContent(aetabin,ptbin);
        else errData = hDataEff->GetBinError(aetabin,ptbin);

        float effMc = hMcEff->GetBinContent(aetabin,ptbin);
        float errMc = 0.0; 
        if (hMcEff) errMc = hMcEff->GetBinContent(aetabin,ptbin);
        else errMc = hMcEff->GetBinError(aetabin,ptbin);

        if (ptbin == hDataEff->GetNbinsY() ) ptmax = 8000.; // highest is open, current recommendation
        if (aetabin == hDataEff->GetNbinsX()  ) aetamax = aetamax+0.0001; // put it slightly larger  to get 2.4 as well
        add(ptmin,ptmax,aetamin,aetamax,effData,effMc,errData,errMc);
    }
    f->Close(); // delete?
    delete f;
    f = NULL;
}

// normalize migmatrix column-wise
TH2* normalizeMigMat(TH2* h)
{
  TH2* hclone = (TH2*) h->Clone();
  const int xbins = hclone->GetNbinsX();
  const int ybins = hclone->GetNbinsY();
    
  for(int x=0; x<xbins; x++)  
  {
    double integ = hclone->Integral(x+1, x+1, 1, ybins);
    for(int y=0; y<ybins; y++)
    {
      hclone->SetBinContent(x+1,y+1, hclone->GetBinContent(x+1, y+1)/integ);
    }
  }
  
  return hclone;
}

void MirrorBorders( TH2& hist ) {
    int numx = hist.GetNbinsX();
    int numy = hist.GetNbinsY();

    Float_t val;
    // corner points
    hist.SetBinContent(0,0,hist.GetBinContent(1,1));
    hist.SetBinContent(numx+1,numy+1,hist.GetBinContent(numx,numy));
    hist.SetBinContent(numx+1,0,hist.GetBinContent(numx,1));
    hist.SetBinContent(0,numy+1,hist.GetBinContent(1,numy));

    for(int i=1; i<=numx; i++){
        hist.SetBinContent(i,0,	   hist.GetBinContent(i,1));
        hist.SetBinContent(i,numy+1, hist.GetBinContent(i,numy));
    }
    for(int i=1; i<=numy; i++) {
        hist.SetBinContent(0,i,      hist.GetBinContent(1,i));
        hist.SetBinContent(numx+1,i, hist.GetBinContent(numx,i));
    }
}

TH2* smoothHoles(const TH2* originalHist)
{
  TH2* hist = (TH2*)originalHist->Clone("_smoothed");
  int xMax = hist->GetNbinsX();
  int yMax = hist->GetNbinsY();

  int xMin = 0;
  
  for(int xBin = 1; xBin <= xMax; xBin++)
    {
      for(int yBin = 1; yBin <= yMax; yBin++)
	{
	  if(hist->GetBinContent(xBin,yBin)>0)
	    {
	      xMin = xBin;
	      yBin = yMax+1;
	      xBin = xMax+1;
	    }
	}
    } 
  //  for(unsigned int i = 0; i< 1000; i++) smoothHistAcross(hist,xMin);
  for(unsigned int i = 0; i< 1000; i++) interpolateHistAcross(hist,xMin);
  return hist;
}

//.........这里部分代码省略.........
      if(iChannel > double(readoutMap->GetEntries())/reductionFactor) break;
      iChannel++;
      TString objName;
      uint32_t fedKey =  SiStripFedKey(*fedId,SiStripFedKey::feUnit(*fedCh),SiStripFedKey::feChan(*fedCh)).key();
      std::stringstream stream;
      stream << std::hex << fedKey;
      string fedKeyStr = stream.str();
      if(fedKeyStr.size() == 4)
	objName = Form("DQMData/SiStrip/ControlView/FecCrate%d/FecSlot%d/FecRing%d/CcuAddr%d/CcuChan%d/ExpertHisto_PedsFullNoise_FedKey0x0000%s_LldChannel%d_Noise2D",*fecCrate,*fecSlot,*fecRing,*ccuAdd,*ccuChan,fedKeyStr.c_str(),*lldChannel);      
      else if(fedKeyStr.size() == 5)
	objName = Form("DQMData/SiStrip/ControlView/FecCrate%d/FecSlot%d/FecRing%d/CcuAddr%d/CcuChan%d/ExpertHisto_PedsFullNoise_FedKey0x000%s_LldChannel%d_Noise2D",*fecCrate,*fecSlot,*fecRing,*ccuAdd,*ccuChan,fedKeyStr.c_str(),*lldChannel);      
      else
	cerr<<"hex number to short "<<fedKeyStr<<" --> please check "<<endl;

      inputFile->GetObject(objName.Data(),histoNoise);
      // extract single strip noise histogram --> loop on the y-axis
      uint16_t apvID = 0;
      uint16_t stripID = 0;       
      if(histoNoiseStrip == 0 or histoNoiseStrip == NULL){
	histoNoiseStrip = new TH1F ("histoNoiseStrip","",histoNoise->GetNbinsX(),histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	histoNoiseStrip->Sumw2();
      }
      for(int iBinY = 0; iBinY < histoNoise->GetNbinsY(); iBinY++){
	histoNoiseStrip->Reset();
	histoNoiseStrip->SetDirectory(0);
	// two multiplexed APV per line
	if(iBinY < histoNoise->GetNbinsY()/2) apvID = 1;
	else apvID = 2;
	// strip id
	stripID++;
	if(stripID > 128) stripID = 1;
	// loop on x-axis bin
	for(int iBinX = 0; iBinX < histoNoise->GetNbinsX(); iBinX++){
	  histoNoiseStrip->SetBinContent(iBinX+1,histoNoise->GetBinContent(iBinX+1,iBinY+1));
	  histoNoiseStrip->SetBinError(iBinX+1,histoNoise->GetBinError(iBinX+1,iBinY+1));	    
	}
     	
	// to initialize branches
	detid_ = 0; fedKey_ = 0; fecCrate_ = 0; fecSlot_ = 0; fecRing_ = 0; ccuAdd_ = 0; ccuChan_ = 0; lldChannel_ = 0; fedId_ = 0; fedCh_ = 0; apvId_ = 0; stripId_ = 0; 
	noiseMean_ = 0.; noiseRMS_ =  0.; noiseSkewness_ = 0.; noiseKurtosis_ = 0.; 
	fitGausMean_ = 0.; fitGausSigma_ = 0.;fitGausNormalization_ = 0.;
	fitGausMeanError_ = 0.; fitGausSigmaError_ = 0.;fitGausNormalizationError_ = 0.;	  	  
	fitChi2_ = 0.; fitChi2Probab_ = 0.; fitStatus_ = -1.; 
	noiseIntegral3Sigma_ = 0.; noiseIntegral3SigmaFromFit_ = 0.; 
	noiseIntegral4Sigma_ = 0.; noiseIntegral4SigmaFromFit_ = 0.; 
	noiseIntegral5Sigma_ = 0.; noiseIntegral5SigmaFromFit_ = 0.; 
	kSProbab_ = 0.; jBProbab_ = 0.;
	kSValue_ = 0.; jBValue_ = 0.; 
	aDValue_= 0.; aDProbab_ = 0.;
	nBin_ = 0.; xMin_ = 0.; xMax_ = 0.;
	
	// basic info
	detid_ = *detid;
	fedKey_ = fedKey;
	fecCrate_ = *fecCrate;
	fecSlot_ = *fecSlot;
	fecRing_ = *fecRing;
	ccuAdd_  = *ccuAdd;
	ccuChan_ = *ccuChan;
	lldChannel_ = *lldChannel;
	fedId_   = *fedId;
	fedCh_   = *fedCh;
	apvId_   = apvID;
	stripId_ = stripID;
	
	// basic info of nioise distribution