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

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");

}

Double_t calc_dNdY(const char *data_file = "pt_RFE", Bool_t show = kFALSE)
{
  TGraphErrors *g = new TGraphErrors(data_file);
  if (g->IsZombie()) return;

  TF1 *flt = new TF1("flt", LevyTsallis, 0., 5., 4);
  flt->SetParameters(2.93483e-02, 2.80382e-01, 8.10224e+00, 1.01944e+00);
  flt->FixParameter(3, 1.019445);

  Double_t fitmin = 0.25, fitmax = 5.25; // !!!
  g->Fit(flt, "Q", "", fitmin, fitmax);
  g->Fit(flt, "Q", "", fitmin, fitmax);
  g->Fit(flt, "Q", "", fitmin, fitmax);

  Double_t first, last;
  // !!! bining sensibility !!!
  Int_t graph_range[2] = {2, 23}; // [2] = 0.05 - 0.01, [23] = 4.70 + 0.01

  // interval 1
  first = 0.0;   // !!!
  last  = g->GetX()[graph_range[0]] - g->GetErrorX(graph_range[0]);
  Double_t fi1 = flt->Integral(first, last);

  // interval 2
  first = last;
  last  = g->GetX()[graph_range[1]] + g->GetErrorX(graph_range[1]);
  Double_t fi2 = flt->Integral(first, last);
  Double_t gi2 = 0.0;
  for (Int_t ip = graph_range[0]; ip <= graph_range[1]; ip++)
    gi2 += g->GetY()[ip]*g->GetErrorX(ip)*2.0;

  // interval 3
  first = last;
  last  = 30.0;   // !!!
  Double_t fi3 = flt->Integral(first, last);

  Double_t result = fi1 + gi2 + fi3;

  if (!show) {
    delete g;
    delete flt;
    return result;
  }

  Printf("function: %f \t %f \t %f", fi1, fi2, fi3);
  Printf("graph: \t \t \t %f", gi2);
  Printf("result => %f", result);
  g->Draw("AP");
  gPad->SetLogy();
  return result;
}

double ll_matrix::matrix_calibration::get_min_mass( TGraphErrors & graph )
{
    double m_min_x = graph.GetX()[0] , m_min_y = graph.GetY()[0];
    for( int i = 0 ; i < graph.GetN() ; i++ )
    {
        if( graph.GetY()[i] < m_min_y )
        {
            if( graph.GetX()[i] < LOWEST_MASS || graph.GetX()[i] > HIGHEST_MASS ) 
                continue;
            m_min_x = graph.GetX()[i];
            m_min_y = graph.GetY()[i];
        }
    }
    return m_min_x;
}

void covertQuan(){
    TFile *fin = TFile::Open("pPb.root");
//    TFile *fout = new TFile("pPb_xq.root","Recreate");
    int mult[4] = {120,150,185,220};
    int n[4] = {2,4,6,8};
    for(int i=0;i<4;i++){
    for(int j=0;j<4;j++){
        TGraphErrors *g = (TGraphErrors*)fin->Get(Form("v2Pt_%d_%d",n[j],mult[i]));
        if(n[j]==4){
        cout<<mult[i]<<endl;
        for(int ibin=0;ibin<g->GetN();ibin++){
            cout<<g->GetX()[ibin]<<",";
        }
        cout<<endl;
        for(int ibin=0;ibin<g->GetN();ibin++){
            cout<<g->GetY()[ibin]<<",";
        }
        cout<<endl;
        for(int ibin=0;ibin<g->GetN();ibin++){
            cout<<g->GetEY()[ibin]<<",";
        }
        cout<<endl;
        }
       // fout->cd();
       // g->Write("AP");
       // g1->Write();
    }
    }
}

void points(TString filename) {

	
	TString cmssw;
	// 167
	
	cmssw = "$1.6.7$";
	
	TFile *f = TFile::Open(filename);

	std::vector< TString > taggers;
	taggers.push_back( "gTC2_udsg" );
	taggers.push_back( "gTC3_udsg" );
	taggers.push_back( "gTP_udsg" );
	taggers.push_back( "gJBP_udsg" );
 	taggers.push_back( "gSSV_udsg" );
 	taggers.push_back( "gCSV_udsg" );

	std::vector< TString > discriminators;
	discriminators.push_back( "discTC2_udsg" );
	discriminators.push_back( "discTC3_udsg" );
	discriminators.push_back( "discTP_udsg" );
	discriminators.push_back( "discJBP_udsg" );
 	discriminators.push_back( "discSSV_udsg" );
 	discriminators.push_back( "discCSV_udsg" );

	//TCanvas *cv_TC = new TCanvas("cv_TC","cv_TC",700,700);
	//TCanvas *cv_TP = new TCanvas("cv_TP","cv_TP",700,700);

	std::cout << "Tagger & Point & Discriminator & light mistag & b-efficiency \\\\ \\hline" << std::endl;
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {

		TString tag = taggers[itagger];
		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("Histograms/MCtruth/"+tag);

		//if (taggers == "gTC2_udsg" || taggers =
		
		TGraph *dgraph = (TGraph*) gDirectory->Get("Histograms/MCtruth/"+discriminators[itagger]);
		dgraph->Sort();

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		
		//cv[itagger] = new TCanvas("cv","cv",600,600);
		//g->Draw("ACP");
		TString se = " & ";
		
		
		std::cout << tag << se << "Loose" << se << std::setprecision(3) << dgraph->Eval(0.1) << se << "0.1" << se << std::setprecision(2) << g->Eval(0.1) << "\\\\" << std::endl;
		std::cout << tag << se << "Medium" << se << std::setprecision(3) << dgraph->Eval(0.01) << se << "0.01" << se << std::setprecision(2) << g->Eval(0.01) << "\\\\" << std::endl;
		std::cout << tag << se << "Tight" << se << std::setprecision(3) << dgraph->Eval(0.001) << se << "0.001" << se << std::setprecision(2) << g->Eval(0.001) << "\\\\ \\hline" << std::endl;
		
	}

	
}

