C++ TH1F::SetYTitle方法代码示例

void gerrors2() {
   TCanvas *c1 = new TCanvas("c1","gerrors2",200,10,700,500);
   c1->SetFillColor(42);
   c1->SetGrid();

      // draw a frame to define the range
   TH1F *hr = c1->DrawFrame(-0.4,0,1.2,12);
   hr->SetXTitle("X title");
   hr->SetYTitle("Y title");
   c1->GetFrame()->SetFillColor(21);
   c1->GetFrame()->SetBorderSize(12);

      // create first graph
   const Int_t n1 = 10;
   Double_t x1[]  = {-0.22, 0.05, 0.25, 0.35, 0.5, 0.61,0.7,0.85,0.89,0.95};
   Double_t y1[]  = {1,2.9,5.6,7.4,9,9.6,8.7,6.3,4.5,1};
   Double_t ex1[] = {.05,.1,.07,.07,.04,.05,.06,.07,.08,.05};
   Double_t ey1[] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8};
   TGraphErrors *gr1 = new TGraphErrors(n1,x1,y1,ex1,ey1);
   gr1->SetMarkerColor(kBlue);
   gr1->SetMarkerStyle(21);
   gr1->Draw("LP");

      // create second graph
   const Int_t n2 = 10;
   Float_t x2[]  = {-0.28, 0.005, 0.19, 0.29, 0.45, 0.56,0.65,0.80,0.90,1.01};
   Float_t y2[]  = {0.82,3.86,7,9,10,10.55,9.64,7.26,5.42,2};
   Float_t ex2[] = {.04,.12,.08,.06,.05,.04,.07,.06,.08,.04};
   Float_t ey2[] = {.6,.8,.7,.4,.3,.3,.4,.5,.6,.7};
   TGraphErrors *gr2 = new TGraphErrors(n2,x2,y2,ex2,ey2);
   gr2->SetMarkerColor(kRed);
   gr2->SetMarkerStyle(20);
   gr2->Draw("LP");
}

void read25cm() {
    
    TFile *f = new TFile("read25c.root");    //root file name
    TTree *ntuple = (TTree*)f->Get("ntuple");
    
    Float_t dRm,dRmbar,dRb,dRbbar;  //Values
    
    ntuple->SetBranchAddress("dRm",&dRm);   //Values
    ntuple->SetBranchAddress("dRmbar",&dRmbar);
    ntuple->SetBranchAddress("dRb",&dRb);
    ntuple->SetBranchAddress("dRbbar",&dRbbar);
    
    //create histogram
    TH1F *mudR     = new TH1F("dRm","dR distribution of mu and jet",100,0.,2 * M_PI);
    TH1F *mubardR   = new TH1F("dRmbar","dR distribution of mubar and jet",100,0.,2 * M_PI);
    TH1F *bdR     = new TH1F("dRb","dR distribution of b and jet",100,0.,2 * M_PI);
    TH1F *bbardR     = new TH1F("dRbbar","dR distribution of bbar and jet",100,0.,2 * M_PI);
    
    //read all entries and fill the histograms
    Long64_t nentries = ntuple->GetEntries();
    for (Long64_t q=0;q<nentries;q++) {
        ntuple->GetEntry(q);
        mudR->Fill(dRm);        //Histogram value
        mubardR->Fill(dRmbar);
        bdR->Fill(dRb);
        bbardR->Fill(dRbbar);
        }
    
    {
        mudR->SetXTitle("dR");
        mudR->SetYTitle("Event");
        mubardR->SetXTitle("dR");
        mubardR->SetYTitle("Event");
        bdR->SetXTitle("dR");
        bdR->SetYTitle("Event");
        bbardR->SetXTitle("dR");
        bbardR->SetYTitle("Event");
    }
    
    if (gROOT->IsBatch()) return;
    f->Write();
    new TBrowser();
    ntuple->StartViewer();
    
}

TCanvas * plot (TH1F* histoDataIn, TString legendData, TH1F* histoSimulationIn, TString legendSimulation,
                TString & canvasName, Float_t maximum = 0.15, TString xAxisTitle = "#eta",
                TString yAxisTitle = "Number of Clusters", TString error = "", bool useLegend = true,
                TString text = "", Float_t textX = 0.7, Float_t textY = 0.4, Float_t rebin = 0 ) {

  TH1F * histoData = (TH1F*)histoDataIn->Clone();
  TH1F * histoSimulation = (TH1F*)histoSimulationIn->Clone();

  // histoData->Sumw2();
  histoData->Scale(1/(histoData->Integral()));
  histoSimulation->Scale(1/(histoSimulation->Integral()));

  // Create also the legend and add the histograms
  TLegend * legend = new TLegend( 0.55, 0.65, 0.76, 0.82 );
  legend->AddEntry( histoData, xAxisTitle );
  legend->AddEntry( histoSimulation, yAxisTitle, "F" );

  cout << "histoData = " << histoData << endl;
  cout << "histoSimulation = " << histoSimulation << endl;

  TCanvas * c = new TCanvas( canvasName, canvasName, 1000, 800 );
  c->Draw();

  histoSimulation->SetMaximum(maximum);

  histoSimulation->SetFillColor(kRed);
  // histoSimulation->SetLineWidth(0);
  histoSimulation->SetLineColor(kRed);
  histoSimulation->SetXTitle(xAxisTitle);
  histoSimulation->SetYTitle(yAxisTitle);
  histoSimulation->SetTitleOffset(1.6,"Y");
  histoSimulation->SetTitle();

  histoData->SetLineStyle(1);
  histoData->SetLineWidth(2.5);

  histoSimulation->Draw();
  histoData->Draw("same");

  legend->SetFillColor(kWhite);
  if (useLegend) legend->Draw("same");

  if ( text != "" ) {
    TPaveText * pt = new TPaveText(textX, textY, textX+0.2, textY+0.17, "NDC" ); // "NDC" option sets coords relative to pad dimensions
    pt->SetFillColor(0); // text is black on white
    pt->SetTextSize(0.08); 
    pt->SetBorderSize(0);
    pt->SetTextAlign(12);
    pt->AddText(text);
    pt->Draw("same");       //to draw your text object
  }

  return c;
};

