C++ TBox::SetLineColor方法代码示例

void Draw_CMS_Y1S_RaaVsRap(TLegend *lgd)
{
  const int nbinsRap=6;

  Double_t RapCMS[nbinsRap]={0.2,0.6,1.0,1.4,1.8,2.2};
  Double_t ErrRapCMS[nbinsRap]={0};
  
  Double_t RaaRapCMS[nbinsRap] = {0.402,0.377,0.452,0.461,0.466,0.35}; 
  Double_t RaaRapStatErrCMS[nbinsRap] = {0.025,0.025,0.030,0.034,0.039,0.053};
  Double_t RaaRapSystErrCMS[nbinsRap] = {0.0404,0.038,0.046,0.0466,0.0484,0.0373};
  
  TGraphErrors *grRaaRapCMS = new TGraphErrors(nbinsRap, RapCMS, RaaRapCMS, ErrRapCMS, RaaRapStatErrCMS);  
  grRaaRapCMS->SetMarkerStyle(20);
  grRaaRapCMS->SetMarkerColor(4);
  grRaaRapCMS->GetYaxis()->SetRangeUser(0,3.0);
  grRaaRapCMS->GetXaxis()->SetTitle("p_{T}(GeV/c)");
  grRaaRapCMS->GetYaxis()->SetTitle("R_{AA}");
   
  TAxis *Xaxis2 = grRaaRapCMS->GetXaxis();
  Xaxis2->SetLimits(0.,2.4);
  
  grRaaRapCMS->Draw("AP");
  
  TLine *lh4 = new TLine(0.0,1.0,20.0,1.0);
  lh4->SetLineColor(1);
  lh4->SetLineStyle(1);
  lh4->SetLineWidth(2);
  lh4->Draw("same");


  TLatex *tb= new TLatex;
  tb->SetNDC(); 
  tb->SetTextAlign(12);
  tb->SetTextColor(1);
  tb->SetTextSize(0.040);
  //tb->DrawLatex(0.55,0.22,"PbPb #sqrt{s_{NN}} = 2.76 TeV");
  //tb->DrawLatex(0.55,0.16,"#varUpsilon #rightarrow #mu^{+} #mu^{-}, p_{T}^{#varUpsilon} > 0.0 GeV/c");  
  
  TBox *RaaRapCMSSys[nbinsRap];
  for(int j=0;j<nbinsRap;j++){
    RaaRapCMSSys[j] = new TBox(RapCMS[j]-0.1,  RaaRapCMS[j]-RaaRapSystErrCMS[j], RapCMS[j]+0.1,  RaaRapCMS[j]+RaaRapSystErrCMS[j]);
  }
  
  for(int j=0;j<nbinsRap;j++){
    RaaRapCMSSys[j]->SetFillStyle(0000);
    RaaRapCMSSys[j]->SetLineColor(4);
    RaaRapCMSSys[j]->Draw("same"); 
  }
    
  TBox *CMSGlobalSysRap;
  CMSGlobalSysRap = new TBox(18-0.2, 1 - 0.083, 18+0.2, 1 + 0.083);
  CMSGlobalSysRap->SetFillStyle(3001);
  CMSGlobalSysRap->SetLineColor(4);
  CMSGlobalSysRap->SetFillColor(4);
  //CMSGlobalSysRap->Draw("same"); 
    
  lgd->AddEntry(grRaaRapCMS,"CMS Data", "P");  

  }

void Draw_ALICEMid_JPsi_RaaVsNpart(TLegend *lgd)
{
  //ALICE MID rapidity
  int nbinsALICEMid=3;
  Double_t NPartALICEMid[3]={357,193,46};
  Double_t ErrNPartALICEMid[3]={0};
  Double_t RaaALICEMid[3] = {0.82,0.65,0.73}; 
  Double_t SystErrALICEMid[3] = {0.15,0.10,0.16};
  
  TGraphErrors *grRaaALICEMid = new TGraphErrors(nbinsALICEMid, NPartALICEMid, RaaALICEMid, ErrNPartALICEMid, SystErrALICEMid);  
  grRaaALICEMid->SetMarkerStyle(20);
  grRaaALICEMid->SetMarkerColor(2);
  //grRaaALICEMid->SetLineColor(2);
  grRaaALICEMid->GetYaxis()->SetRangeUser(0,1.5);
  
  TAxis *Xaxis = grRaaALICEMid->GetXaxis();
  Xaxis->SetLimits(0.,430.0);
  grRaaALICEMid->GetXaxis()->SetTitle("N_{Part}");
  grRaaALICEMid->GetYaxis()->SetTitle("R_{AA}");
  
  grRaaALICEMid->Draw("Psame");
  
  //TLatex *tb= new TLatex;
  //tb->SetNDC(); 
  //tb->SetTextAlign(12);
  //tb->SetTextColor(1);
  //tb->SetTextSize(0.040);
  //tb->DrawLatex(0.55,0.90,"PbPb #sqrt{s_{NN}} = 2.76 TeV");
  //tb->DrawLatex(0.22,0.16,"J/#psi #rightarrow #mu^{+} #mu^{-}, p_{T}^{J/#psi} > 0.0 GeV/c");  
  //tb->DrawLatex(0.55,0.85,"p_{T}^{J/#psi} > 0.0 GeV/c");  
  
  TLine *lh3 = new TLine(0.0,1.0,420,1.0);
  lh3->SetLineColor(1);
  lh3->SetLineStyle(1);
  lh3->SetLineWidth(1.5);
  lh3->Draw("same");
    
  TBox *RaaJPsiALICEMidSys[12];
  for(int j=0;j<3;j++){
    RaaJPsiALICEMidSys[j] = new TBox(NPartALICEMid[j]-3,  RaaALICEMid[j]-SystErrALICEMid[j], NPartALICEMid[j]+3, RaaALICEMid[j]+SystErrALICEMid[j]);
  }
  
  for(int j=0;j<3;j++){
    RaaJPsiALICEMidSys[j]->SetFillStyle(0000);
    RaaJPsiALICEMidSys[j]->SetLineColor(2);
    RaaJPsiALICEMidSys[j]->Draw("same"); 
  }
  
  TBox *ALICEMidGlobalSysJPsi;
  ALICEMidGlobalSysJPsi = new TBox(385-5, 1 - 0.26, 385+5, 1 + 0.26);
  
  ALICEMidGlobalSysJPsi->SetFillStyle(3001);
  ALICEMidGlobalSysJPsi->SetLineColor(2);
  ALICEMidGlobalSysJPsi->SetFillColor(2);
  ALICEMidGlobalSysJPsi->Draw("same"); 

  lgd->AddEntry(grRaaALICEMid,"ALICE Data p_{T}^{J/#psi} > 0.0 GeV/c, |y^{J/#psi}| #leq 1.0","P");
  
}