printgraph(const char* file, const char* graphname, const char* graphname_syst, bool reverseta=false, double scale=1./34.62)
{
   TFile * f = new TFile(file);
   TGraphErrors *gr = (TGraphErrors*) f->Get(graphname);
   TGraphErrors *grsyst = (TGraphErrors*) f->Get(graphname_syst);

   int nbins = gr->GetN();
   double *x = new double[nbins];
   if (reverseta)
      for (int i=0; i<nbins; i++)
         x[i] = gr->GetX()[nbins-i-1];
   else 
      x = gr->GetX();

   double *y = gr->GetY();
   double *ey = gr->GetEY();
   double *ey_syst = grsyst->GetEY();

   for (int i=0; i<nbins; i++)
      cout << x[i] << ": " << y[i]*scale << "+-" << ey[i]*scale << "(stat.) +-" << ey_syst[i]*scale << "(syst.)" << endl;
}

TGraphErrors *ReadMWGraph(const char *name, Int_t flag)
{
  Double_t xreject = 0.49;

  TGraphErrors *g = new TGraphErrors(name);
  if (g->IsZombie()) return 0;
  while (g->GetX()[0] < xreject)
    g->RemovePoint(0);

  TGraphErrors *g2 = new TGraphErrors(name);
  if (g2->IsZombie()) return 0;
  while (g2->GetX()[0] < xreject)
    g2->RemovePoint(0);
  g2->SetMarkerStyle(4);
  g2->SetMarkerSize(1.00);
  g2->SetMarkerColor(kBlack);
  g2->SetLineColor(kBlack);

  TGraphErrors *gsys = new TGraphErrors(name, "%lg %lg %lg %*lg %lg");
  if (gsys->IsZombie()) return 0;
  while (gsys->GetX()[0] < xreject)
    gsys->RemovePoint(0);
  for (Int_t i = 0; i < gsys->GetN(); i++)
    gsys->SetPointError(i, gsys->GetErrorX(i)*0.75, gsys->GetErrorY(i));
  gsys->SetFillColor(kGray+2);
  gsys->SetLineColor(kGray+2);
  gsys->SetFillStyle(3000);

  if (flag == 1 || flag == 3) {
    TGraphErrors *gt = new TGraphErrors(Form("%s_%s", name, "trues"));
    if (gt->IsZombie()) return 0;
    while (gt->GetX()[0] < xreject)
      gt->RemovePoint(0);
    gt->SetMarkerStyle(20);
    gt->SetMarkerSize(0.75);
    gt->SetMarkerColor(kGreen+1);
    gt->SetLineColor(kGreen+1);

    TGraphErrors *gbw = new TGraphErrors(Form("%s_%s", name, "gen"));
    if (gbw->IsZombie()) return 0;
    while (gbw->GetX()[0] < xreject)
      gbw->RemovePoint(0);
    gbw->SetMarkerStyle(20);
    gbw->SetMarkerSize(0.75);
    gbw->SetMarkerColor(kBlue+1);
    gbw->SetLineColor(kBlue+1);

    for (Int_t i = 0; i < g->GetN(); i++) {
      g->SetPointError(i, g->GetEX()[i], 0.);
      gt->SetPointError(i, gt->GetEX()[i], 0.);
      gbw->SetPointError(i, gbw->GetEX()[i], 0.);
    }

    for (Int_t i = 0; i < g2->GetN(); i++) {
      g2->SetPoint(i, g2->GetX()[i], g2->GetY()[i] - gt->GetY()[i] + gbw->GetY()[i]);
      g2->SetPointError(i, g2->GetEX()[i], TMath::Sqrt(g2->GetEY()[i]*g2->GetEY()[i] + gt->GetEY()[i]*gt->GetEY()[i] +
                                                       gbw->GetEY()[i]*gbw->GetEY()[i]));

      //      g2->SetPoint(i, g2->GetX()[i], g2->GetY()[i] - gt->GetY()[i] + 1.01947);
      //      g2->SetPointError(i, g2->GetEX()[i], TMath::Sqrt(g2->GetEY()[i]*g2->GetEY()[i] + gt->GetEY()[i]*gt->GetEY()[i] +
      //                                                       7.78680e-06*7.78680e-06));

      gsys->SetPoint(i, gsys->GetX()[i], g2->GetY()[i]);
    }
  }

  g->SetTitle();
  g->SetName(name);
  g->GetXaxis()->SetTitle("p_{T}, GeV/c");
  g->SetMarkerStyle(20);
  g->SetMarkerSize(0.95);
  g->SetMarkerColor(kRed+1);
  g->SetLineColor(kRed+1);

  const Double_t mass = 1.019455;
  const Double_t mass_delta = 0.000020;
  const Double_t width = 0.00426;
  const Double_t width_delta = 0.00004;
  if (flag == 1) { // mass
    g->GetYaxis()->SetTitleOffset(1.50);
    g->GetYaxis()->SetTitle("mass, GeV/c^{2}");
    g->SetMaximum(mass+0.0015);
    g->SetMinimum(mass-0.0015);
    TBox *box = new TBox(g->GetXaxis()->GetXmin(), mass - mass_delta, g->GetXaxis()->GetXmax(), mass + mass_delta);
    box->SetFillColor(kGray+1);
    box->SetFillStyle(3001);
    g->GetListOfFunctions()->Add(box);
    g->GetListOfFunctions()->Add(g2, "CP");
    g->GetListOfFunctions()->Add(gt, "CP");
    g->GetListOfFunctions()->Add(gbw, "CP");
  }
  else if (flag == 3) { // mass simple
    g2->SetTitle();
    g2->SetName(Form("%s_only", name));
    g2->GetXaxis()->SetTitle("p_{T}, GeV/c");
    g2->SetMarkerStyle(20);
    g2->SetMarkerSize(0.75);
    g2->SetMarkerColor(kBlack);
    g2->SetLineColor(kBlack);
    g2->GetYaxis()->SetTitleOffset(1.50);
//.........这里部分代码省略.........

Int_t VisualizeSurface()
{

    std::string gfname;
    std::cout << "\nFor the graph ";
    ListAllFilesInDirOfType (".",".root");
    std::cout << "\nEnter filename : ";
    std::getline(std::cin, gfname);

    TGraphErrors* mygr = GetGraph(gfname);
    mygr->SetMarkerColor(kRed);
    mygr->SetLineColor(kRed);
    mygr->SetMarkerStyle(20);

    std::string parfname = GetParamFile();
    std::ifstream pfile(parfname.data());
    Parameters params(pfile);
    if (! params.KeysAreSensible()) return -1;
    AngDistC W(params, 1, ELECTRIC);

    std::cout << params << std::endl;

    TGraph* gr = new TGraph(mygr->GetN());
    gr->SetMarkerColor(kBlue);
    gr->SetLineColor(kBlue);
    gr->SetMarkerStyle(20);

    Double_t num, denom;
    Double_t ratio;
    Double_t x[2];
    for (UInt_t i=0; i<mygr->GetN(); i++)
    {
        x[0] = mygr->GetX()[i];
        x[1] = 0;
        num = W(x);

        x[1] = TMath::Pi()/2.;
        denom = W(x);
        ratio = num/denom;
        gr->SetPoint(i, mygr->GetX()[i], TMath::Abs(ratio));
    }

    std::cout << "\n Phase2ChiSqC setup completed " << std::endl;

    TPaveText* pt = new TPaveText(0.5, 0.75, 0.8, 0.95);
    std::ostringstream os("", std::ios::out|std::ios::app);
    os << "a = " << params[0]->GetValue();
    pt->AddText(os.str().data());
    os.str("b = ");
    os << params[1]->GetValue();
    pt->AddText(os.str().data());
    os.str("c = ");
    os << params[2]->GetValue();
    pt->AddText(os.str().data());
    os.str("d = ");
    os << params[3]->GetValue();
    pt->AddText(os.str().data());

    TH1* h = new TH2D("Ratio Comparison", "", 10, -1,1, 500, 0,100);
    h->SetStats(0);

    TCanvas* c = new TCanvas("c");
    h->Draw();
    gr->Draw("LP");
    mygr->Draw("LP");

    gPad->Modified();
    c->Update();


    TLegend *leg = new TLegend(0.5, 0.75, 0.8, 0.95);
    leg->AddEntry(gr, "Fit", "PL");
    leg->AddEntry(mygr, "Data", "PL");
    leg->Draw();

    return 0;
}

void AnalyzeClipping(TString inputWaveName = "sum trigger input ch5 960mV",
TString outputWaveName = "sum trigger output ch5 - 2V clip - 960mV input", 
Double_t inputDelay = 1.1E-8, Double_t lowerCut = 16E-9, Double_t upperCut = 23E-9,
const char *inFile = "Data.root",const char *WaveformsFile = "Waveforms.root") {
	
	//try to access data file and in case of failure return
	if(gSystem->AccessPathName(inFile,kFileExists)) {
		cout << "Error: file " << inFile << " does not exsist. Run .x DataParse.C to create it" << endl;
		return;
	}

	TFile *f = TFile::Open(inFile);
	TFolder *dataSet;
	
	TString dataFolderS = "SumTriggerBoardData";
	dataFolderS.Append(";1");

	dataSet = (TFolder*)f->Get(dataFolderS);
	
	cout << dataSet << endl;

	cout << dataSet->GetName() << endl;
	
	Int_t nScope = 150; // number of measurements done by the scope evey time
	
	//try to access waveforms file and in case of failure return
	if(gSystem->AccessPathName(WaveformsFile,kFileExists)) {
		cout << "Error: file " << WaveformsFile << " does not exsist. Run .x WaveformsFileMaker.C to create it" << endl;
		return;
	}
	
	TFile *f = TFile::Open(WaveformsFile);
	TList *listOfKeys = f->GetListOfKeys();
	Int_t numberOfKeys = listOfKeys->GetEntries();
	TList *listOfGraphs = new TList();
	
	// if the waveform file name begins with the string "comparator" it goes in this list
	TList *listOfCompWaves = new TList();
	// if the waveform file name begins with the string "sum output" it goes in this list
	TList *listOfAdderWaves = new TList();

	for(Int_t i = 0; i < numberOfKeys; i++) {
		TString *keyName = new TString(listOfKeys->At(i)->GetName());
		TTree *tree = (TTree*)f->Get(keyName->Data());
		Float_t x = 0;
		Float_t y = 0;
		tree->SetBranchAddress("x",&x);
		tree->SetBranchAddress("y",&y);
		Int_t nentries = tree->GetEntries();

		TString *gName = new TString(keyName->Data());
		gName->Append(" graph");
		TGraphErrors *gWave = new TGraphErrors(nentries);
		gWave->SetName(gName->Data());
		gWave->SetTitle(gName->Data());
		gWave->GetXaxis()->SetTitle("Time");
		gWave->GetYaxis()->SetTitle("Voltage");

		for (Int_t j = 0; j < nentries; j++) {
			tree->GetEntry(j);
			gWave->SetPoint(j,x,y);
		}

		listOfGraphs->Add(gWave);
	}
	
	// Global variables
	
	Double_t *xInput, *xOutput, *yInput, *yOutput;
	
	// V input 960 mV
	
	TString path = "Clipping/Output width analysis/Channel 5/V input 960mV/";
	
	TGraphErrors *gClip960mV = TBGraphErrors(dataSet,path,"V clip","Output FWHM",1,nScope);
	
	gClip960mV->SetMarkerStyle(20);
	gClip960mV->SetMarkerSize(0.8);
	gClip960mV->GetXaxis()->SetTitle("V clipping (mV)");
	gClip960mV->GetYaxis()->SetTitle("Output FWHM (ns)");
	
	TCanvas *cClip960mV = new TCanvas("cClip960mV","Output FWHM in function of V clipping",800,600);
	gClip960mV->Draw("APEL");
	
	// Expected output FWHM
	
	TGraphErrors *gInput960mV = listOfGraphs->FindObject("sum trigger input ch5 960mV graph");
	
	Double_t *xClip = gClip960mV->GetX();
	Int_t nClip = gClip960mV->GetN();
	cout << "nClip = " << nClip << endl;
	Long64_t graphPoints = gInput960mV->GetN();
	yInput = gInput960mV->GetY();
	xInput = gInput960mV->GetX();
	
	vector<double> xFirst(nClip);
	vector<double> xLast(nClip);
	Double_t half;
	
	Int_t flag = 0;
//.........这里部分代码省略.........

void disceff(TString filename) {
	gROOT->SetStyle("Plain");
	
	TString cmssw;
	// 167
	
	cmssw = "$3.1.0_pre9$";
		
	TFile *f = TFile::Open(filename);

	std::vector< TString > taggers;
	taggers.push_back( "TC2" );
	taggers.push_back( "TC3" );
	taggers.push_back( "TP" );
	taggers.push_back( "BTP" );
	taggers.push_back( "SSV" );
	taggers.push_back( "CSV" );
	taggers.push_back( "MSV" );
	taggers.push_back( "SMT" );
	taggers.push_back( "SETbyIP3d" );
	taggers.push_back( "SETbyPt" );
	taggers.push_back( "SMTbyIP3d" );
	taggers.push_back( "SMTbyPt" );

	std::vector< TString > discriminators;
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
		discriminators.push_back( "disc"+taggers[itagger]+"_udsg" );
	}
//	discriminators.push_back( "discTC3_udsg" );
//	discriminators.push_back( "discTP_udsg" );

	const int dim=taggers.size();
        TCanvas *cv[dim];
	TMultiGraph* mg[dim];
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
		TString tag = "g"+taggers[itagger]+"_udsg";
		TString tagb = "g"+taggers[itagger]+"_b";
		TString tagc = "g"+taggers[itagger]+"_c";
		cv[itagger] = new TCanvas("cv_"+taggers[itagger],"cv_"+taggers[itagger],700,700);
		TLegend *legend0 = new TLegend(0.68,0.70,0.88,0.90);

		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);
		TGraphErrors *bgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagb);
		TGraphErrors *cgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagc);

		TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
		TGraph *bvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),bgraph->GetY());
		TGraph *cvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),cgraph->GetY());
		TGraph *lvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),agraph->GetY());
		TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
		dgraph->Sort();
