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

void QA_Draw_CEMC_G4Hit(
    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")

{
  SetsPhenixStyle();
  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());
  }

  TCanvas *c1 = new TCanvas("QA_Draw_CEMC_G4Hit", "QA_Draw_CEMC_G4Hit", 1800, 900);
  c1->Divide(4, 2);
  int idx = 1;
  TPad *p;

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

  TH2F *h_QAG4Sim_CEMC_G4Hit_XY = (TH2F *) qa_file_new->GetObjectChecked(
      "h_QAG4Sim_CEMC_G4Hit_XY", "TH2F");
  assert(h_QAG4Sim_CEMC_G4Hit_XY);
  h_QAG4Sim_CEMC_G4Hit_XY->GetYaxis()->SetTitleOffset(1.5);
  h_QAG4Sim_CEMC_G4Hit_XY->Draw("COLZ");

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

  TH2F *h_QAG4Sim_CEMC_G4Hit_RZ = (TH2F *) qa_file_new->GetObjectChecked(
      "h_QAG4Sim_CEMC_G4Hit_RZ", "TH2F");
  assert(h_QAG4Sim_CEMC_G4Hit_RZ);
  h_QAG4Sim_CEMC_G4Hit_RZ->GetYaxis()->SetTitleOffset(1.5);
  h_QAG4Sim_CEMC_G4Hit_RZ->Draw("COLZ");

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

  {
    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->ProjectionX(
            "qa_file_new_h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection_px");

    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->ProjectionX(
          "qa_file_ref_h_QAG4Sim_CEMC_G4Hit_LateralTruthProjection_px");
      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->SetLogz();

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

void readroot(char* input="/gpfs/mnt/gpfs02/phenix/mpcex/liankun/Run16/Ana/offline/analysis/mpcexcode/Liankun/macro/offline_monitor/453692/Run16Ana_MinBias_NoCMN_Sub-453692-0.root"){
  char input_file[5000];
  strcpy(input_file,input);
  strtok(input_file, "-");
  int runnumber = atoi(strtok(0, "-"));
  int segment = atoi(strtok(strtok(0, "-"), "."));

  cout<<"Run Number: "<<runnumber<<" segment: "<<segment<<endl;

  gSystem->Load("libMyMpcEx.so");
  TFile* ifile = new TFile(input,"READONLY");
  if(!ifile){
    cout<<"Open "<<input<<" failed !"<<endl;
    return;
  }
  TH2D* hkey_adc_high;
  TH2D* hkey_adc_low;
  TH2D* hkey_rawadc_high;
  TH2D* hkey_rawadc_low;

  TH2D* hHL_sensor[2][8][24];
  TH2D* hHL_sensor_raw[2][8][24];
  Exogram* hgrammy_high[2];
  Exogram* hgrammy_low[2];
  Exogram* hgrammy_combine[2];
  TH2D* hadc_mpc_e[2];
  TH2D* hlayer_adc_high[2];
  TH2D* hlayer_adc_low[2];
  TH2D* htower_e;

  TH2D* hbbc_adc[2];
  TH2D* hbbc_adc_low[2];
  TH2D* hbbc_nhits[2];

  char hname[500];
  vector<TH2D*> h2d_list;
  hkey_adc_high = (TH2D*)ifile->Get("hkey_adc_high");
  h2d_list.push_back(hkey_adc_high);

  hkey_adc_low = (TH2D*)ifile->Get("hkey_adc_low");
  h2d_list.push_back(hkey_adc_low);
  
  hkey_rawadc_high = (TH2D*)ifile->Get("hkey_rawadc_high");
  h2d_list.push_back(hkey_rawadc_high);

  hkey_rawadc_low = (TH2D*)ifile->Get("hkey_rawadc_low");
  h2d_list.push_back(hkey_rawadc_low);

  hbbc_nhits[0] = (TH2D*)ifile->Get("hbbc_nhits_arm0");
  hbbc_nhits[1] = (TH2D*)ifile->Get("hbbc_nhits_arm1");
  h2d_list.push_back(hbbc_nhits[0]);
  h2d_list.push_back(hbbc_nhits[1]);
  
  hbbc_adc[0] = (TH2D*)ifile->Get("hbbc_adc_arm0");
  hbbc_adc[1] = (TH2D*)ifile->Get("hbbc_adc_arm1");
  h2d_list.push_back(hbbc_adc[0]);
  h2d_list.push_back(hbbc_adc[1]);
  
  hbbc_adc_low[0] = (TH2D*)ifile->Get("hbbc_adc_low_arm0");
  hbbc_adc_low[1] = (TH2D*)ifile->Get("hbbc_adc_low_arm1");
  h2d_list.push_back(hbbc_adc_low[0]);
  h2d_list.push_back(hbbc_adc_low[1]);
  
  hlayer_adc_high[0] = (TH2D*)ifile->Get("hlayer_adc_high_arm0");
  hlayer_adc_high[1] = (TH2D*)ifile->Get("hlayer_adc_high_arm1");
  h2d_list.push_back(hlayer_adc_high[0]);
  h2d_list.push_back(hlayer_adc_high[1]);

  hlayer_adc_low[0] = (TH2D*)ifile->Get("hlayer_adc_low_arm0");
  hlayer_adc_low[1] = (TH2D*)ifile->Get("hlayer_adc_low_arm1");
  h2d_list.push_back(hlayer_adc_low[0]);
  h2d_list.push_back(hlayer_adc_low[1]);

  htower_e = (TH2D*)ifile->Get("htower_e");
  h2d_list.push_back(htower_e);

  hadc_mpc_e[0] = (TH2D*)ifile->Get("hadc_mpc_e_arm0");
  hadc_mpc_e[1] = (TH2D*)ifile->Get("hadc_mpc_e_arm1");
  h2d_list.push_back(hadc_mpc_e[0]);
  h2d_list.push_back(hadc_mpc_e[1]);

  for(unsigned int i = 0;i < h2d_list.size();i++){
    string s = h2d_list[i]->GetName();
    stringstream ss("");
    ss<<s<<"_"<<runnumber<<"_"<<segment;
    s=ss.str();
//    cout << s <<endl;
    TCanvas* c = new TCanvas(s.c_str(),s.c_str(),1200,800);
    c->SetLogz();
    h2d_list[i]->Draw("colz");
    ss.str("");
    ss << runnumber <<"/"<<s<<".gif";
    s=ss.str();
    c->Print(s.c_str(),"gif");
    delete c;
  }

  vector<Exogram*> hexo_list;
  hgrammy_high[0] = (Exogram*)ifile->Get("hgrammy_high0");
  hgrammy_high[1] = (Exogram*)ifile->Get("hgrammy_high1");
//.........这里部分代码省略.........

void plot2d(int hid) {
  gStyle->SetOptStat(0);
  char c[50];
  c1->Clear();
  if(hid!=2 && hid!=4){
    c1->Divide(2,3);
    for(int disc=0; disc<kFgtNumDiscs; disc++){
      TPad *pad = c1->cd(disc+1);
      pad->SetLogz(1);
      pad->SetTopMargin(0.01);   pad->SetBottomMargin(0.02);
      sprintf(c,"Disc%1d%s",disc+1,c2dHist[hid]);
      //printf("Getting %s\n",c);
      TH2F *h = hist2[disc][hid] = (TH2F*)file->Get(c);
      h->Draw("COLZ");  
    }
  }else if(hid==2){  // special case for timing per APVboard
    gStyle->SetOptTitle(0);
    gStyle->SetOptFit(0);
    c1->Divide(2,1);
    char txt[100];
    //int disc=3; char name[100]="MaxAdc";
    int disc=4; char name[100]="LandauMPV";
    //int disc=5; char name[100]="LandauMPV-3Sing";
    TVirtualPad *pad2 = c1->cd(1);
    pad2->Divide(1,2); 
    TVirtualPad *pad3=pad2->cd(1);
    pad3->SetLogz(1);  pad2->SetLogz(0); 
    sprintf(c,"Disc%1d%s",disc,c2dHist[hid]);
    TH2F *h = hist2[disc][hid] = (TH2F*)file->Get(c);
    h->Draw("COLZ");
    TText *tt1= new TText(0.05,0.1,"(RDO-1)*12+ARM*2+GRP"); tt1->SetTextAngle(90); tt1->SetNDC(); tt1->Draw();    
    sprintf(txt,"Tbin for %s",name);
    TText *tt2= new TText(0.3,0,txt); tt2->SetNDC(); tt2->Draw();
    
    TVirtualPad* pad4 = c1->cd(2);
    pad4->SetTopMargin(0.01); pad4->SetBottomMargin(0.1);
    int maxid=0;
    float off,max=0;
    TH1D *h1[24];
    float mean[24];
    for(int i=0; i<24; i++){
      char ccc[10]; sprintf(ccc,"_%d_%d",disc,i);
      h1[i] = h->ProjectionX(ccc,i+1,i+1); 
      if(h1[i]->GetMaximum() > max && i!=0) {max=h1[i]->GetMaximum(); maxid=i; }
    }
    off=max/4.0;
    printf("max=%f off=%f\n",max,off);
    for(int i=0; i<24; i++){
      h1[i]->GetXaxis()->SetRangeUser(2,11);
      int res = h1[i]->Fit("gaus","0Q");
      TF1* f=h1[i]->GetFunction("gaus");
      if(h1[i]->GetMaximum()>max/3 && res==0){
	mean[i] = f->GetParameter(1);
	//mean[i]=h1[i]->GetMean();
      }else{mean[i]=0;};
      //printf("%d mean=%f\n",i,mean[i]);
    }
    //h1[maxid]->SetLineColor(maxid+1); h1[maxid]->SetLineWidth(2); h1[maxid]->Draw("PL");
    for(int rdo=1; rdo<=2; rdo++){
      for(int arm=0; arm<6; arm++){
	for(int grp=0; grp<2; grp++){
	  i=(rdo-1)*12+arm*2+grp;
	  int nb=h1[i]->GetNbinsX();
	  for(int t=0; t<nb; t++){ h1[i]->AddBinContent(t+1,off*i); }
	  h1[i]->SetLineColor(i%6+1); h1[i]->SetLineWidth(3);
	  if(i==0) {
	    h1[i]->SetMinimum(0);
	    h1[i]->SetMaximum(max*6.5);
	    h1[i]->Draw("PL");
	  } else {h1[i]->Draw("PL same");}	
	  char name[100]; 
	  sprintf(name,"Rdo%1dArm%1dGrp%1d",rdo,arm,grp);
	  TText *tx = new TText(8.5,(max/4.0)*(i+0.2),name); tx->SetTextColor(i%6+1); tx->SetTextSize(0.03);
	  tx->Draw();
	}
      }
    }
    //    TText *tt3= new TText(0.95,0.1,"offsets added by (RDO-1)*12+ARM*2+GRP"); tt3->SetTextAngle(90); tt3->SetNDC(); tt3->Draw();
    TText *tt4= new TText(0.4,0,txt); tt4->SetNDC(); tt4->Draw();
    
    //correlation
    float t2[24]={-8.47, -5.16, -0.21, -2.23,  1.11, -4.09, 
		  -3.13, -9.08, -5.88, -7.01, -6.22, -9.79,  
		  0.75, -8.91,  0.16,  1.12, -0.99, -4.56,  
		  7.57, -3.68,  7.12, -6.54, -4.08, -8.21};
    TGraph *g= new TGraph(1);
    int j=0;
    for(int i=0; i<24; i++){
      if(mean[i]>0) {g->SetPoint(j,(mean[i]-6.0)*27,t2[i]); j++;}
    }
    TVirtualPad* pad5=pad2->cd(2);
    g->SetMarkerStyle(20+i/6); g->SetMarkerSize(1);
    g->Draw("ap");
    for(int i=0; i<24; i++){
      TGraph *g2= new TGraph(1);
      if(mean[i]>0) {g2->SetPoint(j,(mean[i]-6.0)*27,t2[i]); j++;}
      g2->SetMarkerStyle(20+i/6); g2->SetMarkerSize(2); g2->SetMarkerColor(i%6+1);
      g2->Draw("p");
    }
    
//.........这里部分代码省略.........

// -----------------------------------------------------------------------------
//
void test() {

  time_t start = TTimeStamp().GetSec();

  set_plot_style();

  bool draw = true;
  //bool debug = true;

  // Define analysis configuration
  PSet ps;
  defaultPSet(ps);

  // Response plots
  if (false) {
    //xSectDistr(ps);
    responseProfile();
    return;
  }
  
  // Print configuration
  std::stringstream ss;
  printPSet(ps,ss);
  std::cout << ss.str() << std::endl;
  
  // Params to store
  DoubleVV ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err;
  IntV length;
  clear( ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err, length );
  init( ps, ratio, ratio_errh, ratio_errl, pass, pass_err, fail, fail_err, length );
  
  // Loop through Meff bins 
  int loop = 0;
  int nloops = ps.nmeff;
  for ( int imeff = 0; imeff < ps.nmeff; ++imeff ) {

    // Generate numbers in (x1,x2) plane
    DoubleVV dalitz;
    generateTruth( ps, imeff, dalitz, true );

    // Integrate across dalitz plane
    integrate( ps, imeff, dalitz, ratio, ratio_errh, ratio_errl, 
	       pass, pass_err, fail, fail_err, length );
    
    // Labeling
    std::stringstream ss;
    ss << "Meff" << int( ps.meff_bins[imeff] );
    
    // New canvas for plots
    TCanvas* c1 = 0;
    if (draw) c1 = new TCanvas( TString("Canvas"+ss.str()), "" );
    
    // Pad for cross-section plot
    TPad* pad = 0;
    if (draw) pad = new TPad(TString("Pad"+ss.str()),"",0.,0.,1.,1.);
    if (pad) {
      pad->SetGrid();
      pad->Draw();
      pad->cd();
      pad->SetLogz();
    }
    TH1F* hr = 0;
    if (draw) hr = pad->DrawFrame(ps.min,ps.min,ps.max,ps.max);
    
    // Histo title
    if (hr) {
      std::stringstream sss;
      sss << "M_{eff}=" << ps.meff_bins[imeff] << " GeV"
	  << ", p_{T1}=" << dr(ps.pt1_bins[imeff],1) << " GeV"
	  << ", p_{T2}=" << dr(ps.pt2_bins[imeff],1) << " GeV"
	  << ", p_{T3}=" << dr(ps.pt3_bins[imeff],1) << " GeV";
      hr->SetTitle( sss.str().c_str() );
      hr->GetXaxis()->SetTitle( "x_{2}" );
      hr->GetYaxis()->SetTitle( "x_{1}" );
    }
    
    // Create 2D cross-section plot
    TH2D* his = 0;
    if (draw) his = new TH2D(TString("Histo"+ss.str()),"",
			     ps.nbins,ps.min,ps.max,
			     ps.nbins,ps.min,ps.max);
    
    //double x3 = ( 2. * ps.pt3_bins[imeff] ) / ( ps.meff_bins[imeff] + ps.pt3_bins[imeff] );
    
    // Fill 2D cross-section plot
    for ( int x2_bin = 0; x2_bin < ps.nbins; ++x2_bin ) { 
      for ( int x1_bin = 0; x1_bin < ps.nbins; ++x1_bin ) { 
// 	std::cout << " Fill:"
// 		  << " x2_bin: " << x2_bin 
// 		  << " x2: " << val(x2_bin,nbins) 
// 		  << " x1_bin: " << x1_bin 
// 		  << " x1: " << val(x1_bin,nbins) 
// 		  << " val: " << dalitz[x2_bin][x1_bin]
// 		  << std::endl;
	if (his) his->Fill( val(x2_bin,ps)+ps.width/2.,
			    val(x1_bin,ps)+ps.width/2.,
			    dalitz[x2_bin][x1_bin] ); 
      }
//.........这里部分代码省略.........

void EMCDistribution_ADC(bool log_scale = true)
{
  TString gain = "RAW";

  TText *t;
  TCanvas *c1 = new TCanvas(
      "EMCDistribution_ADC_" + gain + TString(log_scale ? "_Log" : "") + cuts,
      "EMCDistribution_ADC_" + 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();

      if (log_scale)
      {
        p->SetLogz();
      }
      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 TH2F(hname,
                     Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 24, -.5,
                     23.5,
                     //                128+64, 0, 3096);
                     4098, -1, 4097);
      //          else
      //            h = new TH2F(hname,
      //                Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100,
      //                -.050, .5,128,0,2048);

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

      //          if (log_scale)
      //            QAHistManagerDef::useLogBins(h->GetYaxis());

      TString sdraw = "TOWER_" + gain + "_CEMC[].signal_samples[]:fmod(Iteration$,24)>>" + hname;
      TString scut =
          Form(
              "TOWER_%s_CEMC[].get_bineta()==%d && TOWER_%s_CEMC[].get_binphi()==%d",
              gain.Data(), ieta, gain.Data(), iphi);

      cout << "T->Draw(\"" << sdraw << "\",\"" << scut << "\");" << endl;

      T->Draw(sdraw, scut, "colz");

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

// draw the same thing but after reco
void genPlots02(std::string fullPath, int nOverlay = 500, bool custBinning = false)
{
    const int fVerbose(1);
    setTDRStyle();
    gStyle->SetOptStat(112211);
    gStyle->SetPalette(1);
    // Canvas
    c = new TCanvas("c2","c2",1000,600);
    const unsigned int nPadX = 1;
    const unsigned int nPadY = 1;
    c->Divide(nPadX,nPadY);
    const unsigned int nPads=nPadX*nPadY;
    for(unsigned int i=1; i<=nPads; i++)
    {
	TPad* pad= (TPad*)c->cd(i);
	pad->SetTopMargin(0.10);
	pad->SetRightMargin(0.20);
	pad->SetLeftMargin(0.15);
    }
    // Open file
    TFile *f = TFile::Open(fullPath.c_str());
    if (f==0)
    {
	cout << "File " << fullPath << " not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0)
	cout << "Succesfully opened file " << fullPath << endl;
    // Get TTree
    TTree* t = (TTree*) f->Get("events");
    if(fVerbose>0) cout << "Got TTree with " << t->GetEntries() << " entries" << endl;
    // Do a cut, if needed
    //t->Draw(">>lst","chi2lb>.1&&mlb>5.61&&mlb<5.63");
    //t->Draw(">>lst","chi2lb>.1&&isSig==1");
    t->Draw(">>lst","(rid1m&4)==4&&(rid2m&4)==4&&mjp>2.895&&mjp<3.295&&prob1m>0.1&&prob2m>0.1&&ptjp>2&&probjp>0.005&&ml0>1.101&&ml0<1.129&&probpr>0.02&&probpi>0.02&&rptpr>rptpi&&ptl0>3&&rptpr>1&&rptpi>0.5&&probl0>0.02&&alphal0<0.3&&d3l0>1&&d3l0/d3El0>10&&problb>0.001&&alphalb<0.3");
    TEventList *lst;
    lst = (TEventList*)gDirectory->Get("lst");
    t->SetEventList(lst);
    if(fVerbose>0) cout << "Got TTree with " << t->GetEntries() << " entries" << endl;

    // Do plots
    c->cd(1);
    //doPlot2d(t,"hrzL0vtx", "vrl0:TMath::Abs(vzl0)",30,0,300,30,0,120,"Tit","|z|","r","cm","cm");
    if (custBinning)
    {
	double newbinsX[]={0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50};
	const int newbinsX_size = sizeof(newbinsX)/sizeof(double);
	std::vector<double> binvecX(newbinsX,newbinsX+newbinsX_size);
	//double newbinsY[]={0,1,2,3,4,5,6,7,8,9,10};
	double newbinsY[]={0,0.5,1,2,4,8,16,32};
	const int newbinsY_size = sizeof(newbinsY)/sizeof(double);
	std::vector<double> binvecY(newbinsY,newbinsY+newbinsY_size);
	doPlot2d(t,"hrzL0vtxreco", "vrl0:TMath::Abs(vzl0)",binvecX, binvecY,"#Lambda vertices","|z|","r","cm","cm");
    }
    else
    {
	doPlot2d(t,"hrzL0vtxreco", "vrl0:TMath::Abs(vzl0)",30,0,50,30,0,30,"#Lambda vertices","|z|","r","cm","cm");
    }

    // add tracker
    TPad* pad;
    pad = (TPad*)c->cd(1);
    pad->Modified();
    pad->Update();
    repositionPalette("hrzL0vtxreco");
    pad->Update();
    pad->SetLogz();
    drawTracker(pad);

    if (nOverlay<=0) return;
    int maxN = nOverlay;
    if (maxN > t->GetEntries()) maxN = t->GetEntries();
    double vrl0,vzl0,ppr,ppi,etapr,etapi;
    t->SetBranchAddress("vrl0",&vrl0);
    t->SetBranchAddress("vzl0",&vzl0);
    t->SetBranchAddress("ppr",&ppr);
    t->SetBranchAddress("etapr",&etapr);
    t->SetBranchAddress("ppi",&ppi);
    t->SetBranchAddress("etapi",&etapi);
    double scalepr = 4;
    double scalepi = 8;

    { // reference indicator
	const double x1pr = 0; const double y1pr = -3;
	const double x2pr = scalepr; const double y2pr = y1pr;
	const double versatz = 14;
	const double x1pi = x1pr+versatz; const double y1pi = -3;
	const double x2pi = x2pr+versatz+scalepi; const double y2pi = y1pi;
	TArrow *a;
	a = new TArrow(x1pr,y1pr,x2pr,y2pr,.01,">");
	a->SetLineColor(24);
	a->Draw();
	TLatex tl;
	tl.SetTextSize(20);
	tl.SetTextFont(4);
	tl.DrawLatex(x1pr,y2pr-1.2,"p(p) / 1 GeV");
	a = new TArrow(x1pi,y1pi,x2pi,y2pi,.01,">");
	a->SetLineColor(20);
	a->Draw();
//.........这里部分代码省略.........

// -----------------------------------------------------------------------------
//
void talk() {
  
  time_t start = TTimeStamp().GetSec();

  bool draw = true;
  bool debug = true;
  
  // Binning
  int xbins = 100;
  int ybins = 100;
  double xmin = 0.0;
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.0;
  double xbin_centre = ( ( xmax - xmin ) / xbins ) / 2.;
  double ybin_centre = ( ( ymax - ymin ) / ybins ) / 2.;
  
  // AlphaT values
  const int nat = 1;
  double at[nat];
  for ( int ii = 0; ii < nat; ++ii ) { at[ii] = 0.55 + ii * 0.001; } 
  
  // HT regions
  const int nht = 3;
  double ht_min[nht] = { 250., 300., 350. };

  // Jet pT thresholds
  double pt1_min[nht] = { 71.4, 85.7, 100. };
  double pt2_min[nht] = { 71.4, 85.7, 100. };
  double pt3_min[nht] = { 35.7, 42.9., 50. };

  // x fractions
  double x1_min[nht];
  double x2_min[nht];
  double x3_max[nht];
  for ( int ii = 0; ii < nht; ++ii ) { x1_min[ii] = ( 2. * pt1_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  for ( int ii = 0; ii < nht; ++ii ) { x2_min[ii] = ( 2. * pt2_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  for ( int ii = 0; ii < nht; ++ii ) { x3_max[ii] = ( 2. * pt3_min[ii]) / ( ht_min[ii] + pt3_min[ii] ); }
  
  // Loop through bins 
  int loop = 0;
  int nloops = nht;
  for ( int iht = 0; iht < nht; ++iht ) {
    std::cout << "Completed " 
	      << 100.*float(loop)/float(nloops) 
	      << "%..." 
	      << std::endl; 
    loop++;
    
    // Labeling
    std::stringstream ss;
    ss << "HT" << int(ht_min[iht]);
	  
    // New canvas for plots
    TCanvas* c1 = 0;
    if (draw) c1 = new TCanvas( TString("Canvas"+ss.str()), "" );
	  
    // Pad for cross-section plot
    TPad* pad = 0;
    if (draw) pad = new TPad(TString("Pad"+ss.str()),"",0.,0.,1.,1.);
    if (pad) {
      pad->SetGrid();
      pad->Draw();
      pad->cd();
      pad->SetLogz();
    }
    TH1F* hr = 0;
    if (draw) hr = pad->DrawFrame(0.,0.,1.,1.);

    // Histo title
    if (hr) {
      std::stringstream sss;
      sss << "H_{T}=" << ht_min[iht]
	  << "(p_{T1},p_{T2},p_{T3})=" 
	  << pt1_min[iht] << ","
	  << pt2_min[iht] << ","
	  << pt3_min[iht] << ")"
	  << ", (x_{1},x_{2},x_{3})=" 
	  << x1_min[iht] << ","
	  << x2_min[iht] << ","
	  << x3_max[iht] << ")";
      hr->SetTitle( sss.str().c_str() );
      hr->GetXaxis()->SetTitle( "x_{2}" );
      hr->GetYaxis()->SetTitle( "x_{1}" );
    }
	  
    // Create 2D cross-section plot
    TH2D* his = 0;
    if (draw) his = new TH2D(TString("Histo"+ss.str()),"",
			     xbins,xmin,xmax,
			     ybins,ymin,ymax);
	  
    // Fill 2D cross-section plot
    for ( int xbin = 0; xbin < xbins; ++xbin ) { 
      for ( int ybin = 0; ybin < ybins; ++ybin ) { 
	double x2 = ( ( xmax - xmin ) / xbins ) * xbin + xmin;
	double x1 = ( ( ymax - ymin ) / ybins ) * ybin + ymin;
	double val = ( x1*x1 + x2*x2 ) / ( ( 1 - x1 ) * ( 1 - x2 ) ); 
//.........这里部分代码省略.........

void EMC_HodoScope_Calibration(TString Cherenkov_Choice)
{
  TString cut_C(Cherenkov_Choice + " > 1000");
  TString cuts = Cherenkov_Choice;

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);

  TH1F *h1_HCalOut = new TH1F("h1_HCalOut",
                              "(CEMC MIP && HCal_{IN} MIP);EMCal + HCal_{IN} + HCal_{OUT} (GeV)", 100,
                              0, abs(beam_momentum_selection) * 2.);

  TH2 *h2_hodo_v_h =
      new TH2F("h2_hodo_v_h",
               "Hodoscope hit distribution;Average Horizontal Hodoscope Idx;Average Vertical Hodoscope Idx",
               8, -.5, 7.5, 8, -.5, 7.5);

  TH2 *h2_hodoE_v_h =
      new TH2F("h2_hodoE_v_h",
               "Event-Averaged Tower-Summed EMCal Energy [A.U.];Average Horizontal Hodoscope Idx;Average Vertical Hodoscope Idx",
               8, -.5, 7.5, 8, -.5, 7.5);
  TH2 *h2_hodoE_v_h_norm =
      new TH2F("h2_hodoE_v_h_norm",
               "Event-Averaged Tower-Summed EMCal Energy [A.U.];Average Horizontal Hodoscope Idx;Average Vertical Hodoscope Idx",
               8, -.5, 7.5, 8, -.5, 7.5);

  TH2 *h2_hodo_c_h =
      new TH2F("h2_hodo_c_h",
               "Column VS Horizontal Hodoscope;Average Horizontal Hodoscope Idx;Average EMCal Column",
               24, -.5, 7.5, 24, -.5, 7.5);
  TH2 *h2_hodo_r_v =
      new TH2F("h2_hodo_r_v",
               "Row VS Vertical Hodoscope;Average Vertical Hodoscope Idx;Average EMCal Row",
               24, -.5, 7.5, 24, -.5, 7.5);

  TH2 *h2_hodoE_h =
      new TH2F("h2_hodoE_h",
               "Total EMCal Energy VS Horizontal Hodoscope;Average Horizontal Hodoscope Idx;Total EMCal Energy [A.U.]",
               8, -.5, 7.5, 80, -.5, 70);
  TH2 *h2_hodoE_v =
      new TH2F("h2_hodoE_v",
               "Total EMCal Energy VS Vertical Hodoscope;Average Vertical Hodoscope Idx;Total EMCal Energy [A.U.]",
               8, -.5, 7.5, 80, -.5, 70);

  TText *t;
  TCanvas *c1 = new TCanvas("EMC_HodoScope_Calibration" + cuts,
                            "EMC_HodoScope_Calibration" + cuts, 1800, 1100);
  c1->Divide(3, 2);
  int idx = 1;
  TPad *p;

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

  T->Draw("Average_HODO_VERTICAL:Average_HODO_HORIZONTAL>>h2_hodo_v_h", cut_C, "colz");

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

  T->Draw("Average_column:Average_HODO_HORIZONTAL>>h2_hodo_c_h", cut_C, "colz");

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

  T->Draw("Energy_Sum_RAW_CEMC:Average_HODO_HORIZONTAL>>h2_hodoE_h", cut_C, "colz");

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

  T->Draw("Average_HODO_VERTICAL:Average_HODO_HORIZONTAL>>h2_hodoE_v_h", "(" + cut_C + ")*(Energy_Sum_RAW_CEMC) * (Energy_Sum_RAW_CEMC>0.1)", "goff");
  T->Draw("Average_HODO_VERTICAL:Average_HODO_HORIZONTAL>>h2_hodoE_v_h_norm", "(" + cut_C + " ) && (Energy_Sum_RAW_CEMC>0.1)", "goff");
  h2_hodoE_v_h->Divide(h2_hodoE_v_h_norm);
  h2_hodoE_v_h->Draw("colz");

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

  T->Draw("Average_row:Average_HODO_VERTICAL>>h2_hodo_r_v", cut_C, "colz");

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

  T->Draw("Energy_Sum_RAW_CEMC:Average_HODO_VERTICAL>>h2_hodoE_v", cut_C, "colz");

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

//.........这里部分代码省略.........
  //    p->SetLogx();
  p->SetLogy();

    {

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

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

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

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

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

      DrawReference(h_new, h_ref);
    }

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

  TH2F * h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection =
      (TH2F *) qa_file_new->GetObjectChecked(
          "h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection", "TH2F");
  assert(h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection);
  h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection->GetYaxis()->SetTitleOffset(
      1.5);
//  h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection->GetXaxis()->SetRangeUser(-5,
//      5);
//  h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection->GetYaxis()->SetRangeUser(-5,
//      5);
  h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection->Draw("COLZ");

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

    {

      TH2F * h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection =
          (TH2F *) qa_file_new->GetObjectChecked(
              "h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection", "TH2F");
      assert(h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection);

      TH1D * proj_new =
          h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection->ProjectionX(
              "qa_file_new_h_QAG4Sim_HCALIN_Cluster_LateralTruthProjection_px");
      proj_new->Rebin(4);

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

      TH1D * proj_ref = NULL;

// -----------------------------------------------------------------------------
//
void rob() {

  // Binning
  
  int xbins = 100;
  int ybins = 100;
  
  double xmax = 1.0;
  double xmin = 0.0;
//   double xrange = xmax - xmin;
//   xmax += xrange / xbins / 2;
//   xmin -= xrange / xbins / 2;
//   //xbins++;
  
  double ymax = 1.0;
  double ymin = 0.0;
//   double yrange = ymax - ymin;
//   ymax += yrange / ybins / 2;
//   ymin -= yrange / ybins / 2;
//   //ybins++;

  std::cout << " Binning: "
	    << " xbins: " << xbins
	    << " ybins: " << xbins
	    << " xmin: " << xmin
	    << " xmax: " << xmax
	    << " ymin: " << ymin
	    << " ymax: " << ymax
	    << std::endl;
  
  // Examples values of pt1, pt2, mht, x1, x2, x3, sigma and alpha_t
//   double pt1 = 50.;
//   double pt2 = 50.;
//   double mht = 50.;
//   double x1 = ( 2. * pt1 ) / ( pt1 + pt2 + mht );
//   double x2 = ( 2. * pt2 ) / ( pt1 + pt2 + mht );
//   double x3 = 2 - x1 - x2;
//   double sigma =  ( x1*x1 + x2*x2 ) / ( (1-x1) * (1-x2) );
//   double alpha_t =  x2 / ( 2 * sqrt(x1+x2-1) );
//   std::cout << " pt1: " << pt1
// 	    << " pt2: " << pt2
// 	    << " mht: " << mht 
// 	    << " x1: " << x1
// 	    << " x2: " << x2
// 	    << " x3: " << x3
// 	    << " sigma: " << sigma
// 	    << " alpha_t: " << alpha_t
// 	    << std::endl;
  
  // Cross section
  TCanvas* c1 = new TCanvas( "Contours", "" );
  //c1->SetGridx(1);
  //c1->SetGridy(1);

  TPad* pad = new TPad("pad","",0.,0.,1.,1.);
  pad->SetGrid();
  pad->Draw();
  pad->cd();
  pad->SetLogz(1);

  TH1F* hr = pad->DrawFrame(0.,0.,1.,1.);
  
  const int nx = 3;
  const int ny = 10;
  double pt[nx] = { 30., 50., 100. };
  double ht[ny];
  for ( int ii = 0; ii < ny; ++ii ) { ht[ii] = 150. + ii * 50.; }
  double ratio[nx][ny];
  
  double xbin_centre = ( ( xmax - xmin ) / xbins ) / 2.;
  double ybin_centre = ( ( ymax - ymin ) / ybins ) / 2.;
  TH2D* his = new TH2D("Contours","",xbins,xmin,xmax,ybins,ymin,ymax);
  for ( int ii = 0; ii < nx; ++ii ) {
    for ( int jj = 0; jj < ny; ++jj ) {
      double x3 = ( 2. * pt[ii] ) / ( ht[jj] + pt[ii] );
      double n = 0.;
      double d = 0.;
      for ( int xbin = 0; xbin < xbins; ++xbin ) { 
	for ( int ybin = 0; ybin < ybins; ++ybin ) { 
	  double x2 = ( ( xmax - xmin ) / xbins ) * xbin + xmin;
	  double x1 = ( ( ymax - ymin ) / ybins ) * ybin + ymin;
	  double val = ( x1*x1 + x2*x2 ) / ( ( 1 - x1 ) * ( 1 - x2 ) ); 
	  if ( x1 < x2 ||      // jet ordering by Pt
	       x1 + x2 > 2. || // from relation "x1 + x2 + x3 = 2"
	       x1 > 1.0 ||     // from "lost jet" and relation "xmiss = -x1 -x2"
	       x1 + x2 < 1.    // from "lost jet" and relation "xmiss = -x1 -x2"
	       ) { continue; }
	  if ( ( x1 + x2 ) < ( 2 - x3 ) ) { continue; }
	  d += val;
	  double alpha_t =  x2 / ( 2 * sqrt(x1+x2-1) );
	  if ( alpha_t > 0.5 ) n+= val;
	  if ( ii == 2 && jj == 9 ) { his->Fill( x2+xbin_centre, x1+ybin_centre, val ); }
	}
      }
      double r = 0.;
      if ( d > 0. ) { r = n/d; }
      ratio[ii][jj] = r;
      std::cout << " Pt: " << pt[ii]
//.........这里部分代码省略.........