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

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
 
	
}

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

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

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

//.........这里部分代码省略.........
  Double_t ymax    = nChannel + ymin + ylegend;
  
  TH2F* h2 = new TH2F("h2_" + suffix, "", 100, xmin, xmax, 100, ymin, ymax);

  h2->Draw();
  
  
  // NLO WZ cross-section
  //----------------------------------------------------------------------------
  TBox* nlo = new TBox(1. - xs_nlo_left [theCharge] / xs_nlo[theCharge], ymin,
		       1. + xs_nlo_right[theCharge] / xs_nlo[theCharge], ymax - ylegend);

  nlo->SetLineColor(0);
  nlo->SetFillColor(kGray);
  nlo->SetFillStyle(1001);

  nlo->Draw("e2,same");

  TLine* line = new TLine(1., ymin, 1., ymax - ylegend);

  line->SetLineColor(kGray+1);
  line->SetLineWidth(2);

  line->Draw("same");


  // Cross sections
  //----------------------------------------------------------------------------
  gLumi->Draw("p,same");
  gSyst->Draw("p,same");
  gStat->Draw("p,same");


  // Labels
  //----------------------------------------------------------------------------
  for (UInt_t i=0; i<nChannel; i++) {

    Double_t x = gStat->GetX()[i];
    Double_t y = gStat->GetY()[i];

    Double_t gStatError = gStat->GetErrorX(i);
    Double_t gSystError = gSyst->GetErrorX(i);
    Double_t gLumiError = gLumi->GetErrorX(i);

    DrawTLatex(42, xmin+0.06, y+0.15, 0.035, 12,
	       Form("%s %.2f #pm %.2f",
		    lChannel[i].Data(), x, gLumiError), 0);

    gLumiError = sqrt(gLumiError*gLumiError - gSystError*gSystError);
    gSystError = sqrt(gSystError*gSystError - gStatError*gStatError);

    DrawTLatex(42, xmin+0.06, y-0.15, 0.025, 12,
	       Form("%.2f #pm %.2f #pm %.2f #pm %.2f",
		    x, gStatError, gSystError, gLumiError), 0);
  }

  DrawTLatex(42, 0.050, 0.975, _bigLabelSize, 13, "CMS Preliminary");

  DrawTLatex(42, 0.940, 0.983, _bigLabelSize, 33, 
	     Form("#sqrt{s} = 7 TeV, L = %.1f fb^{-1}", luminosity/1e3));

  TString swz = "";

  if      (theCharge == WPlus)  swz = "W^{+}Z";
  else if (theCharge == WMinus) swz = "W^{-}Z";
  else                          swz = "W^{#pm}Z";
  
  h2->GetXaxis()->CenterTitle();
  h2->GetXaxis()->SetTitleOffset(1.4);
  h2->GetXaxis()->SetTitle(Form("#sigma_{%s}^{exp} / #sigma_{%s}^{theory}",
				   swz.Data(),
				   swz.Data()));

  h2->GetYaxis()->SetTitle("");


  // Remove y-axis labels
  //----------------------------------------------------------------------------
  TAxis* yaxis = h2->GetYaxis();
  
  for (Int_t j=1; j<yaxis->GetNbins(); j++) yaxis->SetBinLabel(j, "");


  // Additional legend
  //----------------------------------------------------------------------------
  DrawLegend(0.645, 0.840, gStat, " stat.",  "lp");
  DrawLegend(0.645, 0.795, nlo,   " theory", "f");
  DrawLegend(0.800, 0.840, gSyst, " syst.",  "l");
  DrawLegend(0.800, 0.795, gLumi, " lumi.",  "l");


  // Save
  //----------------------------------------------------------------------------
  canvas->Update();
  canvas->GetFrame()->DrawClone();
  canvas->RedrawAxis();

  canvas->SaveAs(Form("pdf/ratioNLO_%s.pdf", suffix.Data()));
  canvas->SaveAs(Form("png/ratioNLO_%s.png", suffix.Data()));
}

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

//.........这里部分代码省略.........
  ewktau->Add(ewkVV);
  
  //TH1* qcd = (TH1*)f->Get("QCDInv");
  TH1* qcd = (TH1*)f->Get("QCD");
  ci = TColor::GetColor("#ffcc33");
  qcd->SetFillColor(ci);
  qcd->SetLineWidth(0);
  TH1* fakett = (TH1*)f->Get("fake_tt");
  ci = TColor::GetColor("#669900");
  fakett->SetFillColor(ci);
  fakett->SetLineWidth(0);
  TH1* fakeW = (TH1*)f->Get("fake_W");
  ci = TColor::GetColor("#cc3300");
  fakeW->SetFillColor(ci);
  fakeW->SetLineWidth(0);
  TH1* faket = (TH1*)f->Get("fake_t");

  TH1F *hFrame = new TH1F("hFrame","",20,0,400);
  hFrame->SetMinimum(ymin);
  if (log)
    hFrame->SetMaximum(ymax*1.5);
  else
    hFrame->SetMaximum(ymax);
  hFrame->SetDirectory(0);
  hFrame->SetStats(0);
  hFrame->SetLineStyle(0);
  hFrame->SetMarkerStyle(20);
  hFrame->SetXTitle("Transverse mass (#tau jet, E_{T}^{miss}), (GeV/c^{2})");
  if (paperStatus)
    hFrame->SetXTitle("Transverse mass (#tau_{h}, E_{T}^{miss}), (GeV/c^{2})");
  hFrame->SetYTitle("Events / 20 GeV/c^{2}");
  hFrame->GetXaxis()->SetTitleSize(0);
  hFrame->GetXaxis()->SetLabelSize(0);
  hFrame->GetYaxis()->SetTitleFont(43);
  hFrame->GetYaxis()->SetTitleSize(27);
  hFrame->GetYaxis()->SetTitleOffset(1.3);
  

  // signal
  hh->Scale(br*br);
  hw->Scale(2*br*(1.0-br));
  TH1* signal = (TH1*)hh->Clone();
  signal->Add(hw);

  ci = TColor::GetColor("#ff3399");
  signal->SetLineColor(ci);
  signal->SetLineStyle(2);
  signal->SetLineWidth(2);

  // Fakes
  TH1* fakes = (TH1*)(fakett->Clone());
  fakes->Add(fakeW);
  fakes->Add(faket);

  // stacked backgrounds
  THStack *exp = new THStack();
  exp->SetName("exp");
  exp->SetTitle("exp");
  exp->Add(fakes);
  exp->Add(ewktau);
  exp->Add(qcd);
  exp->Add(signal);
  
  TH1* hExpBkg = (TH1*)fakes->Clone();
  hExpBkg->Add(ewktau);
  hExpBkg->Add(qcd);

//.........这里部分代码省略.........
      gg = addGraph(gg,gg2);
    }
    
    TGraphErrors *gf = (TGraphErrors*)ff->Get(Form("%s/eta00-13/%s_%s_a30",ct,cm,cs));
    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;