C++ TGraphErrors::SetPoint方法代码示例

TGraphErrors* ExtractCalCurve(vector<calData>& FitPositions, vector<Float_t>& Energies)
//void ExtractCalCurve(vector<Float_t>& FitPositions, vector<Float_t>& Energies)
{
    //TF1* calFit = new TF1("calFit","[0]*x + [1]",0,800E3);
    TGraphErrors* calCurve = new TGraphErrors(Energies.size());
    
    //Get highest energy peak and estimate [0]
    std::vector<calData>::iterator pit;
    std::vector<Float_t>::iterator eit;
    
    pit = max_element(FitPositions.begin(),FitPositions.end(),CompareByadc);
    eit = max_element(Energies.begin(),Energies.end());
    
    //Extrapolate coefficient and place first point
    int currPoint = 1;
    calCurve->SetPoint(currPoint,(*pit).fadc,*eit);
    calCurve->SetPointError(currPoint,(*pit).efadc,0.0);
    TFitResultPtr r = calCurve->Fit("pol1","SQ");   
    Float_t a = r->Parameter(1);

    Float_t CurrentPeak = 0.0;
    Float_t CurrentEnergy = 0.0;
    Float_t CurrentEnergyEst = 0.0;
    // Loop through found peaks and locate closest estimated energy
    // Assume fitted peaks are already ordered from lowest to highest
    for(std::vector<calData>::iterator i = --(FitPositions.end()); i!=FitPositions.begin(); --i)
    {
        currPoint++;
        CurrentPeak = (*i).fadc;
        CurrentEnergyEst = CurrentPeak*a;
        Float_t CurrentDelta = 800E3;
        for(std::vector<Float_t>::iterator j = Energies.begin(); j!=Energies.end(); j++)
	{
		if( abs(*j - CurrentEnergyEst) < CurrentDelta)
		{
			CurrentDelta = abs(*j - CurrentEnergyEst);
                        CurrentEnergy = *j;
		}
	} 
	
	calCurve->SetPoint(currPoint,CurrentPeak,CurrentEnergy);
        calCurve->SetPointError(currPoint,(*i).efadc,CurrentDelta);
        r = calCurve->Fit("pol1","SQ");
        a = r->Parameter(1);
    }
    r->Print("V");
    return calCurve; 
}

int RampCorrellation_single(std::string file_name, TGraphErrors &graph, TGraphErrors &graph2, vector<double> &vec_errors)
{

	TTree *t = new TTree();
	const std::string file_path = "/home/jlab/github/MagCloak_Analysis/calibration/data-calib/DATA_MegaVIEW/";
//	std::string file_name = "DataFile_2016-12-08_06-59-11.csv";
	std::string file = file_path+file_name;

	t->ReadFile(file.c_str());
//	t->Print();

	TCanvas *c = new TCanvas();
	t->Draw("TMath::Abs(B1/B2):time");
	TGraph *gh = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double ratio_mean = gh->GetMean(2);
	double ratio_std = gh->GetRMS(2);

	t->Draw("Bnom:time");
	TGraph *gh1 = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double nom_mean = gh1->GetMean(2);
	double nom_std = 0;

//	cout << "At " << nom_mean << " mT, B1/B2 is: " << ratio_mean << " +/- " << ratio_std << endl;
	
	c->Close();

	int n = graph.GetN();
	graph.SetPoint(n,nom_mean,ratio_mean);
	graph.SetPointError(n,nom_std,ratio_std);

//-------------------------------------------------------------------------------

	TCanvas *c2 = new TCanvas();
	t->Draw("B3:time");
	TGraph *gh2 = (TGraph*)c2->GetListOfPrimitives()->FindObject("Graph");
	double B3_mean = -1*(gh2->GetMean(2));
	double B3_std = gh2->GetRMS(2);
	
	int n2 = graph2.GetN();
	graph2.SetPoint(n2, nom_mean, B3_mean);
	graph2.SetPointError(n2, nom_std, B3_std);
	
	c2->Close();

	if(nom_mean < 500) vec_errors.push_back(B3_std);

	return 0;
}

