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

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;
}

void pointsed(TString filename = "btagpatanalyzerpy.root", TString check="") {
	TString cmssw;
	// 167
	
	cmssw = "$3.1.0_pre9$";
	
	TFile *f = TFile::Open(filename);
/////////////////////////////////////////////////////////////////////////////////////////
	std::vector< TString > checks;
	checks.push_back("");
	checks.push_back("Corr30_");
	checks.push_back("Corr30t0_");
	checks.push_back("Corr30t1_");
	checks.push_back("Corr30t2_");
	checks.push_back("Corr30t0_nConstituent_");
	checks.push_back("Corr30t1_nConstituent_");
	checks.push_back("Corr30t2_nConstituent_");
	checks.push_back("Uncor10_");
	checks.push_back("Uncor10t0_");
	checks.push_back("Uncor10t1_");
	checks.push_back("Uncor10t2_");
	checks.push_back("Uncor10t0_nConstituent_");
	checks.push_back("Uncor10t1_nConstituent_");
	checks.push_back("Uncor10t2_nConstituent_");
/////////////////////////////////////////////////////////////////////////////////////////
	std::vector< TString > taggers;
	std::vector< TString > labeltag;
        std::vector< int> colorlines;
	taggers.push_back( "TC2" );colorlines.push_back( 1 );labeltag.push_back( "TCHE");
	taggers.push_back( "TC3" );colorlines.push_back( 2 );labeltag.push_back( "TCHP");
	taggers.push_back( "TP"  );colorlines.push_back( 3 );labeltag.push_back( "JP");
	taggers.push_back( "SSV" );colorlines.push_back( 4 );labeltag.push_back( "SSV");
	taggers.push_back( "CSV" );colorlines.push_back( 5 );labeltag.push_back( "CSV");
	taggers.push_back( "MSV" );colorlines.push_back( 6 );labeltag.push_back( "CSVMVAB");
	taggers.push_back( "SMT" );colorlines.push_back( 8 );labeltag.push_back( "SMT");
	taggers.push_back( "BTP" );colorlines.push_back( 9 );labeltag.push_back( "JBP");
	taggers.push_back( "SMTbyIP3d" );colorlines.push_back( 11 );labeltag.push_back( "SMTByIP3d");
	taggers.push_back( "SMTbyPt" );colorlines.push_back( 12 );labeltag.push_back( "SMTByPt");
	taggers.push_back( "SETbyIP3d" );colorlines.push_back( 13 );labeltag.push_back( "SETByIP3d");
	taggers.push_back( "SETbyPt" );colorlines.push_back( 14 );labeltag.push_back( "SETByPt");
//	taggers.push_back( "IPM" );colorlines.push_back( 11 );labeltag.push_back( "IPMVAB");
//	taggers.push_back( "SMNIPT" );colorlines.push_back( 12 );labeltag.push_back( "SMTNoIP");

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
//ed	for ( size_t ichks = 0; ichks < checks.size(); ++ichks ) {
//ed		TString check = checks[ichks];
	
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
        TString PerfTitle= check;
        if (check == "") PerfTitle="uncorr pt>30";
	if (check == "Corr30_") PerfTitle="pt>30";
	if (check == "Corr30t0_") PerfTitle="pt>30, #tracks > 0";
	if (check == "Corr30t1_") PerfTitle="pt>30, #tracks > 1";
	if (check == "Corr30t2_") PerfTitle="pt>30, #tracks > 2";
	if (check == "Corr30t0_nConstituent_") PerfTitle="pt>30, #tracks > 0 & nConstituents >1";
	if (check == "Corr30t1_nConstituent_") PerfTitle="pt>30, #tracks > 1 & nConstituents >1";
	if (check == "Corr30t2_nConstituent_") PerfTitle="pt>30, #tracks > 2 & nConstituents >1";
	if (check == "Uncor10_") PerfTitle="uncorr pt>10";
	if (check == "Uncor10t0_") PerfTitle="uncorr pt>10 #tracks > 0";
	if (check == "Uncor10t1_") PerfTitle="uncorr pt>10 #tracks > 1";
	if (check == "Uncor10t2_") PerfTitle="uncorr pt>10 #tracks > 2";
	if (check == "Uncor10t0_nConstituent_") PerfTitle="uncorr pt>10 #tracks > 0 & nConstituents >1";
	if (check == "Uncor10t1_nConstituent_") PerfTitle="uncorr pt>10 #tracks > 1 & nConstituents >1";
	if (check == "Uncor10t2_nConstituent_") PerfTitle="uncorr pt>10 #tracks > 2 & nConstituents >1";

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

	TCanvas *cv = new TCanvas("cv","cv",900,900);
	TMultiGraph *mg =new TMultiGraph();
	TLegend *legend0 = new TLegend(0.68,0.12,0.88,0.32);
        // inverted axis
	TMultiGraph *mginv =new TMultiGraph();
	TLegend *legend1 = new TLegend(0.65,0.18,0.85,0.48);
	std::ofstream salida("BTagPATop3"+check+".txt");
	
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {

		TString tag    = check+"g"+taggers[itagger]+"_udsg";
		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);

		
		TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
//		TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
		dgraph->Sort();
//		udsgvsdgraph->Sort();

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		g->SetLineColor(itagger+1);
                legend0 -> AddEntry(g,taggers[itagger],"l");
		mg->Add(g);

                //inverted axis
		TGraphErrors *ginv = new TGraphErrors(agraph->GetN(),agraph->GetX(),agraph->GetY(),agraph->GetEX(),agraph->GetEY());
		ginv->Sort();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
  double sol = 299.97;//[mm/ns]
  TLine* solLine = new TLine( -1, -1./sol, 25, -1./sol);
  //can->cd( 6 );


  TCanvas* can1 = new TCanvas("can1","",800,400);
  can1->Divide( 2,1);
  can1->cd(1);
  grSlope->GetXaxis()->SetRangeUser( -1, 25);  
  grSlope->Draw("AP");
  solLine->SetLineStyle(2);
  solLine->Draw("same");
  TCanvas* can2 = new TCanvas("can2","",1200,800);
  can2->Divide(3,2);
  for( int i = 0; i< 5; i++){
    can2->cd( i+1 );
    CosmicP1TrigLarge[i][0]->Draw();
    for( int j = 0; j< 5; j++){
      CosmicP1TrigLarge[i][j]->Draw("same");
    }
  }
  can2->cd(6);
  grSlopeLarge->Draw("AP");

  
  Double_t H_trig = 2800;
  Double_t R      = 1200;
  TGraph* grDelta1 = new TGraph();
  TGraph* grDelta2 = new TGraph();

  TGraphErrors * grVCluster = new TGraphErrors();

  for( int i = 0; i< 5; i++){
    for( int j  =0; j< 5; j++){
      Double_t iv_0 = TMath::Abs(grSlope->GetY()[ i*5 +j ]);
      Double_t ive_0 = grSlope->GetEY()[ i*5 + j];
      

      Double_t deltaZ = (i-j)*100*R/H_trig;//mm            

      Double_t EdeltaZ = TMath::Sqrt(2)*50*R/H_trig;

      Double_t l  = TMath::Sqrt( R*R + deltaZ*deltaZ);
      Double_t El = deltaZ*EdeltaZ/TMath::Sqrt( R*R + deltaZ*deltaZ); 
      Double_t delta_t0 = l/sol;
      Double_t delta_t1 = iv_0*R - delta_t0;
      Double_t Edelta_t1 = QuadSum(ive_0*R, El/sol);
      
      Double_t V_CsI = deltaZ / delta_t1;
      Double_t EV_CsI = QuadSum( EdeltaZ / delta_t1, Edelta_t1/delta_t1);

      Double_t delta_t1_adj = deltaZ/70.;
      Double_t delta_t1_adj1 = deltaZ/80.;
      Double_t delta_t1_adj2 = deltaZ/90.;
      
      grDelta->SetPoint(grDelta->GetN(), (delta_t0+delta_t1_adj)/R , iv_0);
      grDelta1->SetPoint(grDelta1->GetN(), (delta_t0+delta_t1_adj1)/R , iv_0);
      grDelta2->SetPoint(grDelta2->GetN(), (delta_t0+delta_t1_adj2)/R , iv_0);
      
      if( TMath::Abs(deltaZ) <= 2*EdeltaZ ){ continue; }
      grVCluster->SetPoint( grVCluster->GetN(), i*5+j+1 , V_CsI);
      grVCluster->SetPointError( grVCluster->GetN()-1 , 0, EV_CsI);
      
    }
  }
  can1->cd(2);
  grVCluster->SetMarkerStyle( 6);
  grVCluster->Draw("AP");
  
  grVCluster->Fit("pol0","","",0,26);
  /*
  grDelta1->SetMarkerColor(2);
  grDelta2->SetMarkerColor(3);
  grDelta1->SetMarkerStyle(6);
  grDelta2->SetMarkerStyle(6);


  can1->cd(2);
  grDelta->Draw("AP");
  grDelta1->Draw("P");
  grDelta2->Draw("P");
  
  TF1* func = new TF1("func","x",0,0.1);
  func->Draw("same");

  */


  /*
  for( int i = 0; i< 5; i++){
    CosmicP1UpTrig[i]->Write();
    for( int j  =0 ;j < 5; j++){
      CosmicP1Trig[i][j]->Write();
    }
  }

  tfHist->Close();
  */
}