void DrawSmooth(Int_t pad, const char *title, const char *xt, const char *yt)
{
   vC1->cd(pad);
   TH1F *vFrame = gPad->DrawFrame(0,-130,60,70);
   vFrame->SetTitle(title);
   vFrame->SetTitleSize(0.2);
   vFrame->SetXTitle(xt);
   vFrame->SetYTitle(yt);
   grin->Draw("P");
   grout->DrawClone("LPX");
}

//================================================
void toyMC_Cent()
{
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1);
  gStyle->SetStatY(0.9);                
  gStyle->SetStatX(0.9);  
  gStyle->SetStatW(0.2);                
  gStyle->SetStatH(0.2);

  myRandom = new TRandom3();
  TDatime *clock = new TDatime();
  myRandom->SetSeed(clock->GetTime());


  // Get true refMult vs Ncoll distribution
  TH2F *hRefMultVsNcoll = 0x0;

  // Get different efficiency curve vs. refMult
  const int nEff = 2;
  const int nMult = 2;
  TF1 *funcEff[nEff];
  funcEff[0] = new TF1("funcEff_0", "pol1",0,1000);
  funcEff[0]->SetParameters(0.9, -2e-4);
  funcEff[1] = new TF1("funcEff_1", "pol0",0,1000);
  funcEff[1]->SetParameter(0,0.6);
  TCanvas *c = new TCanvas("cEff","cEff",800,600);
  TH1F *hplot = new TH1F("hplot",";Mult;",100,0,1000);
  hplot->SetYTitle("Efficiency");
  hplot->DrawCopy();
  for(int i=0; i<2; i++)
    {
      funcEff[i]->SetLineColor(i+1);
      funcEff[i]->SetLineWidth(1.5);
      funcEff[i]->Draw("sames");
    }
  TPaveText *title = GetTitleText("TPC tracking efficiency");
  title->Draw();
  TLegend *leg = new TLegend(0.4,0.2,0.7,0.4);
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->SetTextSize(0.04);
  leg->AddEntry(funcEff[0], "Efficiency 1", "L");
  leg->AddEntry(funcEff[1], "Efficiency 2", "L");
  leg->Draw();

  // Get ngTrack distribution to mimic efficiency distribution
  
  TH1F *hMcMult[nMult];
  
}

void binomialEfficiency1D(TH1F * numerator,TH1F * denominator){
  TH1F * efficiency = numerator->Clone("efficiency");
  efficiency->SetXTitle(numerator->GetXaxis()->GetTitle());
  efficiency->SetYTitle("#epsilon");
  for(int j=0;j<=numerator->GetXaxis()->GetNbins() ;j++){
    if(denominator->GetBinContent(j)!=0){
      float eff = numerator->GetBinContent(j)/denominator->GetBinContent(j);
      float err = sqrt(eff*(1-eff)/denominator->GetBinContent(j));
      efficiency->SetBinContent(j,eff);
      efficiency->SetBinError(j,err);
      cout<<"1Deff "<<j<<" "<<eff<<" +/- "<<err<<endl;
    }
  }
  efficiency->Draw("E");
  efficiency->SetMarkerColor(kRed);
  efficiency->SetMarkerStyle(23);
  efficiency->SetMarkerSize(2);
}

void DrawSmooth(Int_t pad, const char *title, const char *xt,
   const char *yt)
{
   vC1->cd(pad);
   TH1F *vFrame = gPad->DrawFrame(0,0,15,150);
   vFrame->SetTitle(title);
   vFrame->SetTitleSize(0.2);
   vFrame->SetXTitle(xt);
   vFrame->SetYTitle(yt);
   grxy->SetMarkerColor(kBlue);
   grxy->SetMarkerStyle(21);
   grxy->SetMarkerSize(0.5);
   grxy->Draw("P");
   grin->SetMarkerColor(kRed);
   grin->SetMarkerStyle(5);
   grin->SetMarkerSize(0.7);
   grin->Draw("P");
   grout->DrawClone("LP");
}