//______________________________________________________________________________
TGraphErrors *GraphRatio(TGraphErrors *g1,TGraphErrors *g2)
{
  Int_t n = g1->GetN();
  if (n != g2->GetN()) {
    Printf("!!!!!!!!!!!!!!!!");
    return 0;
  }

  TGraphErrors *g = new TGraphErrors(n);
  g->SetMarkerStyle(20);
  g->GetXaxis()->SetTitle("p_{t}, GeV/c");
  g->SetTitle(Form("#frac{%s}{%s}", g1->GetTitle(), g2->GetTitle()));
  Double_t x, y1, y2;
  Double_t c0, c1, c12, e0, e1;
  Double_t *ex = g1->GetEX();
  Double_t *ey = g1->GetEY();
  for (Int_t i = 0; i < n; i++) {
    g1->GetPoint(i, x, y1);
    g2->GetPoint(i, x, y2);
    if (TMath::Abs(y2) < 1.E-12) y2 = 1.0;
    c0 = y1;
    c1 = y2;
    c12 = c1*c1;
    e0  = g1->GetErrorY(i);
    e1  = g2->GetErrorY(i);
    g->SetPoint(i, x, y1/y2);
    g->SetPointError(i, ex[i], (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12));
  }

  return g;
}

TGraphErrors* addGraph(TGraphErrors *g1, TGraphErrors *g2) {
  
  assert(g1->GetN()==g2->GetN());
  TGraphErrors *g = (TGraphErrors*)g1->Clone();//new TGraphErrors(g1->GetN());
  for (int i = 0; i != g1->GetN(); ++i) {

    double x1 = g1->GetX()[i];
    double x2 = g2->GetX()[i];
    assert(fabs(x1/x2-1)<0.1);

    double ex1 = g1->GetEX()[i];
    double ex2 = g2->GetEX()[i];

    double y1 = g1->GetY()[i];
    double y2 = g2->GetY()[i];

    double ey1 = g1->GetEY()[i];
    double ey2 = g2->GetEY()[i];

    double x = (ey2*ey2*x1 + ey1*ey1*x2) / (ey2*ey2+ey1*ey1);
    double y = (ey2*ey2*y1 + ey1*ey1*y2) / (ey2*ey2+ey1*ey1);
    double ex = (ey2*ey2*ex1 + ey1*ey1*ex2) / (ey2*ey2+ey1*ey1);
    double ey = (ey2*ey2*ey1 + ey1*ey1*ey2) / (ey2*ey2+ey1*ey1);

    g->SetPoint(i, x, y);
    g->SetPointError(i, ex, ey);
  }

  return g;
}

void add_point( double coefover, double sigma, double dsigma ){//============================================================


      TGraphErrors *fisig;
      if ( gROOT->GetListOfSpecials()->FindObject("fisig")!=NULL ){
	fisig=(TGraphErrors*) gROOT->GetListOfSpecials()->FindObject("fisig");
	printf("fisig exists\n%s","");
      }else{

	printf("\n\n NO fisig exists NONONO \n\n\n%s","");

	fisig=new TGraphErrors();
	fisig->SetTitle("fisig");
	fisig->SetName("fisig");
	gROOT->GetListOfSpecials()->Add( fisig );// TRICK TO KEEP IT
	gDirectory->Add( fisig );
	fisig->SetMarkerStyle(7);
      }


      int n=fisig->GetN();
      fisig->SetPoint(n,coefover,sigma);
      fisig->SetPointError(n,0.0,dsigma);
      //      fisig->Print(); // too many data 

     char text[100];
      sprintf(text,"%.3f ... %d",coefover, fisig->GetN()  );

      //      TText *t=new TText(from,0,text ); 
      //      t->SetTextSize(0.08);t->Draw();
 

}//============================================================

