C++ TPad::SetLogy方法代码示例

// Modes: linear or log
plotHistogram(const char* canvas, const char* mode, const char* title,
	      double xlo, double xhi, int nbin,
	      double* sig, double* sigerr, double* bgr, double* bgrerr)
{
  TCanvas *c 
    = new TCanvas(canvas,"SPR Classifier Output",200,10,600,400);
  gStyle->SetPalette(1);
  TLegend *leg = new TLegend(0.1,0.85,0.5,1.,"Classifier Output","NDC");
  double dx = (xhi-xlo) / nbin;
  TH1D* hs = new TH1D("signal",    title,nbin,xlo,xhi);
  TH1D* hb = new TH1D("background",title,nbin,xlo,xhi);
  leg->AddEntry(hs,"Signal","L");
  leg->AddEntry(hb,"Background","L");
  for( int i=0;i<nbin;i++ ) {
    hs->SetBinContent(i+1,sig[i]);
    hs->SetBinError(i+1,sigerr[i]);
    hb->SetBinContent(i+1,bgr[i]);
    hb->SetBinError(i+1,bgrerr[i]);
  }
  TPad* pad = new TPad("pad","pad",0,0,1,1);
  if( strcmp(mode,"log") == 0 ) pad->SetLogy(1);
  pad->Draw();
  pad->cd();
  hs->SetLineColor(2);
  hs->SetLineWidth(3);
  hb->SetLineColor(4);
  hb->SetLineWidth(3);
  hb->Draw();
  hs->Draw("same");
  leg->Draw();
}

void LHCHiggsExample() 
{ 

#ifdef __CINT__
  gROOT->LoadMacro("LHCHiggsUtils.C");
#endif

  SetLHCHiggsStyle();

  TCanvas* c1 = new TCanvas("c1","Higgs Cross Section",50,50,600,600);
  TPad* thePad = (TPad*)c1->cd();
  thePad->SetLogy();

  Double_t ymin=1.e-2;  Double_t ymax=1.e2;
  Double_t xmin=90.00;  Double_t xmax=600.;
  TH1F *h1 = thePad->DrawFrame(xmin,ymin,xmax,ymax);
  h1->SetYTitle("#sigma_{pp #rightarrow H} [pb]");
  h1->SetXTitle("M_{H}  [GeV]");
  h1->GetYaxis()->SetTitleOffset(1.4);
  h1->GetXaxis()->SetTitleOffset(1.4);
  //h1->GetXaxis()->SetNdivisions(5);
  h1->Draw();

  myText(0.2,0.88,1,"#sqrt{s}= 13 TeV");
  myBoxText(0.55,0.67,0.05,5,"NNLO QCD");

  LHCHIGGS_LABEL(0.98,0.725);
  myText(0.2,0.2,1,"Preliminary");

  c1->Print("LHCHiggsExample.eps");
  c1->Print("LHCHiggsExample.png");
  c1->Print("LHCHiggsExample.pdf");

}

/** 
 * Draw the Poisson estimate of the occupancy in a given ring.
 * 
 * @param p            List 
 * @param d            Detector
 * @param r            Ring
 * 
 * @return The occupancy (in percent)
 *
 * @deprecated Use QATrender instead
 * @ingroup pwglf_forward_scripts_qa
 */
Double_t
DrawRingOccupancy(TList* p, UShort_t d, Char_t r)
{
  if (!p) return 0;

  TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
  if (!ring) { 
    Error("DrawOccupancy", "List FMD%d%c not found in %s",d,r,p->GetName());
    return 0;
  }
  
  TH1* corr = static_cast<TH1*>(ring->FindObject("occupancy"));
  if (!corr) { 
    Error("DrawRingOccupancy", "Histogram occupancy not found in FMD%d%c",
	  d, r);
    return 0;
  }
  corr->Rebin(4);

  TPad* pad = static_cast<TPad*>(gPad);
  pad->SetGridy();
  pad->SetGridx();
  pad->SetLogy();
  pad->SetFillColor(0);
    pad->SetRightMargin(0.01);
#if 0
  if (d == 3) { 
    pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99, 
		 pad->GetYlowNDC()+pad->GetHNDC());
    pad->SetRightMargin(0.15);
  }
#endif

  corr->Draw("hist");

  TLatex* ltx = new TLatex(.95, .95, Form("FMD%d%c", d, r));
  ltx->SetNDC();
  ltx->SetTextAlign(33);
  ltx->SetTextSize(.08);
  ltx->Draw();

  return corr->GetMean();
}

void EMCDistribution(TString gain = "CALIB", bool log_scale = false)
{
  TText *t;
  TCanvas *c1 = new TCanvas(
      "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
      "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
      1000);
  c1->Divide(8, 8, 0., 0.01);
  int idx = 1;
  TPad *p;

  for (int iphi = 8 - 1; iphi >= 0; iphi--)
  {
    for (int ieta = 0; ieta < 8; ieta++)
    {
      p = (TPad *) c1->cd(idx++);
      c1->Update();

      p->SetLogy();
      p->SetGridx(0);
      p->SetGridy(0);

      TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi) + TString(log_scale ? "_Log" : "");

      TH1 *h = NULL;

      if (log_scale)
        h = new TH1F(hname,
                     Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300,
                     5e-3, 3096);
      else
        //            h = new TH1F(hname,
        //                Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 196,
        //                1900, 2096);
        h = new TH1F(hname,
                     Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 596,
                     -96, 500);

      h->SetLineWidth(0);
      h->SetLineColor(kBlue + 3);
      h->SetFillColor(kBlue + 3);
      h->GetXaxis()->SetTitleSize(.09);
      h->GetXaxis()->SetLabelSize(.08);
      h->GetYaxis()->SetLabelSize(.08);

      if (log_scale)
        QAHistManagerDef::useLogBins(h->GetXaxis());

      T->Draw(
          "TOWER_" + gain + "_CEMC[].get_energy_power_law_exp()>>" + hname,
          Form(
              "TOWER_%s_CEMC[].get_bineta()==%d && TOWER_%s_CEMC[].get_binphi()==%d",
              gain.Data(), ieta, gain.Data(), iphi),
          "");

      TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
      t->SetTextAlign(33);
      t->SetTextSize(.15);
      t->SetNDC();
      t->Draw();

      //          return;
    }
  }

  SaveCanvas(c1,
             TString(_file0->GetName()) + TString("_DrawPrototype3EMCalTower_") + TString(c1->GetName()), false);
}