//.........这里部分代码省略.........
    if ( line.find("#") != string::npos )
      {
        cout << "Skip line " << line << endl;
        continue;
      }

    if ( line == "" )
      continue;

    //cout << "Processing file " << line << endl;


    TString infilename("data_calib/");
    infilename.Append(line);

    TFile *fin = new TFile( infilename );
    TTree *tin = (TTree*)fin->Get("t");

    ctemp->cd();
    tin->Draw("Bi:time");
    TGraph *gtime = new TGraph(tin->GetEntries(), &(tin->GetV2()[0]), &(tin->GetV1()[0]));
    gtime->SetLineColor(colorcounter);
    colorcounter++;

    TH1F* hBext = new TH1F("hBext","",100,0,1000);
    tin->Draw("Bo >> hBext");

    cres->cd();
    gtime->Draw("lsame");

    double Bext_i = hBext->GetMean();
    double BextErr_i = hBext->GetRMS();

    double Bint_i = gtime->Eval(590);
    double BintErr_i = 0;

    /* add legend entry */
    TString legname("B_ext ~ ");
    legname += (int)Bext_i;
    leg->AddEntry(gtime,legname,"l");

    cout << "B_ext: " << Bext_i << " \t B_int: " << Bint_i << endl;

    v_Bint.push_back(Bint_i);
    v_BintErr.push_back(BintErr_i);
    v_Bext.push_back(Bext_i);
    v_BextErr.push_back(BextErr_i);

  }

  cres->cd();
  leg->Draw();

  TGraphErrors *gfinal = new TGraphErrors(v_Bext.size(), &(v_Bext[0]), &(v_Bint[0]), &(v_BextErr[0]), &(v_BintErr[0]));
  gfinal->Sort();
  gfinal->SetName("Bint_Vs_Bext");
  gfinal->SetTitle("");
  gfinal->GetXaxis()->SetTitle("B_{ext} (mT)");
  gfinal->GetYaxis()->SetTitle("B_{int} (mT)");

  TCanvas *cfinal = new TCanvas();
  gfinal->Draw("APL");

  /* Save output graph */
  TString outfilename("output/");
  outfilename.Append(outfile);
  TFile *fout = new TFile(outfilename,"RECREATE");

  cres->Write();
  gfinal->Write();

  fout->Close();

  /* Write result to txt output file */
  TString outfilenametxt = outfilename;
  outfilenametxt.ReplaceAll(".root",".txt");

  ofstream foutxt;
  foutxt.open( outfilenametxt );
  foutxt <<  "# Bo sig_Bo Bi sig_Bi shield sig_shield sf sig_sf time_dependent" << endl;

  for ( int i = 0; i < gfinal->GetN(); i++ )
  {
    double Bo = gfinal->GetX()[i];
    double sig_Bo = gfinal->GetEX()[i];
    double Bi = gfinal->GetY()[i];
    double sig_Bi = gfinal->GetEY()[i];
    double shield = 0;
    double sig_shield = 0;
    double sf = 0;
    double sig_sf = 0;
    double time_dependent = 0;

    foutxt <<  Bo << " " << sig_Bo << " " << Bi << " " << sig_Bi << " "
         << shield << " " << sig_shield << " " << sf << " " << sig_sf
         << " " << time_dependent << endl;
  }

  return 0;
}

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);
//.........这里部分代码省略.........

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
 
	
}