//		udsgvsdgraph->Sort();
//		udsgvsdgraph->SetLineColor(1);
//              legend0 -> AddEntry(udsgvsdgraph,"control","l");
		lvsdgraph->Sort();
		lvsdgraph->SetLineColor(2);
                legend0 -> AddEntry(lvsdgraph,tag,"l");
		cvsdgraph->Sort();
		cvsdgraph->SetLineColor(3);
                legend0 -> AddEntry(cvsdgraph,tagc,"l");
		bvsdgraph->Sort();
		bvsdgraph->SetLineColor(4);
                legend0 -> AddEntry(bvsdgraph,tagb,"l");
		mg[itagger]= new TMultiGraph();
//		mg[itagger]->Add(udsgvsdgraph);
		mg[itagger]->Add(lvsdgraph);
		mg[itagger]->Add(cvsdgraph);
		mg[itagger]->Add(bvsdgraph);
//		mg[itagger]->Add(dgraph);

		cv[itagger]->cd(1);
		mg[itagger]->Draw("ALP");
		mg[itagger]->GetYaxis()->SetTitle("eff");
		mg[itagger]->GetXaxis()->SetTitle("discriminant");
		legend0 -> Draw();
		cv[itagger]->Update();
	        cv[itagger]->cd(0);
		cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".eps");
		cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".ps");

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		g->SetLineColor(itagger+1);

		std::cout << " Tagger: " << tag << std::endl;
		std::cout << " Loose(10%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.1) << " b-eff = " << g->Eval(0.1) << std::endl;
		std::cout << " Medium(1%):  " << " cut > " << std::setprecision(4) << dgraph->Eval(0.01) << " b-eff = " << g->Eval(0.01) << std::endl;
		std::cout << " Tight(0.1%) = " << " cut > " << std::setprecision(4) << dgraph->Eval(0.001) << " b-eff = " << g->Eval(0.001) << std::endl;

	}//end for
 
	
}