void EMCDistribution_PeakSample_Fast(bool full_gain = false)
{
  const TString gain = "RAW";

  TString hname = "EMCDistribution_" + gain + TString(full_gain ? "_FullGain" : "") + cuts;

  TH2 *h2 = NULL;
  {
    if (full_gain)
    {
      h2 = new TH2F(hname,
                    Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 100,
                    .05 * 100, 25 * 100, 64, -.5, 63.5);
      QAHistManagerDef::useLogBins(h2->GetXaxis());
    }
    else
    {
      h2 = new TH2F(hname,
                    Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 260,
                    -.2 * 100, 5 * 100, 64, -.5, 63.5);
    }
    T->Draw(
        "TOWER_" + gain + "_CEMC[].get_bineta() + 8* TOWER_" + gain + "_CEMC[].get_binphi():(TOWER_RAW_CEMC[].signal_samples[10] - TOWER_RAW_CEMC[].signal_samples[0])*(-1)>>" + hname, "", "goff");
  }

  TText *t;
  TCanvas *c1 = new TCanvas(
      "EMCDistribution_PeakSample_Fast_" + TString(full_gain ? "_FullGain" : "") + cuts,
      "EMCDistribution_PeakSample_Fast_" + TString(full_gain ? "_FullGain" : "") + cuts, 1800, 950);
  c1->Divide(8, 8, 0., 0.01);
  int idx = 1;
  TPad *p;

  for (int iphi = 8 - 1; iphi >= 0; iphi--)
  {
    for (int ieta = 0; ieta < 8; ieta++)
    {
      p = (TPad *) c1->cd(idx++);
      c1->Update();

      p->SetLogy();
      if (full_gain)
      {
        p->SetLogx();
      }
      p->SetGridx(0);
      p->SetGridy(0);

      TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi) + TString(full_gain ? "_FullGain" : "");

      TH1 *h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
                               ieta + 8 * iphi + 1);  // axis bin number is encoded as ieta+8*iphi+1

      h->SetLineWidth(0);
      h->SetLineColor(kBlue + 3);
      h->SetFillColor(kBlue + 3);

      h->GetXaxis()->SetTitleSize(.09);
      h->GetXaxis()->SetLabelSize(.08);
      h->GetYaxis()->SetLabelSize(.08);

      h->Draw();

      if (full_gain)
        h->Fit("x*gaus", "M");
      else
        h->Fit("landau", "M");

      double peak = -1;

      TF1 *fit = ((TF1 *) (h->GetListOfFunctions()->At(0)));
      if (fit)
      {
        fit->SetLineColor(kRed);
        peak = fit->GetParameter(1);
      }

      cout << Form("Finished <Col%d Row%d> = %.1f", ieta, iphi, peak)
           << endl;

      TText *t = new TText(.9, .9,
                           Form("<Col%d Row%d> = %.1f", ieta, iphi, peak));
      t->SetTextAlign(33);
      t->SetTextSize(.15);
      t->SetNDC();
      t->Draw();
    }
  }

  SaveCanvas(c1,
             TString(_file0->GetName()) + TString("_DrawPrototype3EMCalTower_") + TString(c1->GetName()), false);
}

//.........这里部分代码省略.........
  TH1F *hptFOR = (TH1F*) hptFO->Rebin(nPt,hptFO->GetTitle(),binPt);
  if (hptFOR->GetSumw2N() == 0)
    hptFOR->Sumw2();

  TH1F* hptSU = (TH1F*) f1->Get( sample+FOtag+"/PreselCut_"+sysu+"photonPt" );
  TH1F* hptSD = (TH1F*) f1->Get( sample+FOtag+"/PreselCut_"+sysd+"photonPt" );

  TH1F *hptSUR = (TH1F*) hptSU->Rebin(nPt,hptSU->GetTitle(),binPt);
  if (hptSUR->GetSumw2N() == 0)
    hptSUR->Sumw2();

  TH1F *hptSDR = (TH1F*) hptSD->Rebin(nPt,hptSD->GetTitle(),binPt);
  if (hptSDR->GetSumw2N() == 0)
    hptSDR->Sumw2();

  if ( hptFOR != 0 ) {
    for (int b = 0; b < hptFOR->GetNbinsX(); b++) {

      Double_t mainHistoContent   = hptFOR->GetBinContent(b);
      Double_t systUpHistoContent = hptSUR->GetBinContent(b);
      Double_t systDnHistoContent = hptSDR->GetBinContent(b);
      Double_t systDiffUp = fabs( (double) systUpHistoContent - mainHistoContent );
      Double_t systDiffDn = fabs( (double) mainHistoContent - systDnHistoContent );

      // use average error for histogram
      Double_t systDiff = ( systDiffUp + systDiffDn ) / 2.;
      
      Double_t statErr   = hptFOR->GetBinError(b);
      Double_t combError = sqrt( systDiff * systDiff + statErr * statErr );
      hptFOR->SetBinError(b, combError);
    } //for
  }//if
    
  cPt->SetLogy(1);
  gPad->Update();
  cPt->Update();

  //  hptR->Scale(1,"width");
  //  hptR->SetMinimum(0.02);
  //  hptR->SetMaximum(1000); 
  hptR->GetXaxis()->SetTitle("1^{st} photon p_{T} [GeV]");
  hptR->GetYaxis()->SetTitle("Number of Events / GeV");

  if(isMC) hptR->SetMarkerSize(0);
  if(isMC) hptR->Draw("histE");
  else     hptR->Draw("E X0");


  hptFOR->SetMarkerSize(0);
  hptFOR->SetLineColor(2);
  hptFOR->SetFillColor(2);
  hptFOR->SetFillStyle(3004);
  hptFOR->Draw("same hist ][ E2");

  legend->Clear();
  if(isMC) legend->SetHeader("#gamma/QCD (Sim)");
  else     legend->SetHeader("#gamma/QCD (Data)");
  if(isMC) legend->AddEntry(hptR, "#gamma", "l");
  else     legend->AddEntry(hptR, "#gamma", "p");
  legend->AddEntry(hptFOR, "Pred (from #gamma_{jet})", "f");
  legend->Draw();

  as->DrawLatex(0.17, 0.93, outLumi.c_str() );
  cPt->Update();
  cPt->SetBottomMargin(0.2 + 0.8 * cPt->GetBottomMargin() - 0.2 * cPt->GetTopMargin());

