C++ TH2D::DrawCopy方法代码示例

TCanvas* pHitSpecPos( )
{
  TCanvas* c = new TCanvas("cHitSpecPos","cHitSpecPos",1000,1100);
  c->Divide(2,2);
  TVirtualPad* p; TH2D *h;
  TText t; t.SetTextColor(4);

  p =c->cd(1); p->SetLogy(); p->SetGrid(1,0);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhi_BarMinus"); 
  h->GetYaxis()->SetRange(4,33);
  h->DrawCopy("box");
  t.DrawTextNDC(0.17,0.15, "BARREL MU MINUS");

  p =c->cd(2); p->SetLogy(); p->SetGrid(1,0);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhi_BarPlus"); 
  h->GetYaxis()->SetRange(4,32);
  h->DrawCopy("box");
  t.DrawTextNDC(0.17,0.15, "BARREL MU PLUS");

  p =c->cd(3); p->SetLogy(); p->SetGrid(1,0);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhi_EndMinus"); 
  h->GetYaxis()->SetRange(4,32);
  h->DrawCopy("box");
  t.DrawTextNDC(0.17,0.15, "ENDCAP MU MINUS");

  p =c->cd(4); p->SetLogy(); p->SetGrid(1,0);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhi_EndPlus"); 
  h->GetYaxis()->SetRange(4,32);
  h->DrawCopy("box");
  t.DrawTextNDC(0.17,0.15, "ENDCAP MU PLUS");

  return c;
}

void PlotImages()
{
gStyle->SetOptStat(0);
gStyle->SetPadRightMargin(0.2);
string filename[1];
char buffer[100];
char buffertext[100];

for (int seed=55; seed<56; seed++)
  {
    int n=sprintf(buffer, "/home/les67/locust_faketrack_waterfall.root");

    const char *file = buffer;

    TFile* f = TFile::Open(file);
    if (!((!f)||f->IsZombie()))
    {

    TH2D* hspectrogram = GetSpectrogram(file);
    TH2D* hlabels = GetLabels(hspectrogram, 40e-21);  // threshold for labeling goes here.
//    TGraph* grlabels = GetLabelGraph(hlabels);
//    PrintLabels(grlabels);  // print labels to terminal.

    hspectrogram->GetXaxis()->SetRangeUser(0.,0.02);
    hspectrogram->GetYaxis()->SetRangeUser(149.e6, 170.e6);

    TCanvas *c = new TCanvas;
    n=sprintf(buffertext, "hspectrogram_%d.png", seed);
    const char *pngname = buffertext;
    hspectrogram->GetYaxis()->SetTitleOffset(1.25);
    hspectrogram->DrawCopy("colz");

/*
    grlabels->SetMarkerColor(2);
    grlabels->SetMarkerStyle(8);
    grlabels->SetMarkerSize(0.6);
    grlabels->SetLineWidth(3.);
    grlabels->SetLineColor(2);
    grlabels->Draw("psame");
*/


 /*  
    TImage *img = TImage::Create();
    img->FromPad(c);
    img->WriteImage(pngname);
    delete img;
    delete hspectrogram;
    delete c;
*/
    f->Close();
    }

  }

}

TCanvas* pHitSpecDet(const std::string& s)
{
  TCanvas* c = new TCanvas( ("HitSpecDet"+s).c_str(), ("HitSpecDet"+s).c_str(), -2);
  c->SetLogx();
  TH2D* h = (TH2D*)gROOT->FindObject(s.c_str());
  h->GetXaxis()->SetRange(4,32);
  std::stringstream str;
  str<<h->GetTitle();
  TText t;
  t.SetTextColor(4);
  t.SetTextAlign(11);
  t.SetTextSize(0.035);
  h->DrawCopy("box");
  t.DrawTextNDC(0.17,0.2, str.str().c_str());
  return c;
}