//.........这里部分代码省略.........
    leg0->SetFillColor(10);
    leg0->SetBorderSize(0);
    leg0->SetTextFont(42);
    leg0->SetTextSize(0.047);
    leg0->AddEntry(graph,"From Ancestor fitting","p");
    leg0->AddEntry(graph1,"From One Comp fitting","p");
//    leg0->AddEntry(Gri055_graph,"Gribov #Omega=0.55","p");
//    leg0->AddEntry(Gri101_graph,"Gribov #Omega=1.01","p");
 if(str=="Npart")   leg0->AddEntry(graphNpartDil,"Npart from run I","p");
 if(str=="Ncoll")   leg0->AddEntry(graphNcollDil,"Ncoll from run I","p");
	leg0->Draw();	
c1->SaveAs(Form("%sGri.png",str.Data()));
c1->SaveAs(Form("%sGri.pdf",str.Data()));


 TCanvas *c2 = new TCanvas("c2","c2",1,1,550,460);
  c2->SetFillColor(10);
  c2->SetFrameFillColor(0);
  c2->SetFrameBorderSize(0);
  c2->SetFrameBorderMode(0);
  c2->SetLeftMargin(0.15);
  c2->SetBottomMargin(0.15);
  c2->SetTopMargin(0.02);
  c2->SetRightMargin(0.02);
  c2->SetTicks(-1);