//.........这里部分代码省略.........
    TH2F *h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection =
        (TH2F *) qa_file_new->GetObjectChecked(
            "h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection", "TH2F");
    assert(h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection);

    TH1D *proj_new =
        h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection->ProjectionY(
            "qa_file_new_h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection_py");

    proj_new->Scale(1. / proj_new->GetSum());

    TH1D *proj_ref = NULL;
    if (qa_file_ref)
    {
      TH2F *h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection =
          (TH2F *) qa_file_ref->GetObjectChecked(
              "h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection", "TH2F");
      assert(h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection);

      proj_ref = h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection->ProjectionY(
          "qa_file_ref_h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection_py");
      proj_ref->Scale(1. / proj_ref->GetSum());
    }

    proj_new->GetYaxis()->SetTitle("Normalized energy distribution");
    proj_new->GetXaxis()->SetRangeUser(-10, 10);

    DrawReference(proj_new, proj_ref);
  }

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogx();
  p->SetLogy();

  {
    TH1F *h_new = (TH1F *) qa_file_new->GetObjectChecked(
        "h_QAG4Sim_CEMC_G4Hit_HitTime", "TH1F");
    assert(h_new);

    h_new->Rebin(5);
    h_new->Scale(1. / h_new->GetSum());

    TH1F *h_ref = NULL;
    if (qa_file_ref)
    {
      h_ref = (TH1F *) qa_file_ref->GetObjectChecked(
          "h_QAG4Sim_CEMC_G4Hit_HitTime", "TH1F");
      assert(h_ref);

      h_ref->Rebin(5);
      h_ref->Scale(1. / h_ref->GetSum());
    }

    h_new->GetYaxis()->SetTitleOffset(1.5);
    h_new->GetYaxis()->SetTitle("Normalized energy per bin");
    //      h_new->GetXaxis()->SetRangeUser(-0, .1);

    DrawReference(h_new, h_ref);
  }

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  //    p->SetLogx();
  p->SetLogy();

void
QA_Draw_Calorimeter_Sum_TrackProjEP(
    const char * qa_file_name_new =
        "/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap/G4Hits_sPHENIX_pi-_eta0.30_32GeV-0000.root_qa.root",
    const char * qa_file_name_ref =
        "/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap/G4Hits_sPHENIX_pi+_eta0.30_32GeV-0000.root_qa.root")
