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

void spectrum(char *filename = "example.root")
{
   TFile *infile = new TFile(filename);
   TTree *data = (TTree*) infile->FindObjectAny("data");
   
   TH1F *h[10];
   
   TCanvas *c = new TCanvas("c","",960,480);
   TPad *p = new TPad("p","",0,0,1,1);
   p->Divide(4,2,0.01,0.01);
   p->Draw();

   double norm;
   for (int i=0;i<8;i++)
   { 
      p->cd(i+1);
      p->GetPad(i+1)->SetLogy();
      h[i] = new TH1F(Form("h%d",i),Form("Ch %d",i),60,-0.1,1.1);
//      data->Draw(Form("-min[%d]>>h%d",i,i),"min[0]<min[1]&&min[0]<min[2]&&min[0]<min[3]&&min[0]<min[4]&&min[0]<min[5]&&min[0]<min[6]&&min[0]<min[7]");
      data->Draw(Form("-min[%d]/256./0.86>>h%d",i,i));
      if (i==0) norm = h[i]->GetEntries();
      h[i]->SetTitle(Form("Ch %d",i));
      h[i]->SetXTitle("Pulse height (Volt)");
      h[i]->SetYTitle("#");
      h[i]->SetMarkerSize(0.5);
      h[i]->SetMarkerColor(2);
      h[i]->Sumw2();
      h[i]->Scale(1./norm);
      h[i]->SetAxisRange(1e-7,10,"Y");
      h[i]->Draw("e");
   }
}