void Draw_ALICEFor_JPsi_RaaVsNpart(TLegend *lgd)
{
  int nbinsALICE=9;
  Double_t NPartALICE[10]={357,262,187,128,86,53,30,16,8};
  Double_t ErrNPartALICE[10]={0};
  Double_t RaaALICE[10] = {0.47,0.48,0.51,0.51,0.52,0.61,0.70,0.74,0.94}; 
  Double_t SystErrALICE[10] = {0.03,0.02,0.02,0.02,0.03,0.05,0.06,0.09,0.1};
  
  TGraphErrors *grRaaALICE = new TGraphErrors(nbinsALICE, NPartALICE, RaaALICE, ErrNPartALICE, SystErrALICE);  
  grRaaALICE->SetMarkerStyle(20);
  grRaaALICE->SetMarkerColor(4);
  grRaaALICE->GetYaxis()->SetRangeUser(0,1.7);
  
  TAxis *Xaxis2 = grRaaALICE->GetXaxis();
  Xaxis2->SetLimits(0.,430.0);
  
  grRaaALICE->GetXaxis()->SetTitle("N_{Part}");
  grRaaALICE->GetYaxis()->SetTitle("R_{AA}");

  grRaaALICE->Draw("AP");
 
  TLine *lh3 = new TLine(0.0,1.0,420,1.0);
  lh3->SetLineColor(1);
  lh3->SetLineStyle(1);
  lh3->SetLineWidth(1.5);
  lh3->Draw("same");
 
  //TLatex *tb= new TLatex;
  //tb->SetNDC(); 
  //tb->SetTextAlign(12);
  //tb->SetTextColor(1);
  //tb->SetTextSize(0.040);
  //tb->DrawLatex(0.22,0.22,"PbPb #sqrt{s_{NN}} = 2.76 TeV");
  //tb->DrawLatex(0.22,0.16,"J/#psi #rightarrow #mu^{+} #mu^{-}, p_{T}^{J/#psi} > 0.0 GeV/c");  
   
  
  TBox *RaaJPsiALICESys[12];
  for(int j=0;j<9;j++){
    RaaJPsiALICESys[j] = new TBox(NPartALICE[j]-3,  RaaALICE[j]-SystErrALICE[j], NPartALICE[j]+3, RaaALICE[j]+SystErrALICE[j]);
  }
  
  for(int j=0;j<9;j++){
    RaaJPsiALICESys[j]->SetFillStyle(0000);
    RaaJPsiALICESys[j]->SetLineColor(4);
    RaaJPsiALICESys[j]->Draw("same"); 
  }
  
  TBox *ALICEGlobalSysJPsi;
  ALICEGlobalSysJPsi = new TBox(400-5, 1 - 0.14, 400+5, 1 + 0.14);
  ALICEGlobalSysJPsi->SetFillStyle(3001);
  ALICEGlobalSysJPsi->SetLineColor(4);
  ALICEGlobalSysJPsi->SetFillColor(4);
  ALICEGlobalSysJPsi->Draw("same"); 

  Draw_ALICEMid_JPsi_RaaVsNpart(lgd);

  lgd->AddEntry(grRaaALICE,"ALICE Data p_{T}^{J/#psi} > 0.0 GeV/c, 2.5 #leq y^{J/#psi} #leq 4.0","P");

}

void Draw_CMS_JPsi_RaaVsRap(TLegend *lgd)
{
  //=============== CMS Raa Vs Rap Data ===============================================================//
  //AvpT 10.92,9.65,8.92
  int nbinsRapCMS=3;
  Double_t RapCMSD[3]={0.6,1.4,2.0};
  Double_t ErrRapCMS[3]={0.6,0.2,0.4};
  Double_t RaaRapCMS[3] = {0.31,0.33,0.36}; 
  Double_t RaaRapStatErrCMS[3] = {0.02,0.03,0.03};
  Double_t RaaRapSystErrCMS[3] = {0.03,0.04,0.04};

  TGraphErrors *grRaaRapCMS = new TGraphErrors(nbinsRapCMS, RapCMSD, RaaRapCMS, ErrRapCMS, RaaRapStatErrCMS);  
  grRaaRapCMS->SetMarkerStyle(20);
  grRaaRapCMS->SetMarkerColor(2);
  grRaaRapCMS->GetYaxis()->SetRangeUser(0,1.5);
  grRaaRapCMS->GetXaxis()->SetTitle("|y|");
  grRaaRapCMS->GetYaxis()->SetTitle("R_{AA}");
  
  TAxis *XaxisgrRaaRapCMS = grRaaRapCMS->GetXaxis();
  XaxisgrRaaRapCMS->SetLimits(0.0,2.4);
  
  grRaaRapCMS->Draw("AP"); 

  TLine *lh_grRaaRapCMS = new TLine(0.0,1.0,2.4,1.0);
  lh_grRaaRapCMS->SetLineColor(1);
  lh_grRaaRapCMS->SetLineStyle(1);
  lh_grRaaRapCMS->SetLineWidth(1.5);
  lh_grRaaRapCMS->Draw("same");

  TLatex *tb= new TLatex;
  tb->SetNDC(); 
  tb->SetTextAlign(12);
  tb->SetTextColor(1);
  tb->SetTextSize(0.040);
  tb->DrawLatex(0.20,0.20,"Pb+Pb #sqrt{s_{NN}} = 2.76 TeV");
  tb->DrawLatex(0.20,0.15,"J/#psi #rightarrow #mu^{+} #mu^{-}, p_{T}^{J/#psi} > 6.5 GeV/c");  
  
  TBox *RaaRapJPsiCMSSys[4];
  for(int j=0;j<3;j++){
    RaaRapJPsiCMSSys[j] = new TBox(RapCMSD[j]-0.05,  RaaRapCMS[j]-RaaRapSystErrCMS[j], RapCMSD[j]+0.05,  RaaRapCMS[j]+RaaRapSystErrCMS[j]);
  }
  
  for(int j=0;j<3;j++){
    RaaRapJPsiCMSSys[j]->SetFillStyle(0000);
    RaaRapJPsiCMSSys[j]->SetLineColor(2);
    RaaRapJPsiCMSSys[j]->Draw("same"); 
  }
  
  TBox *CMSGlobalSysJPsiRap;
  CMSGlobalSysJPsiRap = new TBox(0.2-0.05, 1 - 0.05, 0.2+0.05, 1 + 0.05);
  CMSGlobalSysJPsiRap->SetFillStyle(3001);
  CMSGlobalSysJPsiRap->SetLineColor(2);
  CMSGlobalSysJPsiRap->SetFillColor(2);
  CMSGlobalSysJPsiRap->Draw("same"); 


  lgd->AddEntry(grRaaRapCMS,"CMS Data", "P");
}