TH1F * DrawOverflowBin(TH1F * h)
{
  // This function paint the histogram h with an extra bin for overflows                                                                               
  
  UInt_t nx = h->GetNbinsX()+1;
  Double_t *xbins= new Double_t[nx+1];
  for (UInt_t i=0;i<nx;i++){
    xbins[i]=h->GetBinLowEdge(i+1);
  }
  xbins[nx]=xbins[nx-1]+h->GetBinWidth(nx);
  
  //char *tempName= new char[strlen(h->GetName())+10];
  //sprintf(tempName,"%swtOverFlow",h->GetName());
  
  stringstream tempName;
  tempName << h->GetName() << "wtOverFlow";
  // Book a temporary histogram having ab extra bin for overflows                                                                                      
  TH1F *htmp = new TH1F(tempName.str().c_str(), h->GetTitle(), nx, xbins);
  
  // Reset the axis labels                                                                                                                             
  htmp->SetXTitle(h->GetXaxis()->GetTitle());
  htmp->SetYTitle(h->GetYaxis()->GetTitle());
  
  // Fill the new hitogram including the extra bin for overflows                                                                                       
  for (UInt_t i=1; i<=nx; i++){
    htmp->Fill(htmp->GetBinCenter(i), h->GetBinContent(i));
    htmp->SetBinError(i, h->GetBinError(i));
  }
  //htmp->Draw();
  // Fill the underflows                                                                                                                               
  //htmp->Fill(h->GetBinLowEdge(1)-1, h->GetBinContent(0));                                                                                            
  // Restore the number of entries                                                                                                                       
 
  htmp->SetEntries(h->GetEntries());
   /*
  // FillStyle and color                                                                                                                               
  htmp->SetFillStyle(h->GetFillStyle());
  htmp->SetFillColor(h->GetFillColor());
  */
  return htmp;                          
  //return h;
}

//.........这里部分代码省略.........
//     grthSM->SetPoint(19, 1750, 4.1121E-05);
//     grthSM->SetPoint(20, 1800, 3.2378E-05);
//     grthSM->SetPoint(21, 1850, 2.5507E-05);
//     grthSM->SetPoint(22, 1900, 2.0215E-05);
//     grthSM->SetPoint(23, 1950, 1.6020E-05);
//     grthSM->SetPoint(24, 2000, 1.2714E-05);
//     grthSM->SetPoint(25, 2050, 1.0133E-05);
//     grthSM->SetPoint(26, 2100, 8.0785E-06);
//     grthSM->SetPoint(27, 2150, 6.4583E-06);
//     grthSM->SetPoint(28, 2200, 5.1774E-06);
//     grthSM->SetPoint(29, 2250, 4.1620E-06);
//     grthSM->SetPoint(30, 2300, 3.3440E-06);
//     grthSM->SetPoint(31, 2350, 2.7018E-06);
//     grthSM->SetPoint(32, 2400, 2.1753E-06);
//     grthSM->SetPoint(33, 2450, 1.7626E-06);
//     grthSM->SetPoint(34, 2500, 1.4225E-06);
//   }
  grthSM10->SetName("SMXSection_2nd");

  // double fr_left = 590.0, fr_down = 1E-5, fr_right = 2000.0, fr_up = 0.5; 
//   double fr_left = 590.0, fr_down = 5E-5, fr_right = 2000.0, fr_up = 5;
   double fr_left = 500.0, fr_down = 5E-7, fr_right = 4500.0, fr_up = 5E-1;

  TCanvas *cMCMC = new TCanvas("c_lim_Asymptotic", "canvas with limits for Asymptotic CLs", 630, 600);
  cMCMC->cd();
  cMCMC->SetGridx(1);
  cMCMC->SetGridy(1);
  // draw a frame to define the range

  TH1F *hr = cMCMC->DrawFrame(fr_left, fr_down, fr_right, fr_up, "");
  TString VV = "ZH";
  
  hr->SetXTitle("M_{X} [GeV]");
  hr->SetYTitle("#sigma_{95%} [pb]"); // #rightarrow 2l2q
  

  gr95_cls->SetFillColor(kYellow);
  gr95_cls->SetFillStyle(1001);//solid
  gr95_cls->SetLineStyle(kDashed);
  gr95_cls->SetLineWidth(3);
  gr95_cls->GetXaxis()->SetTitle("M_{V'} [GeV]");
  gr95_cls->GetYaxis()->SetTitle("#sigma_{95%} #times BR(V' #rightarrow " + VV + ") [pb]"); // #rightarrow 2l2q
  gr95_cls->GetXaxis()->SetRangeUser(fr_left, fr_right);

//  gr95_cls->Draw("3");

  gr68_cls->SetFillColor(kGreen);
  gr68_cls->SetFillStyle(1001);//solid
  gr68_cls->SetLineStyle(kDashed);
  gr68_cls->SetLineWidth(3);
