C++ TGraph::Set方法代码示例

void DoEvolutions( const TString &sim, Int_t time, Int_t Nbins=1, const TString &options="") { 
  
#ifdef __CINT__  
  gSystem->Load("libptools.so");
#endif

  PGlobals::Initialize();
    
  // Palettes!
  gROOT->Macro("PPalettes.C");

  TString opt = options;
  // cout << "options = " << opt << endl;

  // Load PData
  PData *pData = PData::Get(sim.Data());
  pData->LoadFileNames(time);
  if(!pData->IsInit()) return;

  Bool_t CYL = kFALSE;
  if(sim.Contains("cyl")) { CYL = kTRUE; opt += "cyl"; } 
    
  Bool_t ThreeD = kFALSE;
  if(sim.Contains("3D")) ThreeD = kTRUE; 
  
  // Some plasma constants
  Double_t n0 = pData->GetPlasmaDensity();
  Double_t kp = pData->GetPlasmaK();
  Double_t skindepth = 1.0;
  if(kp!=0.0) skindepth = 1/kp;
  Double_t E0 = pData->GetPlasmaE0();

  // Some initial beam properties:
  Float_t Ebeam = pData->GetBeamEnergy() * PUnits::MeV;
  Float_t gamma = pData->GetBeamGamma();
  Float_t vbeam = pData->GetBeamVelocity();
  
  Double_t rms0 = pData->GetBeamRmsY() * kp;
  if(CYL)  rms0 = pData->GetBeamRmsR() * kp;
  
  // Time in OU
  Float_t Time = pData->GetRealTime();
  // z start of the plasma in normalized units.
  Float_t zStartPlasma = pData->GetPlasmaStart() * kp;
  // z start of the beam in normalized units.
  Float_t zStartBeam = pData->GetBeamStart() * kp;
  
  if(opt.Contains("center")) {
    Time -= zStartPlasma;
    if(opt.Contains("comov"))      // Centers on the head of the beam.
      Time += zStartBeam;
  }

  // Beam charge 2D and 1D histogram (on-axis)
  // ------------------------------------------------------------------
  TH2F *hDen2D = NULL;
  if(pData->GetChargeFileName(1)) {
    char hName[24];
    sprintf(hName,"hDen2D");
    hDen2D = (TH2F*) gROOT->FindObject(hName);
    if(hDen2D) { delete hDen2D; hDen2D = NULL; }

    if(!ThreeD)
      hDen2D = pData->GetCharge(1,opt);
    else
      hDen2D = pData->GetCharge2DSliceZY(1,-1,1,opt+"avg");

    hDen2D->SetName(hName);
    hDen2D->GetXaxis()->CenterTitle();
    hDen2D->GetYaxis()->CenterTitle();
    hDen2D->GetZaxis()->CenterTitle();
    
    if(opt.Contains("comov"))
      hDen2D->GetXaxis()->SetTitle("k_{p}#zeta");
    else
      hDen2D->GetXaxis()->SetTitle("k_{p}z");
    
    if(CYL) 
      hDen2D->GetYaxis()->SetTitle("k_{p}r");
    else
      hDen2D->GetYaxis()->SetTitle("k_{p}y");

    hDen2D->GetZaxis()->SetTitle("n_{b}/n_{0}");

 
  }
  
  // Define ranges from the charge 2D histogram:
  // Binning for 2D histograms:
  // We get this values from the 2D density histogram.
  Int_t   x1Nbin    = hDen2D->GetNbinsX();
  Float_t x1Range   = (hDen2D->GetXaxis()->GetXmax() - hDen2D->GetXaxis()->GetXmin());
  Float_t x1Mid     = (hDen2D->GetXaxis()->GetXmax() + hDen2D->GetXaxis()->GetXmin())/2.;
  Float_t x1Min     = hDen2D->GetXaxis()->GetXmin();
  Float_t x1Max     = hDen2D->GetXaxis()->GetXmax();
  
  Int_t   x2Nbin    = hDen2D->GetNbinsY();      
  Float_t x2Range   = (hDen2D->GetYaxis()->GetXmax() - hDen2D->GetYaxis()->GetXmin());
  Float_t x2Mid     = (hDen2D->GetYaxis()->GetXmax() + hDen2D->GetYaxis()->GetXmin())/2.;
  Float_t x2Min     = x2Mid - x2Range/2;
//.........这里部分代码省略.........

///
/// Find an interpolated x value near a certain bin position of a histogram that is the
/// best estimate for h(x)=y. Interpolates by means of fitting a second grade polynomial
/// to up to five adjacent points. Because that's giving us two solutions, we use the central
/// value and knowledge about if it is supposed to be an upper or lower boundary to pick
/// one.
///
/// \param h - the histogram to be interpolated
/// \param i - interpolate around this bin. Must be a bin such that i and i+1 are above and below val
/// \param y - the y position we want to find the interpolated x for
/// \param central - Central value of the solution we're trying to get the CL interval for.
/// \param upper - Set to true if we're computing an upper interval boundary.
/// \param val - Return value: interpolated x position
/// \param err - Return value: estimated interpolation error
/// \return true, if inpterpolation was performed, false, if conditions were not met
///
bool CLIntervalMaker::interpolatePol2fit(const TH1F* h, int i, float y, float central, bool upper,
		float &val, float &err) const
{
	// cout << "CLIntervalMaker::interpolatePol2fit(): i=" << i << " y=" << y << " central=" << central << endl;
	// check if too close to border so we don't have enough bins to fit
	if ( !( 2 <= i && i <= h->GetNbinsX()-1 ) ) return false;
	// check if all necessary bins are on the same side. They are not always,
	// if e.g. in a plugin there's statistical fluctuations
	if ( binsOnSameSide(i-1, y) && binsOnSameSide(i, y) ){
		//cout << "CLIntervalMaker::interpolatePol2fit() : ERROR : bins i-1, i, and i+1 on same side of y" << endl;
		return false;
	}

	// create a TGraph that we can fit
	TGraph *g = new TGraph(3);
	g->SetPoint(0, h->GetBinCenter(i-1), h->GetBinContent(i-1));
	g->SetPoint(1, h->GetBinCenter(i),   h->GetBinContent(i));
	g->SetPoint(2, h->GetBinCenter(i+1), h->GetBinContent(i+1));

	// see if we can add a 4th and 5th point:
	if ( 3 <= i && i <= h->GetNbinsX()-2 ){
		// add a point to the beginning
		if ( (h->GetBinContent(i-2) < h->GetBinContent(i-1) && h->GetBinContent(i-1) < h->GetBinContent(i))
		  || (h->GetBinContent(i-2) > h->GetBinContent(i-1) && h->GetBinContent(i-1) > h->GetBinContent(i)) ){
			TGraph *gNew = new TGraph(g->GetN()+1);
			gNew->SetPoint(0, h->GetBinCenter(i-2), h->GetBinContent(i-2));
			Double_t pointx, pointy;
			for ( int i=0; i<g->GetN(); i++ ){
				g->GetPoint(i, pointx, pointy);
				gNew->SetPoint(i+1, pointx, pointy);
			}
			delete g;
			g = gNew;
		}
		// add a point to the end
		if ( (h->GetBinContent(i+2) < h->GetBinContent(i+1) && h->GetBinContent(i+1) < h->GetBinContent(i))
		  || (h->GetBinContent(i+2) > h->GetBinContent(i+1) && h->GetBinContent(i+1) > h->GetBinContent(i)) ){
			g->Set(g->GetN()+1);
			g->SetPoint(g->GetN()-1, h->GetBinCenter(i+2), h->GetBinContent(i+2));
		}
	}

	// debug: show fitted 1-CL histogram
	// if ( y<0.1 )
	// // if ( methodName == TString("Plugin") && y<0.1 )
	// {
	//   // TString debugTitle = methodName + Form(" y=%.2f ",y);
	//   // debugTitle += upper?Form("%f upper",central):Form("%f lower",central);
	// 		TString debugTitle = "honk";
	//   TCanvas *c = newNoWarnTCanvas(getUniqueRootName(), debugTitle);
	//   g->SetMarkerStyle(3);
	//   g->SetHistogram(const_cast<TH1F*>(h));
	//   const_cast<TH1F*>(h)->Draw();
	//   g->DrawClone("p");
	// }

	// fit
	TF1 *f1 = new TF1("f1", "pol2", h->GetBinCenter(i-2), h->GetBinCenter(i+2));
	g->Fit("f1", "q");    // fit linear to get decent start parameters
	g->Fit("f1", "qf+");  // refit with minuit to get more correct errors (TGraph fit errors bug)
	float p[3], e[3];
	for ( int ii=0; ii<3; ii++ ){
		p[ii] = f1->GetParameter(ii);
		e[ii] = f1->GetParError(ii);
	}

	// get solution by solving the pol2 for x
	float sol0 = pq(p[0], p[1], p[2], y, 0);
	float sol1 = pq(p[0], p[1], p[2], y, 1);

	// decide which of both solutions to use based on the position of
	// the central value
	// \todo we probably don't need the central value to decide which solution
	// to use. Just take the one closest to the original bin! Can't think of a
	// situation where this should fail. So we don't need the central value in
	// this method at all, to be removed.
	int useSol = 0;
	//if ( (sol0<central && sol1>central) || (sol1<central && sol0>central) ){
		//if ( upper ){
			//if ( sol0<sol1 ) useSol = 1;
			//else             useSol = 0;
		//}
		//else{
			//if ( sol0<sol1 ) useSol = 0;
//.........这里部分代码省略.........