void HiggsPlot::PlotBF(int iDecay, double tBmH_max, double BF_max){
  if(iDecay<0 || iDecay>2) {cout<<"iDecay must be 0, 1 or 2"<<endl; return;}
  styles style; style.setPadsStyle(-1); style.setDefaultStyle();
  int nBins = 1000;
  TString hName, epsName = "public_html/Higgs_BF_TEMP_BaBar.eps", label, Llabel, Rlabel;
  TString yTitle[] = {"BF(B#rightarrow#tau#nu) (10^{-5})", "R(D)", "R(D*)"};
  TString TagDecay[] = {"BF", "R(D)", "R(D*)"};
  TCanvas can;
  TH1F *hBF[2];
  for(int his=0; his<1; his++) {
    hName = "hBF"; hName += his;
    hBF[his] = new TH1F(hName,"",nBins,0,tBmH_max);
  }
  double tBmH, BF[2];
  for(int bin=1; bin<=nBins; bin++){
    tBmH = hBF[0]->GetBinCenter(bin);
    Compute(tBmH,BF,iDecay);
    hBF[0]->SetBinContent(bin, BF[0]);
    hBF[0]->SetBinError(bin, BF[1]);
  }
  hBF[0]->SetFillColor(2); hBF[0]->SetLineColor(2);
  hName += "1";
  hBF[1] = static_cast<TH1F*>(hBF[0]->Clone(hName));
  for(int bin=1; bin<=nBins; bin++)hBF[1]->SetBinError(bin,0);
  hBF[1]->SetFillColor(0);hBF[1]->SetLineColor(1); hBF[1]->SetLineWidth(2);

  TBox box; box.SetLineColor(4);box.SetFillColor(4);
  TLine line; line.SetLineStyle(1); line.SetLineColor(4); line.SetLineWidth(3);
  if(BF_max>0) hBF[0]->SetMaximum(BF_max);
  hBF[0]->Draw("e3");
  box.SetFillStyle(3002); 
  box.DrawBox(0,Measurement[iDecay][0]-Measurement[iDecay][1],
	      tBmH_max,Measurement[iDecay][0]+Measurement[iDecay][1]);
  line.DrawLine(0,Measurement[iDecay][0],tBmH_max,Measurement[iDecay][0]);
  hBF[0]->Draw("e3 same");
  hBF[1]->Draw("c same");

  Compute(0,BF,iDecay);
  label = "#splitline{"; label += TagDecay[iDecay]; label += "_{SM} = ";
  if(iDecay==0){
    label+="(";label+=RoundNumber(BF[0],1); label+=" #pm ";
    label+=RoundNumber(BF[1],1); label+=")#times10^{-5}}{BF_{exp} = (";
    label+=RoundNumber(Measurement[iDecay][0],1); label+=" #pm ";
    label+=RoundNumber(Measurement[iDecay][1],1); label+=")#times10^{-5}}";
    Llabel = ""; Rlabel = label;
  } else {
    label+=RoundNumber(BF[0],3); label+=" #pm ";
    label+=RoundNumber(BF[1],3); label+="}{"; label += TagDecay[iDecay]; label += "_{exp} = ";
    label+=RoundNumber(Measurement[iDecay][0],3); label+=" #pm ";
    label+=RoundNumber(Measurement[iDecay][1],3); label+="}";
    Rlabel = ""; Llabel = label;
  }
  style.setTitles(hBF[0],"tan#beta/m_{H^{+}} (GeV^{-1})",yTitle[iDecay],Llabel,Rlabel);
  epsName.ReplaceAll("TEMP",DecayName[iDecay]);
  can.SaveAs(epsName);
  for(int his=0; his<2; his++) hBF[his]->Delete();
}

void Draw_CMS_JPsi_RaaVsNpart(TLegend *lgd)
{

int nbins = 12;
  Double_t RaaCMS[20]={0.23,0.24,0.29,0.33,0.38,0.40,0.45,0.41,0.52,0.52,0.64,0.64};
  Double_t StatErrCMS[20]={0.01,0.01,0.02,0.02,0.02,0.03,0.03,0.03,0.04,0.04,0.05,0.05};
  Double_t SystErrCMS[20]={0.02,0.03,0.03,0.03,0.04,0.04,0.05,0.04,0.05,0.05,0.07,0.13};
  Double_t NPartCMS[20]={Npart(0,2),Npart(2,4),Npart(4,6),Npart(6,8),Npart(8,10),Npart(10,12),
			 Npart(12,14),Npart(14,16),Npart(16,18),Npart(18,20),Npart(20,24),
			 Npart(24,40)};
  Double_t ErrNPartCMS[20]={0};
  
  TGraphErrors *grRaaCMS = new TGraphErrors(nbins, NPartCMS, RaaCMS, ErrNPartCMS, StatErrCMS);  
  grRaaCMS->SetMarkerStyle(21);
  grRaaCMS->SetMarkerColor(2);
  grRaaCMS->GetYaxis()->SetRangeUser(0,1.7);
  grRaaCMS->GetXaxis()->SetTitle("N_{Part}");
  grRaaCMS->GetYaxis()->SetTitle("R_{AA}");
  
  grRaaCMS->Draw("AP");
  TBox *RaaJPsiSys[12];
  for(int j=0;j<12;j++){
    RaaJPsiSys[j] = new TBox(NPartCMS[j]-3,  RaaCMS[j]-SystErrCMS[j], NPartCMS[j]+3, RaaCMS[j]+SystErrCMS[j]);
  }
  
  for(int j=0;j<12;j++){
    RaaJPsiSys[j]->SetFillStyle(0000);
    RaaJPsiSys[j]->SetLineColor(6);
    RaaJPsiSys[j]->Draw("same"); 
  }
  
  TLine *lh2 = new TLine(0.0,1.0,400,1.0);
  lh2->SetLineColor(1);
  lh2->SetLineStyle(1);
  lh2->SetLineWidth(2);
  lh2->Draw("same");
  
  TBox *GlobalSysJPsi;
  GlobalSysJPsi = new TBox(400-5, 1 - 0.06, 400+5, 1 + 0.06);
  GlobalSysJPsi->SetFillStyle(3001);
  GlobalSysJPsi->SetLineColor(6);
  GlobalSysJPsi->SetFillColor(6);
  GlobalSysJPsi->Draw("same"); 
  
  TLatex *tb= new TLatex;
  tb->SetNDC(); 
  tb->SetTextAlign(12);
  tb->SetTextColor(1);
  tb->SetTextSize(0.040);
  tb->DrawLatex(0.17,0.23,"PbPb #sqrt{s_{NN}} = 2.76 TeV");
  tb->DrawLatex(0.17,0.17,"p_{T}^{J/#psi} > 6.5 GeV/c");

  lgd->AddEntry(grRaaCMS,"CMS Data", "P");

}

//--------------------------
// DrawTargetBoundaries
//--------------------------
void DrawTargetBoundaries(int color=target_color)
{
  double Zlo = 50.0;
  double Zhi = Zlo + 30.0;
  double Rlo = 0.0;
  double Rhi = Rlo + 1.5;
  if (full) { Rlo=-1.5; Rhi=Rlo+3; }

  TBox *box = new TBox(Zlo, Rlo, Zhi, Rhi);
  box->SetLineWidth(2.0);
  if (fill_detectors) { box->SetFillColor(color); box->SetFillStyle(1001); }
  else box->SetFillStyle(0);
  box->SetLineColor(color);
  box->Draw();

  TLatex *lab = new TLatex(Zlo-5, 0, "target");
  if (full) lab->SetTextAlign(32); else lab->SetTextAlign(31);
  lab->SetTextSize(0.02);
  lab->SetTextColor(color);
  lab->Draw();
}