//  gr68_cls->Draw("3same");
  grmedian_cls->GetXaxis()->SetTitle("M_{V'} [GeV]");
  grmedian_cls->GetYaxis()->SetTitle("#sigma_{95%} #times BR(V' #rightarrow " + VV + ") [pb]"); // #rightarrow 2l2q
  grmedian_cls->SetMarkerStyle(24);//25=hollow squre
  grmedian_cls->SetMarkerColor(kBlack);
  grmedian_cls->SetLineStyle(2);
  grmedian_cls->SetLineWidth(3);
  grmedian_cls->SetMinimum(0.0);
  grmedian_cls->SetMaximum(8.0);

  grobslim_cls->SetMarkerColor(kBlack);
  grobslim_cls->SetMarkerStyle(21);//24=hollow circle
  grobslim_cls->SetMarkerSize(1.0);
  grobslim_cls->SetLineStyle(1);
  grobslim_cls->SetLineWidth(3);

void unfoldPt(int mode=0)
{

   // Matched Tracklets
   TFile *inf = new TFile("match-10TeV-12.root");
   TNtuple *nt = (TNtuple*)inf->FindObjectAny("nt");

   // Test sample
   TFile *infTest = new TFile("./TrackletTree-Run123596.root");
   TNtuple *ntTest = (TNtuple*)infTest->FindObjectAny("TrackletTree12");
   
   TFile *pdfFile;
   if (mode==0) pdfFile = new TFile("pdf.root","recreate");
           else pdfFile = new TFile("pdf.root");

   double nPtBin=15;
   double minPt=log(0.05);
   double maxPt=log(10);
   double nDphiBin=600;
   double maxDphi=0.4;
   
   char* mycut = Form("abs(eta)<2&&log(pt)>%f&&log(pt)<%f",minPt,maxPt);
   char* mycut1=Form("abs(eta)<2&&log(pt)>%f&&log(pt)<%f&&abs(eta-eta1)<0.01&&abs(deta)<0.01",minPt,maxPt);

   TH2F *h;
   TH1F *hdphi = new TH1F("hdphi","",nDphiBin,0,maxDphi);
   TH1F *hdphi2;
   TH1F *hpt;
   TH1F *hptH = new TH1F("hptH","",nPtBin,minPt,maxPt);
   
   TH1F *hptUnfold = new TH1F("hptUnfold","",nPtBin,minPt,maxPt);
   TH1F *hptMC = new TH1F("hptMC","",nPtBin,minPt,maxPt);
   TH1F *hptTemp = new TH1F("hptTemp","",nPtBin,minPt,maxPt);

   // Delta phi as a function of matched genparticle transverse momentum
   TCanvas *c = new TCanvas("c","",600,600);
   
   if (mode == 0) {
      h = new TH2F("h","",nPtBin,minPt,maxPt,nDphiBin,0,maxDphi);
      hdphi2 = new TH1F("hdphiMC","",nDphiBin,0,maxDphi);
      hpt = new TH1F("hpt","",nPtBin,minPt,maxPt);      
      h->SetXTitle("ln(P_{T}) GeV/c");
      h->SetYTitle("|#Delta#phi|");
      nt->Draw("abs(dphi):log(pt)>>h",mycut1,"col");
      // used to generate pdf
      nt->Draw("abs(dphi)>>hdphiMC",mycut,"");
      nt->Draw("log(pt)>>hpt",mycut,"");
      h->Write();      
      hpt->Write();      
      hdphi2->Write();      
   } else {
      h = (TH2F*) pdfFile->FindObjectAny("h");
      hdphi2 = (TH1F*) pdfFile->FindObjectAny("hdphiMC");
      hpt = (TH1F*) pdfFile->FindObjectAny("hpt");
   }
   // Delta phi fit
   TCanvas *c2 = new TCanvas("c2","",600,600);
   c2->SetLogy();
   c2->SetLogx();

 
   // dphi for unfolding and MC truth:  
   ntTest->Draw("abs(dphi)>>hdphi","abs(eta1)<2&&abs(deta)<0.1","",200000);
   ntTest->Draw("log(pt)>>hptH",mycut,"",200000);
   
   histFunction2D *myfun = new histFunction2D(h);
   
   TF1 *test = new TF1("histFun",myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);
   TF1 *test2 = new TF1("histFunMC",myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);

   for (int i=0;i<nPtBin+1;i++)
   {  
      test->SetParameter(i,1);   
   }


   hdphi2->SetXTitle("|#Delta#phi|");
   hdphi2->SetYTitle("Arbitrary Normalization");
   hdphi2->Fit("histFunMC","M");

   hdphi->SetXTitle("|#Delta#phi|");
   hdphi->SetYTitle("Arbitrary Normalization");
   hdphi->Fit("histFun","M");
   hdphi->SetStats(0);
   hdphi->Draw();

   
   for (int i=0;i<nPtBin+1;i++) {
      TF1 *testPlot = new TF1(Form("histFun%d",i),myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);

      testPlot->SetParameter(i,test->GetParameter(i));
      testPlot->SetLineColor(i+2);
      testPlot->Draw("same");
   }
   
   int total=0,totalMC=0;


   for (int i=0;i<nPtBin;i++){
      if (test->GetParameter(i)==0) continue;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  fit_status = hsig->Fit(f1,"","",-1., 5.);

  hsig->Draw();
  f1->Draw("same");
  if ( fit_status > 0 ) {
     printf("fit signal template failed. QUIT \n");
     return fitted;
  }
  if(para_index>0 && para_index<4){
    double tmppar = f1->GetParameter(para_index);
    f1->SetParameter(para_index, tmppar+para_sigma*f1->GetParError(para_index));
  }

  TF1 *fmcsig = (TF1*)f1->Clone();
  TF1 *fmcsigcorr = (TF1*)f1->Clone();
  fmcsig->SetNpx(10000);
  fmcsigcorr->SetNpx(10000);
  fmcsigfit->SetNpx(10000);
  
  TCanvas *c101 = new TCanvas("c101","c101",1000,500);
  c101->Divide(2,1);
  c101->cd(1);

  fmcsig->SetLineColor(1);
//   fmcsig->Draw();
//   f1->Draw("same");
  TH1F *htmp1 = (TH1F*)fmcsig->GetHistogram();
//   TH1F *htmp2 = (TH1F*)fmcsigcorr->GetHistogram();
  
  TH2F *htmp2 = new TH2F("htmp2","",210, -1., 20., 100, 0., htmp1->GetMaximum()*1.25);

  htmp2->SetNdivisions(505,"XY");
  htmp2->SetXTitle("Iso");
  htmp2->SetYTitle("A.U.");
  htmp2->SetLineColor(1);

//   htmp2->Draw();
//   htmp1->Draw("same");
//   htmp2->Add(htmp1,-1);
//   htmp2->Divide(htmp1);
  htmp2->GetXaxis()->SetRangeUser(-1., 10.);
  htmp2->SetMinimum(-1.);
  //htmp2->SetMaximum(1.5);
  htmp2->Draw();
  fmcsig->Draw("same");
//   fmcsigcorr->Draw("same");
  
  TLegend *tleg1 = new TLegend(0.5, 0.7, 0.93, 0.92);
  tleg1->SetHeader("");
  tleg1->SetFillColor(0);
  tleg1->SetShadowColor(0);
  tleg1->SetBorderSize(0);
  tleg1->AddEntry(fmcsig,"Zee data","l");
  //tleg1->AddEntry(fmcsigcorr,"corrected shape","l");
  tleg1->AddEntry(fmcsigfit,"shape from data","l");
  tleg1->Draw();

  //return fitted;
       
  SigPDFnorm = f1->Integral(-1.,20.);
  printf("status %d, sig area %3.3f \n", fit_status,f1->Integral(-1., 20.));


  f1->SetParameter(2,f1->GetParameter(2)+0.2);
  f1->SetParameter(3,f1->GetParameter(3)+0.1);

//.........这里部分代码省略.........
	} else {
	  x3[np3] = radiusH[j]*cos(phiH[j]);
	  y3[np3] = radiusH[j]*sin(phiH[j]);
	  np3++;
	}
      }
    }
    //    std::cout << "i " << i << " " <<np1 << " " <<np2 << " " <<np3 <<"\n";
    if (np1 > 0) {
      gr[0][i] = new TGraph(np1,x1,y1); gr[0][i]->SetTitle(""); 
      gr[0][i]->SetMarkerStyle(symbol[i]);  gr[0][i]->SetMarkerColor(colr[0]);
    } else 
      gr[0][i] = 0;
    if (np2 > 0) {
      gr[1][i] = new TGraph(np2,x2,y2); gr[1][i]->SetTitle(""); 
      gr[1][i]->SetMarkerStyle(symbol[i]);  gr[1][i]->SetMarkerColor(colr[1]);
    } else 
      gr[1][i] = 0;
    if (np3 > 0) {
      gr[2][i] = new TGraph(np3,x3,y3); gr[2][i]->SetTitle(""); 
      gr[2][i]->SetMarkerStyle(symbol[i]);  gr[2][i]->SetMarkerColor(colr[2]);
    } else 
      gr[2][i] = 0;
  }
  np1 = np2 = np3 = 0;
  for (int j=0; j<np; j++) {
    if (etaIndex[j] < 29 || etaIndex[j] > 41) {
      int k = 2;
      if (depths[j] == 3) k = 0;
      else if (depths[j] == 4) k = 1;
      if (k == 0) {
	x1[np1] = radiusH[j]*cos(phiH[j]);
	y1[np1] = radiusH[j]*sin(phiH[j]);
	np1++;
	if (np1 == 1) std::cout << i << " 0 " <<x1[0] << " " <<y1[0] << "\n";
      } else if (k==1) {
	x2[np2] = radiusH[j]*cos(phiH[j]);
	y2[np2] = radiusH[j]*sin(phiH[j]);
	np2++;
	if (np2 == 1) std::cout << i << " 1 " <<x2[0] << " " <<y2[0] << "\n";
      } else {
	x3[np3] = radiusH[j]*cos(phiH[j]);
	y3[np3] = radiusH[j]*sin(phiH[j]);
	np3++;
      }
    }
  }
  //    std::cout << "i " << i << " " <<np1 << " " <<np2 << " " <<np3 <<"\n";
  if (np1 > 0) {
    gr[0][13] = new TGraph(np1,x1,y1); gr[0][13]->SetTitle(""); 
    gr[0][13]->SetMarkerStyle(symbol[13]);  gr[0][13]->SetMarkerColor(colr[0]);
  } else 
    gr[0][13] = 0;
  if (np2 > 0) {
    gr[1][13] = new TGraph(np2,x2,y2); gr[1][13]->SetTitle(""); 
    gr[1][13]->SetMarkerStyle(symbol[13]);  gr[1][13]->SetMarkerColor(colr[1]);
  } else 
    gr[1][13] = 0;
  if (np3 > 0) {
    gr[2][13] = new TGraph(np3,x3,y3); gr[2][13]->SetTitle(""); 
    gr[2][13]->SetMarkerStyle(symbol[13]);  gr[2][13]->SetMarkerColor(colr[2]);
  } else 
    gr[2][13] = 0;

  TCanvas *c0  = new TCanvas("c0","PMT Hits",800,600); 
  TH1F *vFrame = c0->DrawFrame(1000.0,-250.0,1500.0,250.0);
  vFrame->SetXTitle("x (mm)");
  vFrame->SetYTitle("y (mm)");
  for (int i=0; i<=13; i++) {
    for (int j=0; j<3; j++) {
      if (gr[j][i] != 0) {
	gr[j][i]->Draw("p");
	gr[j][i]->SetLineColor(colr[j]); gr[j][i]->SetLineWidth(2);
	//	std::cout << "Next " << i << " " << j << "\n";
      }
    }
  }
  TLegend *leg1 = new TLegend(0.75,0.55,0.90,0.90);
  char list[40];
  for (i=0; i<= 13; i++) {
    if (i < 13) sprintf (list, "#eta = %d", i+29);
    else        sprintf (list, "Unknown #eta");
    if      (gr[0][i] != 0) leg1->AddEntry(gr[0][i], list, "P");
    else if (gr[1][i] != 0) leg1->AddEntry(gr[1][i], list, "P");
  }
  for (i=0; i<2; i++) {
    if (i == 0) sprintf (list, "Long Fibre");
    else        sprintf (list, "Short Fibre");
    if      (gr[i][0] != 0) leg1->AddEntry(gr[i][0], list, "L");
    else if (gr[i][1] != 0) leg1->AddEntry(gr[i][1], list, "L");
    else if (gr[i][2] != 0) leg1->AddEntry(gr[i][2], list, "L");
  }
  leg1->SetFillColor(0); leg1->SetTextSize(0.03); 
  leg1->SetBorderSize(1); leg1->Draw();

  if (save != 0) {
    if (save > 0) c0->SaveAs("PMTHits.eps");
    else          c0->SaveAs("PMTHits.gif");
  }
}