//.........这里部分代码省略.........
				nTest++;
			}while(!hit);
			nEvent++;
		}
		MCnTry->Fill(1.0*nTry/EVENTS);
	}
	clock0->Stop();
	
	delete fRandomEngine;
	cout<<nEvent<<" events have been generated."<<endl;
	cout<<"There is a total of "<<nInside<<"/"<<nTest<<" muons that passed close enough to the detector."<<endl;
	cout<<"Including the given weight this corresponds to ";
	cout<<kmax*weighttest/xdist/zdist*10000/123.3044<<" spills (1 spill = "<<xdist*zdist*123.3044/10000;
	cout<<" real cosmic muons = "<<EVENTS<<" simulated events)."<<endl;
   
   cout<<weighttest<<endl;
   clock0->Print();
   
   Double_t meanflux = 0;
   Int_t binsum = 0;
   for (Int_t ix = 2; ix<29;ix++){
	   for (Int_t iz = 2; iz<89;iz++){
		 	binsum++;
		   meanflux += MCXZ->GetBinContent(ix,iz);
		}
	}
	cout<< "meanflux: "<<meanflux/binsum<<"  "<< meanflux<<endl<<endl;
   printf("MCnTry: %.9f +- %.9f",MCnTry->GetMean(),MCnTry->GetMeanError());
   cout<<endl<<endl;
     
    TCanvas *c1 = new TCanvas("c1","c1",400,400);
	 c1->Divide(1,1);
	 c1->cd(1);
	 MCnTry->DrawCopy();
	
	 TCanvas *c4 = new TCanvas("c4","c4",400,400);
	 c4->Divide(1,1);
	 c4->cd(1);
	 StartPTheta->DrawCopy("SURF2");
	 gPad->SetLogx();
	
	 //TCanvas *c2 = new TCanvas("c2","c2",400,400);
	 //c2->Divide(1,1);
	 //c2->cd(1);
	 //gPad->SetLogy();
	 //wei->DrawCopy();
	
    TCanvas *c3 = new TCanvas("c3","c3",1600,800);
	 c3->Divide(4,2);
	 c3->cd(1);
	 StartXZ->DrawCopy("COLZ");
	 c3->cd(2);
	 //MCP->SetLineColor(kGreen);
	// MCP->DrawCopy();
	 //StartP->DrawCopy();
	 //TF1 *fs = new TF1("fs",NEWs,1,100,2);
	 //fs->FixParameter(0, 0);
	 //fs->FixParameter(1, 500);
	 //StartP->Add(fs,-1);
	 StartP->DrawCopy();
	 //StartP->Fit(fs,"I");
	 
	 
	// 
	 //fs->DrawCopy("SAME");
	 