void RDK2AnalysisPlotter::makeEPGPlot(CoDet detType)
{
     gROOT->cd();
    TCanvas* theCanvas;
    TPad* mainPad;
    TPad* titlePad;
    TPad* gEPad;
    TPad* detPad;
    TPad* gEPadSubs[3];


    int numTitleLines=getTitleBoxLines( detType);
    int canvasNumPixelsYPlot=600;
    int canvasNumPixelsYPlots=2*canvasNumPixelsYPlot;
    int canvasNumPixelsYTitle=40*numTitleLines;
    int canvasNumPixelsY=canvasNumPixelsYPlots+canvasNumPixelsYTitle;

    theCanvas = new TCanvas("EPGExpMCAnalysisComparisonPlot","EPGExpMCAnalysisComparisonPlot",10,10,1200,canvasNumPixelsY);
     mainPad=(TPad*) theCanvas->GetPad(0);


    double ylow,yhigh;

    yhigh=1;
    ylow=1.-canvasNumPixelsYTitle/(double)canvasNumPixelsY;
    titlePad=new TPad("titlePad", "titlePad", 0., ylow, 1., yhigh, -1, 1, 1);

    yhigh=ylow;
    ylow=ylow-canvasNumPixelsYPlot/(double)canvasNumPixelsY;
    gEPad=new TPad("gEPad", "gEPad", 0., ylow, 1., yhigh, -1, 1, 1);

    yhigh=ylow;
    ylow=0;
    detPad=new TPad("detPad", "detPad", 0., ylow, 1., yhigh, -1, 1, 1);

    titlePad->Draw();
    gEPad->Draw();
    detPad->Draw();



//    titlePad=(TPad*) theCanvas->GetPad(1);
//    gEPad=(TPad*) theCanvas->GetPad(2);
//    detPad=(TPad*) theCanvas->GetPad(3);
//    egTPad=(TPad*) theCanvas->GetPad(4);

    titlePad->SetMargin(0.05,0.0,0.0,0.0);
    gEPad->SetMargin(0.05,0.0,0.0,0.05);

    if(numExp>1 && numMC>1)
    {
        gEPad->Divide(2);

        gEPadSubs[0]=(TPad*) gEPad->GetPad(1);

        TPad* tempPad;
        tempPad=(TPad*) gEPad->GetPad(2);
        tempPad->Divide(1,2);
        gEPadSubs[1]=(TPad*) tempPad->GetPad(1);
        gEPadSubs[2]=(TPad*) tempPad->GetPad(2);


        TPad* allBasePads[5]={titlePad,gEPadSubs[0],gEPadSubs[1],gEPadSubs[2],detPad};
         for (int i = 1;i< 5;i++)
        {
            allBasePads[i]->SetGrid(1,1);
            allBasePads[i]->SetTickx(1);
            allBasePads[i]->SetTicky(1);
        }

        gEPadSubs[2]->SetLogx();

        gEPadSubs[0]->SetMargin(0.15,0.1,.1,0.1);
        gEPadSubs[1]->SetMargin(0.1,0.1,.1,0.1);
        gEPadSubs[2]->SetMargin(0.1,0.1,.1,0.1);

        //Plot gE

        gEPadSubs[0]->cd(); drawPlot(detType, PLOTVAR_GE,PLOT_COMP);
        gEPadSubs[1]->cd(); drawPlot(detType, PLOTVAR_GE,PLOT_RESID);
        gEPadSubs[2]->cd(); drawPlot(detType, PLOTVAR_GEVAR,PLOT_NORMRESID);

    }
    else
    {

        TPad* allBasePads[3]={titlePad,gEPad,detPad};
        for (int i = 1;i< 3;i++)
        {
            allBasePads[i]->SetGrid(1,1);
            allBasePads[i]->SetTickx(1);
            allBasePads[i]->SetTicky(1);
        }

        gEPad->SetMargin(0.13,0.1,.1,0.1);


        //Plot gE

        gEPad->cd(); drawPlot(detType, PLOTVAR_GE,PLOT_COMP);
//.........这里部分代码省略.........

void advancedNoiseAnalysis( unsigned int runNumber, unsigned int loop = 1) {
  
  string inputFileName = "./histo/run00" + toString( runNumber ) + "-ped-histo.root";
  string outputFileName = "./histo/run00" + toString( runNumber ) + "-adv-noise.root";
  

  // before opening the input and the output files, try to see if they
  // are not opened yet and in case close them before continue   
  TList * listOfOpenedFile = (TList*) gROOT->GetListOfFiles();
  for ( int i = 0; i < listOfOpenedFile->GetSize() ; ++i ) {
    TFile * file = (TFile*) listOfOpenedFile->At( i ) ;
    TString fileName(file->GetName());
    TString inputFileName1( inputFileName.c_str() );
    TString outputFileName1( outputFileName.c_str() );

    if (  ( fileName.Contains( inputFileName1 ) ) ||
	  ( inputFileName1.Contains( fileName ) ) ||
	  ( fileName.Contains( outputFileName1 ) ) ||
	  ( outputFileName1.Contains( fileName ) ) ) {
      cout << "Closing " << fileName << " before reopen " << endl;
      file->Close();
    }
  }


  // close also all the previously opened canvas
  TList * listOfOpenedCanvas = (TList*) gROOT->GetListOfCanvases();
  for ( int i = 0 ; i < listOfOpenedCanvas->GetSize() ; ++i ) {
    TCanvas * canvas = (TCanvas*) listOfOpenedCanvas->At( i );
    TString canvasName2 = canvas->GetName();
    if ( canvasName2.Contains( "det" ) ) {
      canvas->Close();
    }
  }

	 
  // now safely open the file
  TFile * inputFile = TFile::Open( inputFileName.c_str() ) ;
  TFile * outputFile = TFile::Open( outputFileName.c_str(), "RECREATE") ;
  TList * outputHistoList = new TList;

  // look into the inputFile for a folder named
  string pedeProcessorFolderName = "PedestalAndNoiseCalculator";
  TDirectoryFile * pedeProcessorFolder = (TDirectoryFile*) inputFile->Get( pedeProcessorFolderName.c_str() );
  
  if ( pedeProcessorFolder == 0 ) { 
    cerr << "No pedestal processor folder found in file " << inputFileName << endl;
    return ;
  }

  // this folder should contain one folder for each loop.
  string loopFolderName = "loop-" + toString( loop );
  TDirectoryFile * loopFolder = (TDirectoryFile *) pedeProcessorFolder->Get( loopFolderName.c_str() );
  
  if ( loopFolder == 0 ) {
    cerr << "No " << loopFolderName << " found in file " << inputFileName << endl;
    return ;
  }

  // guess the number of sensors from the number of subfolder in the loopfolder
  size_t nDetector = loopFolder->GetListOfKeys()->GetSize();
  cout << "This file contains " << nDetector << " detectors" << endl;

  // prepare arrays to store the mean and the rms of the noise distribution
  if ( noiseMean == NULL ) {
    delete [] noiseMean;
    noiseMean = NULL;
  }
  if ( noiseRMS == NULL ) {
    delete [] noiseRMS;
    noiseRMS = NULL;
  }
  if ( channel == NULL ) {
    delete [] channel;
    channel = NULL;
  }

  noiseMean = new double[ nDetector * kNChan ];
  noiseRMS  = new double[ nDetector * kNChan ];
  channel   = new double[ kNChan ];

  string canvasName = "comparison";
  string canvasTitle = "Noise comparison";

  TCanvas * comparisonCanvas = new TCanvas( canvasName.c_str(), canvasTitle.c_str(), 1000, 500 );
  comparisonCanvas->Divide(1,2);
  
  TPad * topPad = (TPad*) comparisonCanvas->cd(1);
  topPad->Divide( nDetector );
  
  TPad * middlePad = (TPad *) comparisonCanvas->cd(2);
  middlePad->Divide( kNChan );


  // for each detector we have to get the noise map and to prepare 4
  // separe histos and maps
  for ( unsigned int iDetector = 0; iDetector < nDetector; iDetector++ ) {

    // get the noise map.
    string noiseMapName = "detector-" + toString( iDetector ) ;
//.........这里部分代码省略.........

/*#include "TTree.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TRint.h"
#include <time.h>
#include <stdio.h>

const short MaxNN = 6;
*/
void make_scatterplot(TString data_file, TString plot_name, TString header)
{
	//TApplication program = new TRInt(); 
	//pull in data
	TTree *t = new TTree(); 
	t->ReadFile(data_file);
	t->SetName("t");
	
	//TString name = "run2.png";
	TCanvas *BG = new TCanvas("c1", "Read Velocity on Local Disk for CMS3 Files", 1920, 1080);
	BG->cd();
	TPad *c = new TPad("MainPad", "My main pad", 0, 0, 1, 0.9);
	c->Draw();
	c->Divide(2,2);

	//In first slot have Time vs. Buffer for 1 Concurrent Read
	c->cd(1);
	int n1 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==1", "goff");
	TGraph *ghist1 = new TGraph(n1, t->GetV2(), t->GetV1());
	ghist1->SetName("ghist1");
	ghist1->SetMarkerStyle(3);
	ghist1->SetMarkerColor(1);
	ghist1->SetTitle("Single File Read");
	ghist1->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist1->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist1->Draw("ap");
	
	//In second slot have Time vs. Buffer for 3 Concurrent Reads
	c->cd(2);
	int n2 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==3", "goff");
	TGraph *ghist2 = new TGraph(n2, t->GetV2(), t->GetV1());
	Double_t DP2x[n2], DP2y[n2];
	ghist2->SetName("ghist2");
	ghist2->SetMarkerStyle(3);
	ghist2->SetMarkerColor(1);
	ghist2->SetTitle("3 Concurrent Reads");
	ghist2->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist2->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist2->Draw("ap");
	
	//In third slot have Time vs. Buffer for 3 Concurrent Reads
	c->cd(3);
	int n3 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==6", "goff");
	TGraph *ghist3 = new TGraph(n3, t->GetV2(), t->GetV1());
	Double_t DP3x[n3], DP3y[n3];
	ghist3->SetName("ghist3");
	ghist3->SetMarkerStyle(3);
	ghist3->SetMarkerColor(1);
	ghist3->SetTitle("6 Concurrent Reads");
	ghist3->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist3->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist3->Draw("ap");
	
	//In fourth slot have Time vs. Buffer for 10 Concurrent Reads
	c->cd(4);
	int n4 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==10", "goff");
	TGraph *ghist4 = new TGraph(n4, t->GetV2(), t->GetV1());
	ghist4->SetMarkerStyle(3);
	ghist4->SetName("ghist4");
	ghist4->SetMarkerColor(1);
	ghist4->SetTitle("10 Concurrent Reads");
	ghist4->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist4->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist4->Draw("ap");
	
	//Draw to screen
	//Initialize Canvas
	
	c->cd(0);
	c->Draw();
	
	BG->cd();
	TText *title = new TText(.5,.95, header);
	title->SetTextAlign(22);
	title->Draw();
	
	gDirectory->Add(ghist1);
	gDirectory->Add(ghist2);
	gDirectory->Add(ghist3);
	gDirectory->Add(ghist4);
	gDirectory->Add(t);
	gPad->SaveAs(plot_name);
}

//.........这里部分代码省略.........
      pt->Draw();
      cv->cd();
      pD = new TPad("dis","dis", 0, 0.0, 1, 0.93);
      pD->Draw();
      pD->cd();
    }
    
    pD->Clear();
    pD->cd();
    std::cout<<"Save : "<<dirName.c_str()<<"/"<<h1->GetName()<<std::endl;

    if (! isH2){
      h1->SetLineWidth(2);
      h1->SetLineColor(1);
      h1->SetMarkerColor(1);
      h2->SetLineColor(2);
      h2->SetMarkerColor(2);
      //      if (h1->GetNbinsX() > 25) h1 = h1->Rebin();
 //     if (h2->GetNbinsX() > 25) h2 = h2->Rebin();
      //      if (h1->GetNbinsX() > 50) h1 = h1->Rebin(5);
      //      if (h2->GetNbinsX() > 50) h2 = h2->Rebin(5);
      double max1 = h1->GetMaximum();
      double max2 = h2->GetMaximum();
      double min1 = h1->GetMinimum();
      double min2 = h2->GetMinimum();
      if (max2> max1) h1->SetMaximum(max2+0.15*fabs(max2));
      if (min2 < min1) h1->SetMinimum(min2-0.15*fabs(min2));
      //      pD->SetLogy();
      if ((logmod&1)) pD->SetLogx();
      if ((logmod&2)) pD->SetLogy();
      h1->Draw();
      h2->Draw("sames");   

      if (std::string(h1->GetName())==std::string("reconstruction_step_module_total")
	  || std::string(h1->GetName())==std::string("validation_step_module_total")){
	TPaveText ksPt(0,0, 0.35, 0.04, "NDC"); ksPt.SetBorderSize(0); ksPt.SetFillColor(0);
	ksPt.AddText(Form("P(KS)=%g, diffBins=%g, eblk %g ered %g",ksProb, bDiff, h1->GetEntries(), h2->GetEntries()));
	//	ksPt.AddText(h1->GetName());
	ksPt.Draw();
	cv->Print("diff.ps");

	int nX = h1->GetNbinsX();
	TAxis* h1Ax = h1->GetXaxis();
	int nRanges = nX/20 + 1;
	double h1Int = h1->Integral();
	float curBMargin = pD->GetBottomMargin(); 
	pD->SetBottomMargin(0.3);
	for (int iR = 0; iR < nRanges; iR++){
	  h1Ax->SetRange(iR*20+1, iR*20+20);
	  double bDiffL = 0;
	  double max1L = -1;
	  double max2L = -1;
	  double min1L = h1->GetMaximum();
	  double min2L = h2->GetMaximum();
	  for (int iBL = iR*20+1; iBL<= iR*20+20; ++iBL){
	    double h1L = h1->GetBinContent(iBL);
	    double h2L = h2->GetBinContent(iBL);
	    bDiffL += std::abs(h1L-h2L);
	    if (max1L < h1L) max1L = h1L;
	    if (max2L < h2L) max2L = h2L;
	    if (min1L > h1L) min1L = h1L;
	    if (min2L > h2L) min2L = h2L;
	  }
	  if (max2L> max1L)  h1->SetMaximum(max2L+0.15*std::abs(max2L));
	  else h1->SetMaximum(max1L+0.15*std::abs(max1L));
	  if (min2L < min1L) h1->SetMinimum(min2L-0.15*std::abs(min2L));
	  else h1->SetMinimum(min1L-0.15*std::abs(min1L));
	  h1->Draw();
	  h2->Draw("sames");
	  TPaveText ksPtL(0,0, 0.35, 0.04, "NDC"); ksPtL.SetBorderSize(0); ksPtL.SetFillColor(0);
	  ksPtL.AddText(Form("P(KS)=%g, diffBinsL=%g(%g), eblk %g ered %g",ksProb, bDiffL, bDiffL/h1Int, h1->GetEntries(), h2->GetEntries()));
	  ksPtL.Draw();
	  cv->Print("diff.ps");
	}
	pD->SetBottomMargin(curBMargin);

      }
    }
    if (isH2){
      pD->Divide(2);
      pD->cd(1);
      h1->Draw("colz");
      pD->cd(2);
      h2->Draw("colz");
    }
    TPaveText ksPt(0,0, 0.55, 0.06, "NDC"); ksPt.SetBorderSize(0); ksPt.SetFillColor(0);
    ksPt.AddText(Form("P(KS)=%g, diffBins=%g, eblk %g ered %g",ksProb, bDiff, h1->GetEntries(), h2->GetEntries()));
    ksPt.AddText(h1->GetName());
    ksPt.Draw();
    cv->Print("diff.ps");
    cv->Print("diff.pdf");


  }

  //  std::cout<<"Done in "<<dirName.c_str()<<std::endl;
  //  delete pH; delete pD;
  if (cv) delete cv;
  //  cv->Print("diff.ps]");
}

void PlotFinal(TGraphAsymmErrors* g_data_final, TGraphAsymmErrors* g_mg_final, TGraphAsymmErrors* g_ph_final, TGraphAsymmErrors* g_ratio_phistar, TGraphAsymmErrors* g_ratio_mg_phistar, TGraphAsymmErrors* g_ratio_ph_phistar, bool isPlot2 = 0, TGraphAsymmErrors* g_re_final = 0, TGraphAsymmErrors* g_ratio_re_phistar = 0) {

    vector<TGraphAsymmErrors*> g_data = SplitGraph(g_data_final);
    vector<TGraphAsymmErrors*> g_mg = SplitGraph(g_mg_final);
    vector<TGraphAsymmErrors*> g_ph = SplitGraph(g_ph_final);
    vector<TGraphAsymmErrors*> r_data = SplitGraph(g_ratio_phistar, 1);
    vector<TGraphAsymmErrors*> r_mg = SplitGraph(g_ratio_mg_phistar);
    vector<TGraphAsymmErrors*> r_ph = SplitGraph(g_ratio_ph_phistar);
    vector<TGraphAsymmErrors*> g_re = SplitGraph(g_re_final);
    vector<TGraphAsymmErrors*> r_re = SplitGraph(g_ratio_re_phistar);
    //New Graphs
    TGraphAsymmErrors* g_ANlo_final = ResbosFromRaj(2);
    TGraphAsymmErrors* r_ANlo_Ratio_final = CreateRatio(g_data_final, g_ANlo_final, 0);
    TGraphAsymmErrors* g_PowPyth8_final = ResbosFromRaj(1);
    TGraphAsymmErrors* r_PowPyth8_Ratio_final = CreateRatio(g_data_final, g_PowPyth8_final, 0);

    vector<TGraphAsymmErrors*> g_ANlo = SplitGraph(g_ANlo_final);
    vector<TGraphAsymmErrors*> r_ANlo = SplitGraph(r_ANlo_Ratio_final);
    vector<TGraphAsymmErrors*> g_Pyth8 = SplitGraph(g_PowPyth8_final);
    vector<TGraphAsymmErrors*> r_Pyth8 = SplitGraph(r_PowPyth8_Ratio_final);

    vector<TGraphAsymmErrors*> g_dummy = CreateDummy(g_data);
    OneDYPlot(g_data, g_mg, g_ph, g_re, g_ANlo, g_Pyth8);
    for (uint i = 0; i < ny; i++) {
        //  for (uint i=0; i<5; i++){ 
        std::ostringstream strs;
        strs << i;
        //TO HERE
        std::string Canvasname = "EventEff_Bin" + strs.str();

        TCanvas* FinalPhiTot = new TCanvas(Canvasname.c_str(), Canvasname.c_str(), 800, 900);
        FinalPhiTot->Divide(1, 2);
        FinalPhiTot->cd(1);
        gPad->SetPad("p1", "p1", 0, 2.5 / 9.0, 1, 1, kWhite, 0, 0);
        gPad->SetBottomMargin(0.01);
        gPad->SetTopMargin(0.06);
        gPad->SetLeftMargin(0.15);
        gPad->SetRightMargin(0.06);
        gPad->SetLogx(1);
        gPad->SetLogy(1);
        g_dummy[i]->GetXaxis()->SetRangeUser(0.001, 10.0);
        if (doNorm) g_dummy[i]->GetYaxis()->SetRangeUser(0.00000005, 1000.0);
        else g_dummy[i]->GetYaxis()->SetRangeUser(0.00005, 1000000.0);
        g_dummy[i]->GetXaxis()->CenterTitle();
        g_dummy[i]->GetYaxis()->CenterTitle();
        g_dummy[i]->Draw("A2");
        g_mg[i]->SetMarkerColor(kBlue - 7);
        g_mg[i]->SetLineColor(kBlue - 7);
        g_mg[i]->SetMarkerSize(1);
        g_mg[i]->SetLineWidth(2);
        g_mg[i]->SetMarkerStyle(21);
        g_mg[i]->Draw("PEsame");
        g_ph[i]->SetMarkerColor(kRed);
        g_ph[i]->SetLineColor(kRed);
        g_ph[i]->SetMarkerSize(1);
        g_ph[i]->SetLineWidth(2);
        g_ph[i]->SetMarkerStyle(22);
        g_ph[i]->Draw("PEsame");
        if (!isPlot2 && elec == 1 && g_re_final) {
            g_re[i]->SetMarkerColor(kGreen + 1);
            g_re[i]->SetLineColor(kGreen + 1);
            g_re[i]->SetMarkerSize(1);
            g_re[i]->SetLineWidth(2);
            g_re[i]->SetMarkerStyle(23);
            g_re[i]->Draw("PEsame");
        }
        g_data[i]->SetFillColor(kYellow);
        g_data[i]->SetMarkerSize(1);
        g_data[i]->SetLineWidth(2);
        g_data[i]->SetMarkerStyle(20);
        g_data[i]->Draw("PEsame");
        g_data[i]->SetFillColor(kYellow);

        TLegend* leg;
        if (isPlot2) leg = new TLegend(0.15, 0.06, 0.80, 0.27);
        else leg = new TLegend(0.15, 0.06, 0.80, 0.31); //TLegend(0.45,0.73,0.94,0.93);//.19 0.06
        leg->SetFillStyle(0);
        leg->SetBorderSize(0);
        leg->SetLineWidth(1);
        leg->SetNColumns(1);
        leg->SetTextFont(42);

        if (!isPlot2) {
            leg->AddEntry(g_data[i], "2012 data", "PEF");
            if (Type == "elec") {
                leg->AddEntry(g_mg[i], "#gamma*/Z #rightarrow ee (MadGraph+Pythia6 Z2*)", "P");
                leg->AddEntry(g_ph[i], "#gamma*/Z #rightarrow ee (Powheg+Pythia6 Z2*)", "P");
                //ToDo AMCAt decisions
                //if (elec == 1 && g_re_final) leg->AddEntry(g_re[i], "#gamma*/Z #rightarrow ee (Resbos)", "P");
                if (elec == 1 && g_re_final) leg->AddEntry(g_re[i], "#gamma*/Z #rightarrow ee (POWHEG+Pythia8)", "P");
            }
            if (Type == "muon") {
                leg->AddEntry(g_mg[i], "#gamma*/Z #rightarrow #mu#mu (MadGraph+Pythia6 Z2*)", "P");
                leg->AddEntry(g_ph[i], "#gamma*/Z #rightarrow #mu#mu (Powheg+Pythia6 Z2*)", "P");
                if (elec == 1 && g_re_final) leg->AddEntry(g_re[i], "#gamma*/Z #rightarrow #mu#mu (Resbos)", "P");
                //if (elec == 1 && g_re_final) leg->AddEntry(g_re[i], "#gamma*/Z #rightarrow #mu#mu (POWHEG+Pythia8 CUETP8M1)", "P");
            }
            if (Type == "combined") {
                leg->AddEntry(g_mg[i], "#gamma*/Z #rightarrow ll (MadGraph+Pythia6 Z2*)", "P");
                leg->AddEntry(g_ph[i], "#gamma*/Z #rightarrow ll (Powheg+Pythia6 Z2*)", "P");
//.........这里部分代码省略.........

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

    // saving histos
    TFile f("dist.root", "RECREATE");
    for (int i = 0; i < MNP; i++) {
      perfec[i]->Write();
      doNoth[i]->Write();
      fixBad[i]->Write();
      fixAll[i]->Write();
    }
    f.Close();
  }

  // histos loaded, now plotting
  for (int i = 0; i < MNP; i++) {
    doNoth[i]->SetLineColor(kRed);
    fixBad[i]->SetLineColor(kBlue);
    fixAll[i]->SetLineColor(kGreen);

    perfec[i]->GetXaxis()->SetLabelSize(0.08);
    perfec[i]->GetYaxis()->SetLabelSize(0.08);
    perfec[i]->GetXaxis()->SetLabelOffset(-0.02);
    perfec[i]->GetYaxis()->SetLabelOffset(0.02);
  }

  double xmom[10];
  double ymomA[10];
  double ymomB[10];
  double ymomC[10];

  TCanvas* res = new TCanvas("res", "Photon Yields", 1500, 1200);

  TPad* pres = new TPad("pres", "pres", 0.01, 0.01, 0.98, 0.9);
  pres->Divide(5, 2);
  pres->Draw();

  TLatex lab;
  lab.SetTextFont(102);
  lab.SetTextColor(kBlue + 2);
  lab.SetTextSize(0.06);
  lab.SetNDC();

  for (int i = 2; i < 12; i++) {
    pres->cd(i - 1);

    perfec[i]->Draw();
    fixAll[i]->Draw("same");
    doNoth[i]->Draw("same");
    fixBad[i]->Draw("same");

    double perfectCounts = perfec[i]->Integral(40, 250);
    double donothCounts = doNoth[i]->Integral(40, 250);
    double fixBadCounts = fixBad[i]->Integral(40, 250);
    double fixAllCounts = fixAll[i]->Integral(40, 250);

    double dmom = (max_m - min_m) / MNP;
    double mom = min_m + i * dmom;

    xmom[i - 2] = mom;
    ymomA[i - 2] = 100 * donothCounts / perfectCounts;
    ymomB[i - 2] = 100 * fixBadCounts / perfectCounts;
    ymomC[i - 2] = 100 * fixAllCounts / perfectCounts;

    cout << " momentum: " << mom << " perfect: " << perfectCounts
         << "   nothing: " << donothCounts << "    fix bad: " << fixBadCounts
         << "   fix all " << fixAllCounts << endl;

void plotFinalEnergy()  {

    
    gStyle->SetStatH(0.3);
    gStyle->SetStatY(0.88);
    gStyle->SetStatW(0.3);
    char canTitle[180];
    sprintf(canTitle,"Energy at Terminal Interaction (10000 tries)"); 
    TCanvas *canFinalEnergy = new TCanvas("canFinalEnergy","Final Energy",800,800);
    TPaveLabel *pl = new TPaveLabel(0.1,0.96,0.9,0.99,canTitle,"br NDC");
    pl->SetBorderSize(0);
    pl->SetFillColor(0);
    pl->SetFillStyle(0);
    pl->Draw();
    TPad *subCanFinalEnergy = new TPad("subCanFinalEnergy","",0,0,1,0.95);
    subCanFinalEnergy->Draw();
    subCanFinalEnergy->cd();
    subCanFinalEnergy->Divide(2,4);
    subCanFinalEnergy->Update();
    


//    gStyle->SetOptStat(0);
    char fileName[80];   
    char histTitle[80];

    char theEnergies[4][5]={"1e9","1e10","1e11","1e12"};

    int theColours[2][3]={{50,42,46},{40,30,38}};
    
    for(int isATau=0;isATau<=1;isATau++) {
	for(int i=0;i<4;i++) {		 
	    sprintf(fileName,
		    "newest%sFile%sIce.root",getParticleNameCaps(isATau),theEnergies[i]);
	    sprintf(histTitle,"%s -- %s GeV",getParticleNameCaps(isATau),theEnergies[i]);
	    TFile *fp = new TFile(fileName);
	    TH1F *histEnergy = 
		new TH1F("histEnergy",histTitle,100,7.5,12.5);
	    TH1F *histEnergy2 = 
		new TH1F("histEnergy2","Last Energy",100,7.5,12.5);
	    TH1F *histEnergy3 = 
		new TH1F("histEnergy3","Last Energy",100,7.5,12.5);

	    TTree *theTree = (TTree*) fp->Get("theTree");
	    cout << theTree->GetEntries() << endl;
	    subCanFinalEnergy->cd((2*i)+isATau+1);
	    gPad->SetTopMargin(0.12);
//	gPad->SetBottomMargin(0.2);
	    histEnergy->SetLineWidth(3);
	    histEnergy->SetLineColor(theColours[isATau][0]);
	    histEnergy2->SetLineWidth(3);
	    histEnergy2->SetLineColor(theColours[isATau][1]);
	    histEnergy3->SetLineWidth(3);
	    histEnergy3->SetLineColor(theColours[isATau][2]);
	    theTree->Draw("log10(stepIntEnergy)>>histEnergy","stepIntType>=4");
//	histEnergy->GetXaxis()->SetTitle("IntType (#mus)");
// 	histEnergy->GetXaxis()->SetBinLabel(2,"Bremsstrahlung");
// 	histEnergy->GetXaxis()->SetBinLabel(3,"Pair Production");
// 	histEnergy->GetXaxis()->SetBinLabel(4,"Photonuclear");
// 	histEnergy->GetXaxis()->SetLabelSize(0.09);
// 	histEnergy->GetXaxis()->SetLabelOffset(0.02);
 	    theTree->Draw("log10(stepIntEnergy)>>histEnergy2","stepIntType==4");
 	    theTree->Draw("log10(stepIntEnergy)>>histEnergy3","stepIntType==5");

	    histEnergy->DrawCopy();
	    if(histEnergy2->GetEntries())
		histEnergy2->DrawCopy("same");
	    if(histEnergy3->GetEntries())
		histEnergy3->DrawCopy("same");
	    
	    if(histEnergy->GetEntries())
		gPad->SetLogy();

	    if(i==0) {
		TLegend *leg = new TLegend(0.7,0.2,0.9,0.6);
		leg->SetFillColor(0);
		leg->SetFillStyle(0);
		leg->SetBorderSize(0);
		leg->AddEntry(histEnergy2,"Decay","l");
		leg->AddEntry(histEnergy3,"Weak (CC)","l");
		leg->AddEntry(histEnergy,"Either","l");
		leg->Draw("same");
	    }
	}
    } 
    for(int i=0;i<8;i++) {
//	cout << "Doing title "<< i+1 << endl;
	subCanFinalEnergy->cd(i+1);
	sortOutTitle(0.07);
    }
//    gStyle->SetOptStat(1110);

}

void quarks () {
   TCanvas *c1 = new TCanvas("c1", "c1",10,10,630,760);
   c1->SetFillColor(kBlack);
   Int_t quarkColor  = 50;
   Int_t leptonColor = 16;
   Int_t forceColor  = 38;
   Int_t titleColor  = kYellow;
   Int_t border = 8;

   TLatex *texf = new TLatex(0.90,0.455,"Force Carriers");
   texf->SetTextColor(forceColor);
   texf->SetTextAlign(22); texf->SetTextSize(0.07);
   texf->SetTextAngle(90);
   texf->Draw();

   TLatex *texl = new TLatex(0.11,0.288,"Leptons");
   texl->SetTextColor(leptonColor);
   texl->SetTextAlign(22); texl->SetTextSize(0.07);
   texl->SetTextAngle(90);
   texl->Draw();

   TLatex *texq = new TLatex(0.11,0.624,"Quarks");
   texq->SetTextColor(quarkColor);
   texq->SetTextAlign(22); texq->SetTextSize(0.07);
   texq->SetTextAngle(90);
   texq->Draw();

   TLatex tex(0.5,0.5,"u");
   tex.SetTextColor(titleColor); tex.SetTextFont(32);
   tex.SetTextAlign(22);
   tex.SetTextSize(0.14);
   tex.DrawLatex(0.5,0.93,"Elementary");
   tex.SetTextSize(0.12);
   tex.DrawLatex(0.5,0.84,"Particles");
   tex.SetTextSize(0.05);
   tex.DrawLatex(0.5,0.067,"Three Generations of Matter");

   tex.SetTextColor(kBlack); tex.SetTextSize(0.8);

// ------------>Create main pad and its subdivisions
   TPad *pad = new TPad("pad", "pad",0.15,0.11,0.85,0.79);
   pad->Draw();
   pad->cd();
   pad->Divide(4,4,0.0003,0.0003);

   pad->cd(1); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"u");

   pad->cd(2); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"c");

   pad->cd(3); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"t");

   pad->cd(4); gPad->SetFillColor(forceColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.55,"#gamma");

   pad->cd(5); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"d");

   pad->cd(6); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"s");

   pad->cd(7); gPad->SetFillColor(quarkColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"b");

   pad->cd(8); gPad->SetFillColor(forceColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.55,"g");

   pad->cd(9); gPad->SetFillColor(leptonColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"#nu_{e}");

   pad->cd(10); gPad->SetFillColor(leptonColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"#nu_{#mu}");

   pad->cd(11); gPad->SetFillColor(leptonColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"#nu_{#tau}");

   pad->cd(12); gPad->SetFillColor(forceColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"Z");

   pad->cd(13); gPad->SetFillColor(leptonColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.5,"e");

   pad->cd(14); gPad->SetFillColor(leptonColor);
   gPad->SetBorderSize(border);
   tex.DrawLatex(.5,.56,"#mu");
//.........这里部分代码省略.........