//.........这里部分代码省略.........
  if(cut_value[i]==515 || cut_value[i]==535 || cut_value[i]==555) {nXm = 10;}// lost M800
  if(cut_value[i]==575 || cut_value[i]==595 || cut_value[i]==615) {nXm = 10;}// lost M800
  if(cut_value[i]==635 ) {nXm = 10;}// lost M800
  if(cut_value[i]==835) {nXm = 8;}// lost M800, M1000, M1200


  cout<<endl;
  cout<<"dir_path[i]: "<< dir_path[i] <<" cut_value[i]: "<< cut_value[i]<<" nXm: "<< nXm <<endl;
  cout<<endl;


  grmedian_cls_new[i] =  plot_Asymptotic( dir_path[i] ,grmedian_cls, nXm );

//  leg_name[i] = Form("CL_{S} Expected Limit for cut pt>%d",cut_value[i] ) ;
  leg_name[i] = Form("cut pt>%d",cut_value[i] ) ;
  }





  TCanvas *c1 = new TCanvas("c1", "canvas with limits for Asymptotic CLs", 630, 600);


  TLegend *leg = new TLegend(.75, .35, 1., .85);
  leg->SetFillColor(0);
  leg->SetShadowColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);

  c1->SetGridx(1);
  c1->SetGridy(1);
    gPad->SetLogy();


  
  c1->cd();