void ExpManager::GetExp1DGraphX(TString NameTitle, double zmin, double zmax, double ymin, double ymax, TString grid){   // do the same thing as for the polar interpolation

    TGraphErrors *fGraph = new TGraphErrors(); //= new TGraph2DErrors(np, x_array, y_array, bz_array, ex, ey, ez);
    fGraph->SetTitle("Experimental Data;X (mm);Magnetic Field (mT)");
    fGraph->SetMarkerSize(1.2);
    fGraph->SetMarkerStyle(20);
    fGraph->SetMarkerColor(kBlue);
    fGraph->SetLineColor(kBlue);
    fGraph->SetLineWidth(2);
       
    int graph_counter = 0 ; 
   for (unsigned i=0; i< fExpY.size(); i++)   {
    //cout <<  " X "  << fExpX.at(i) ;  
        if( (fExpY.at(i) >= ymin && fExpY.at(i) <= ymax) && (fExpZ.at(i) >= zmin && fExpZ.at(i) <= zmax)  && fGrid.at(i)== grid ){
            cout << "  < ----- " ; 
            fGraph->SetPoint(graph_counter,fExpX.at(i),fExpB.at(i));    // CHECK, add new the stuff for emag
            fGraph->SetPointError(graph_counter,fExpXErr.at(i),fExpBErr.at(i));   // CHECK, add new the stuff for emag
            graph_counter++;
        }
        cout << endl ; 
   }

    fGraph->SetTitle(NameTitle+Form(" Experimental Data : %.2f < Depth < %.2f mm  __  %.2f < Y < %.2f mm;X (mm);Magnetic Field (mT)",zmin,zmax,ymin,ymax));
    fGraph->SetName(NameTitle+Form("_Exp_Depth_%.2f_%.2fmm_Y_%.2f_%.2fmm",zmin,zmax,ymin,ymax));
    fGraph->Write();

}

TGraphErrors bestFit(TTree *tree, TString alphaName, TString constTermName){
  Double_t alpha, constTerm;
  alphaName.ReplaceAll("-","_");
  constTermName.ReplaceAll("-","_");

  tree->SetBranchAddress(alphaName,&alpha);
  tree->SetBranchAddress(constTermName,&constTerm);

  //Long64_t nEntries=genTree->GetEntries();

  TGraphErrors graph;
  Int_t iPoint=0;

  tree->GetEntry(0);
  graph.SetPoint(0, constTerm, alpha);
  graph.SetPointError(0,0, 0);
  iPoint++;

  graph.Set(iPoint);
  tree->ResetBranchAddresses();
  graph.Draw("A P");
//   graph.SetFillColor(kBlue);
//   graph.SetLineColor(kYellow);
//   graph.GetXaxis()->SetTitle("Energy [GeV]");
//   graph.GetYaxis()->SetTitle("Additional smearing [%]");
  return graph;
}

void draw() {

  gStyle->SetOptStat(0000);

  TGraphErrors* gr = new TGraphErrors("papel_um.txt", "%lg %lg %lg %lg");
  gr->Draw("A");

  gr->GetXaxis()->SetTitle("Diametro (mm)");
  gr->GetYaxis()->SetTitle("Massa (u.m.)");
  gr->SetTitle("Dimensao Fractal - Papel");
  gr->SetMarkerStyle(20);
  gr->SetMarkerSize(gr->GetMarkerSize()/2.);

  TCanvas* c1 = new TCanvas();
  gr->Draw("AP");
  c1->Print("gr1.png");


  TGraphErrors* grlog = new TGraphErrors();
  for (int i=0; i<gr->GetN(); i++) {
    grlog->SetPoint(i, TMath::Log(gr->GetX()[i]), TMath::Log(gr->GetY()[i]));
    grlog->SetPointError(i, gr->GetEX()[i]/gr->GetX()[i], 0.0);
  }
  grlog->Draw("A");
  grlog->GetXaxis()->SetTitle("Logaritmo do diametro (diametro em mm)");
  grlog->GetYaxis()->SetTitle("Logaritmo da massa (massa em unidade arbitraria)");
  grlog->SetTitle("Dimensao Fractal - Papel");
  grlog->SetMarkerStyle(20);
  //  grlog->SetMarkerSize(gr->GetMarkerSize()/2.);
  grlog->Fit("pol1");
  TCanvas* c2 = new TCanvas();
  grlog->Draw("AP");
  c2->Print("gr2.png");

}