void scan2() {

    /*  TFile f1("hist.root");
      TFile f2("SysTot.root");

      TH1D *datahist = (TH1D*)f1.Get("datahist");
      TH1D *fithist = (TH1D*)f1.Get("fithist");

      TH1D *totlow = (TH1D*)f2.Get("totlow");
      TH1D *tothigh = (TH1D*)f2.Get("tothigh");
    */

    TF1 *fun1 = new TF1("fun",fun,0,1000,3);

    TFile f("Output.root");
    TH1F *datahist = dynamic_cast<TH1F*> (f.Get("InvMass"));
    TH1F *fithist = dynamic_cast<TH1F*> (f.Get("hist_err"));
    TH1F *tothigh = dynamic_cast<TH1F*> (f.Get("SystPlus"));
    TH1F *totlow = dynamic_cast<TH1F*> (f.Get("SystMinus"));
    new TCanvas();
    datahist->DrawCopy();
    fithist->DrawCopy("same");
    gPad->SetEditable(0);

    /*	if (datahist->IsZombie()) { std::cout << "Err! datahist not found " << std::endl; return; }
    	if (fithist->IsZombie()) { std::cout << "Err! err_hist not found " << std::endl; return; }
    	if (tothigh->IsZombie()) { std::cout << "Err! syst_plus hist not found " << std::endl; return; }
    	if (totlow->IsZombie()) { std::cout << "Err! syst_minus hist not found " << std::endl; return; }
    */

    TH1D *p = new TH1D("p","Probability",100,150,900);

    for (int i = 1; i<=fithist->GetNbinsX(); i++)
    {
        double mass = fithist->GetBinCenter(i);

        if (mass>440 && mass<450)
        {
            double sigma = sqrt(2.0)*sqrt(pow(0.135*sqrt(mass/2.0),2)+
                                          pow(0.02*(mass/2.0),2));
            cout<< " ===================== mass +/- sigma = " << mass<<"+/-"<<sigma<<endl;
            int bin1 = fithist->FindBin(mass-sigma/2);
            int bin2 = fithist->FindBin(mass+sigma/2);
            cout<<mass<<" "<<bin1<<" "<<bin2<<endl;
            double data = 0;
            double bg = 0;
            double err = 0;

            for (int j = bin1; j<=bin2; j++)
            {
                data+=datahist->GetBinContent(j);
                bg+=fithist->GetBinContent(j);
                double err1 = -totlow->GetBinContent(j);
                double err2 = tothigh->GetBinContent(j);
                err+=TMath::Max(err1,err2)*bg; //why multiply by bg???
            }
            cout << "Total Data/Bg+/-err in mass window[" << mass << "] = "<< data <<"/ "<< bg << "+/-" << err <<endl;
            double prob = 0;
            fun1->SetParameter(0,bg);
            fun1->SetParameter(1,err);

            for (int j = int(data+0.001); j<100; j++) {
                fun1->SetParameter(2,j);
                //fun1->Print();
                //cout << "Evaluating Intrgral for j = " << j << " from x0= " << TMath::Max(0.0,bg-10*err) << " to x1 = " << bg+10*err << endl;
                double val = fun1->Integral(TMath::Max(0.0,bg-10*err),bg+10*err);
                //double val = fun1->Integral(TMath::Max(0.0,bg-2*err),bg+2*err);
                /*for (int z=TMath::Max(0.0,bg-2*err); z < bg+2*err; ++z)
                {
                	if (c<4)
                	{
                	std::cout << "func at [" << z << "]=" << 	fun1->Eval(z) << std::endl;
                	}
                }
                */
                prob += val;
            }
            cout<< "Prob for mass[" << mass<<"]="<< prob <<endl;
            p->SetBinContent(p->FindBin(mass),prob);
        }
    }
    /*
    delete gRandom;
    gRandom = (TRandom*) new TRandom3;
    gRandom->SetSeed(3);

    TH1D *minp = new TH1D("minp","Minimum Probability of Each PseudoExpt",100,0,0.2);

    //int nexp = 50000;
    int nexp = 10;

    TH1D *htemp = (TH1D*)datahist->Clone("htemp");

    for (int iexp = 0; iexp<nexp; iexp++){
      //if (iexp%10==0) cout<<iexp<<endl;
      //generate pseudo-experiments
      htemp->Reset();
      for (int i = 1; i<=htemp->GetNbinsX(); i++){
        double mass = htemp->GetBinCenter(i);
        if (mass>150&&mass<650){
//.........这里部分代码省略.........

int genbod_diff(const char* filename1, const char* filename2, bool ifSave = false)
{
  TFile* f1 = new TFile(filename1,"READ");
  TFile* f2 = new TFile(filename2,"READ");

  TH2D* numF1 = ((TH2D*) f1->Get("cnumepNonIdEP"));
  TH2D* numF2 = ((TH2D*) f2->Get("cnumepNonIdEP"));
  TH2D* denF1 = ((TH2D*) f1->Get("cdenepNonIdEP"));
  TH2D* denF2 = ((TH2D*) f2->Get("cdenepNonIdEP"));
  TH2D* numB1 = ((TH2D*) f1->Get("cnumepNonIdEPTrue"));
  TH2D* numB2 = ((TH2D*) f2->Get("cnumepNonIdEPTrue"));
  TH2D* denB1 = ((TH2D*) numF1->Clone());
  TH2D* denB2 = ((TH2D*) numF2->Clone());
  TH2D* numC1 = ((TH2D*) numB1->Clone());
  TH2D* numC2 = ((TH2D*) numB2->Clone());
  TH2D* denC1 = ((TH2D*) f1->Get("cdenepNonIdEPTrue"));
  TH2D* denC2 = ((TH2D*) f2->Get("cdenepNonIdEPTrue"));

  double eta = 1.1;
  gStyle->SetOptStat(111);
  //TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,800,600);
   TCanvas* canv = new TCanvas("canv", "GENBOD results", 10,10,1600,1200);
  canv->Divide(3,3);
  //TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,800,600);
  TCanvas* canvND = new TCanvas("canvND", "liczniki i mianowniki ", 10,10,1600,1200);
  canvND->Divide(2,2);

  //////   Pure Corr fun
  TH2D* numdiff = (TH2D*)numF1->Clone();
  numdiff->Divide(numF2);
  canvND->cd(1); numdiff->DrawCopy("colz");
  TH2D* dendiff = (TH2D*)denF1->Clone();
  dendiff->Divide(denF2);
  canvND->cd(2); dendiff->DrawCopy("colz");
  double sF1 = numF1->GetEntries()/denF1->GetEntries();
  numF1->Divide(denF1);
  numF1->Scale(1./sF1);
  numF1->GetYaxis()->SetRangeUser(-eta,eta);
  double sF2 = numF2->GetEntries()/denF2->GetEntries();
  numF2->Divide(denF2);
  numF2->Scale(1./sF2);
  numF2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(1);
  TH2D* corrF1 = (TH2D*) numF1->Clone();
  corrF1->Draw("surf1");
  DrawNicely(corrF1, 0,0, filename1);
  canv->cd(2);
  TH2D* corrF2 = (TH2D*) numF2->Clone();
  corrF2->Draw("surf1");
  DrawNicely(corrF2, 0,0, filename2);
  canv->cd(3);
  numF1->Divide(numF2);
  DrawNicely(numF1, 0,0, "Femtoscopic component");
  numF1->Draw("colz");
  
  ////// Pure Bkg
  TH2D* numdiffTrue = (TH2D*)numB1->Clone();
  numdiffTrue->Divide(numB2);
  canvND->cd(3); numdiffTrue->DrawCopy("colz");
  double sB1 = numB1->GetEntries()/denB1->GetEntries();
  numB1->Divide(denB1);
  numB1->Scale(1./sB1);
  numB1->GetYaxis()->SetRangeUser(-eta,eta);
  double sB2 = numB2->GetEntries()/denB2->GetEntries();
  numB2->Divide(denB2);
  numB2->Scale(1./sB2);
  numB2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(4);
  TH2D* corrB1 = (TH2D*) numB1->Clone();
  corrB1->Draw("surf1");
  DrawNicely(corrB1, 0,0, filename1);
  canv->cd(5);
  TH2D* corrB2 = (TH2D*) numB2->Clone();
  corrB2->Draw("surf1");
  DrawNicely(corrB2, 0,0, filename2);
  canv->cd(6);
  numB1->Divide(numB2);
  DrawNicely(numB1, 0,0, "Background component");
  numB1->DrawCopy("colz");
  
  //////   Norm Corr fun
  TH2D* dendiffTrue = (TH2D*)denC1->Clone();
  dendiffTrue->Divide(denC2);
  canvND->cd(4); dendiffTrue->DrawCopy("colz");
  double sC1 = numC1->GetEntries()/denC1->GetEntries();
  numC1->Divide(denC1);
  numC1->Scale(1./sC1);
  numC1->GetYaxis()->SetRangeUser(-eta,eta);
  double sC2 = numC2->GetEntries()/denC2->GetEntries();
  numC2->Divide(denC2);
  numC2->Scale(1./sC2);
  numC2->GetYaxis()->SetRangeUser(-eta,eta);
  canv->cd(7);
  TH2D* corrC1 = (TH2D*) numC1->Clone();
  corrC1->Draw("surf1");
  DrawNicely(corrC1, 0,0, filename1);
  canv->cd(8);
  TH2D* corrC2 = (TH2D*) numC2->Clone();
  corrC2->Draw("surf1");
  DrawNicely(corrC2, 0,0, filename2);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
 
  for (int i = startChip; i < startChip+nChips; i++) { addVcalThreshold(dirName, i, mThreshold); }
  for (int i = startChip; i < startChip+nChips; i++) { addChip("vcals_xtalk", i, mBumps); }

  TFile *f1 = new TFile(Form("%s/%s", dirName, adFileName));
  for (int i = startChip; i < startChip+nChips; i++) { addChip("AddressDecoding", i, mAddr);}

  if ( nChips < 16 && startChip == 0 ) { 
    for (int i = 8; i < nChips+8; i++) { removeChip(i, mThreshold, -99); }
    for (int i = 8; i < nChips+8; i++) { removeChip(i, mBumps, -99); }
    for (int i = 8; i < nChips+8; i++) { removeChip(i, mAddr, -99); }
  }

  if ( nChips < 16 && startChip == 8 ) { 
    for (int i = 0; i < nChips; i++) { removeChip(i, mThreshold, -99); }
    for (int i = 0; i < nChips; i++) { removeChip(i, mBumps, -99); }
    for (int i = 0; i < nChips; i++) { removeChip(i, mAddr, -99); }
  }

  TString noslash(dirName);
  noslash.ReplaceAll("/", "");
  noslash.ReplaceAll("..", "");
  
  c1->cd(1);

  if ( mThreshold->GetMaximum() < mThresholdmax ) { 
    mThresholdmax = mThreshold->GetMaximum();
  }
  if ( mThreshold->GetMinimum() > mThresholdmin ) {
    mThresholdmin = mThreshold->GetMinimum();
  }
  mThreshold->GetZaxis()->SetRangeUser(mThresholdmin,mThresholdmax);
  mThreshold->DrawCopy("colz");
  tl->DrawLatex(0.1, 0.92, "Vcal threshold");
  tl->DrawLatex(0.75, 0.92, Form("%s",noslash.Data()));

  if ( nChips < 16 && startChip == 0 ) { 
   
    box->SetFillColor(29);
    box->DrawBox( 0, 0,  416,  80);
  }

  if ( nChips < 16  && startChip == 8 ) { 
   
    box->SetFillColor(29);
    box->DrawBox( 0, 80,  416,  160);
  }



  c1->cd(2);
  mBumps->SetMaximum(2.);
  mBumps->SetMinimum(-2.);  
  mBumps->DrawCopy("colz");
  tl->DrawLatex(0.1, 0.92, "Bump bonding map");

  if ( nChips < 16 && startChip == 0 ) { 
   
    box->SetFillColor(29);
    box->DrawBox( 0, 0,  416,  80);
  }

  if ( nChips < 16 && startChip == 8 ) { 
   
    box->SetFillColor(29);

//.........这里部分代码省略.........
    if (i == nThetaHistos)  {
      hTemp = hPull;
      thetaString = "tan(#Theta) integrated";
    }
    else {
      hTemp = hPullTheta[i];
      thetaString = Form("%.2f #leq |tan(#Theta)| < %.2f", tThetaLow[i], tThetaHigh[i]);
    }
    
    normaliseHisto(hTemp);
    hTemp->FitSlicesY();
    hTemp->GetYaxis()->SetNdivisions(12);
    hTemp->GetXaxis()->SetMoreLogLabels(kTRUE);
    TH1D* hTempMean = (TH1D*)gDirectory->Get(Form("%s_1", hTemp->GetName()));
    hTempMean->SetTitle(Form("mean(pull), %s", thetaString.Data()));
    hTempMean->GetXaxis()->SetMoreLogLabels(kTRUE);
    hTempMean->SetLineWidth(2);
    hTempMean->SetMarkerStyle(20);
    TH1D* hTempSigma = (TH1D*)gDirectory->Get(Form("%s_2", hTemp->GetName()));
    hTempSigma->SetTitle(Form("#sigma(pull), %s", thetaString.Data()));
    hTempSigma->GetXaxis()->SetMoreLogLabels(kTRUE);
    hTempSigma->SetLineColor(kMagenta);
    hTempSigma->SetMarkerStyle(20);
    hTempSigma->SetMarkerColor(kMagenta);
    hTempSigma->SetLineWidth(2);
    TH1D* hTempChi2 = (TH1D*)gDirectory->Get(Form("%s_chi2", hTemp->GetName()));
    hTempChi2->SetTitle(Form("#chi^{2} / NDF (pull), %s", thetaString.Data()));
    hTempChi2->GetXaxis()->SetMoreLogLabels(kTRUE);
    hTempChi2->SetLineColor(kMagenta + 2);
    hTempChi2->SetMarkerStyle(20);
    hTempChi2->SetMarkerColor(kMagenta + 2);
    hTempChi2->SetLineWidth(2);

    hTemp->DrawCopy("colz");
    hTempMean->DrawCopy("same");
    hTempSigma->DrawCopy("same");
    hTempChi2->Scale(-1./10.);
    hTempChi2->DrawCopy("same");
    hTempChi2->Scale(-10.);

    canvPullMean->cd();
    hTempMean->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempMean->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempMean->DrawCopy((i == 0 ? "" : "same"));

    canvPullSigma->cd();
    hTempSigma->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempSigma->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempSigma->DrawCopy((i == 0 ? "" : "same"));

    canvPullChi2->cd();
    hTempChi2->SetLineColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempChi2->SetMarkerColor(1 + ((j >= 9) ? (39 + 2 * (j - 9)) : j));
    hTempChi2->DrawCopy((i == 0 ? "" : "same"));
  }

  canvPullMean->BuildLegend();
  canvPullSigma->BuildLegend();
  canvPullChi2->BuildLegend();
*/
  // Histograms with additional correction
  TCanvas* canvPullMeanCorr = 0x0;
  TCanvas* canvPullSigmaCorr = 0x0;
  TCanvas* canvPullChi2Corr = 0x0;
  TCanvas* canvPullCorr[nThetaHistos + 1];
  for (Int_t i = 0; i < nThetaHistos + 1; i++) 

TCanvas* pHitSpecPosGen( )
{
  TCanvas* c = new TCanvas("cHitSpecPosGen","cHitSpecPosGen",1200,600);
  c->Divide(2,1);
  TVirtualPad* p; TH2D *h; TH2D *hAtPhi0;
  TText t; t.SetTextColor(4);

  p =c->cd(1); p->SetLogy(); p->SetGrid(0,1);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGen_Minus"); 
  hAtPhi0 = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGenAtPhi0_Minus"); 
  h->GetYaxis()->SetRange(4,33);
  h->GetXaxis()->SetNdivisions(505); h->GetXaxis()->SetLabelSize(0.04);
  h->DrawCopy("box");
  hAtPhi0->SetLineColor(2);
  hAtPhi0->DrawCopy("box same");
  t.DrawTextNDC(0.17,0.15, "BARREL MU MINUS");
  std::cout <<h->GetTitle() << std::endl;
  for (int iy = 1; iy <=h->GetNbinsY(); iy++) {
  std::cout <<" pt: " << h->GetYaxis()->GetBinLowEdge(iy);
  double xmin=100.;
  double xAtMax = 0.; double valAtMax = 0.;
  double xmax=-100.;
  for (int ix = 1; ix <=h->GetNbinsX(); ix++) {
     double val = h->GetBinContent(ix,iy); 
     double xbmin = h->GetXaxis()->GetBinLowEdge(ix);
     double xcent = h->GetXaxis()->GetBinCenter(ix);
     double xbmax = h->GetXaxis()->GetBinUpEdge(ix);
     if (val > valAtMax) { valAtMax = val; xAtMax = xcent; }
     if (val > 1 && (xbmin < xmin)) xmin = xbmin;    
     if (val > 1 && (xbmax > xmax)) xmax = xbmax;    
  }
  std::cout <<"  set DPHI0 = "<<1.025-xAtMax 
            <<"; set DPHI_MARGIN = " <<std::max( (xAtMax-xmin), (xmax-xAtMax))
<<";" << std::endl;
  }

  p =c->cd(2); p->SetLogy(); p->SetGrid(1,1);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGen_Plus"); 
  hAtPhi0 = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGenAtPhi0_Plus"); 
  h->GetYaxis()->SetRange(4,32);
  h->GetXaxis()->SetNdivisions(505); h->GetXaxis()->SetLabelSize(0.04);
  h->DrawCopy("box");
  hAtPhi0->SetLineColor(2);
  hAtPhi0->DrawCopy("box same");
  t.DrawTextNDC(0.17,0.15, "BARREL MU PLUS");
  std::cout <<h->GetTitle() << std::endl;
  for (int iy = 1; iy <=h->GetNbinsY(); iy++) {
  std::cout <<" pt: " << h->GetYaxis()->GetBinLowEdge(iy);
  double xmin=100.;
  double xAtMax = 0.; double valAtMax = 0.;
  double xmax=-100.;
  for (int ix = 1; ix <=h->GetNbinsX(); ix++) {
     double val = h->GetBinContent(ix,iy); 
     double xbmin = h->GetXaxis()->GetBinLowEdge(ix);
     double xcent = h->GetXaxis()->GetBinCenter(ix);
     double xbmax = h->GetXaxis()->GetBinUpEdge(ix);
     if (val > valAtMax) { valAtMax = val; xAtMax = xcent; }
     if (val > 1 && (xbmin < xmin)) xmin = xbmin;    
     if (val > 1 && (xbmax > xmax)) xmax = xbmax;    
  }
  std::cout <<"  set DPHI0 = "<<1.025-xAtMax 
            <<"; set DPHI_MARGIN = " <<std::max( (xAtMax-xmin), (xmax-xAtMax))
<<";" << std::endl;
  }

  return c;
}

//.........这里部分代码省略.........
				// Find next regular term only
				bool foundRegTerm = false;
				for (int jTerm = iTerm + 1; jTerm < nTerms && !foundRegTerm; ++jTerm) {
					const Student::Enrollment jEnrollment = student->Enrollments()[jTerm];
					if (!MyFunctions::regularSemester(jEnrollment.term)) continue;
					foundRegTerm = true;
					for (Student::Grade jGrade : jEnrollment.grades) {
						if (jGrade.course == iGrade.course) continue;
						if (!MyFunctions::ValidGrade(jGrade.grade)) continue;
						
						double prediction_j = student->CourseGradePrediction(jGrade, Student::DISTRIBUTION);
						double delta_j = jGrade.quality - prediction_j;
						
						corrMap[std::make_pair(iGrade.course, jGrade.course)].Add(delta_i, delta_j);
						++nPairAll;
						
					}
				}
			} 
		}
		
	}
	myBenchmark->Stop("Main Loop");
	
	std::cout << "nPairAll     = " << nPairAll << std::endl;
	std::cout << "Unique Pairs = " << corrMap.size() << std::endl;
	
	TH1D* rHist = new TH1D("rHist", "Correlation Coefficient, #rho", 120, -1.2, 1.2);
	TH1D* pHist = new TH1D("pHist", "Probablity Distribution", 100, 0., 1.);
	TH1D* nHist = new TH1D("nHist", "Number of entries", 100, 0., 2000.);
	TH2D* pVrHist = new TH2D("pVrHist", "Prob vs. #rho", 100, -1., 1., 100, 0., 1.);
	
	myBenchmark->Start("Prune");
	for (auto iter = corrMap.begin(); iter != corrMap.end();) {
		if (iter->second.n() < nCut) {
			corrMap.erase(iter++);
			continue;
		}
		double p = iter->second.p();
		double r = iter->second.r();

		// Test for nan?
		if (p != p) {
			std::cout << "Found p = nan: n = " << iter->second.n() << std::endl;
			corrMap.erase(iter++);
			continue;
		}
		if (p < 0.) {
			corrMap.erase(iter++);
			continue;
		}
		rHist->Fill(r);
		pHist->Fill(p);
		nHist->Fill(iter->second.n());
		pVrHist->Fill(r, p);
		if (p < 1. - prob && p > prob) {
			corrMap.erase(iter++);
		}
		else {
			std::cout << "r = " << iter->second.r() << ", p = " << p << std::endl;
			++iter;
		}
	}
	myBenchmark->Stop("Prune");
	std::cout << "Post Cut     = " << corrMap.size() << std::endl;
	
	myBenchmark->Start("Sort");
	std::vector<std::pair<std::pair<TString, TString>, CorrelationCalculator>> corrVec(corrMap.begin(), corrMap.end());
	std::sort(corrVec.begin(), corrVec.end(), &sortFunc);
	myBenchmark->Stop("Sort");
	
	int printTop = 50;
	int printed = 0;
	for (auto const& entry : corrVec) {
		std::cout << entry.first.first << " : " << entry.first.second << "\t, n = " << entry.second.n() << "\t, r = " << entry.second.r() 
			<< "\t, p = " << entry.second.p() << std::endl;
		++printed;
		if (printed >= printTop) break;
	}

	TCanvas* c1 = new TCanvas("c1", "Grade Correlation", 1600, 1200);
	c1->Divide(2,2);
	c1->cd(1);
	TF1* myGaus = new TF1("myGaus", "gaus", -1., 1.);
	myGaus->SetParameters(600., 0., 0.2);
	myGaus->FixParameter(1, 0.);
	rHist->Fit(myGaus, "0B", "", -1., 0.);
	rHist->DrawCopy();
	myGaus->DrawCopy("SAME");
	c1->cd(2);
	pHist->DrawCopy();
	c1->cd(3);
	nHist->DrawCopy();
	c1->cd(4);
	pVrHist->DrawCopy();
	
	float rt, cp;
	myBenchmark->Summary(rt, cp);
	
}