//   double fr_left = 500.0, fr_down = 5E-2, fr_right = 4500.0, fr_up = 3E2;
//   double fr_left = 500.0, fr_down = 2E-1, fr_right = 1500.0, fr_up = 7E-1;
//   double fr_left = 1800.0, fr_down = 1.1E-1, fr_right = 2200.0, fr_up = 1.5E-1;
   double fr_left = 2800.0, fr_down = 1E-1, fr_right = 3200.0, fr_up = 1.1E-1;
  TH1F *hr = c1->DrawFrame(fr_left, fr_down, fr_right, fr_up, "");

  int colorN = 1;

// loop

  for(int i=0;i<n_cut_index;i++){

        if (i==0){colorN = kRed;}
        if (i==2){colorN = kOrange;}
        if (i==4){colorN = kCyan;}
        if (i==6){colorN = kBlue;}
        if (i==8){colorN = kGreen;}
	
  	grmedian_cls_new[i]->SetLineColor(colorN);
        grmedian_cls_new[i]->SetLineStyle(1);
  	grmedian_cls_new[i]->Draw();

//  	leg->AddEntry(grmedian_cls_new[i], "CL_{S} Expected Limit", "L");
        leg->AddEntry(grmedian_cls_new[i],leg_name[i] , "L");
  	leg->Draw();

	colorN = colorN+2;
  }