TGraphErrors g(TTree *genTree, TString constTermName=""){
  Double_t alpha, constTerm;
  genTree->SetBranchAddress("alpha",&alpha);
  genTree->SetBranchAddress("constTerm",&constTerm);

  Long64_t nEntries=genTree->GetEntries();

  TH1F h("smearingHist","",10000,0,0.1);
  TGraphErrors graph;
  Int_t iPoint=0;
  for(Double_t energy=20; energy<150; energy+=10){
    h.Reset();
    for(Long64_t jentry=0; jentry<nEntries; jentry++){
      genTree->GetEntry(jentry);
      h.Fill(sqrt(alpha*alpha/energy+constTerm*constTerm));
    }
    graph.SetPoint(iPoint, energy, h.GetMean());
    graph.SetPointError(iPoint,0, h.GetRMS());
    iPoint++;
  }
  h.SetTitle(constTermName);
  h.SaveAs("tmp/h-"+constTermName+".root");
  graph.Set(iPoint);
  genTree->ResetBranchAddresses();
  graph.Draw("A L3");
  graph.SetFillColor(kBlue);
  graph.SetLineColor(kYellow);
  graph.GetXaxis()->SetTitle("Energy [GeV]");
  graph.GetYaxis()->SetTitle("Additional smearing [%]");
  return graph;
}

TGraphErrors* getMoliereGraph( std::vector<DataFile> dataFiles ) {

  TGraphErrors* graph = new TGraphErrors(0);

  for( unsigned i=0; i<dataFiles.size(); ++i ) {

    if( dataFiles[i].file==0 ) continue;

    Double_t x = (i==0) ? -1 : dataFiles[i].absThickness;

    Double_t y = ((TH1D*)(dataFiles[i].file->Get("Rmoliere")))->GetMean();
    Double_t yerr = ((TH1D*)(dataFiles[i].file->Get("Rmoliere")))->GetMeanError();

    graph->SetPoint(i, x, y);
    graph->SetPointError(i, 0, yerr);

    if( i==0 ) std::cout << "LYSO Moliere radius: " << y << std::endl;

  }

  graph->SetName( "mol" );
  graph->SetTitle( "Moliere Radius [mm]" );

  return graph;

}

TGraphErrors *TruePoint(TH1 *hT, TH1 *hp, int t, double y)
{
  double ymin = hp->GetMinimum();
  double ymax = hp->GetMaximum();
  double yDrawData = y>0? y : ymin + 0.2*(ymax-ymin);
  TGraphErrors *truePoint = new TGraphErrors();
  truePoint->SetPoint(0, hT->GetBinContent(t+1), yDrawData);
  SetGraphProps(truePoint, kRed+2, kGray, kRed+2, kOpenCircle, 1.5);
  return truePoint;
}

int main(int narg,char **arg)
{
  int i;
  double x,y;
  ifstream file;
  TApplication myapp("App",NULL,NULL);
  TCanvas tela;
  TGraphErrors graf;
  TF1 fun("fun","[0]*exp(-x*x/(2*[1]*[1]))");

  if(narg<2)
    {
      cerr<<"usa: "<<arg[0]<<" file"<<endl;
      exit(0);
    }

  file.open(arg[1]);

  if(file.good()==0)
    {
      cout<<"File: "<<arg[1]<<" inesistente!"<<endl;
      exit(0);
    }

  i=0;
  double dum,N=1;
  while(file>>x>>dum>>y)
    {
      if(i==0) N=y;
      graf.SetPoint(i,x,y);
      graf.SetPointError(i,0,0.0001);
      i++;
    }
  
  fun.SetParameter(0,N);
  fun.SetParameter(1,1);

  fun.SetParName(0,"N");
  fun.SetParName(1,"S");

  graf.Draw("AP");

  graf.Fit("fun");

  //graf.GetXaxis()->SetLimits(0,x);

  tela.Modified();
  tela.Update();

  myapp.Run(true);

  return 0;
}

///
/// Make an evolution graph for one parameter.
///
TGraphErrors* ParameterEvolutionPlotter::makeEvolutionGraphErrors(vector<RooSlimFitResult*> results, TString parName)
{
	TGraphErrors *g = new TGraphErrors(results.size());
	int iGraph = 0;
	for ( int i=0; i<results.size(); i++ ){
		if ( results[i] ){
			g->SetPoint(iGraph, results[i]->getParVal(scanVar1), results[i]->getParVal(parName));
			g->SetPointError(iGraph, 0, results[i]->getParErr(parName));
			iGraph++;
		}
	}
	return g;
}