void Draw_CMS_JPsi_RaaVspT(TLegend *legend_ratio)
{
  
  int nbinsPtCMS=4;
  Double_t PtCMSD[4]={7.31,8.97,11.32,16.52};
  Double_t ErrPtCMS[4]={0};
  Double_t RaaPtCMS[4] = {0.35,0.34,0.34,0.29}; 
  Double_t RaaPtStatErrCMS[4] = {0.09,0.03,0.03,0.04};
  Double_t RaaPtSystErrCMS[4] = {0.04,0.04,0.03,0.03};
  
  TGraphErrors *grRaaPtCMS = new TGraphErrors(nbinsPtCMS, PtCMSD, RaaPtCMS, ErrPtCMS, RaaPtStatErrCMS);  
  grRaaPtCMS->SetMarkerStyle(21);
  grRaaPtCMS->SetMarkerColor(kRed+1);
  grRaaPtCMS->SetMarkerSize(1.3);
  //grRaaPtCMS->SetFillColor(kRed);
  grRaaPtCMS->SetLineColor(1);
  grRaaPtCMS->GetYaxis()->SetRangeUser(0,1.5);
  grRaaPtCMS->GetXaxis()->SetTitle("p_{T}(GeV/c)");
  grRaaPtCMS->GetYaxis()->SetTitle("R_{AA}");
  
  TLine *lh5 = new TLine(0.0,1.0,20.0,1.0);
  lh5->SetLineColor(1);
  lh5->SetLineStyle(1);
  lh5->SetLineWidth(1.5);
  
  TLatex *tb= new TLatex;
  tb->SetNDC(); 
  tb->SetTextAlign(12);
  tb->SetTextColor(1);
  tb->SetTextSize(0.040);
  
  TBox *RaaPtJPsiCMSSys[4];
  for(int j=0;j<4;j++){
    RaaPtJPsiCMSSys[j] = new TBox(PtCMSD[j]-0.4,  RaaPtCMS[j]-RaaPtSystErrCMS[j], PtCMSD[j]+0.4,  RaaPtCMS[j]+RaaPtSystErrCMS[j]);
  }
  
  for(int j=0;j<4;j++){
    RaaPtJPsiCMSSys[j]->SetFillStyle(1001);
    RaaPtJPsiCMSSys[j]->SetFillColor(kPink+1);
    RaaPtJPsiCMSSys[j]->SetLineColor(1);
  }
  
  TBox *CMSGlobalSysJPsiPt;
  //CMSGlobalSysJPsiPt = new TBox(27.0-1.0, 1 - 0.075, 27.0+1.0, 1 + 0.075);
  CMSGlobalSysJPsiPt = new TBox(19.165,0.923,19.948,1.067);
  CMSGlobalSysJPsiPt->SetFillStyle(3001);
  CMSGlobalSysJPsiPt->SetLineColor(2);
  CMSGlobalSysJPsiPt->SetFillColor(2);
  
  TAxis *XaxisgrRaaPtCMS = grRaaPtCMS->GetXaxis();
  XaxisgrRaaPtCMS->SetLimits(0.0,20.0);
  grRaaPtCMS->Draw("AP");
  for(int j=0;j<4;j++) RaaPtJPsiCMSSys[j]->Draw("same"); 
  grRaaPtCMS->Draw("Psame");
  CMSGlobalSysJPsiPt->Draw("same"); 
  tb->SetTextSize(0.030);
  //tb->DrawLatex(0.20,0.20,"PbPb #sqrt{s_{NN}} = 2.76 TeV");
  tb->DrawLatex(0.71476,0.91289,"J/#psi #rightarrow #mu^{+} #mu^{-}");  
  lh5->Draw("same");
  //legend_ratio->AddEntry(grRaaPtCMS,"PbPb #sqrt{s_{NN}} = 2.76 TeV and Au-Au #sqrt{s_{NN}} = 0.2 TeV", "");
  legend_ratio->AddEntry(grRaaPtCMS,"Pb+Pb #sqrt{s_{NN}} = 2.76 TeV", "");
  legend_ratio->AddEntry(grRaaPtCMS,"Prompt J/#psi CMS: |y| < 2.4", "P");

  //legend_ratio->AddEntry(grRaaPtCMS,"PbPb #sqrt{s_{NN}} = 2.76 TeV and Au-Au #sqrt{s_{NN}} = 0.2 TeV", "");
  //legend_ratio->AddEntry(grRaaPtCMS,"CMS: J/#psi #rightarrow #mu^{+} #mu^{-}, |y| < 2.4 & p_{T}^{J/#psi} > 6.5 GeV/c", "P");
  //legend_ratio->AddEntry(grRaaPtCMS,"CMS: J/#psi #rightarrow #mu^{+} #mu^{-}, |y| < 2.4 & p_{T}^{J/#psi} > 6.5 GeV/c", "P");
}