c2->cd();
TGraphErrors* fdvf1 = (TGraphErrors*)graph->Clone("fdvf1");
TGraphErrors* f1dvDil = (TGraphErrors*)graph->Clone("f1dvDil");
TGraphErrors* fdvDil = (TGraphErrors*)graph->Clone("fdvDil");
TGraphErrors* Dilsys = (TGraphErrors*)graph->Clone("Dilsys");
TGraphErrors* Dilsys2 = (TGraphErrors*)graph->Clone("Dilsys2");
for(int ip = 0;ip<fdvf1->GetN();ip++){
    double x = graph->GetX()[ip];
    double ey = graph->GetEY()[ip];
    double y = graph->GetY()[ip];
    double ey1 = graph1->GetEY()[ip];
    double y1 = graph1->GetY()[ip];
    double eyNcoll = graphNcollDil->GetEY()[ip];
    double yNcoll = graphNcollDil->GetY()[ip];
    double eyNpart = graphNpartDil->GetEY()[ip];
    double yNpart = graphNpartDil->GetY()[ip];
    if(y1!=0 && y!=0){
        fdvf1->SetPoint(ip,x,y/y1);
        fdvf1->SetPointError(ip,0,y/y1*sqrt((ey/y)**2+(ey1/y1)**2));
        if(str=="Ncoll"){
        fdvDil->SetPoint(ip,x,y/yNcoll);
        double yerr = y/yNcoll*sqrt((ey/y)**2+(eyNcoll/yNcoll)**2);
        fdvDil->SetPointError(ip,0,0);
        f1dvDil->SetPoint(ip,x,y1/yNcoll);
        double y1err = y1/yNcoll*sqrt((ey1/y1)**2+(eyNcoll/yNcoll)**2);
        f1dvDil->SetPointError(ip,0,0);
        Dilsys->SetPoint(ip,x,1.);
        Dilsys->SetPointError(ip,0,TMath::Max(fabs(y1/yNcoll-1),fabs(y/yNcoll-1)));
        Dilsys2->SetPoint(ip,x,1.);
        Dilsys2->SetPointError(ip,0,TMath::Max(yerr,y1err));
        }
        else if(str=="Npart"){
        fdvDil->SetPoint(ip,x,y/yNpart);
        double yerr = y/yNpart*sqrt((ey/y)**2+(eyNcoll/yNpart)**2);
        fdvDil->SetPointError(ip,0,0);
        f1dvDil->SetPoint(ip,x,y1/yNpart);
        double y1err = y1/yNpart*sqrt((ey1/y1)**2+(eyNcoll/yNcoll)**2);
        f1dvDil->SetPointError(ip,0,0);
        Dilsys->SetPoint(ip,x,1.);
        Dilsys->SetPointError(ip,0,TMath::Max(fabs(y1/yNpart-1),fabs(y/yNpart-1)));
        Dilsys2->SetPoint(ip,x,1.);

int postprocessingSysError(){

  cout<<endl<<endl<<endl<<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Postproccess all systematic uncertainties! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<endl<<endl;
  gErrorIgnoreLevel = 1001;

  const TString method  = "RMS99";
  const TString type    = "PFCHS";

  const int nEta =4;
  double eta_bins[5] = {0., 0.5, 1.1, 1.7, 2.3};
  // For looking at different systematic uncertainties independently
  const bool QCD    = true;
  const bool JEC    = true;
  const bool flavor = true;
  const bool PU     = true;
  const bool MC     = true;
  

  TString etaString, filename;   

  TString rootFiles, AuxString;  
  TString JetType = "PFCHS";
  TString Method  = "RMS99";   

  double *ratioEtaBinnedX  = new double[nEta];
  double *ratioEtaBinnedY  = new double[nEta];
  double *ratioEtaBinnedEX = new double[nEta];
  double *ratioEtaBinnedEY = new double[nEta];

  double *ratioEtaBinnedQCDUpY   = new double[nEta];
  double *ratioEtaBinnedQCDDownY = new double[nEta];

  TF1 *QCDuncertainty;

  if(QCD){
    rootFiles   = (TString) "scripts/plotsQCD/FinalErrorsQCD_" + type + (TString) "_" + method + (TString) ".root";
    TFile *_file = TFile::Open(rootFiles);    
    _file->GetObject("function",QCDuncertainty);
  }
  
  for(int eta = 0; eta < nEta; eta++){
    
    //cout<< endl<<endl<<endl<<eta+1<<". eta Bin!!"<<endl;

    // Read the MC and data results 
    rootFiles = (TString) "root_files_FINAL_data/Resolution_for_" + (long) (eta+1) + (TString) "_eta_bin_" + JetType + (TString) "_data_" + Method + (TString) ".root";
    TGraphErrors* JERData = readTGraphErrors(rootFiles,"Graph;1","Graph");
    rootFiles = (TString) "root_files_FINAL_mc/Resolution_for_" + (long) (eta+1) + (TString) "_eta_bin_" + JetType + (TString) "_mc_" + Method + (TString) ".root";
    TGraphErrors* JERMC = readTGraphErrors(rootFiles,"Graph","Graph");
    
    if(eta+1 == 1) etaString = Form("JER for |#eta| < %4.1f",etaBins[eta+1]);
    else           etaString = Form("JER for %4.1f <|#eta|< %4.1f",etaBins[eta+1],etaBins[eta+2]);
 
    //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // 1.) Calculate the ratio w/o systematic Uncertainties  


    int nData    = JERData->GetN();

    double *dataX  = JERData->GetX();
    double *dataY  = JERData->GetY();
    double *dataEX = JERData->GetEX();
    double *dataEY = JERData->GetEY();
    
    double *mcX  = new double[nData];
    double *mcY  = new double[nData];
    double *mcEX = new double[nData];
    double *mcEY = new double[nData];
    
    double *ratioX  = new double[nData];
    double *ratioY  = new double[nData];
    double *ratioEX = new double[nData];
    double *ratioEY = new double[nData];


    //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // Initialize some stuff for QCD uncertainty
    double *ratioQCDUpY    = new double[nData];
    double *ratioQCDDownY  = new double[nData];    
    //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
    int idx = 0;
    for(int i=0; i<nData; i++){

      JERMC   -> GetPoint(idx,mcX[i],mcY[i]);
      mcEX[i] = JERMC -> GetErrorX(idx);
      mcEY[i] = JERMC -> GetErrorY(idx);

      idx += 1;

      if(TMath::Abs(dataX[i]/mcX[i] - 1.) > 0.1){
	i -= 1;
	continue;
      }
     
      ratioX[i]  = 1./2.*(dataX[i] + mcX[i]);
      ratioY[i]  = dataY[i]/mcY[i];
      ratioEX[i] = 1./2.*TMath::Sqrt(TMath::Power(dataEX[i],2)+TMath::Power(mcEX[i],2));
      ratioEY[i] = TMath::Sqrt(TMath::Power((1./mcY[i]),2)*TMath::Power(dataEY[i],2)+TMath::Power((dataY[i]/(TMath::Power(mcY[i],2))),2)*TMath::Power(mcEY[i],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

//.........这里部分代码省略.........
    assert(gf);
    if (ff2) {
      TGraphErrors *gf2 = (TGraphErrors*)ff2->Get(Form("%s/eta00-13/%s_%s_a30",ct,cm,cs));
      assert(gf2);
      gf = addGraph(gf,gf2);
    }
    
    if (!(gf->GetN()==gg->GetN())) {

      // Remove highest pT point is that is the offender (BCD vs GH)
      if (gg->GetN()>gf->GetN() &&
	  fabs(gg->GetX()[gg->GetN()-1]/gf->GetX()[gf->GetN()-1]-1)>0.1 &&
	  fabs(gg->GetX()[gg->GetN()-2]/gf->GetX()[gf->GetN()-1]-1)<0.1) {
	cout << "Remove point B(N-1)" << endl;
	gg->RemovePoint(gg->GetN()-1);
      }
      else {
	cout << "sample " << samples[is] << " method " << methods[im]
	     << " gf->N: " << gf->GetN() << " gg->N: " << gg->GetN() << endl;
	cout << " x_gf(N-1)=" << gf->GetX()[gf->GetN()-1]
	     << " x_gg(N-1)=" << gg->GetX()[gg->GetN()-1]
	     << " x_gg(N-2)=" << gg->GetX()[gg->GetN()-2] << endl;
      }

      assert(gf->GetN()==gg->GetN());
    }

    TGraphErrors *g = (TGraphErrors*)gg->Clone(Form("ge_%s_%s",cm,cs));
    for (int i = 0; i != g->GetN(); ++i) {
      double yg = gg->GetY()[i];
      double yf = gf->GetY()[i];
      g->SetPoint(i, gg->GetX()[i], yf / yg);
      double ex = gg->GetEX()[i];
      double eg = gg->GetEY()[i];
      double ef = gf->GetEY()[i];
      g->SetPointError(i, ex, yf/yg*sqrt(pow(eg/yg,2)+pow(ef/yf,2)));
    }
    //g->Draw(is==0 ? "AP" : "SAMEP");
    g->SetLineWidth(1+is);
    g->Draw("SAMEPZ");

    if (samples[is]=="gamjet" && methods[im]=="mpfchs1" && nogjmpf) {
      tex->SetTextColor(kBlue);
      tex->DrawLatex(0.20,0.63,"#gamma+jet MPF excl. from fit");
      tex->SetTextColor(kBlack);
    }
    else if (samples[is]=="gamjet" && methods[im]=="ptchs" && nogjptb) {
      tex->SetTextColor(kBlue);
      tex->DrawLatex(0.20,0.63,"#gamma+jet p_{T}^{bal} excl. from fit");
      tex->SetTextColor(kBlack);
    }
    else if ((samples[is]=="zmmjet" || samples[is]=="zeejet") &&
	     methods[im]=="ptchs" && nozjptb) {
      tex->SetTextColor(kRed);
      tex->DrawLatex(0.20,0.63,"Z+jet p_{T}^{bal} excl. from fit");
      tex->SetTextColor(kBlack);
    }
    else if (samples[is]=="multijet") {
      g->SetMarkerColor(kGray+1);
      g->SetLineColor(kGray+1);
      if (methods[im]=="ptchs") gmjb = g;
      if (methods[im]=="mpfchs1") gmpf = g;
    }
    else
      mg->Add(g);
  } // for is