void OneBin(Int_t bin, 
	    TH1* h0900,
	    TH1* h2760, 
	    TH1* h7000, 
	    TH1* h8000, 
	    TMultiGraph* mg,
	    TNtuple* tuple,
	    Double_t sysErr=0.076)
{
  Info("OneBin", "Getting one bin %d,%p,%p,%p,%p,%p",
       bin, h0900,h2760,h7000,h8000,tuple);
  
  Double_t eta = h0900->GetXaxis()->GetBinCenter(bin);
  Double_t w   = h0900->GetXaxis()->GetBinWidth(bin);
  Info("", "Eta=%f +/- %f", eta, w);
  Double_t e[] = { 900.,  2760., 7000., 8000., 0 };
  TH1*     h[] = { h0900, h2760, h7000, h8000, 0 };
  Float_t  x[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  x[0] = eta;
  x[1] = w;

  TGraphErrors* g = new TGraphErrors(0);
  g->SetName(Form("eta%03d", bin));
  g->SetTitle(Form("%f", eta));
  g->SetMarkerStyle(bin % 10 + 20);
  g->SetMarkerColor(bin % 6 + 2);
  
  Double_t* pe = e;
  TH1**     ph = h;
  Int_t     i  = 0;
  Int_t     j  = 1;
  while (*pe && *ph) { 
    Double_t c = (*ph)->GetBinContent(bin);
    Double_t v = sysErr*c;
    if (c > 1e-6){
      g->SetPoint(i, *pe, c);
      g->SetPointError(i, w, v);
      x[Int_t(2*j+0)] = c;
      x[Int_t(2*j+1)] = v;
      i++;
    }
    j++;
    pe++;
    ph++;
  }
  if (tuple) tuple->Fill(x);
  if (i > 0) mg->Add(g);
  else delete g;
}

///
/// Plot the discrepancy between the observable and the predicted
/// observable when making predictions about observables by scanning
/// them. This checks if the chi2 term of the observable is tight enough.
/// This only works for 1D scans for now.
///
void ParameterEvolutionPlotter::plotObsScanCheck()
{
	vector<RooSlimFitResult*> results = curveResults;

	cout << "ParameterEvolutionPlotter::plotObsScanCheck() : plotting ..." << endl;
	TCanvas *c2 = newNoWarnTCanvas("plotObsScanCheck"+getUniqueRootName(), title, 800,600);
	c2->SetLeftMargin(0.2);

	// get observable
	TGraphErrors *g = new TGraphErrors(results.size());
	int iGraph = 0;

	for ( int i=0; i<results.size(); i++ ){
		assert(results[i]);
		// get value of observable
		float obsValue = results[i]->getParVal(scanVar1);
		float obsError = w->var(scanVar1)->getError();

		// get value of theory prediction
		setParameters(w,parsName,results[i]);
		TString thName = scanVar1;
		thName.ReplaceAll("_obs","_th");
		if ( !w->function(thName) ){
			cout << "ParameterEvolutionPlotter::plotObsScanCheck() : ERROR : theory value not found: " << thName << endl;
			continue;
		}
		float thValue = w->function(thName)->getVal();
		g->SetPoint(iGraph, iGraph, obsValue-thValue);
		g->SetPointError(iGraph, 0., obsError);
		iGraph++;
	}

	g->SetTitle(scanVar1);
	g->GetXaxis()->SetTitle("scan step");
	g->GetYaxis()->SetTitleSize(0.06);
	g->GetYaxis()->SetLabelSize(0.04);
	g->GetYaxis()->SetTitleOffset(1.5);
	g->GetYaxis()->SetTitle(scanVar1);
	Int_t ci = 927;
	TColor *col = new TColor(ci, 0, 0, 1, " ", 0.5);
	g->SetFillColor(ci);
	g->SetFillStyle(1001);
	g->Draw("a3");
	g->Draw("lxsame");
	c2->Update();

	savePlot(c2, "parEvolutionObsSanCheck_"+name+"_"+scanVar1);
}