// end loop

  hr->SetXTitle("M_{Z'} [GeV]");
  hr->SetYTitle("95% #sigma_{Z'}xBR(Z'->ZH) [pb]");

    TLatex * latex = new TLatex();
    latex->SetNDC();
    latex->SetTextSize(0.04);
    latex->SetTextAlign(31);
    latex->SetTextAlign(11);

    latex->DrawLatex(0.25, 0.91, "CL_{S} Expected Limit, cut electron leading pt");
    latex->DrawLatex(0.18, 0.96, "CMS preliminary 2015");
    latex->DrawLatex(0.60, 0.96, Form("%.1f fb^{-1} at #sqrt{s} = 13 TeV", intLumi));


  TString save_name = "expected_limit_cut_";
  save_name = save_name + variable_name ;

  TString png_name = save_name + ".png";
  TString pdf_name = save_name + ".pdf";

  c1->SaveAs(png_name);
  c1->SaveAs(pdf_name);


//TString path_to_save = "/afs/cern.ch/user/y/yuchang/www/set_limit_optimize_leading_pt/" ;
TString path_to_save = "/afs/cern.ch/user/y/yuchang/www/optimize_cut_with_limit_use_eleIsPassLoose/leading_pt/" ;
pdf_name = path_to_save + pdf_name;

c1->SaveAs( pdf_name );


}// end codes

//.........这里部分代码省略.........
    // pp (average Run-2/3 parameterization)
    // weight = 2.0*pi*pt * (546.2/42.2) / pow(exp(-0.2664*pt) + pt/0.7703,8.718);
    //    weight = weight*(1+0.6362/exp(7.736*pt));

    ptpion->AddBinContent(ibin,weight);
    intpion = intpion + weight*binwidth;

    // if ( pt<1.5 ) weight = 2.0 * pi * pt * otherhagedorn;
    // if ( pt>=1.5 && pt<=5.5 ) 
    //   weight = 2.0 * pi * pt * (trans*otherhagedorn + (1.-trans)*otherpower);
    // if ( pt>5.5 ) weight = 2.0*pi*pt*otherpower;

    weight = 2.0*pi*pt * (450.8/42.2) / 
      pow(exp(-0.5553*sqrt(mtother*mtother-pimass*pimass)) + 
	  sqrt(mtother*mtother-pimass*pimass)/0.6467,8.046);

    ptothermt->AddBinContent(ibin,weight);
    intothermt = intothermt + weight*binwidth;

    weight = norm*InitializePtHydro(omass,t_fo,beta,pt);
    ptother->AddBinContent(ibin,weight);
    intother = intother + weight*binwidth;
  }

  cout << "Integral pion    (dN/dy): " << intpion << endl;
  cout << "Integral other   (dN/dy): " << intother << endl;
  cout << "Integral othermt (dN/dy): " << intothermt << endl;

  norm = intpion/intothermt;
  ptothermt->Scale(norm);
 
  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    weight = ptothermt->GetBinContent(ibin)/ptpion->GetBinContent(ibin);
    ratiomt->AddBinContent(ibin,weight);
  }

  norm = expected/ratiomt->GetMaximum();
  ratiomt->Scale(norm);
  ptothermt->Scale(norm);
  cout << "meson/pion ratio should be set to: " << norm << endl;

  norm = expected*intothermt/intother;
  ptother->Scale(norm);
 
  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    weight = ptother->GetBinContent(ibin)/ptpion->GetBinContent(ibin);
    ratio->AddBinContent(ibin,weight);
  }

  ptpion->SetXTitle("p_{t} [GeV/c]");
  ptpion->SetYTitle("(1/2#pi p_{t})dN/dp_{t}dy [(c/GeV)^{2}]");
  ptpion->SetLineWidth(5);
  ptpion->SetLineColor(1);
  ptother->SetLineWidth(5);
  ptother->SetLineColor(2);
  ptothermt->SetLineWidth(5);
  ptothermt->SetLineColor(4);
  ptpion->Draw("");
  ptothermt->Draw("same");
  ptother->Draw("same");

  TLatex *text_pion0 = new TLatex(3.2,1.5,"#pi^{0}");
  text_pion0->SetTextColor(1);
  text_pion0->Draw();
  TLatex *text_pionm = new TLatex(3.2,1.1,"#eta m_{t} scaled");
  text_pionm->SetTextColor(4);
  text_pionm->Draw();
  TLatex *text_pionp = new TLatex(3.2,0.7,"#eta hydro");
  text_pionp->SetTextColor(2);
  text_pionp->Draw();

   c2 = new TCanvas("c2","A Simple Graph with error bars",200,10,700,500);

   c2->SetFillColor(10);
   c2->SetGrid();
   c2->SetLogy(0);
   c2->SetTicks();
 
   c2->cd();

  ratiomt->SetXTitle("p_{t} [GeV/c]");
  ratiomt->SetYTitle("#eta / #pi^{0}");
  ratio->SetLineWidth(5);
  ratio->SetLineColor(2);
  ratiomt->SetLineWidth(5);
  ratiomt->SetLineColor(4);
  ratiomt->Draw();
  ratio->Draw("same");

  TLatex *text_pionm = new TLatex(3.2,0.42,"#eta m_{t} scaled");
  text_pionm->SetTextColor(4);
  text_pionm->Draw();
  TLatex *text_pionp = new TLatex(3.2,0.35,"#eta hydro");
  text_pionp->SetTextColor(2);
  text_pionp->Draw();

  return;
}