Int_t VisualizeSurface()
{

    std::string gfname;
    std::cout << "\nFor the graph ";
    ListAllFilesInDirOfType (".",".root");
    std::cout << "\nEnter filename : ";
    std::getline(std::cin, gfname);

    TGraphErrors* mygr = GetGraph(gfname);
    mygr->SetMarkerColor(kRed);
    mygr->SetLineColor(kRed);
    mygr->SetMarkerStyle(20);

    std::cout << "\nBegin setup of eta function" << std::endl;
    std::string rfname = GetROOTFile();
    TFile *fle = new TFile(rfname.data());
    if (!fle->IsOpen()) return -1;

    Eta2C eta(fle);
    AngMarginalEta2C marg_eta(eta);
    std::cout << "\neta function set up complete" << std::endl;

    std::cout << "\nBegin setup of denominators" << std::endl;
    std::string parfname = GetParamFile();
    std::ifstream pfile(parfname.data());
    Parameters params(pfile);
    if (! params.KeysAreSensible()) return -1;
    AngDistC W(params, 0, ELECTRIC);

    std::cout << params << std::endl;
    std::cout << "\ndenominator set up complete" << std::endl;

    std::cout << "\nSetting up numerator " << std::endl;
    PolPartAngDistC polpart(1, ELECTRIC);
    std::cout << "\nNumerator setup complete" << std::endl;

    std::cout << "\nBegin setup of h function" << std::endl;

    rfname = GetROOTFile();
    TFile *file = new TFile(rfname.data());
    if (!file->IsOpen()) return -1;

    TDirectory* dir = file->GetDirectory("SolidAnglePDFGeneratorCOutput");
    if (dir==0) return -3;
    std::vector<TH2*> solang_hists = ROOTUtils::GetAllTH2InDirectory(dir);

    std::cout << "h function set up" << std::endl;

    std::vector<UInt_t> indices(32);
    for (UInt_t i=0; i<indices.size(); i++)
    {
        indices[i] = i;
    }

    Phase2ChiSqC FUNC(mygr, W, polpart, marg_eta, solang_hists , indices);
    std::cout << "\nang dist set up" << std::endl;

    Double_t par[] = {params.at(0)->GetValue(),
                       params.at(1)->GetValue(),
                       params.at(2)->GetValue(),
                       params.at(3)->GetValue()};
//    FUNC.SetAllIntegrandParameters(par);
    Double_t chisq_pdf = FUNC(par);

    std::vector<Double_t> ratios = FUNC.GetComputedRatios();

    TGraph* gr = new TGraph(ratios.size());
    gr->SetMarkerColor(kBlue);
    gr->SetLineColor(kBlue);
    gr->SetMarkerStyle(20);

    for (UInt_t i=0; i<ratios.size(); i++)
    {
        gr->SetPoint(i, mygr->GetX()[i], ratios[i]);
    }

    std::cout << "\n Phase2ChiSqC setup completed " << std::endl;

//    Int_t strip;
//    std::cout << "\nEnter desired strip : ";
//    std::cin >> strip;

//    std::shared_ptr<ROOT::Math::IBaseFunctionMultiDim> pi = FUNC.GetDenominator(strip);
//    ParamFunctorAdapter *pfa = new ParamFunctorAdapter(*pi);
//
//    std::cout << "\nSetting up the function" << std::endl;
//    ROOT::Math::ParamFunctor pf;
//    pf.SetFunction(pfa);
//    TF2* f = new TF2("f", pf, -1, 1, -1.0*CLHEP::pi,CLHEP::pi,0);
//    Double_t pars[] = {params[0]->GetValue(),
//                       params[1]->GetValue(),
//                       params[2]->GetValue(),
//                       params[3]->GetValue()};
//    //f->SetParameters(pars);
//    f->SetNpx(50);
//    f->SetNpy(50);

    TPaveText* pt = new TPaveText(0.5, 0.75, 0.8, 0.95);
    std::ostringstream os("", std::ios::out|std::ios::app);
//.........这里部分代码省略.........

void Drawmethodcomp(){ 
 const int nDil = 8;
 const double centDil[nDil+1] = {1.0,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0};
 double centDilbin[nDil];
for(int i=0;i<nDil;i++)
  centDilbin[i]=0.29+i+N-1-nDil;
 const double NpartDil[nDil] = {8.75,30.51,53.30,86.23,130.06,187.35,261.49,355.45};
 const double NpartDilerr[nDil] = {1.13,3.02,3.95,4.35,4.60,4.44,3.96,2.83};
 const double NcollDil[nDil] = {8.01,38.86,86.85,175.76,326.06,563.21,926.79,1484.49};
 const double NcollDilerr[nDil] = {1.41,6.41,12.48,21.13,34.27,52.66,81.37,120.0};
 TGraphErrors *graphNpartDil = new TGraphErrors(nDil,centDilbin,NpartDil,0,NpartDilerr);
 TGraphErrors *graphNcollDil = new TGraphErrors(nDil,centDilbin,NcollDil,0,NcollDilerr);
 TCanvas *c1 = new TCanvas("c1","c1",1,1,550,460);
 c1->SetLogy();
  c1->SetFillColor(10);
  c1->SetFrameFillColor(0);
  c1->SetFrameBorderSize(0);
  c1->SetFrameBorderMode(0);
  c1->SetLeftMargin(0.15);
  c1->SetBottomMargin(0.15);
  c1->SetTopMargin(0.02);
  c1->SetRightMargin(0.02);
  gStyle->SetOptStat(0);
  c1->SetTicks(-1);

	N=N-1;
 TString str="Ncoll";
 TH1D* hist = new TH1D("","",N,0,N);
 hist->GetXaxis()->SetNdivisions(502);
if(method==0)
 hist->SetXTitle("Centrality");
else
 hist->SetXTitle("HF #Sigma E_{T} |#eta|>3");
 hist->SetYTitle(Form("<%s> and systematic errors",str.Data()));
 hist->SetMinimum(1);
 hist->SetMaximum(3999.99);
 hist->GetXaxis()->CenterTitle(0);
 hist->GetYaxis()->CenterTitle(1);
 hist->GetYaxis()->SetTitleOffset(1.1);
 hist->GetXaxis()->SetTitleOffset(1.1);
 hist->GetXaxis()->SetTitleSize(0.056);
 hist->GetYaxis()->SetTitleSize(0.056);
 hist->GetXaxis()->SetLabelSize(0.05);
 hist->GetYaxis()->SetLabelSize(0.05);
hist->GetXaxis()->SetLabelOffset(99);
hist->Draw();
	
	TFile *f=TFile::Open(outG);
	TGraphErrors* graph = (TGraphErrors*)f->Get(Form("std/%s_graph",str.Data()));
//	TGraphErrors* Gri055_graph = (TGraphErrors*)f->Get(Form("Gri055/%s_graph",str.Data()));
//	TGraphErrors* Gri101_graph = (TGraphErrors*)f->Get(Form("Gri101/%s_graph",str.Data()));
        TVectorD *centbin = (TVectorD*)f->Get(Form("std/G0/centbin"));
        TVectorD *kpoint = (TVectorD*)f->Get(Form("std/G0/kpoint"));

	TFile *f1=TFile::Open(Form("../../OneComp/double_side/%s",outG.Data()));
	TGraphErrors* graph1 = (TGraphErrors*)f1->Get(Form("std/%s_graph",str.Data()));
//	TGraphErrors* Gri055_graph = (TGraphErrors*)f->Get(Form("Gri055/%s_graph",str.Data()));
//	TGraphErrors* Gri101_graph = (TGraphErrors*)f->Get(Form("Gri101/%s_graph",str.Data()));
        TVectorD *centbin1 = (TVectorD*)f1->Get(Form("std/G0/centbin"));
        TVectorD *kpoint1 = (TVectorD*)f1->Get(Form("std/G0/kpoint"));

graph->SetTitle("g1");
graph->SetMarkerStyle(20);
graph->SetMarkerColor(1);
graph->SetLineColor(1);
graph->SetLineWidth(2);
graph->SetMarkerSize(1.2);
graph->Draw("Psameez");

graph1->SetMarkerStyle(24);
graph1->SetMarkerColor(4);
graph1->SetLineColor(4);
graph1->SetLineWidth(2);
graph1->SetMarkerSize(1.2);
graph1->Draw("Psameez");
/*
Gri055_graph->SetTitle("g2");
Gri055_graph->SetMarkerStyle(33);
Gri055_graph->SetMarkerColor(2);
Gri055_graph->SetLineColor(2);
Gri055_graph->SetLineWidth(2);
Gri055_graph->SetMarkerSize(1.2);
Gri055_graph->Draw("Psameez");

Gri101_graph->SetTitle("g3");
Gri101_graph->SetMarkerStyle(34);
Gri101_graph->SetMarkerColor(4);
Gri101_graph->SetLineColor(4);
Gri101_graph->SetLineWidth(2);
Gri101_graph->SetMarkerSize(1.2);
Gri101_graph->Draw("Psameez");
*/
 graphNpartDil->SetMarkerSize(1.2);
 graphNpartDil->SetLineColor(2);
 graphNpartDil->SetLineWidth(2);
 graphNpartDil->SetMarkerStyle(33);
 graphNpartDil->SetMarkerColor(2);
 graphNcollDil->SetMarkerSize(1.2);
 graphNcollDil->SetLineColor(2);
 graphNcollDil->SetLineWidth(2);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TH1D *h3 = new TH1D("h3",";p_{T,ref} (GeV);FSR sensitivity: -dR/d#alpha [%]",
		      1270,30,1300);

  cout << "Reading in data" << endl << flush;
  // Read in plots vs pT (and alpha)
  map<string, map<string, map<string, map<int, TGraphErrors*> > > > gemap;
  map<string, map<string, map<string, map<int, TGraphErrors*> > > > gamap;
  for (int itype = 0; itype != ntypes; ++itype) {
    for (int  imethod = 0; imethod != nmethods; ++imethod) {
      for (int  isample = 0; isample != nsamples; ++isample) {

	for (int  ialpha = 0; ialpha != nalphas; ++ialpha) {

	  fin->cd();
	  assert(gDirectory->cd(types[itype]));
	  assert(gDirectory->cd(bin));
	  TDirectory *d = gDirectory;

	  const char *ct = types[itype];
	  const char *cm = methods[imethod];
	  const char *cs = samples[isample];
	  const int a = alphas[ialpha];
	  // Get graph made vs pT
	  string s = Form("%s/%s/%s_%s_a%d",types[itype],bin,cm,cs,a);
	  TGraphErrors *g = (TGraphErrors*)fin->Get(s.c_str());
	  if (!g) cout << "Missing " << s << endl << flush;
	  assert(g);

	  // Clean out empty points
	  // as well as trigger-biased ones for dijets
	  // as well as weird gamma+jet high pT point
	  for (int i = g->GetN()-1; i != -1; --i) {
	    if (g->GetY()[i]==0 || g->GetEY()[i]==0 ||
		(string(cs)=="dijet" && g->GetX()[i]<70.) ||
		(string(cs)=="gamjet" && g->GetX()[i]>600. && etamin!=0))
	      g->RemovePoint(i);
	  }

	  gemap[ct][cm][cs][a] = g;
	  
	  // Sort points into new graphs vs alpha
	  TH1D *hpt = (isample==0 ? hpt2 : hpt1);
	  TProfile *ppt = (isample==0 ? ppt2 : ppt1);
	  if (isample==3) { hpt = hpt4; ppt = ppt4; } // pas-v6
	  for (int i = 0; i != g->GetN(); ++i) {
	    
	    double pt = g->GetX()[i];
	    ppt->Fill(pt, pt);
	    int ipt = int(hpt->GetBinLowEdge(hpt->FindBin(pt))+0.5);
	    //int ipt = int(pt+0.5);
	    TGraphErrors *ga = gamap[ct][cm][cs][ipt];
	    if (!ga) {
	      ga = new TGraphErrors(0);
	      ga->SetMarkerStyle(g->GetMarkerStyle());
	      ga->SetMarkerColor(g->GetMarkerColor());
	      ga->SetLineColor(g->GetLineColor());
	      gamap[ct][cm][cs][ipt] = ga;
	    }
	    int n = ga->GetN();
	    ga->SetPoint(n, 0.01*a, g->GetY()[i]);
	    ga->SetPointError(n, 0, g->GetEY()[i]);
	  } // for i 

	} // for ialpha

      } // for isample

//.........这里部分代码省略.........
  Double_t FitP0[2];
  Double_t FitP1[2];
  Double_t FitChisq[2];
  Double_t CSIDigiDeltaT0[nCSI];//nCSIDigi
  Double_t CSIDigiDeltaT1[nCSI];//nCSIDigi
  Int_t    CosmicTrigUp;
  Int_t    CosmicTrigDn;
  Double_t Roh;
  Double_t Theta;

  trin->SetBranchAddress( "RunNumber"     , &RunNumber      );
  trin->SetBranchAddress( "EventNumber"   , &EventNumber    );
  trin->SetBranchAddress( "ScintiSignal"  , &ScintiSignal   );
  trin->SetBranchAddress( "ScintiHHTimne" , &ScintiHHTime   );
  trin->SetBranchAddress( "ScintiTime"    , &ScintiTime     );
  trin->SetBranchAddress( "nCSIDigi"      , &nCSIDigi       );
  trin->SetBranchAddress( "CSIDigiE"      , CSIDigiE        );
  trin->SetBranchAddress( "CSIDigiTime"   , CSIDigiTime     );
  trin->SetBranchAddress( "CSIDigiHHTime" , CSIDigiHHTime   );
  trin->SetBranchAddress( "CSIDigiID"     , CSIDigiID       );
  trin->SetBranchAddress( "CSIDigiSignal" , CSIDigiSignal   );
  trin->SetBranchAddress( "CSIDigiDeltaT0" , CSIDigiDeltaT0 );
  trin->SetBranchAddress( "CSIDigiDeltaT1" , CSIDigiDeltaT1 );
  trin->SetBranchAddress( "FitP0"         , FitP0           );
  trin->SetBranchAddress( "FitP1"         , FitP1           );
  trin->SetBranchAddress( "FitChisq"      , FitChisq        );
  trin->SetBranchAddress( "CosmicTrigUp"  , &CosmicTrigUp   );
  trin->SetBranchAddress( "CosmicTrigDn"  , &CosmicTrigDn   );
  trin->SetBranchAddress( "Roh"           , &Roh            );
  trin->SetBranchAddress( "Theta"         , &Theta          );

  TFile* tfout = new TFile("CosmicOut_hist3.root", "recreate");
  TH2D* hisDeltaChannel = new TH2D("hisDeltaChannel","hisDeltaChannel",2716,0,2716,100,-10,10);
  TH1D* hisDelta[2716];
  TGraphErrors* grDelta = new TGraphErrors();
  TGraphErrors* grRES   = new TGraphErrors(); 
  TCanvas *can = new TCanvas("can","",800,800);


  for( int i = 0; i< 2716; i++){
    hisDelta[i] = new TH1D(Form("hisDelta%d",i ),Form("hisDelta%d",i),100,-10,10);
  }

  for( int ievent = 0; ievent < trin->GetEntries(); ievent++){

    trin->GetEntry(ievent);
    
    for( int idigi = 0; idigi < nCSIDigi ; idigi++){
      hisDelta[ CSIDigiID[ idigi ] ]->Fill( CSIDigiDeltaT1[ idigi ] );
      hisDeltaChannel->Fill( CSIDigiID[ idigi  ] , CSIDigiDeltaT1[ idigi ] );
      //std::cout  << CSIDigiID[ idigi ] << std::endl;
    }
  }
  for( int i = 0; i< 2716; i++){
    //std::cout << hisDelta[i]->GetEntries() << std::endl;
    if( hisDelta[i]->GetEntries() > 10){
      int rst = hisDelta[i]->Fit("gaus","Q","",hisDelta[i]->GetBinCenter( hisDelta[i]->GetMaximumBin() ) - 3, hisDelta[i]->GetBinCenter( hisDelta[i]->GetMaximumBin() )+3);
      TF1* func = NULL;
      func = hisDelta[i]->GetFunction("gaus");
      if( func != NULL ){
	grDelta->SetPoint( grDelta->GetN(), i, func->GetParameter(1));
	grDelta->SetPointError( grDelta->GetN()-1, 0, func->GetParError(2));
	grRES->SetPoint( grRES->GetN() , i , func->GetParameter(2));
      }
    }
    /*
    hisDelta[i]->Draw();
    can->Modified();
    can->Update();
    getchar();
    */
    hisDelta[ i ]  ->Write(); 
  }

  std::ofstream ofs("TimeResolutionCosmic3.dat");
  int    ID[2716];
  double Delta[2716];
  double Resolution[2716];
  for( int i = 0; i< 2716; i++){
    Resolution[i] = 0xFFFF;
    Delta[i]      = 0xFFFF;
  }
  for( int i = 0; i< grRES->GetN(); i++){
    Delta[(int)(grDelta->GetX()[i])] = grDelta->GetY()[i];
    Resolution[(int)(grRES->GetX()[i])] = grRES->GetY()[i] ;
  }
  for( int i = 0; i< 2716; i++){
    ofs << i             << "\t" 
	<< Delta[i]      << "\t"
	<< Resolution[i] << "\n";
  }
  
  grDelta->SetNameTitle("grDelta","grDelta");
  grRES->SetNameTitle("grRES","grRES");
  grDelta->Write();
  grRES->Write();
  hisDeltaChannel->Write();
  tfout->Close();
  ofs.close();
}

// Determine sensitivity to tracker dynamic inefficiency
// by studying ratio of jet responses in Runs G and F (and BCD / F, E / F)
void drawAvsB() {

  setTDRStyle();

  string epocha = "BCD";//"BCD";//"H";//"F";//"BCD";//"F";//"E";//"BCD";//"F";
  string epochb = "GH";//"G";//"BCD";//"G";//"E";//"E";//"F";//"G";

  // Add the rest as well
  string epocha2 = "";//"EF";
  string epochb2 = "";//"G";

  string type = "data";

  vector<string> methods;
  methods.push_back("mpfchs1");
  methods.push_back("ptchs");
  bool nozjptb = false;
  bool nogjmpf = false;
  bool nogjptb = true;
  bool mjvsjes = false;
  
  vector<string> samples;
  samples.push_back("zeejet");
  samples.push_back("zmmjet");
  samples.push_back("gamjet");
  //samples.push_back("multijet");

  cout << "draw"<<epocha<<"vs"<<epochb<<endl;
  const char *ct = type.c_str();
  const char *pa = epocha.c_str();
  const char *pb = epochb.c_str();

  const char *pa2 = epocha2.c_str();
  const char *pb2 = epochb2.c_str();

  TFile *fg = new TFile(Form("rootfiles/jecdata%s.root",pb),"READ");
  assert(fg && !fg->IsZombie());

  TFile *ff = new TFile(Form("rootfiles/jecdata%s.root",pa),"READ");
  assert(ff && !ff->IsZombie());

  TFile *fg2(0), *ff2(0);
  if (epochb2!="") fg2 = new TFile(Form("rootfiles/jecdata%s.root",pb2),"READ");
  if (epocha2!="") ff2 = new TFile(Form("rootfiles/jecdata%s.root",pa2),"READ");

  TH1D *h = new TH1D("h",
		     Form(";p_{T,ref} (GeV);%s ratio (%s / %s)",
			  (type=="ratio" ? "Data/MC" :
			   type=="data" ? "Data/data" : "MC/MC"),
			  (epocha + (epocha2!="" ? "+"+epocha2 : "")).c_str(),
			  (epochb + (epochb2!="" ? "+"+epochb2 : "")).c_str()),
		     3470,30,3500);
  h->SetMinimum(0.90);
  h->SetMaximum(1.15);
  h->GetXaxis()->SetMoreLogLabels();
  h->GetXaxis()->SetNoExponent();

  if (epocha=="F" && epochb=="G")
    lumi_13TeV = "Run2016F+G, 3.1+7.1 fb^{-1}";
  if (epocha=="BCD" && epochb=="G")
    lumi_13TeV = "Run2016BCD+H, 12.9+8.8 fb^{-1}";
  if (epocha=="BCD" && epochb=="G")
    lumi_13TeV = "Run2016BCD+FearlyGH, 12.9+16.8 fb^{-1}";
  if (epocha=="BCD" && epochb=="F")
    lumi_13TeV = "Run2016BCD+F, 13+3.1 fb^{-1}";
  if (epocha=="BCD" && epochb=="E")
    lumi_13TeV = "Run2016BCD+E, 13+4.0 fb^{-1}";
  if (epocha=="E" && epochb=="F")
    lumi_13TeV = "Run2016E+F, 4.0+3.1 fb^{-1}";
  if (epocha=="F" && epochb=="E")
    lumi_13TeV = "Run2016E+F, 4.0+3.1 fb^{-1}";

  if ((epocha=="BCDEF" && epochb=="GH") ||
      (epocha=="BCD" && epocha2=="EF" && epochb=="H" && epochb2=="G")) 
    lumi_13TeV = "Run2016BCDEF+GH, 19.7+16.8 fb^{-1}";
  if (epocha=="EF" && epochb=="BCD")
    lumi_13TeV = "Run2016BCD+EF, 12.9+6.8 fb^{-1}";
  if (epocha=="H" && epochb=="G")
    lumi_13TeV = "Run2016G+H, 8.0+8.8 fb^{-1}";

  if ((epocha=="BCD" && epocha2=="EF" && epochb=="G" && epochb2=="H")) 
    lumi_13TeV = "Run2016BCDFearly+FlateGH, 19.7+16.8 fb^{-1}";

  if ((epocha=="BCD" && epocha2=="" && ((epochb=="GH" && epochb2=="") ||
					(epochb=="G" && epochb2=="H"))))
    lumi_13TeV = "Run2016BCD+FlateGH, 12.9+16.8 fb^{-1}";
  if ((epocha=="EF" && epocha2=="" && ((epochb=="GH" && epochb2=="") ||
				       (epochb=="G" && epochb2=="H"))))
    lumi_13TeV = "Run2016EF+FlateGH, 6.8+16.8 fb^{-1}";

  if ((epocha=="EF" && epocha2=="" && epochb=="G" && epochb2=="H")) 
    lumi_13TeV = "Run2016EFearly+FlateGH, 6.8+16.8 fb^{-1}";


  TCanvas *c1 = tdrCanvas("c1",h,4,11,true);
  c1->SetLogx();

  TLatex *tex = new TLatex();
//.........这里部分代码省略.........