//.........这里部分代码省略.........
  h0->GetXaxis()->SetLabelOffset(0.005);
  h0->GetXaxis()->SetLabelSize(0.04);
  h0->GetXaxis()->SetLabelFont(42);
  h0->GetXaxis()->SetTitleFont(42);
  h0->GetYaxis()->SetNdivisions(505);
  h0->GetYaxis()->CenterTitle();
  h0->GetYaxis()->SetTitle("Anisotropy Parameter, v_{2}");
  h0->GetYaxis()->SetTitleOffset(.9);
  h0->GetYaxis()->SetTitleSize(0.05);
  h0->GetYaxis()->SetLabelOffset(0.005);
  h0->GetYaxis()->SetLabelSize(0.04);
  h0->GetYaxis()->SetLabelFont(42);
  h0->GetYaxis()->SetTitleFont(42);
  h0->Draw("c");

  
  double pt[2] = {1.766, 4.459};
  double pt1Bin[1] = {1.866};
  double v2[2] = {0.08329931, 0.126231279};
  double eptLow[2] = {.766, .459};
  double eptHigh[2] = {2.234, 2.341};
  double eptSyst[2] = { 0.0, 0.0};
  double eptSyst1Bin[1] = { 0.0};
  double ev2Low[2] = { 0.00736501, 0.012241649};
  double ev2High[2] = { 0.00736501, 0.012241649};
  double ev2Syst[2] = {.0124948965, .02145931743 };
  
  
  TGraphAsymmErrors *AlexData = new TGraphAsymmErrors(2, pt, v2, eptLow, eptHigh, ev2Low, ev2High);
  
  AlexData->SetMarkerColor(2);
  AlexData->SetMarkerStyle(21);
  AlexData->SetMarkerSize(2);
  AlexData->SetLineColor(2);
  AlexData->SetLineWidth(3);
  
  TGraphErrors *AlexDataSyst = new TGraphErrors(2, pt, v2, eptSyst, ev2Syst);
  
  AlexDataSyst->SetMarkerColor(2);
  AlexDataSyst->SetMarkerStyle(21);
  AlexDataSyst->SetMarkerSize(1.1);
  AlexDataSyst->SetLineColor(2);
  AlexDataSyst->SetLineWidth(33);
  
  
  
  // 8 data points
   
 
  
  double weights[6] = {0.4034, .3058, .1696, .0778, .0316, .0119};

  double avgV2[1] = {0};
  double avgV2StatErrorSquared[1] = {0};
  double avgV2SystErrorSquared[1] = {0};
  double avgV2StatError[1] = {0};
  double avgV2SystError[1] = {0};
  
  for(int i = 0; i < 6; i++){

    cout << "V2: " << y_data_cen[1][i] << endl;
    
    avgV2[0] = avgV2[0] + (weights[i]*y_data_cen[1][i]);
    
    avgV2StatErrorSquared[0] = avgV2StatErrorSquared[0] + (weights[i]*weights[i])*(ye_data_cen[1][i]*ye_data_cen[1][i]);
    avgV2SystErrorSquared[0] = avgV2SystErrorSquared[0] + (weights[i]*weights[i])*(yes_data_cen[1][i]*yes_data_cen[1][i]);

void CorrPtCut(char *prodname = "LHC10d4 PYTHIA D6T 7 TeV p+p", char *shortprodname = "LHC10d4", char *filename="Et.ESD.new.sim.LHC10d4.pp.merged.root"){

  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  TCanvas *c = new TCanvas("c","c",500,400);
  c->SetTopMargin(0.04);
  c->SetRightMargin(0.04);
  c->SetLeftMargin(0.181452);
  c->SetBottomMargin(0.134409);
  c->SetBorderSize(0);
  c->SetFillColor(0);
  c->SetFillColor(0);
  c->SetBorderMode(0);
  c->SetFrameFillColor(0);
  c->SetFrameBorderMode(0);
  float etacut = 0.7;
  cout<<"Pt cut = 150 MeV/c"<<endl;
  TH1D *High = GetHisto(0.15-.001,"High",filename,etacut);
  float tpcHigh=highbound;
  float tpcLow=lowbound;
  float tpcsyserr = syserr;
  float tpcmean = mean;
  float x1 = High->GetXaxis()->GetBinLowEdge(1);
  //TBox *tpcBox = new TBox(-x1*.99,1.0-tpcLow,x1*.99,1.0-tpcHigh);
  TBox *tpcBox = new TBox(-x1*.99,1.0-(mean-syserr),x1*.99,1.0-(mean+syserr));
  tpcBox->SetFillColor(5);
  tpcBox->SetLineColor(0);
  tpcBox->SetFillStyle(1001);
  cout<<"Pt cut = 100 MeV/c"<<endl;
  TH1D *Low = GetHisto(0.1-.001,"Low",filename,etacut);
  float itsHigh=highbound;
  float itsLow=lowbound;
  float itssyserr = syserr;
  float itsmean = mean;

  cout<<Form("dataset & %2.4f \\pm %2.4f &  %2.4f \\pm %2.4f \\",itsmean,itssyserr,tpcmean,tpcsyserr)<<endl;
  float x = Low->GetXaxis()->GetBinLowEdge(1);
  //TBox *itsBox = new TBox(-x*.99,1.0-itsLow,x*.99,1.0-itsHigh);
  TBox *itsBox = new TBox(-x1*.99,1.0-(mean-syserr),x1*.99,1.0-(mean+syserr));
  itsBox->SetFillColor(5);
  itsBox->SetLineColor(0);
  itsBox->SetFillStyle(1001);
  cout<<"Pt cut = 50 MeV/c"<<endl;
  TH1D *Lowest = GetHisto(0.05-.001,"Lowest",filename,etacut);
  TF1 *funcLow = new TF1("funcLow","[0]",-.7,.7);
  funcLow->SetParameter(0,0.01);
  Low->Fit(funcLow);
  TF1 *funcHigh = new TF1("funcHigh","[0]",-.7,.7);
  funcHigh->SetParameter(0,0.02);
  High->Fit(funcLow);
  High->SetMaximum(0.06);
  High->SetMinimum(0.0);
  High->SetMarkerColor(2);
  Low->SetMarkerColor(4);
  High->SetLineColor(2);
  Low->SetLineColor(4);
  High->SetMinimum(0.0);
  High->SetMarkerStyle(20);
  Low->SetMarkerStyle(21);
  Lowest->SetMarkerStyle(22);
  High->Draw();
  tpcBox->Draw("f");
  High->Draw("same");
  itsBox->Draw("f");
  //return;
  Low->Draw("same");
  //Lowest->Draw("same");
  TLatex *tex = new TLatex(-0.723444,0.0373593+0.019,prodname);
  tex->SetTextSize(0.0537634);
  tex->Draw();
  TLegend *leg = new TLegend(0.217742,0.696237,0.477823,0.873656);
  leg->AddEntry(High,"p_{T} cut-off = 0.15 GeV/c");
  leg->AddEntry(Low,"p_{T} cut-off = 0.1 GeV/c");
  //leg->AddEntry(Lowest,"p_{T} cut-off = 0.05 GeV/c");
  leg->SetFillStyle(0);
  leg->SetFillColor(0);
  leg->SetBorderSize(0);
  leg->SetTextSize(0.0537634);
  leg->Draw();
  c->SaveAs(Form("pics/%s/fptcut.eps",shortprodname));
  c->SaveAs(Form("pics/%s/fptcut.png",shortprodname));
  c->SaveAs(Form("pics/%s/fptcut.pdf",shortprodname));
}

void plotXY(char* fname){

  // input: root file from GBL

  gStyle->SetOptStat(0);
  Double_t w = 1.2; // number of RMS for gaussian fit
  TString filename = fname;
  TCanvas *cTop = new TCanvas("cTop","Residuals in lab frame, TOP", 800, 1000);
  TCanvas *cBot = new TCanvas("cBot","Residuals in lab frame, BOTTOM",800, 1000);
  cTop->SetFillColor(kBlue-10); cTop->SetFillStyle(3001);
  cBot->SetFillColor(kMagenta-10); cBot->SetFillStyle(3001);
  Double_t halfWid = 38.34/2.;
  Double_t halfLen = 98.33/2.;
  TBox *module = new TBox(-halfLen, -halfWid, halfLen, halfWid);
  module->SetLineColor(kMagenta);
  module->SetFillStyle(0);

  Int_t ic=0;
  TFile *f = new TFile(filename.Data());
  // extract string 
  Int_t idx2 = filename.Index("/");
  Int_t idx1 = filename.Index("_");
  TString outn = filename(idx1+1,idx2-(idx1+1));
  TString outname ="predXY_"; outname += outn.Data(); outname += ".txt";
  TString topName ="predXY_"; topName += outn.Data(); topName += "_topO.gif";
  TString botName ="predXY_"; botName += outn.Data(); botName += "_botO.gif";
  //  ofstream outf(outname.Data());
  // top
  Int_t nrow = 13;
  Double_t space = 0.007;
  Double_t space1 = 0.03;
  Double_t y2 = 1./13*nrow-space1;
  Double_t y1 = 1./13*(--nrow)-space1;
  TPad *tp1 = new TPad("tp1","tp1",0.33,y1,0.66,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp2 = new TPad("tp2","tp2",0.33,y1,0.66,y2);
  y2 = y1-space; y1 = 1./13.*(--nrow)-space1;
  TPad *tp3 = new TPad("tp3","tp3",0.33,y1,0.66,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp4 = new TPad("tp4","tp4",0.33,y1,0.66,y2);
  y2 = y1-space; y1 = 1./13.*(--nrow)-space1;
  TPad *tp5 = new TPad("tp5","tp5",0.33,y1,0.66,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp6 = new TPad("tp6","tp6",0.33,y1,0.66,y2);
  y2 = y1-space; y1 = 1./13.*(--nrow)-space1;
  Double_t off = 0.12;
  off = 0.165;
  y1 -= space1; y2 -= space1;
  space1 *= 2;
  TPad *tp7 = new TPad("tp7","tp7",0.+off,y1,0.33+off,y2);
  TPad *tp8 = new TPad("tp8","tp8",0.66-off,y1,1.-off,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp9 = new TPad("tp9","tp9",0.+off,y1,0.33+off,y2);
  TPad *tp10 = new TPad("tp10","tp10",0.66-off,y1,1.-off,y2);
  y2 = y1-space; y1 = 1./13.*(--nrow)-space1;
  TPad *tp11 = new TPad("tp11","tp11",0.+off,y1,0.33+off,y2);
  TPad *tp12 = new TPad("tp12","tp12",0.66-off,y1,1.-off,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp13 = new TPad("tp13","tp13",0.+off,y1,0.33+off,y2);
  TPad *tp14 = new TPad("tp14","tp14",0.66-off,y1,1.-off,y2);
  y2 = y1-space; y1 = 1./13.*(--nrow)-space1;
  TPad *tp15 = new TPad("tp15","tp15",0.+off,y1,0.33+off,y2);
  TPad *tp16 = new TPad("tp16","tp16",0.66-off,y1,1.-off,y2);
  y2 = y1; y1 = 1./13.*(--nrow)-space1;
  TPad *tp17 = new TPad("tp17","tp17",0.+off,y1,0.33+off,y2);
  TPad *tp18 = new TPad("tp18","tp18",0.66-off,y1,1.-off,y2);
  cout << " TOP " << endl;

  cTop->cd();
  tp1->Draw();
  tp2->Draw();
  tp3->Draw();
  tp4->Draw();
  tp5->Draw();
  tp6->Draw();
  tp7->Draw();
  tp8->Draw();
  tp9->Draw();
  tp10->Draw();
  tp11->Draw();
  tp12->Draw();
  tp13->Draw();
  tp14->Draw();
  tp15->Draw();
  tp16->Draw();
  tp17->Draw();
  tp18->Draw();

  ic = 2;
  Int_t ipad=0;
  for(Int_t i=1; i<4; i++){
    TString hisname = "h_xy_module_L"; hisname  += i;
    hisname += "t_halfmodule_axial_sensor0";
    TString layer = "L"; layer+=i; layer+= "TA"; 
    TH2D *his2 = (TH2D*) f->Get(hisname.Data());
    TH1D *his = (TH1D*)his2->ProjectionY();
    ipad++;
    TString pd = "tp"; pd+=ipad;
    TPad *pad = (TPad*) cTop->GetListOfPrimitives()->FindObject(pd.Data());
    pad->cd();
//.........这里部分代码省略.........

void FakePhotonSystematic(){
  
  TFile FMCRSSig("DsTolnu/SP1005Skim/DsReco/Final/HistosForFit.root");
  TFile FMCWSSig("DsTolnu/SP1005Skim/DsReco/FinalWS/HistosForFit.root");

  //TFile FMCRS("DsTolnu/MC/Final/HistosForFit.root");
  //TFile FMCWS("DsTolnu/MC/FinalWS/HistosForFit.root");
  TFile FMCRS("DsTolnu/SP1005Skim/DsReco/Final/HistosForFit.root");
  TFile FMCWS("DsTolnu/SP1005Skim/DsReco/FinalWS/HistosForFit.root");

  TFile FDataRS("DsTolnu/Data/Final/HistosForFit.root");
  TFile FDataWS("DsTolnu/Data/FinalWS/HistosForFit.root");

  

  TCanvas C;
  C.Clear();
  C.Print("FakePhotonSystematic.ps[");

  C.Clear();
  C.Divide(2,2);
  C.cd(1);  
//   TH1F* HMCRS=(TH1F*)FMCRS.Get("HEvtMultiplicity");
//   HMCRS->SetTitle(TString("MC RS =")+(long)HMCRS->GetMean()+"."+(long)(1000*(HMCRS->GetMean()-1.)));
//   HMCRS->Draw();
//   C.cd(2);
//   TH1F* HMCWS=(TH1F*)FMCWS.Get("HEvtMultiplicity");
//   HMCWS->SetTitle(TString("MC WS =")+(long)HMCWS->GetMean()+"."+(long)(1000*(HMCWS->GetMean()-1.)));
//   HMCWS->Draw();
  C.cd(3);  
  TH1F* HDataRS=(TH1F*)FDataRS.Get("HEvtMultiplicity");
  HDataRS->SetTitle(TString("Data RS =")+(long)HDataRS->GetMean()+"."+(long)(1000*(HDataRS->GetMean()-1.)));
  HDataRS->Draw();
  C.cd(4);
  TH1F* HDataWS=(TH1F*)FDataWS.Get("HEvtMultiplicity");
  HDataWS->SetTitle(TString("Data WS =")+(long)HDataWS->GetMean()+"."+(long)(1000*(HDataWS->GetMean()-1.)));
  HDataWS->Draw();
  C.Print("FakePhotonSystematic.ps");



  C.Clear();
  C.Divide(2,2);
  C.cd(1);  
  TH1F* HSigMCRS=(TH1F*)FMCRSSig.Get("HEvtMultiplicitySig");
  HSigMCRS->SetTitle(TString("MC Sig RS =")+(long)HSigMCRS->GetMean()+"."+(long)(1000*(HSigMCRS->GetMean()-1.)));
  HSigMCRS->Draw();
  C.cd(2);
  TH1F* HSigMCWS=(TH1F*)FMCWSSig.Get("HEvtMultiplicitySig");
  HSigMCWS->SetTitle(TString("MC Sig WS =")+(long)HSigMCWS->GetMean()+"."+(long)(1000*(HSigMCWS->GetMean()-1.)));
  HSigMCWS->Draw();
  C.cd(3);  
  TH1F* HBkgMCRS=(TH1F*)FMCRS.Get("HEvtMultiplicityBkg");
  HBkgMCRS->SetTitle(TString("MC Bkg RS =")+(long)HBkgMCRS->GetMean()+"."+(long)(1000*(HBkgMCRS->GetMean()-1.)));
  HBkgMCRS->Draw();
  C.cd(4);
  TH1F* HBkgMCWS=(TH1F*)FMCWS.Get("HEvtMultiplicityBkg");
  HBkgMCWS->SetTitle(TString("MC Bkg WS =")+(long)HBkgMCWS->GetMean()+"."+(long)(1000*(HBkgMCWS->GetMean()-1.)));
  HBkgMCWS->Draw();
  C.Print("FakePhotonSystematic.ps");
  
  C.Clear();
  THStack HMCStack("HMCStack","");
  HBkgMCRS->SetFillColor(1);HMCStack.Add(HBkgMCRS);
  HSigMCRS->SetFillColor(2);HSigMCRS->SetLineColor(2);HMCStack.Add(HSigMCRS);
  HMCStack.Draw();
  C.Print("FakePhotonSystematic.ps");

  //---------------------------------------------
  //subtract background from Data
  //---------------------------------------------
  TH1F*HDataRSBkg=(TH1F*)HDataWS->Clone("HDataRSBkg");
  TH1F*HBkgMCWSCl=(TH1F*)HBkgMCWS->Clone("HBkgMCWSCl");
  //need to fix the WS first by subtracting the signal component
  HSigMCWS->Scale(11350./HSigMCWS->Integral());//scale to number of WS signal events in Data
  HDataRSBkg->Add(HSigMCWS,-1);
  //construct correction ratio from MC to correct the WS distribution
  TH1F*HBkgMCRSCl=(TH1F*)HBkgMCRS->Clone("HBkgMCRSCl");
  HBkgMCWSCl->Scale(1./HBkgMCWSCl->Integral());
  HBkgMCRSCl->Scale(1./HBkgMCRSCl->Integral());
  TH1F*HBkgMCRatio=(TH1F*)HBkgMCRSCl->Clone("HBkgMCRatio");
  HBkgMCRatio->Divide(HBkgMCWSCl);
  C.Clear();
  delete HBkgMCRatio->FindObject("stats");
  HBkgMCRatio->SetFillColor(0);HBkgMCRatio->SetLineColor(1);
  HBkgMCRatio->SetStats(0);
  HBkgMCRatio->GetYaxis()->SetRangeUser(0,1.5);
  HBkgMCRatio->Draw("hist");
  C.Print("FakePhotonSystematic.ps");
  //correct the WS distribution
  HDataRSBkg->Multiply(HBkgMCRatio);
  C.Clear();
  HDataWS->Scale(1./HDataWS->Integral());
  HDataRSBkg->Scale(1./HDataRSBkg->Integral());
  delete HDataWS->FindObject("stats");
  HDataWS->SetTitle("");
  HDataWS->GetXaxis()->SetTitle("Reconstruction Multiplicity");
  HDataWS->Draw("p");
  HDataRSBkg->SetLineColor(0);
  HDataRSBkg->SetFillColor(4);
//.........这里部分代码省略.........

void Test_UnfoldQ2(int Type = 0){
  TString TypeName[] = {"Efficiency", "Unfold"};
  if(Type==-1) {
    for(int type=0; type<nType; type++)
      cout<<type<<" "<<TypeName[type]<<", "; 
    cout<<endl; 
    return;
  }
  BToDtaunu Dtaunu; BToDstaunu Dstaunu;
  int nBins = 18, nValues = 50, nBinsRes = 80;
  double limQ2[] = {3.5, 12.5}, limRes = 2, ml[] = {Dtaunu.mTau, Dtaunu.mTau, Dtaunu.mMu};
  if(Type==0){limQ2[0] = 4; nBins = 17;}

  TString mcCuts[3][6] = {{"weight*((MCType==5&&candType==1||MCType==11&&candType==3)&&candM2>1.5)",
			   "weight*((MCType==6&&candType==2||MCType==12&&candType==4)&&candM2>1.5)" ,
			   "weight*((MCType==5&&candType==1||MCType==11&&candType==3))",
			   "weight*((MCType==6&&candType==2||MCType==12&&candType==4))",
			   "weight*(((MCType==1||MCType==3)&&candType==1||(MCType==7||MCType==9)&&candType==3))",
			   "weight*(((MCType==2||MCType==4)&&candType==2||(MCType==8||MCType==10)&&candType==4))"},
			  {"weight*((MCType==11&&candType==3)&&candM2>1.5)",
			   "weight*((MCType==12&&candType==4)&&candM2>1.5)" ,
			   "weight*((MCType==11&&candType==3))",
			   "weight*((MCType==12&&candType==4))",
			   "weight*(((MCType==7||MCType==9)&&candType==3))",
			   "weight*(((MCType==8||MCType==10)&&candType==4))"},
			  {"weight*((MCType==5&&candType==1)&&candM2>1.5)",
			   "weight*((MCType==6&&candType==2)&&candM2>1.5)" ,
			   "weight*((MCType==5&&candType==1))",
			   "weight*((MCType==6&&candType==2))",
			   "weight*(((MCType==1||MCType==3)&&candType==1))",
			   "weight*(((MCType==2||MCType==4)&&candType==2))"}};

  TString PadLabel[] = {"a", "b", "x", "d", "e", "f", "(g)", "(h)", "(i)"};
  TString Units[] = {" GeV^{2})", "0.35 GeV^{2})"}, hName;
  TString xTitle[] = {"Q", "R"};
  TBox box; box.SetLineColor(10);box.SetFillColor(10);

  gStyle->SetOptStat(0);
  int nRows = 2, nCols = 3;
  int Colors[] = {8,4,2}, ColorP[] = {4,1,2,8};
  TLatex label; label.SetTextFont(132); label.SetNDC(kTRUE);
  TCanvas can("dataMC","data Vs MC",700,150*nRows); 
  TPad *Pads[nPads][2];
  TH1D *hbini[nPads], *hxini[nPads][3], *hQ2[nPads][3], *hEff[nPads][3], *hRes[nPads], *hRes2[nPads];
  TH1D *hP[2][4], hTemp("hTemp","",100,0,3.4);
  TChain MC("ntp1");
  MC.Add("AWG82/ntuples/small/FitRAllNewx100_RunAll.root");

  TLegend *leg1[3];
  double scaleEff = 1;
  double dRows = (double)nRows, dCols = (double) nCols;
  double bMargin = 0.12, padH = (1-bMargin)/dRows, padW = 1/dCols, LeftMargin = 0.2;
  for(int col=0; col<nCols; col++){
    for(int row=0; row<nRows; row++){
      can.cd(0);
      int pad = nRows*col+row;
      double TextSize = 0.08;
      double RightMargin = 0.04, TopMargin=0, BottomMargin=bMargin/(bMargin+padH)*row;
      double PadXY[2][2] = {{padW*col, padW*(col+1)},{(padH+bMargin)*(dRows-1-row), bMargin+padH*(2-row)}};
      if(row==0) TextSize *= (bMargin+padH)/padH;
      hName = "Pad0_"; hName += pad;
      //cout<<hName<<": "<<PadXY[0][0]<<", "<<PadXY[1][0]<<", "<<PadXY[0][1]<<", "<<PadXY[1][1]<<endl;
      Pads[pad][0] = new TPad(hName,"",PadXY[0][0], PadXY[1][0], PadXY[0][1], PadXY[1][1]);
      Pads[pad][0]->SetLeftMargin(LeftMargin);     Pads[pad][0]->SetRightMargin(RightMargin); 
      Pads[pad][0]->SetBottomMargin(BottomMargin); Pads[pad][0]->SetTopMargin(TopMargin); 
      Pads[pad][0]->Draw(); Pads[pad][0]->cd();

      if(Type==0){
	for(int isNorm=2; isNorm>=0; isNorm--){
	  hName = "xini"; hName += pad; hName += isNorm;
	  hxini[pad][isNorm] = new TH1D(hName,"",nBins, limQ2[0], limQ2[1]);
	  formatHisto(hxini[pad][isNorm], Colors[isNorm], TextSize);
	  MC.Project(hName,"candQ2",mcCuts[col][row+2*isNorm]);
	  if(row==1) hxini[pad][isNorm]->SetXTitle(xTitle[0]);

	  //cout<<pad<<" Calculating q2"<<isNorm<<endl;
	  hName = "hQ2"; hName += pad; hName += isNorm;
	  hQ2[pad][isNorm] = new TH1D(hName,"",nBins, limQ2[0], limQ2[1]);
	  hQ2[pad][isNorm]->SetLineColor(Colors[isNorm]);
	  for(int bin=1; bin<=nBins; bin++){
	    double q2 = hQ2[pad][isNorm]->GetBinLowEdge(bin), val=0;
	    double dq2 = (hQ2[pad][isNorm]->GetBinLowEdge(bin)-q2)/(double)nValues;
	    q2 += dq2/2.;
	    for(int ival=0; ival<nValues; ival++) {
	      if(row==0) val += 1e16*Dtaunu.Compute(q2,1, ml[isNorm]);
	      else       val += 1e16*Dstaunu.Compute(q2,1, ml[isNorm]);
	      q2 += dq2;
	    }
	    hQ2[pad][isNorm]->SetBinContent(bin,val); hQ2[pad][isNorm]->SetBinError(bin,0);
	  }
	  hEff[pad][isNorm] = (TH1D*)hxini[pad][isNorm]->Clone(); 
	  hEff[pad][isNorm]->Divide(hQ2[pad][isNorm]);
	  if(isNorm==2) scaleEff = 1/hEff[pad][isNorm]->GetMaximum();
	  if(isNorm==1) scaleEff = hEff[pad][2]->Integral(nBins/3,nBins/3*2)/
	    hEff[pad][1]->Integral(nBins/3,nBins/3*2);
	  hEff[pad][isNorm]->Scale(scaleEff);
	  hEff[pad][isNorm]->SetMinimum(0);
	  if(row==1) hEff[pad][isNorm]->SetMaximum(1.72);
	}
	//cout<<pad<<" Drawing"<<endl;
//.........这里部分代码省略.........

void NuclearModificationFactor(TString inputPP="ROOTfiles/CrossSectionPP.root", TString inputPbPb="ROOTfiles/CrossSectionPbPb.root",TString label="PbPb",TString outputfile="RAAfile.root", Float_t centMin=0., Float_t centMax=100.)
{
//  drawDRAA = true;
//  drawChHad = true;
//  drawThm = true;

  float pti = ptBins[0]-2.;
  float pte = ptBins[nBins]+5.;
  if(drawDRAA){
    pti = 1;
    pte = 400.;
  }
  if(drawThm){
    pti = 5;
    pte = 110.;
  }


  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetEndErrorSize(0);
  gStyle->SetMarkerStyle(20);
  gStyle->SetPadRightMargin(0.03);//###0.020
  gStyle->SetPadLeftMargin(0.12);
  gStyle->SetPadTopMargin(0.075);
  gStyle->SetPadBottomMargin(0.12);

  TFile *fpp=new TFile(inputPP.Data());
  TFile *fPbPb=new TFile(inputPbPb.Data());

  TH1D*hSigmaPPStat=(TH1D*)fpp->Get("hPtSigma");
  hSigmaPPStat->SetName("hSigmaPPStat");
  TH1D*hNuclearModification=(TH1D*)fPbPb->Get("hPtSigma");
  hNuclearModification->SetName("hNuclearModification");
  hNuclearModification->Divide(hSigmaPPStat);

  double apt[nBins];
  //bin half width
  double aptl[nBins];
  //number of every rebined bin
  double bin_num[nBins];

  for (int ibin=0; ibin<nBins; ibin++){
    apt[ibin]=(ptBins[ibin+1]+ptBins[ibin])/2.;
    aptl[ibin] = (ptBins[ibin+1]-ptBins[ibin])/2;
    bin_num[ibin]=aptl[ibin]/binsize*2;
  }

  Double_t xr[nBins], yr[nBins], xrlow[nBins], yrlow[nBins],xrhigh[nBins],yrhigh[nBins];
  for(int i=0;i<nBins;i++)
  {
    yr[i] = hNuclearModification->GetBinContent(i+1);
    double systematic=0.01*systematicsForRAA(hNuclearModification->GetBinCenter(i+1),centMin,centMax,0.,0.);
	//double systematic=0.;
    yrlow[i] = hNuclearModification->GetBinContent(i+1)*systematic;
    yrhigh[i] =hNuclearModification->GetBinContent(i+1)*systematic;
  }

  TGraphAsymmErrors* gNuclearModification = new TGraphAsymmErrors(nBins,apt,yr,aptl,aptl,yrlow,yrhigh);
  gNuclearModification->SetName("gNuclearModification");

  TCanvas*canvasRAA=new TCanvas("canvasRAA","canvasRAA",600,600);
  canvasRAA->cd();
  canvasRAA->SetLogx();

  TH2F* hemptyEff=new TH2F("hemptyEff","",50,pti,pte,10.,0,1.55);  
  hemptyEff->GetXaxis()->CenterTitle();
  hemptyEff->GetYaxis()->CenterTitle();
  hemptyEff->GetYaxis()->SetTitle("R_{AA}");
  hemptyEff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
  hemptyEff->GetXaxis()->SetTitleOffset(1.0);
  hemptyEff->GetYaxis()->SetTitleOffset(1.1);
  hemptyEff->GetXaxis()->SetTitleSize(0.05);
  hemptyEff->GetYaxis()->SetTitleSize(0.05);
  hemptyEff->GetXaxis()->SetTitleFont(42);
  hemptyEff->GetYaxis()->SetTitleFont(42);
  hemptyEff->GetXaxis()->SetLabelFont(42);
  hemptyEff->GetYaxis()->SetLabelFont(42);
  hemptyEff->GetXaxis()->SetLabelSize(0.04);
  hemptyEff->GetYaxis()->SetLabelSize(0.04);  
  hemptyEff->SetMaximum(2);
  hemptyEff->SetMinimum(0.);
  hemptyEff->Draw();

  TLine *line = new TLine(pti,1,pte,1);
  line->SetLineStyle(2);
  line->SetLineWidth(2);
  line->Draw();

  //gNuclearModification->SetFillColor(kOrange);//1
  //gNuclearModification->SetFillColor(kPink+7);//1
  gNuclearModification->SetFillColor(kAzure+7);//1
  gNuclearModification->SetFillStyle(3001);//0
  gNuclearModification->SetLineWidth(1);//3
  gNuclearModification->SetMarkerSize(1);
  gNuclearModification->SetMarkerStyle(21);
  //gNuclearModification->SetLineColor(kOrange);//kGreen+4
  //gNuclearModification->SetMarkerColor(kRed);//kGreen+4
  gNuclearModification->SetLineColor(kAzure+7);//kGreen+4
  gNuclearModification->SetMarkerColor(kViolet+3);//kGreen+4
//.........这里部分代码省略.........

//.........这里部分代码省略.........

  // Do the work
  //----------------------------------------------------------------------------
  TGraphErrors* gStat = new TGraphErrors(nChannel);
  TGraphErrors* gSyst = new TGraphErrors(nChannel);
  TGraphErrors* gLumi = new TGraphErrors(nChannel);

  for (UInt_t i=0; i<nChannel; i++)
    {
      Double_t errorSquared = (xs_stat[theCharge][i] * xs_stat[theCharge][i]);

      gStat->SetPointError(i, sqrt(errorSquared) / xs_nlo[theCharge], 0.0);
      
      errorSquared += (xs_syst[theCharge][i] * xs_syst[theCharge][i]);

      gSyst->SetPointError(i, sqrt(errorSquared) / xs_nlo[theCharge], 0.0);

      errorSquared += (xs_lumi[theCharge][i] * xs_lumi[theCharge][i]);

      gLumi->SetPointError(i, sqrt(errorSquared) / xs_nlo[theCharge], 0.0);

      gStat->SetPoint(i, xs_value[theCharge][i] / xs_nlo[theCharge], nChannel-i-1);
      gSyst->SetPoint(i, xs_value[theCharge][i] / xs_nlo[theCharge], nChannel-i-1);
      gLumi->SetPoint(i, xs_value[theCharge][i] / xs_nlo[theCharge], nChannel-i-1);
    }


  // Cosmetics
  //----------------------------------------------------------------------------
  gStat->SetLineWidth  (2);
  gStat->SetMarkerSize (1.3);
  gStat->SetMarkerStyle(kFullCircle);

  gSyst->SetLineColor  (kRed);
  gSyst->SetLineWidth  (2);
  gSyst->SetMarkerSize (1.3);
  gSyst->SetMarkerStyle(kFullCircle);

  gLumi->SetLineColor  (kBlue);
  gLumi->SetLineWidth  (2);
  gLumi->SetMarkerSize (1.3);
  gLumi->SetMarkerStyle(kFullCircle);


  // Draw
  //----------------------------------------------------------------------------
  TString suffix = "7TeV_" + sCharge[theCharge];

  TCanvas* canvas = new TCanvas("ratioNLO_" + suffix, "ratioNLO_" + suffix);

  canvas->SetLeftMargin(canvas->GetRightMargin());

  Double_t xmin    =  0.0;
  Double_t xmax    =  2.0;
  Double_t ylegend =  1.2;
  Double_t ymin    = -0.6;
  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,