void plotHitEnergy(string inputDir, int hittype=1, float radius=0.4) {
    setNCUStyle(true);

    string decversion;
    if(inputDir.find("rfull009")!=std::string::npos)decversion="rfull009";
    else if(inputDir.find("rfull012")!=std::string::npos)decversion="rfull012";

    TH1F* hecal;
    TH1F* hhcal;
    std::string ecalhitname = hittype==1? "EM_BARREL":"EcalBarrelHits";
    std::string hcalhitname = hittype==1? "HAD_BARREL":"HcalBarrelHits";

    TString energy=gSystem->GetFromPipe(Form("file=%s; test=${file##*/}; test2=${test%%mumu*}; echo \"${test2##*tev}\"",inputDir.data()));
    cout << "energy = " << energy.Data() << endl;

    string inputFile = inputDir + "/radius" + Form("%0.1f",radius)+ (hittype==1? "_rawhit.root": "_rawhit_new.root");
    cout << "opening " << inputFile.data() << endl;



    TFile *f = TFile::Open(inputFile.data());
    hecal = (TH1F*)f->FindObjectAny("hecalhit_energy");
    hecal->SetLineWidth(3);
    hecal->SetFillStyle(1001);
    hecal->SetFillColor(4);
    hecal->SetLineColor(4);
    hecal->SetXTitle(Form("%s hit energy [GeV]",ecalhitname.data()));
    hecal->SetYTitle(Form("Number of hits per %.1f GeV",hecal->GetBinWidth(1)));
    hecal->SetTitleOffset(1.2,"X");
    hecal->SetTitleOffset(1.4,"Y");
    hhcal = (TH1F*)f->FindObjectAny("hhcalhit_energy");
    hhcal->SetLineWidth(3);
    hhcal->SetFillStyle(1001);
    hhcal->SetFillColor(2);
    hhcal->SetLineColor(2);
    hhcal->SetXTitle(Form("%s hit energy [GeV]",hcalhitname.data()));
    hhcal->SetYTitle(Form("Number of hits per %.1f GeV",hhcal->GetBinWidth(1)));
    hhcal->SetTitleOffset(1.2,"X");
    hhcal->SetTitleOffset(1.4,"Y");


    TLegend* leg = new TLegend(0.444,0.690,0.990,0.903);
    leg->SetFillColor(0);
    leg->SetFillStyle(0);

    TCanvas* c1 = new TCanvas("c1","",500,500);
    int lastbin=hecal->FindLastBinAbove(0)+20;
    float xmax=hecal->GetBinLowEdge(lastbin);
    hecal->GetXaxis()->SetRangeUser(0,xmax);
    hecal->Draw("hist");
    c1->SetLogy(1);

    leg->SetHeader(decversion.data());
    leg->AddEntry(hecal,Form("%s-TeV Z'#rightarrow qq",energy.Data()),"f");
    leg->Draw("same");

    std::string outputname = decversion + "/" + decversion + "_" + ecalhitname + "hit_energy_" + Form("%s",energy.Data()) + "TeVZp";
    c1->Print(Form("%s.pdf",outputname.data()));
    c1->Print(Form("%s.eps",outputname.data()));
    leg->Clear();


    lastbin=hhcal->FindLastBinAbove(0)+20;
    xmax=hhcal->GetBinLowEdge(lastbin);
    hhcal->GetXaxis()->SetRangeUser(0,xmax);
    hhcal->Draw("hist");
    c1->SetLogy(1);

    leg->SetHeader(decversion.data());
    leg->AddEntry(hhcal,Form("%s-TeV Z'#rightarrow qq",energy.Data()),"f");
    leg->Draw("same");

    outputname = decversion + "/" + decversion + "_" + hcalhitname + "hit_energy_" + Form("%s",energy.Data()) + "TeVZp";

    c1->Print(Form("%s.pdf",outputname.data()));
    c1->Print(Form("%s.eps",outputname.data()));

}