//QA_Draw_Calorimeter_Sum_TrackProjEP(const char * qa_file_name_new =
//    "data/G4sPHENIXCells_100e24GeV.root_qa.root",
//    const char * qa_file_name_ref =
//        "data/G4Hits_sPHENIX_e-_eta0_24GeV-0000.root_qa.root")
{

  SetOKStyle();
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);
  TVirtualFitter::SetDefaultFitter("Minuit2");

  TFile * qa_file_new = new TFile(qa_file_name_new);
  assert(qa_file_new->IsOpen());

  TFile * qa_file_ref = NULL;
  if (qa_file_name_ref)
    {
      qa_file_ref = new TFile(qa_file_name_ref);
      assert(qa_file_ref->IsOpen());
    }

  // obtain normalization
  double Nevent_new = 1;
  double Nevent_ref = 1;

  // obtain normalization
  double Ntrack_new = 0;
  double Ntrack_ref = 0;

  if (qa_file_new)
    {
      TH1D * h_norm = (TH1D *) qa_file_new->GetObjectChecked(
          TString("h_QAG4Sim_CalorimeterSum_Normalization"), "TH1D");
      assert(h_norm);

      Nevent_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
      Ntrack_new = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Track"));
    }
  if (qa_file_ref)
    {
      TH1D * h_norm = (TH1D *) qa_file_ref->GetObjectChecked(
          TString("h_QAG4Sim_CalorimeterSum_Normalization"), "TH1D");
      assert(h_norm);

      Nevent_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Event"));
      Ntrack_ref = h_norm->GetBinContent(h_norm->GetXaxis()->FindBin("Track"));
    }

  TCanvas *c1 = new TCanvas("QA_Draw_Calorimeter_Sum_TrackProjEP",
      "QA_Draw_Calorimeter_Sum_TrackProjEP", 1800, 600);
  c1->Divide(3, 1);
  int idx = 1;
  TPad * p;

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();

  if (Ntrack_new>0)
    {
      TH1F * h_new = (TH1F *) qa_file_new->GetObjectChecked(
          "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP", "TH1F");
      assert(h_new);

      h_new->Sumw2();
      h_new->Scale(1. / Ntrack_new);

      TH1F * h_ref = NULL;
      if (qa_file_ref)
        {
          TH1F * h_ref = (TH1F *) qa_file_ref->GetObjectChecked(
              "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP", "TH1F");
          assert(h_ref);

          h_ref->Scale(1. / Ntrack_ref);
        }

      h_new->GetYaxis()->SetTitleOffset(1.5);
      h_new->GetYaxis()->SetTitle("Count / track / bin");
      //      h_new->GetXaxis()->SetRangeUser(-0, .1);

      DrawReference(h_new, h_ref);
    }


  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();

  if (Ntrack_new>0)
    {
      TH1F * h_new = (TH1F *) qa_file_new->GetObjectChecked(
          "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP", "TH1F");
      assert(h_new);
//.........这里部分代码省略.........

void determineWorkingPoint(TString algo="csv",TString baseURL="~/scratch0/top-nosyst/plotter.root")
{
  float workPoint(0.244);
  float sfb(1.020),    sfberr(0.04);      //from BTV-11-003
  float sflight(1.08), sflighterr(0.13);  //from BTV-11-002
  //  float sfb(0.99),        sfberr(0.099);      //from BTV-11-001
  //  float sflight(1.07882), sflighterr(0.244);  //from BTV-11-001
  if(algo=="tche")
    {
      workPoint=1.7;
      sfb=0.95; sfberr=0.095;
      sflight=1.08018; sflighterr=0.1125;
    }

  TObjArray bjets=getDistributionFromPlotter(algo+"b",baseURL);
  TObjArray lightjets=getDistributionFromPlotter(algo+"light",baseURL);
  
  TH1F *bDisc=(TH1F *) bjets.At(1);
  TH1F *lightDisc=(TH1F *) lightjets.At(1);
  
  setStyle();
  gStyle->SetOptFit(0);

  TCanvas *cnv = getNewCanvas("c","c",false);
  cnv->Clear();
  cnv->SetCanvasSize(1200,1200);
  cnv->SetWindowSize(1200,1200);
  cnv->Divide(2,2);

  //draw the discriminator
  TPad *p = (TPad *)cnv->cd(1); 
  bDisc->SetTitle("b");
  bDisc->SetLineColor(1);
  bDisc->SetMarkerColor(1);
  bDisc->SetMarkerStyle(20);
  bDisc->SetFillStyle(0);
  bDisc->DrawNormalized("histe1");
  lightDisc->SetLineColor(1);
  lightDisc->SetMarkerColor(1);
  lightDisc->SetMarkerStyle(24);
  lightDisc->SetFillStyle(0);
  lightDisc->SetTitle("udcsg");
  lightDisc->DrawNormalized("histe1same");
  TLegend *leg=p->BuildLegend();
  formatForCmsPublic(p,leg,"CMS simulation",2);

  //draw the b/light efficiencies
  p=(TPad *)cnv->cd(2);
  p->SetLogy();
  TGraphAsymmErrors *bEff     = getEfficiencyFrom(bDisc);
  TGraphAsymmErrors *lightEff = getEfficiencyFrom(lightDisc);
  bEff->SetMarkerStyle(20);
  bEff->SetFillStyle(0);
  bEff->Draw("ap");
  bEff->GetXaxis()->SetTitle(bDisc->GetXaxis()->GetTitle());
  bEff->GetYaxis()->SetTitle("Efficiency");
  lightEff->SetMarkerStyle(24);
  lightEff->SetFillStyle(0);
  lightEff->Draw("p");

  //draw relatively to a given working point
  p=(TPad *)cnv->cd(3);
  Double_t baseBEff=bEff->Eval(workPoint);
  TGraphAsymmErrors *relBEff=new TGraphAsymmErrors;
  relBEff->SetMarkerStyle(20);
  relBEff->SetFillStyle(0);
  Double_t baseLightEff=lightEff->Eval(workPoint);
  TGraphAsymmErrors *relLightEff=new TGraphAsymmErrors;
  relLightEff->SetMarkerStyle(24);
  relLightEff->SetFillStyle(0);
  for(int ip=0; ip<bEff->GetN(); ip++)
    {
      Double_t cut, y,ey; 
      bEff->GetPoint(ip,cut,y);     ey = bEff->GetErrorY(ip);
      Double_t relEff(y/baseBEff);
      if(relEff<sfb+2*sfberr && relEff>sfb-2*sfberr)
	{
	  int ipt=relBEff->GetN();
	  relBEff->SetPoint(ipt,cut,relEff);
	  relBEff->SetPointError(ipt,0,0,ey/baseBEff,ey/baseBEff);
	}

      lightEff->GetPoint(ip,cut,y);     ey = lightEff->GetErrorY(ip);
      relEff=y/baseLightEff;
      //      if(relEff<sflight+7*sflighterr && relEff>sflight-7*sflighterr)
      if(relEff<sflight+3*sflighterr && relEff>sflight-3*sflighterr)
	{
	  int ipt=relLightEff->GetN();
	  relLightEff->SetPoint(ipt,cut,relEff);
	  relLightEff->SetPointError(ipt,0,0,ey/baseLightEff,ey/baseLightEff);
	}
    }

  relLightEff->Draw("ap");
  relLightEff->GetXaxis()->SetTitle( bDisc->GetXaxis()->GetTitle() );
  relLightEff->GetYaxis()->SetTitle( "#varepsilon/#varepsilon_{0}" );
  relLightEff->Fit("expo","Q+");
  TF1 *ffunc=relLightEff->GetFunction("expo");
  float newLightCut=(TMath::Log(sflight)-ffunc->GetParameter(0))/ffunc->GetParameter(1);
  float newLightCutErrPlus=(TMath::Log(sflight+sflighterr)-ffunc->GetParameter(0))/ffunc->GetParameter(1)-newLightCut;
//.........这里部分代码省略.........

void limit()
{
//=========Macro generated from canvas: limit/limit
//=========  (Thu Apr 27 14:38:33 2017) by ROOT version6.02/05
   TCanvas *limit = new TCanvas("limit", "limit",0,0,600,600);
   gStyle->SetOptFit(1);
   gStyle->SetOptStat(0);
   gStyle->SetOptTitle(0);
   limit->SetHighLightColor(2);
   limit->Range(0,0,1,1);
   limit->SetFillColor(0);
   limit->SetBorderMode(0);
   limit->SetBorderSize(2);
   limit->SetLeftMargin(0.16);
   limit->SetRightMargin(0.04);
   limit->SetTopMargin(0.06);
   limit->SetBottomMargin(0.12);
   limit->SetFrameFillStyle(0);
   limit->SetFrameBorderMode(0);
  
// ------------>Primitives in pad: pad
   TPad *pad = new TPad("pad", "pad",0,0,1,1);
   pad->Draw();
   pad->cd();
   pad->Range(-59.99999,-2.046818,1565,15.01);
   pad->SetFillColor(0);
   pad->SetBorderMode(0);
    pad->SetLogy();
   pad->SetBorderSize(2);
   pad->SetTickx(1);
   pad->SetTicky(1);
   pad->SetLeftMargin(0.16);
   pad->SetRightMargin(0.04);
   pad->SetTopMargin(0.06);
   pad->SetBottomMargin(0.12);
   pad->SetFrameFillStyle(0);
   pad->SetFrameBorderMode(0);
   pad->SetFrameFillStyle(0);
   pad->SetFrameBorderMode(0);
   
   TH1F *tmp01 = new TH1F("tmp01","Graph",100,200,1500);
   tmp01->SetMinimum(0.);
   tmp01->SetMaximum(13.98659);
   tmp01->SetStats(0);
   tmp01->SetLineStyle(0);
   tmp01->SetMarkerStyle(20);
   tmp01->GetXaxis()->SetTitle("m_{H} (GeV)");
   tmp01->GetXaxis()->SetLabelFont(42);
   tmp01->GetXaxis()->SetLabelOffset(0.01);
   tmp01->GetXaxis()->SetTitleSize(0.05);
   tmp01->GetXaxis()->SetTickLength(0.02);
   tmp01->GetXaxis()->SetTitleOffset(1.08);
   tmp01->GetXaxis()->SetTitleFont(42);
   tmp01->GetYaxis()->SetTitle("95% CL limit on #sigma/#sigma_{SM}");
   tmp01->GetYaxis()->SetLabelFont(42);
   tmp01->GetYaxis()->SetLabelOffset(0.007);
   tmp01->GetYaxis()->SetTitleSize(0.05);
   tmp01->GetYaxis()->SetTickLength(0.02);
   tmp01->GetYaxis()->SetTitleOffset(1.56);
   tmp01->GetYaxis()->SetTitleFont(42);
   tmp01->GetZaxis()->SetLabelFont(42);
   tmp01->GetZaxis()->SetLabelOffset(0.007);
   tmp01->GetZaxis()->SetTitleSize(0.05);
   tmp01->GetZaxis()->SetTickLength(0.02);
   tmp01->GetZaxis()->SetTitleFont(42);
   tmp01->Draw("AXIS");
   
   TH1F *tmp02 = new TH1F("tmp02","Graph",100,200,1500);
   tmp02->SetMinimum(0);
   tmp02->SetMaximum(13.98659);
   tmp02->SetStats(0);
   tmp02->SetLineStyle(0);
   tmp02->SetMarkerStyle(20);
   tmp02->GetXaxis()->SetTitle("m_{H} (GeV)");
   tmp02->GetXaxis()->SetLabelFont(42);
   tmp02->GetXaxis()->SetLabelOffset(0.01);
   tmp02->GetXaxis()->SetTitleSize(0.05);
   tmp02->GetXaxis()->SetTickLength(0.02);
   tmp02->GetXaxis()->SetTitleOffset(1.08);
   tmp02->GetXaxis()->SetTitleFont(42);
   tmp02->GetYaxis()->SetTitle("95% CL limit on #sigma/#sigma_{SM}");
   tmp02->GetYaxis()->SetLabelFont(42);
   tmp02->GetYaxis()->SetLabelOffset(0.007);
   tmp02->GetYaxis()->SetTitleSize(0.05);
   tmp02->GetYaxis()->SetTickLength(0.02);
   tmp02->GetYaxis()->SetTitleOffset(1.56);
   tmp02->GetYaxis()->SetTitleFont(42);
   tmp02->GetZaxis()->SetLabelFont(42);
   tmp02->GetZaxis()->SetLabelOffset(0.007);
   tmp02->GetZaxis()->SetTitleSize(0.05);
   tmp02->GetZaxis()->SetTickLength(0.02);
   tmp02->GetZaxis()->SetTitleFont(42);
   tmp02->Draw("AXIGSAME");
   
   Double_t Graph0_fx3001[27] = {
   200,
   250,
   300,
   350,
   400,
//.........这里部分代码省略.........

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

  //sensor adc for each layer
  double max_sensor_x[2][8][24] = {{{0.}}};
  double max_sensor_y[2][8][24] = {{{0.}}};
  double min_sensor_x[2][8][24] = {{{0.}}};
  double min_sensor_y[2][8][24] = {{{0.}}};
  
  TH1D* hsensor_high[2][8][24];
  TH1D* hsensor_low[2][8][24];
 
  for(int iarm = 0;iarm < 2;iarm++){
    for(int ilayer = 0;ilayer < 8;ilayer++){
      for(int isen = 0;isen < 24;isen++){
        char name[100];
        sprintf(name,"hsensor_arm%d_layer%d_index%d_high",iarm,ilayer,isen);
        hsensor_high[iarm][ilayer][isen] = new TH1D(name,name,300,-40.5,259.5);
        sprintf(name,"hsensor_arm%d_layer%d_index%d_low",iarm,ilayer,isen);
        hsensor_low[iarm][ilayer][isen] = new TH1D(name,name,300,-40.5,259.5);
        max_sensor_x[iarm][ilayer][isen] = -9999;
        max_sensor_y[iarm][ilayer][isen] = -9999;
        min_sensor_x[iarm][ilayer][isen] = 9999;
        min_sensor_y[iarm][ilayer][isen] = 9999;
      }
    }
  }
  
  ifstream sensor_pos("sensor_position.txt");
  string s;
  while(getline(sensor_pos,s)){
    stringstream ss(s);
    int arm = 0;
    int layer = 0;
    int sensor = 0;
    double x0=0;
    double x1=0;
    double y0=0;
    double y1=0;
    ss>>arm>>layer>>sensor>>x0>>x1>>y0>>y1;
//    cout <<arm <<" "<<layer<<" "<<sensor<<" "<<x0<<" "<<x1<<" "<<y0<<" "<<y1<<endl;
    max_sensor_x[arm][layer][sensor] = x1;
    min_sensor_x[arm][layer][sensor] = x0;
    max_sensor_y[arm][layer][sensor] = y1;
    min_sensor_y[arm][layer][sensor] = y0;
  }
 
  MpcExMapper* mapper = MpcExMapper::instance(); 
  for(unsigned int i = 0;i < 50000;i++){
    hkey_adc_high->SetAxisRange(i,i,"X");
    hkey_adc_low->SetAxisRange(i,i,"X");
    TH1D* htemp0 = hkey_adc_high->ProjectionY();
    TH1D* htemp1 = hkey_adc_low->ProjectionY();
    if(htemp0->GetEntries() < 10) continue;
    int arm = mapper->get_arm(i);
    int quadrant = mapper->get_quadrant(i);
    int sensor = mapper->get_sensor_in_quadrant(i);
    int index = 6*quadrant+sensor;
    int layer = mapper->get_layer(i);
    hsensor_high[arm][layer][index]->Add(htemp0);
    hsensor_low[arm][layer][index]->Add(htemp1);
  }
  
  for(int iarm = 0;iarm < 2;iarm++){
    for(int ilayer = 0;ilayer < 8;ilayer++){
      char cname[100];
      sprintf(cname,"csensor_arm%d_layer%d_%d_%d",iarm,ilayer,runnumber,segment);
      TCanvas* c = new TCanvas(cname,cname,1400,800);
      for(int index = 0;index < 24;index++){
        char pname[100];
        sprintf(pname,"sensor_arm%d_layer%d_index%d",iarm,ilayer,index);
        double x0 = min_sensor_x[iarm][ilayer][index];
        double x1 = max_sensor_x[iarm][ilayer][index];
        double y0 = min_sensor_y[iarm][ilayer][index];
        double y1 = max_sensor_y[iarm][ilayer][index];
        if(ilayer%2 == 0){
          y0 = y0 - 0.85;
          y1 = y1 + 0.85;
        }
        if(ilayer%2 == 1){
          x0 = x0 - 0.85;
          x1 = x1 + 0.85;
        }
        TPad* pad = new TPad(pname,pname,0.5+x0/40.,0.5+y0/40.,0.5+x1/40.,0.5+y1/40.);
        pad->cd();
        pad->SetLogy();
	hsensor_high[iarm][ilayer][index]->Draw("");
	hsensor_low[iarm][ilayer][index]->SetLineColor(kRed);
        hsensor_low[iarm][ilayer][index]->Draw("");
        c->cd();
        pad->DrawClone("same");
        delete pad;
      }
      sprintf(cname,"%d/csensor_arm%d_layer%d_%d_%d.gif",runnumber,iarm,ilayer,runnumber,segment);
      c->Print(cname,"gif"); 
      delete c;
    }
  }
}

void FastVsSlowSimPtRes() {

  Int_t plusTPC =0;

  gROOT->LoadMacro("~/fig_template.C"); // figure style
  myOptions(0);
  gROOT->ForceStyle();

  TCanvas *myCan = new TCanvas("myCan");
  myCan->Draw();
  myCan->cd();
  
  TPad *myPad = new TPad("myPad", "The pad",0,0,1,1);
  myPadSetUp(myPad,0.15,0.04,0.04,0.15);
  myPad->Draw();   myPad->cd();
  myPad->SetGridx();   myPad->SetGridy();  myPad->SetLogx(); myPad->SetLogy();


  TLegend *leg = new TLegend(0.44,0.13,0.1.7,0.9,"","brCDN"); 
  leg->SetFillColor(0);

  TGraph *c[6];
  
  // Current ITS +++++++++++++++++++++++++++++++++++++++++
  Int_t color=1; Int_t linewidth=2;

  Int_t pi =0;
 
  DetectorK its("ALICE","ITS");
  its.MakeAliceCurrent(0,plusTPC);
  its.SetMaxRadiusOfSlowDetectors(0.1);
  its.SolveViaBilloir(0);
  Int_t color=1; Int_t linewidth=2;

  c[pi] = its.GetGraphMomentumResolution(color,linewidth);
  c[pi]->Draw("AC");

  leg->AddEntry(c[pi],"FastTool: Current ITS","l");


  // Current ITS +++++++++++++++++++++++++++++++++++++++++

  Int_t color=3; Int_t linewidth=2;
  Int_t pi =2;

  DetectorK its("ALICE","ITS");
  its.MakeAliceCurrent(0,plusTPC);
  
  its.SetRadius("bpipe",2.0);
  its.AddLayer("spd0", 2.2,1,1,1);  

  its.SetRadius("spd0",2.2); its.SetRadiationLength("spd0",X0); its.SetResolution("spd0",resRPhi,resZ);
  its.SetRadius("spd1",4.8);   its.SetRadiationLength("spd1",X0); its.SetResolution("spd1",resRPhi,resZ);
  its.SetRadius("spd2",9.1);   its.SetRadiationLength("spd2",X0); its.SetResolution("spd2",resRPhi,resZ);

  its.SetMaxRadiusOfSlowDetectors(0.1);
  its.SolveViaBilloir(0);

  c[pi] = its.GetGraphMomentumResolution(color,linewidth);
  c[pi]->Draw("C");

  leg->AddEntry(c[pi],"FastTool: \"New SPDs\"","l");


  // ALL NEW +++++++++++++++++++++++++++++++++++++++++++

  color=2; Int_t linewidth=2;
  Int_t pi =1; 


  // for a 0.8,0.2 weight configuration
  
  DetectorK *itsU = new DetectorK((char*)"ALICE",(char*)"ITS");
  
  itsU->AddLayer((char*)"bpipe", 2.0,0.0022); // beam pipe
  itsU->AddLayer((char*)"vertex",  0,     0); // dummy vertex for matrix calculation
  
  itsU->AddLayer("ddd1",  2.2 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd2",  3.8 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd3",  6.8 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd4", 12.4 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd5", 23.5 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd6", 39.6 ,  X0, resRPhi, resZ); 
  itsU->AddLayer("ddd7", 43.0 ,  X0, resRPhi, resZ); 
 
  if(plusTPC) itsU->AddTPC(0.1,0.1);
  itsU->SetMaxRadiusOfSlowDetectors(0.1);
  itsU->SolveViaBilloir(0);
  itsU->PrintLayout();

  
  c[pi] = itsU->GetGraphMomentumResolution(color,linewidth);
  c[pi]->Draw("C");

  leg->AddEntry(c[pi],"FastTool: \"All New\" ","l");

  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



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

void AtlasExample() 
{ 

#ifdef __CINT__
  gROOT->LoadMacro("AtlasUtils.C");
#endif

  SetAtlasStyle();

  Int_t icol1=5;
  Int_t icol2=5;

  TCanvas* c1 = new TCanvas("c1","single inclusive jets",50,50,600,600);
  // TCanvas* c1 = new TCanvas("c1","single inclusive jets");
  TPad* thePad = (TPad*)c1->cd();
  thePad->SetLogy();

  Double_t ymin=1.e-3;  Double_t ymax=2e7;
  Double_t xmin=60.00;  Double_t xmax=3500.;
  TH1F *h1 = thePad->DrawFrame(xmin,ymin,xmax,ymax);
  h1->SetYTitle("d#sigma_{jet}/dE_{T,jet} [fb/GeV]");
  h1->SetXTitle("E_{T,jet}  [GeV]");
  h1->GetYaxis()->SetTitleOffset(1.4);
  h1->GetXaxis()->SetTitleOffset(1.4);
  //h1->GetXaxis()->SetNdivisions(5);
  h1->Draw();

  const Int_t ncut=1;
  TGraphErrors *data[ncut];

  for (Int_t icut=0; icut<ncut; icut++) { // loop over cuts

    TGraphErrors *g1[nren][ncut];
    for (Int_t ir=0; ir<nren; ir++) { // loop over ren scale
      g1[ir][icut]= GetGraph(ir,ir,icut,0);
      if (g1[ir][icut]) 
	cout << g1[ir][icut]->GetTitle() << " found "  << g1[ir][icut]->GetName()  << endl;
      else { 
	cout << " g1 not  found " << endl; 
	return; 
      } 
      g1[ir][icut]->SetLineColor(1);
      g1[ir][icut]->SetMarkerStyle(0);
      //g1[ir][icut]->Draw("C");
    }

    char daname[100];
    sprintf(daname,"data_%d",icut); 
    data[icut]=(TGraphErrors*)g1[0][icut]->Clone(daname); 
    data[icut]->SetMarkerStyle(20);
    data[icut]->SetMarkerColor(1);

    // just invent some data
    for (Int_t i=0; i< data[icut]->GetN(); i++) {
      Double_t x1,y1,e,dx1=0.;
      data[icut]->GetPoint(i,x1,y1);
      Double_t r1 = 0.4*(gRandom->Rndm(1)+2);
      Double_t r2 = 0.4*(gRandom->Rndm(1)+2);
      //cout << " i= " << i << " x1= " << x1 << " y1= " << y1 << " r= " << r << endl;
      Double_t y;
      if (icut==0) y=r1*y1+r1*r2*r2*x1/50000.;
      else         y=r1*y1;
      e=sqrt(y*1000)/200;
      data[icut]->SetPoint(i, x1,y);
      data[icut]->SetPointError(i,dx1,e);
    }
    //data[icut]->Print();
  
    TGraphAsymmErrors* scale[ncut];
    TGraphAsymmErrors* scalepdf[ncut];

    scale[icut]=  myMakeBand(g1[0][icut],g1[1][icut],g1[2][icut]);
    //printf(" band1: \n");
    //scale->Print();

    scalepdf[icut]=(TGraphAsymmErrors* ) scale[icut]->Clone("scalepdf");

    TGraphErrors *gpdf[NUMPDF][ncut];
    for (Int_t ipdf=0; ipdf<NUMPDF; ipdf++) {
      gpdf[ipdf][icut]= GetGraph(0,0,icut,ipdf);
      if (gpdf[ipdf][icut]) 
	cout << gpdf[ipdf][icut]->GetTitle() << " found "  << gpdf[ipdf][icut]->GetName() << endl;
      else { 
	cout << " gpdf not  found " << endl; 
	return; 
      } 
      gpdf[ipdf][icut]->SetLineColor(2);
      gpdf[ipdf][icut]->SetLineStyle(1);
      gpdf[ipdf][icut]->SetMarkerStyle(0);
      myAddtoBand(gpdf[ipdf][icut],scalepdf[icut]); 
    }

    scalepdf[icut]->SetFillColor(icol2);
    scalepdf[icut]->Draw("zE2"); 
    scale[icut]->SetFillColor(icol1);
    scale[icut]->Draw("zE2");
    g1[0][icut]->SetLineWidth(3);
    g1[0][icut]->Draw("z");
    data[icut]->Draw("P");
    
//.........这里部分代码省略.........

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

   TH1D* backgroundsForRatio = (TH1D*)allBackgrounds->Clone("bkgds_for_ratio");

   // create main pad
   const double mainPad_ylow  = 0.3;
   const double mainPad_yhigh = 0.95;
   const double mainPad_xlow = 0.;
   const double mainPad_xhigh = 0.95;
   const double pad_margin_left = 0.2;
   const double pad_margin_right  = 0.02;
   double main_y_max = -99;
   double main_y_min = -99;

   TPad* mainPad = new TPad("main","main",mainPad_xlow,mainPad_ylow,mainPad_xhigh,mainPad_yhigh);
   mainPad->SetMargin(pad_margin_left,pad_margin_right,0.,.05); // left, right, bottom, top
   mainPad->Draw();
   mainPad->cd();

   // find max y
   if(allBackgrounds){
     int maxBinBkg = allBackgrounds->GetMaximumBin();
     double bkgYmax = allBackgrounds->GetBinContent(maxBinBkg);
     main_y_max = bkgYmax;
   }
   if(data){
     int maxBinData = data->GetMaximumBin();
     double dataYmax = data->GetBinContent(maxBinData);
     if(dataYmax > main_y_max) main_y_max = dataYmax;
   }
   if(main_y_max < .002) main_y_max = .1;

   if(doLogAxis){
     mainPad->SetLogy();
     main_y_min = 0.07;
     main_y_max*=60;
   }else{
     main_y_min = 0.001;
     main_y_max*=1.6;
   }

   // draw axis
   allBackgrounds->SetMaximum(main_y_max);
   allBackgrounds->SetMinimum(main_y_min);
   std::string title = std::string(";") + axisName + ";Events";
   std::stringstream binw; 
   
   binw<<allBackgrounds->GetBinWidth(1);
   std::string width = binw.str();
   if (axisName.find("GeV") != std::string::npos){
     width += " GeV";
     title += " / "+width;
   }
   allBackgrounds->SetTitle(title.c_str());   
   allBackgrounds->GetYaxis()->SetTitleOffset(1.4);
   allBackgrounds->GetYaxis()->SetTitleSize(0.06);
   allBackgrounds->GetXaxis()->SetTitleOffset(1.2);
   allBackgrounds->GetYaxis()->SetLabelSize(0.04);
   allBackgrounds->GetXaxis()->SetLabelSize(0.04);
   allBackgrounds->Draw("hist");

   // draw stack with error
   Stack->Draw("hist same");
   if(err){
     err->Draw("E2 same");
   }

//.........这里部分代码省略.........
     }
   } 
   TH1D* backgroundsForRatio = (TH1D*)allBackgrounds->Clone("bkgds_for_ratio");

   // create main pad
   const double mainPad_ylow  = 0.3;
   const double mainPad_yhigh = 0.95;
   const double mainPad_xlow = 0.;
   const double mainPad_xhigh = 0.95;
   const double pad_margin_left = 0.2;
   const double pad_margin_right  = 0.02;
   double main_y_max = -99;
   double main_y_min = -99;

   TPad* mainPad = new TPad("main","main",mainPad_xlow,mainPad_ylow,mainPad_xhigh,mainPad_yhigh);
   mainPad->SetMargin(pad_margin_left,pad_margin_right,0.,.05); // left, right, bottom, top
   mainPad->Draw();
   mainPad->cd();

   // find max y
   if(allBackgrounds){
     int maxBinBkg = allBackgrounds->GetMaximumBin();
     double bkgYmax = allBackgrounds->GetBinContent(maxBinBkg);
     main_y_max = bkgYmax;
   }
   if(data){
     int maxBinData = data->GetMaximumBin();
     double dataYmax = data->GetBinContent(maxBinData);
     if(dataYmax > main_y_max) main_y_max = dataYmax;
   }
   if(main_y_max < .002) main_y_max = .1;

   if(doLogAxis){
     mainPad->SetLogy();
     main_y_min = 0.007;
     main_y_max*=600;
   }else{
     main_y_min = 0.001;
     main_y_max*=1.6;
   }

   // draw axis
   allBackgrounds->SetMaximum(main_y_max);
   allBackgrounds->SetMinimum(main_y_min);
   std::string title = std::string(";") + axisName + ";events/bin";
   allBackgrounds->SetTitle(title.c_str());   
   allBackgrounds->GetYaxis()->SetTitleOffset(1.4);
   allBackgrounds->GetYaxis()->SetTitleSize(0.06);
   allBackgrounds->GetXaxis()->SetTitleOffset(1.2);
   allBackgrounds->GetYaxis()->SetLabelSize(0.04);
   allBackgrounds->GetXaxis()->SetLabelSize(0.04);
   allBackgrounds->Draw("hist");

   // draw stack with error
   Stack->Draw("hist same");
   if(err){
     err->Draw("E2 same");
   }

   // make error histograms for the ratio plot
   TGraphAsymmErrors* ratioErr = tools::getRatio(err, backgroundsForRatio);

   // draw data
   if(!isSignalRegion && data)
      data->Draw("e same");