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

TCanvas* example_plot( int iPeriod, int iPos )
{ 
  //  if( iPos==0 ) relPosX = 0.12;

  int W = 800;
  int H = 600;

  // 
  // Simple example of macro: plot with CMS name and lumi text
  //  (this script does not pretend to work in all configurations)
  // iPeriod = 1*(0/1 7 TeV) + 2*(0/1 8 TeV)  + 4*(0/1 13 TeV) 
  // For instance: 
  //               iPeriod = 3 means: 7 TeV + 8 TeV
  //               iPeriod = 7 means: 7 TeV + 8 TeV + 13 TeV 
  // Initiated by: Gautier Hamel de Monchenault (Saclay)
  //
  int H_ref = 600; 
  int W_ref = 800; 

  // references for T, B, L, R
  float T = 0.08*H_ref;
  float B = 0.12*H_ref; 
  float L = 0.12*W_ref;
  float R = 0.04*W_ref;

  TString canvName = "FigExample_";
  canvName += W;
  canvName += "-";
  canvName += H;
  canvName += "_";  
  canvName += iPeriod;
  if( writeExtraText ) canvName += "-prelim";
  if( iPos%10==0 ) canvName += "-out";
  else if( iPos%10==1 ) canvName += "-left";
  else if( iPos%10==2 )  canvName += "-center";
  else if( iPos%10==3 )  canvName += "-right";

  TCanvas* canv = new TCanvas(canvName,canvName,50,50,W,H);
  canv->SetFillColor(0);
  canv->SetBorderMode(0);
  canv->SetFrameFillStyle(0);
  canv->SetFrameBorderMode(0);
  canv->SetLeftMargin( L/W );
  canv->SetRightMargin( R/W );
  canv->SetTopMargin( T/H );
  canv->SetBottomMargin( B/H );
  canv->SetTickx(0);
  canv->SetTicky(0);

  TH1* h = new TH1F("h","h",40,70,110);
  h->GetXaxis()->SetNdivisions(6,5,0);
  h->GetXaxis()->SetTitle("m_{e^{+}e^{-}} (GeV)");  
  h->GetYaxis()->SetNdivisions(6,5,0);
  h->GetYaxis()->SetTitleOffset(1);
  h->GetYaxis()->SetTitle("Events / 0.5 GeV");  

  h->SetMaximum( 260 );
  if( iPos==1 ) h->SetMaximum( 300 );
  h->Draw();

  int histLineColor = kOrange+7;
  int histFillColor = kOrange-2;
  float markerSize  = 1.0;

  {
    TLatex latex;
				
    int n_ = 2;

    float x1_l = 0.92;
    float y1_l = 0.60;

    float dx_l = 0.30;
    float dy_l = 0.18;
    float x0_l = x1_l-dx_l;
    float y0_l = y1_l-dy_l;

    TPad* legend = new TPad("legend_0","legend_0",x0_l,y0_l,x1_l, y1_l );
    //    legend->SetFillColor( kGray );
    legend->Draw();
    legend->cd();
		
    float ar_l = dy_l/dx_l;
		
    float x_l[1];
    float ex_l[1];
    float y_l[1];
    float ey_l[1];
		
    //    float gap_ = 0.09/ar_l;
    float gap_ = 1./(n_+1);
		
    float bwx_ = 0.12;
    float bwy_ = gap_/1.5;
		
    x_l[0] = 1.2*bwx_;
    //    y_l[0] = 1-(1-0.10)/ar_l;
    y_l[0] = 1-gap_;
    ex_l[0] = 0;
    ey_l[0] = 0.04/ar_l;
//.........这里部分代码省略.........

void PlotShapeSystematics() {

  TFile *file = 0;
  TH1F *DefaultShape = 0;
  TH1F *UpShape = 0;
  TH1F *DownShape = 0;
  TCanvas *cv = 0;
  TLegend *legend = 0;

  //*********************************************************
  //0 Jet Bin - OF
  //*********************************************************

  file = new TFile("/data/smurf/sixie/data/Thesis/cards/130/hwwof_0j.input.root","READ");

  DefaultShape = (TH1F*)file->Get("histo_qqWW");
  UpShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwUp");
  DownShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwDown");

  assert(DefaultShape);
  assert(UpShape);
  assert(DownShape);

  NormalizeHist(DefaultShape);
  NormalizeHist(UpShape);
  NormalizeHist(DownShape);


  cv = new TCanvas("cv","cv",800,600);

  legend = new TLegend(0.2, 0.6, 0.5, 0.9);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);
  legend->SetFillColor(kWhite);
  legend->AddEntry(DefaultShape, "Default Shape", "L");
  legend->AddEntry(UpShape, "Bounding Shape (Up)", "L");
  legend->AddEntry(DownShape, "Bounding Shape (Down)", "L");

  DefaultShape->SetLineColor(kBlack);
  UpShape->SetLineColor(kBlue);
  DownShape->SetLineColor(kRed);
  DefaultShape->SetLineWidth(2);
  UpShape->SetLineWidth(2);
  DownShape->SetLineWidth(2);
  DefaultShape->SetTitle("");
  DefaultShape->GetXaxis()->SetTitle("MVA discriminator");
  DefaultShape->GetYaxis()->SetTitle("Fraction of Events");
  DefaultShape->GetYaxis()->SetTitleOffset(1.4);
  DefaultShape->GetXaxis()->SetRangeUser(-1.0,1.0);
  DefaultShape->GetYaxis()->SetRangeUser(0.0,0.25);


  DefaultShape->Draw("hist");
  UpShape->Draw("same,hist");
  DownShape->Draw("same,hist");
  legend->Draw();

  cv->SaveAs("WWBkgShapeVariation_MCAtNLOScaleVariation_OF0Jet.png");
  cv->SaveAs("WWBkgShapeVariation_MCAtNLOScaleVariation_OF0Jet.eps");



  //*********************************************************
  //0 Jet Bin - SF
  //*********************************************************

  file = new TFile("/data/smurf/sixie/data/Thesis/cards/130/hwwsf_0j.input.root","READ");

  DefaultShape = (TH1F*)file->Get("histo_qqWW");
  UpShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwUp");
  DownShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwDown");

  assert(DefaultShape);
  assert(UpShape);
  assert(DownShape);

  NormalizeHist(DefaultShape);
  NormalizeHist(UpShape);
  NormalizeHist(DownShape);


  cv = new TCanvas("cv","cv",800,600);

  legend = new TLegend(0.2, 0.6, 0.5, 0.9);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);
  legend->SetFillColor(kWhite);
  legend->AddEntry(DefaultShape, "Default Shape", "L");
  legend->AddEntry(UpShape, "Bounding Shape (Up)", "L");
  legend->AddEntry(DownShape, "Bounding Shape (Down)", "L");

  DefaultShape->SetLineColor(kBlack);
  UpShape->SetLineColor(kBlue);
  DownShape->SetLineColor(kRed);
  DefaultShape->SetLineWidth(2);
  UpShape->SetLineWidth(2);
  DownShape->SetLineWidth(2);
  DefaultShape->SetTitle("");
  DefaultShape->GetXaxis()->SetTitle("MVA discriminator");
  DefaultShape->GetYaxis()->SetTitle("Fraction of Events");
//.........这里部分代码省略.........

void dileptonMassFit(const char* pInFileName="PanchoSkim4JanAll.root",
                     // "PromtRecoV2V3V3H_DiMuonPlot_TightSTACutsAll15Dec.root",
                     // "Z0_DataMixPt50_PatDiMuonPlots_NewCutAll14Dec.root",
                     const char* pHistNameOpCh="diMuonsGlobalInvMassVsPt",//diMuonsGlobalInvMassVsPtW",
                     const char* pHistNameSameCh="diMuonsGlobalSameChargeInvMassVsPt",
                     const char* pSpectra="pt",  // pt, y, centr
                     bool doMc=false,
                     int nFitFunction = 3,
                     int getYield = 1)
{
    gROOT->Macro("setStyle.C+");

    //gROOT->Macro("/Users/eusmartass/Software/utilities/setStyle.C+");
    char szBuf[256];

    ////////  definitions of Switches   ///////////
    //  nFitFunction  = 1  RBW + Pol2
    //  nFitFunction  = 2  Gaus + Pol2
    //  nFitFunction  = 3  RBWGaus + Pol2

    //  getYield = 1  Bin counting
    //  getYield = 2  Integral
    ////////////////////////////////////////////////////////////

    // make some choices
    float MassZ0         = 91.1876;
    float WidthZ0        = 2.4952;
    float massFit_low    = 60;
    float massFit_high   = 120;   // Fit ranges
    float massDraw_low   = 30.0;  // 0.
    float massDraw_high  = 130.0; // 200/
    int nrebin           = 80;
    bool isLog           = 0;
    bool isFit           = 1; // draw ranges

    float massCount_low  = 60.0; //78.0
    float massCount_high = 120.0; //102.0

    //___________________________________________________________________________________
    // ------- Open input file
    sprintf(szBuf,"%s",pInFileName);
    TString inFileName(szBuf);
    TFile *pfInFile = new TFile(inFileName);

    // ------- get histograms:
    sprintf(szBuf,"%s",pHistNameOpCh);
    TH2D *phDimuMass_1  = (TH2D*)pfInFile->Get(szBuf)->Clone("phDimuMass_1");

    sprintf(szBuf,"%s",pHistNameSameCh);
    TH2D *phDimuMass_1S = (TH2D*)pfInFile->Get(szBuf)->Clone("phDimuMass_1S");

    phDimuMass_1->SetDirectory(0);
    phDimuMass_1S->SetDirectory(0);

    // Open pp data file

    TFile *ppFile = new TFile("Zmumu_40-200_35pb.root");
    TH1F *Zmumu  = (TH1F*)ppFile->Get("hdata");


    //___________________________________________________________________________________
    // bins definition:
    const char* Xname[] = {" ", "p_{T}^{Dimuon} (GeV/c)", "rapidity", "centrality"};
    bool doPt   = false;
    bool doY    =  false;
    bool doCent = false;

    int GenRange, nBins;
    double binEdge[10];
    char* label;
    sprintf(szBuf,"%s",pSpectra);
    TString wichSpectra(szBuf);

    if ( wichSpectra.CompareTo("pt") == 0) {
        doPt               = true;
        label              = (char*)Xname[1];
        GenRange           = 20;
        nBins              = 1;
        binEdge[0] = 0.0;
        binEdge[1]= 100.0;
        //    double binEdge[10] = {0.0, 10., 20., 100.0};
        if(doMc) {
            nBins              = 1;
            binEdge[0] = 0.0;
            binEdge[1]= 50.0;
            //	  nBins          = 7;
            // binEdge[0] =  0.0;  binEdge[1] =  2.0;  binEdge[2] =  4.0;  binEdge[3] = 8.0;
            //binEdge[4] = 12.0;  binEdge[5] = 16.0;  binEdge[6] = 22.0;  binEdge[7] = 50.0;
        }
    } else {
        if ( wichSpectra.CompareTo("y") == 0) {
            doY              = true;
            label            = (char*)Xname[2];
            nBins            = 3;
            GenRange         = 4.8;
            binEdge[0] = -2.4;
            binEdge[1] = -0.8;
            binEdge[2] =  0.8;
            binEdge[3] =  2.4;
        } else {
//.........这里部分代码省略.........

//------------------------------------------------------------------------------
//Subtraction 
//------------------------------------------------------------------------------
void Subtraction(TString  hname,
		   TString  xtitle,
		   Int_t    ngroup       = -1,
		   Int_t    precision    = 1,
		   TString  units        = "NULL",
		   Double_t xmin         = -999,
		   Double_t xmax         =  999,
		   Bool_t   moveOverflow = true)
{
  TCanvas* canvas = new TCanvas(hname, hname, 800, 800);

  TPad* pad1 = new TPad("pad1", "pad1", 0, 0.0, 1, 1.0);

  pad1->SetTopMargin   (0.08);
  //pad1->SetBottomMargin(0.02);
  pad1->Draw();
      

  //----------------------------------------------------------------------------
  // pad1
  //----------------------------------------------------------------------------
  pad1->cd();

  pad1->SetLogy(_setLogy);

  TH1F* hist[nProcesses];

  for (UInt_t ip=0; ip<nProcesses; ip++) {

    hist[ip] = (TH1F*)input[ip]->Get(hname);
    hist[ip]->SetName(hname + process[ip]);

    if (moveOverflow) MoveOverflowBins  (hist[ip], xmin, xmax);
    else              ZeroOutOfRangeBins(hist[ip], xmin, xmax);

	if (ngroup > 0) hist[ip]->Rebin(ngroup);

	if (_dataDriven && ip == iWW)    hist[ip]->Scale(WWScale[_njet]);
	if (_dataDriven && ip == iDY)    hist[ip]->Scale(ZjScale[_njet]);
	if (_dataDriven && ip == iDYtau) hist[ip]->Scale(ZjScale[_njet]);

  }

  // Data subtraction for Top background estimation
  //----------------------------------------------------------------------------
  TH1F* subData = (TH1F*)hist[iData]->Clone("subData");
  for (UInt_t ip=0; ip<nProcesses; ip++) {
	  if (ip == itt) continue;
	  if (ip == itW) continue;
	  if (ip == iData ) continue;
	  subData->Add(hist[ip],-1);
  }
  subData->SetLineColor(kRed+1);
  Double_t subData_Yield = subData->Integral();
  //subData->SetLineColor();

  // Top background
  //----------------------------------------------------------------------------
  TH1F* Top = (TH1F*)hist[itt]->Clone("Top");
  Top->Add(hist[itW]);
  Top->SetLineColor(kBlue+1);
  Double_t Top_Yield = Top->Integral();

  // Axis labels
  //----------------------------------------------------------------------------
  TAxis* xaxis = subData->GetXaxis();
  TAxis* yaxis = subData->GetYaxis();

  TString ytitle = Form("entries / %s.%df", "%", precision);

  xaxis->SetTitle(xtitle);
  yaxis->SetTitle(Form(ytitle.Data(), subData->GetBinWidth(0)));
  yaxis->SetTitleOffset(1.6);

  if (!units.Contains("NULL")) {
    
    xaxis->SetTitle(Form("%s [%s]", xaxis->GetTitle(), units.Data()));
    yaxis->SetTitle(Form("%s %s",   yaxis->GetTitle(), units.Data()));
  }


  // Draw
  //----------------------------------------------------------------------------
  xaxis->SetRangeUser(xmin, xmax);

  subData->Draw("hist");
  Top->Draw("hist same");

  // Adjust scale
  //----------------------------------------------------------------------------
  subData->SetMinimum(0.0);
  
  Float_t theMax   = GetMaximumIncludingErrors(subData, xmin, xmax);
  Float_t theMaxMC = GetMaximumIncludingErrors(Top,       xmin, xmax);

  if (theMaxMC > theMax) theMax = theMaxMC;

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

void MakePlotDATAMC(std::string var, std::string xlabel){

  TCanvas* c1 = new TCanvas("c1", "c1", 1);
  c1->cd();
  c1->Divide(1,2);
  c1->cd(1);
  TFile* fin_data = TFile::Open("/afs/cern.ch/work/s/soffi/CMSSW_5_3_26-MonoJet/src/GJet/histos_GJets_DATA.root");
  TFile* fin_bkg = TFile::Open("/afs/cern.ch/work/s/soffi/CMSSW_5_3_26-MonoJet/src/GJet/histos_GJets_MC.root");
 
  TH1F* h_data  =  (TH1F*) fin_data->Get(("h_"+var).c_str());
  TH1F* h_bkg  =  (TH1F*) fin_bkg->Get(("h_"+var).c_str());
 
 
  h_data->Sumw2();
  h_bkg->Sumw2();
  double data_norm = h_data->Integral();
  double bkg_norm = h_bkg->Integral()*19.7*1.33;
  std::cout<<data_norm<< " "<<bkg_norm<< std::endl;
  //h_bkg->Scale( h_data->Integral()/h_bkg->Integral());
  h_bkg->Scale(bkg_norm/h_bkg->Integral());
 
  if(var ==  "phmet" || var ==  "phmetW"){
    hdensity(h_data);
    hdensity(h_bkg);
  std::cout<<"\\begin{table}[bthp]"<<std::endl;
  std::cout<<"\\begin{tabular}{||cc||}"<<std::endl;
  std::cout<<"\\hline"<<std::endl;
  // std::cout<<"MET range  \& $\\gamma$+jets events \\\\"<<td::endl;
  std::cout<<"MET range  \& $Z\\rightarrow \\nu\\nu$ events \\\\"<<td::endl;
  std::cout<<"\\hline"<<std::endl;
  std::cout<<"250--300 \&"<<h_data->Integral(h_data->FindBin(250),h_data->FindBin(300))<<" \\\\"<<std::endl;
  std::cout<<"300--350 \&"<<h_data->Integral(h_data->FindBin(300)+1,h_data->FindBin(350))<<" \\\\"<<std::endl;
  std::cout<<"350--400 \&"<<h_data->Integral(h_data->FindBin(350)+1,h_data->FindBin(400))<<" \\\\"<<std::endl;
  std::cout<<"400--500 \&"<<h_data->Integral(h_data->FindBin(400)+1,h_data->FindBin(500))<<" \\\\"<<std::endl;
    std::cout<<"500--1000 \&"<<h_data->Integral(h_data->FindBin(500)+1,h_data->FindBin(1000))<<" \\\\"<<std::endl;
  std::cout<<"\\hline"<<std::endl;
  std::cout<<"\\end{tabular}\\end{table}"<<std::endl;
  }
  h_bkg->SetLineColor(kAzure+7);
  h_bkg->SetFillColor(kAzure+6);
  h_bkg->SetLineWidth(2);
 
  h_data->Draw("PE");
  h_bkg->Draw("histsame");
  h_bkg->GetXaxis()->SetTitle(xlabel.c_str());
  h_data->SetMarkerSize(0.7);
  h_data->Draw("PEsame");

  TLegend* legmc;
  legmc = new TLegend(0.6, 0.55, 0.82, 0.89, "", "brNDC");
 
  legmc->SetTextFont(42);
  legmc->SetBorderSize(0);
  legmc->SetFillStyle(0);
 
  legmc->AddEntry(h_data, "Data", "L");
  legmc->AddEntry(h_bkg, "#gamma + jet", "L");
  legmc->Draw("same");

  CMS_lumi( (TPad*)c1->cd(1),true,0);

  c1->cd(2);

  TH1F* hs = (TH1F*) h_data->Clone();
  hs->Divide(h_bkg);
  hs->Draw("pe");
  hs->GetXaxis()->SetTitle(xlabel.c_str()); 
  hs->GetYaxis()->SetTitle("Data/MC");
  hs->GetYaxis()->SetRangeUser(0., 2.);
 
  TLine line(hs->GetBinCenter(1),1,hs->GetBinCenter(hs->GetNbinsX()+1),1);
  line.SetLineColor(kRed);
  line.Draw("same");
  
  int iPos=0 ;
  CMS_lumi( (TPad*)c1->cd(2),true,iPos );
 
  c1->SaveAs(("~/www/MonoJet/MC-DATA_"+var+"_cat.png").c_str());
  c1->SaveAs(("~/www/MonoJet/MC-DATA_"+var+"_cat.pdf").c_str());
  c1->cd(1)->SetLogy();
  c1->SaveAs(("~/www/MonoJet/MC-DATA_"+var+"_LOG_cat.png").c_str());
  c1->SaveAs(("~/www/MonoJet/MC-DATA_"+var+"_LOG_cat.pdf").c_str());

}

//------------------------------------------------------------------------
void DrawLEDminCFD()
{
  Int_t npeaks = 10;
  Int_t sigma=10.;
  Bool_t down=false;
  Int_t index[20];
  
  TCanvas *c1 = new TCanvas("c1", " LED-CFD C side",0,48,1280,951);
  c1->Divide(4,3);
  
  Char_t buf1[20];
  for (Int_t i=0; i<12; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"LEDminCFD%i",i+1);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      cfd->Draw();
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.2);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0) {
	Float_t *xpeak = s->GetPositionX();
	TMath::Sort(nfound, xpeak, index,down);
	Float_t xp = xpeak[index[0]];
        Int_t xbin = cfd->GetXaxis()->FindBin(xp);
        Float_t yp = cfd->GetBinContent(xbin);
        cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[index[0]]<<"\typeak = "<<yp<<endl;
        Float_t hmin=xp-3*sigma;
        Float_t hmax =xp+3*sigma;
        cfd->GetXaxis()->SetRange(hmin,hmax);
	TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
        cfd->Fit("g1","RQ");

	
      }
      
    }
 TCanvas *c2 = new TCanvas("c2", "LED-CFD A side",0,48,1280,951);
  c2->Divide(4,3);
  
  Char_t buf1[20];
  for (Int_t i=12; i<24; i++)
    {
      c2->cd(i+1-12);
      sprintf(buf1,"LEDminCFD%i",i+1);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      cfd->Draw();
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.2);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0) {
	Float_t *xpeak = s->GetPositionX();
	TMath::Sort(nfound, xpeak, index,down);
	Float_t xp = xpeak[index[0]];
        Int_t xbin = cfd->GetXaxis()->FindBin(xp);
        Float_t yp = cfd->GetBinContent(xbin);
        cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[index[0]]<<"\typeak = "<<yp<<endl;
        Float_t hmin=xp-3*sigma;
        Float_t hmax =xp+3*sigma;
        cfd->GetXaxis()->SetRange(hmin,hmax);
	TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
        cfd->Fit("g1","RQ");

	
      }
      
    }
  /*
  TCanvas *c1 = new TCanvas("c1", "c1",0,48,1280,951);
  c1->Divide(2,2);
  Char_t buf1[10];
  for (Int_t i=0; i<4; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"LED-CFD%i",i+1);
      TH1F *cfd = (TH1F*) file->Get(buf1);
      //  cout<<buf1<<" "<<cfd<<endl;
      //     cfd->Draw();
      //   cfd->GetXaxis()->SetRange(0,100);
      Float_t mean = cfd->GetMean();
      Float_t rms = cfd->GetRMS();
      Float_t hmin=mean - 3*rms;
      Float_t hmax =mean + 3*rms;
      cfd->GetXaxis()->SetRange(hmin-10,hmax+10);
      cout<<" cfd range "<<mean<<" rms "<<rms<<" "<<hmin<<" "<<hmax<<endl;
      //     TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
      //  cfd->Fit("g1","RQ");
      cfd->Draw();
    }
  */

}

void plotTree(TTree *tree_, std::string whichfit, std::string selectString){

	// Create a map for plotting the pullsummaries:
	std::map < const char*, std::pair <double,double> > pullSummaryMap;
	int nPulls=0;

	TObjArray *l_branches = tree_->GetListOfBranches();
	int nBranches = l_branches->GetEntries();

	gStyle->SetPadTopMargin(0.01);

	TCanvas *c = new TCanvas("c","",960,800);

	std::string treename = tree_->GetName();
	c->SaveAs(Form("%s.pdf[",treename.c_str()));
	// File to store plots in 
	TFile *fOut = new TFile(Form("%s.root",treename.c_str()),"RECREATE");

	for (int iobj=0;iobj<nBranches;iobj++){

		TBranch *br =(TBranch*) l_branches->At(iobj);

		// Draw the normal histogram
		const char* name = br->GetName();

                // names with - are not allowed
                string namestr(name);
                if(namestr.find("-")!=string::npos) {
                  std::cout << "Variable " << name << " contains a bad character: -. Skipping. " << std::endl;
                  continue;
                }
		bool fitPull=false;
		bool fitPullf=false;

		bool plotLH=false;

		TGraph *gr=NULL;
		double p_mean =0;
		double p_err  =0;

		int nToysInTree = tree_->GetEntries();
		// Find out if paramter is fitted value or constraint term
		bool isFitted = findNuisancePre(name);
		if (doPull && isFitted){
			
			p_mean = bfvals_[name].first;	// toy constrainits thrown about best fit to data
			if(namestr.find("n_exp")==string::npos) p_err  = prevals_[name].second; // uncertainties taken from card
			std::cout << "******* "<< name << " *******"<<std::endl;
			std::cout << p_mean <<  " " << p_err << std::endl;
			std::cout << "******************************" <<std::endl;

			const char* drawInput;
                        // if the parameter is a normalization, the error is not available. Do the residual instead of the pull
                        if(namestr.find("n_exp")!=string::npos) drawInput = Form("(%s-%f)/%f",name,p_mean,p_mean);
                        else drawInput = Form("(%s-%f)/%f",name,p_mean,p_err);
			tree_->Draw(Form("%s>>%s",drawInput,name),"");
			tree_->Draw(Form("%s>>%s_fail",drawInput,name),selectString.c_str(),"same");
			fitPull  = true;
			fitPullf = true;
			if (doLH) {
			  gr = graphLH(name,p_err,whichfit);
			  if (gr) plotLH=true;
			}
			
		}

		else{
			tree_->Draw(Form("%s>>%s",name,name),"");
			tree_->Draw(Form("%s>>%s_fail",name,name),selectString.c_str(),"same");
		}
		

		TH1F* bH  = (TH1F*) gROOT->FindObject(Form("%s",name))->Clone();
		TH1F* bHf = (TH1F*) gROOT->FindObject(Form("%s_fail",name))->Clone();
		bHf->SetLineColor(2);
		bH->GetXaxis()->SetTitle(bH->GetTitle());
		bH->GetYaxis()->SetTitle(Form("no toys (%d total)",nToysInTree));
		bH->GetYaxis()->SetTitleOffset(1.05);
		bH->GetXaxis()->SetTitleOffset(0.9);
		bH->GetYaxis()->SetTitleSize(0.05);
		bH->GetXaxis()->SetTitleSize(0.05);
		if (isFitted) {bH->GetXaxis()->SetTitle(Form("(%s-#theta_{B})/#sigma_{#theta}",name));}
		else {bH->GetXaxis()->SetTitle(Form("%s",name));}
		
		bH->SetTitle("");	

		if ( bH->Integral() <=0 )  fitPull = false;
		if (fitPull) {bH->Fit("gaus"); bH->GetFunction("gaus")->SetLineColor(4);}
		
		if ( bHf->Integral() <=0 )  fitPullf = false;
		if (fitPullf) {bHf->Fit("gaus"); bHf->GetFunction("gaus")->SetLineColor(2);}

		c->Clear();
		//TPad pad1("t1","",0.01,0.02,0.59,0.98);
		// Pad 1 sizes depend on the parameter type ...
		double pad1_x1,pad1_x2,pad1_y1,pad1_y2;
		if ( !isFitted ) {
			 pad1_x1 = 0.01; 
			 pad1_x2 = 0.98; 
			 pad1_y1 = 0.045; 
//.........这里部分代码省略.........

int
TrackParametrization( TString csvfile="fitslices_out.csv" )
{

  /* Read data from input file */
  TTree *tres = new TTree();
  tres->ReadFile( csvfile, "ptrue:etatrue:psig:psig_err:pmean:pmean_err:norm", ',' );

  /* Print read-in tree */
  tres->Print();

  /* colors array */
  unsigned colors[8] = {1,2,3,4,6,7,14,16};

  /* Create vector of theta values to include for visualization*/
  vector< double > etas_vis;
  etas_vis.push_back(-2.75);
  etas_vis.push_back(-2.25);
  etas_vis.push_back(-1.75);
  etas_vis.push_back(-0.25);
  etas_vis.push_back( 0.25);
  etas_vis.push_back( 1.75);
  etas_vis.push_back( 2.25);

//  etas_vis.push_back(-3.25);
//  etas_vis.push_back(-2.25);
//  etas_vis.push_back(-1.25);
//  etas_vis.push_back(-0.25);
//  etas_vis.push_back( 0.25);
//  etas_vis.push_back( 1.25);
//  etas_vis.push_back( 2.25);
//  etas_vis.push_back( 3.25);

  /* Create vector of theta values to include for fitting*/
  vector< double > etas_fit;
  for ( double eta = -4.45; eta < 4.5; eta += 0.1 )
    etas_fit.push_back( eta );

  /* Create fit function */
  TF1* f_momres = new TF1("f_momres", "sqrt( [0]*[0] + [1]*[1]*x*x )" );

  cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;

  /* Create scratch canvas */
  TCanvas *cscratch = new TCanvas("cscratch");

  /* Create framehistogram */
  TH1F* hframe = new TH1F("hframe","",100,0,40);
  hframe->GetYaxis()->SetRangeUser(0,0.15);
  hframe->GetYaxis()->SetNdivisions(505);
  hframe->GetXaxis()->SetTitle("Momentum (GeV/c)");
  hframe->GetYaxis()->SetTitle("#sigma_{p}/p");

  /* create combined canvas plot */
  TCanvas *c1 = new TCanvas();
  hframe->Draw();

  /* Create legend */
  TLegend* leg_eta = new TLegend( 0.2, 0.6, 0.5, 0.9);
  leg_eta->SetNColumns(2);

  /* Create ofstream to write fit parameter results */
  ofstream ofsfit("track_momres_new.csv");
  ofsfit<<"eta,par1,par1err,par2,par2err"<<endl;

  /* Create resolution-vs-momentum plot with fits for each selected theta value */
  for ( int i = 0; i < etas_fit.size(); i++ )
    {
      /* Switch to scratch canvas */
      cscratch->cd();

      double eta = etas_fit.at(i);

      /* No tracking outside -4 < eta < 4 */
      if ( eta < -4 || eta > 4 )
	continue;

      cout << "\n***Eta = " << eta << endl;

      /* Define range of theta because float comparison with fixed value doesn't work
	 too well for cuts in ROOT trees */
      double eta_min = eta * 0.999;
      double eta_max = eta * 1.001;

      /* Cut for tree */
      TCut cutx( Form("ptrue > 1 && ( (etatrue > 0 && (etatrue > %f && etatrue < %f)) || (etatrue < 0 && (etatrue < %f && etatrue > %f)) )", eta_min, eta_max, eta_min, eta_max) );

      /* "Draw" tree on scratch canvas to fill V1...V4 arrays */
      tres->Draw("psig:ptrue:psig_err:0", cutx );

      /* Create TGraphErrors with selected data from tree */
      TGraphErrors *gres = new TGraphErrors( tres->GetEntries(cutx),
					     &(tres->GetV2())[0],
					     &(tres->GetV1())[0],
					     &(tres->GetV4())[0],
					     &(tres->GetV3())[0] );

      /* reset function parameters before fit */
      f_momres->SetParameter(0,0.1);
      f_momres->SetParameter(1,0.1);
//.........这里部分代码省略.........

void DrawBoosted(TString VAR, float XMIN, float XMAX, int REBIN, TString XTITLE)
{
  gROOT->ForceStyle();
  const int N = 7;
  float XSEC[N] = {1.74e+6,3.67e+5,2.94e+4,6.524e+03,1.064e+03,121.5,2.542e+01};

  TString SAMPLE[N] = {
    "QCD_HT200to300",
    "QCD_HT300to500",
    "QCD_HT500to700", 
    "QCD_HT700to1000",
    "QCD_HT1000to1500",
    "QCD_HT1500to2000",
    "QCD_HT2000toInf" 
  };
 
  TFile *inf[N];
  TH1F  *h0[N],*h1[N];
  TH1F  *hQCD0,*hQCD1; 

  TCanvas *can = new TCanvas("Boosted","Boosted",900,600);
  can->cd(1);
  can->SetBottomMargin(0.3);
  //can->SetRightMargin(0.15);
  for(int k=0;k<N;k++) {
    inf[k] = TFile::Open("Histo_"+SAMPLE[k]+".root");
    h0[k]  = (TH1F*)inf[k]->Get("hadtopBoost/h_"+VAR+"_Cut_ctl");
    h0[k]->Sumw2();
    h0[k]->Rebin(REBIN);
    h1[k]  = (TH1F*)inf[k]->Get("hadtopBoost/h_"+VAR+"_Cut_sig");
    h1[k]->Sumw2();
    h1[k]->Rebin(REBIN);
    h0[k]->Scale(XSEC[k]/((TH1F*)inf[k]->Get("hadtopBoost/TriggerPass"))->GetBinContent(1));
    h1[k]->Scale(XSEC[k]/((TH1F*)inf[k]->Get("hadtopBoost/TriggerPass"))->GetBinContent(1)); 
    cout<<SAMPLE[k]<<": "<<h0[k]->GetEntries()<<" "<<h0[k]->Integral()<<endl;
    if (k == 0) {
      hQCD0 = (TH1F*)h0[k]->Clone();
      hQCD1 = (TH1F*)h1[k]->Clone();
    }
    if (k > 0) {
      hQCD0->Add(h0[k]);
      hQCD1->Add(h1[k]);
    }
  } 
  hQCD0->Scale(1/hQCD0->Integral());
  hQCD1->Scale(1/hQCD1->Integral());
  double max1 = TMath::Max(hQCD0->GetBinContent(hQCD0->GetMaximumBin()),hQCD1->GetBinContent(hQCD1->GetMaximumBin()));
  hQCD0->SetMaximum(1.2*max1);
  hQCD0->SetMinimum(1e-4);
  hQCD0->GetXaxis()->SetLabelSize(0.0);
  hQCD0->GetXaxis()->SetRangeUser(XMIN,XMAX); 
  hQCD0->Draw("HIST");
  hQCD1->Draw("sameE");

  TLegend *leg = new TLegend(0.6,0.7,0.9,0.9);
  leg->SetHeader("QCD Closure");
  leg->SetBorderSize(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.05);
  leg->AddEntry(hQCD1,"Signal sample","P");
  leg->AddEntry(hQCD0,"Control sample","F");
  leg->Draw();

  TH1F *hRatio = (TH1F*)hQCD1->Clone("Ratio");
  hRatio->Divide(hQCD0);

  hRatio->SetLineColor(kBlack);
  hRatio->SetMarkerColor(kBlack);

  TPad* pad = new TPad("pad", "pad", 0., 0., 1., 1.);
  pad->SetTopMargin(0.7);
  //pad->SetRightMargin(0.15);
  pad->SetFillColor(0);
  pad->SetFillStyle(0);
  pad->Draw();
  pad->cd(0);
  gPad->SetGridy();
  hRatio->GetXaxis()->SetTitle(XTITLE);
  hRatio->GetXaxis()->SetRangeUser(XMIN,XMAX); 
  hRatio->GetYaxis()->SetNdivisions(505);
  hRatio->GetYaxis()->SetRangeUser(0,2);
  hRatio->GetYaxis()->SetLabelSize(0.04);
  hRatio->Draw();
  hRatio->Draw("same");
}

void fractionFit()
{
   char name[1000];
  sprintf(name,"/Users/zach/Research/pythia/npeTemplate/outputs/currentB.root");
  TFile *fB = new TFile(name,"READ");
  sprintf(name,"/Users/zach/Research/pythia/npeTemplate/outputs/currentC.root");
  TFile *fC = new TFile(name,"READ");
   sprintf(name,"/Users/zach/Research/rootFiles/run12NPEhPhi/currentData.root");
  TFile *fD = new TFile(name,"READ");
  if (fB->IsOpen()==kFALSE || fC->IsOpen()==kFALSE)
    { std::cout << "!!!!!! Either B,C, or Data File not found !!!!!!" << std::endl
		<< "Looking for currentB.root, currentC.root, and currentData.root" << std::endl;
      exit(1); }
  
  // Set constants and projection bins
  const Int_t numPtBins = 10;
  Float_t lowpt[14] ={2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.5,10.,14.0};
  Float_t highpt[14]={3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.5,10.,14.,200.};
  Float_t hptCut=0.5;

  // Make Canvases
  TCanvas* deltaPhi  = new TCanvas("deltaPhi","Pythia Delta Phi",150,0,1150,1000);
  TCanvas* fitResult0 = new TCanvas("fitResult0","RB Extraction HT0",150,0,1150,1000);
  //  TCanvas* fitResult2 = new TCanvas("fitResult2","RB Extraction HT2",150,0,1150,1000);
  deltaPhi->Divide(2,2);
  fitResult0->Divide(4,3);
  //fitResult2->Divide(4,3);

  // Make histos
  TH1D* projB[numPtBins];
  TH1D* projC[numPtBins];
  TH1D* projData0[numPtBins];
  TH1D* projData2[numPtBins];
  
  // Get and Draw histos
  TPaveText* lbl[numPtBins];
  char textLabel[100];
  Int_t plotbin;

  for(Int_t ptbin=0; ptbin<numPtBins; ptbin++)
    {
      if(ptbin!=0 && ptbin!=2 && ptbin!=4 && ptbin!=6)
	continue;
      if(ptbin==0)plotbin=0;
      if(ptbin==2)plotbin=1;
      if(ptbin==4)plotbin=2;
      if(ptbin==6)plotbin=3;
      // Init necessary plotting tools
      lbl[ptbin] = new TPaveText(.2,.76,.5,.82,Form("NB NDC%i",ptbin));
      sprintf(textLabel,"%.1f < P_{T,e} < %.1f",lowpt[ptbin],highpt[ptbin]);
      lbl[ptbin]->AddText(textLabel);
      lbl[ptbin]->SetFillColor(kWhite);

      projB[ptbin] = (TH1D*)fB->Get(Form("projDelPhi_%i",ptbin));
      projC[ptbin] = (TH1D*)fC->Get(Form("projDelPhi_%i",ptbin));
      projData0[ptbin]= (TH1D*)fD->Get(Form("NPEhDelPhi_0_%i",ptbin));
      projData2[ptbin]= (TH1D*)fD->Get(Form("NPEhDelPhi_2_%i",ptbin));
      Int_t RB = 8;
      projB[ptbin]->Rebin(RB);
      projC[ptbin]->Rebin(RB);
      
      // Draw Templates on own plots
      deltaPhi->cd(plotbin+1);
      projData0[ptbin]->SetLineColor(kBlue);
      projData2[ptbin]->SetLineColor(kGreen+3);
      projB[ptbin]->SetLineColor(kRed);
      projC[ptbin]->SetLineColor(kBlack);
      //      projC[ptbin]->GetYaxis()->SetRangeUser(0.,0.5);
      projC[ptbin]    -> Draw();
      projB[ptbin]    -> Draw("same");
      projData0[ptbin]-> Draw("same");
      projData2[ptbin]-> Draw("same");
      lbl[ptbin]      -> Draw("same");

      TLegend* leg = new TLegend(0.5,0.73,0.85,0.85);
      leg->AddEntry(projB[ptbin],"b#bar{b}->NPE","lpe");
      leg->AddEntry(projC[ptbin],"c#bar{c}->NPE","lpe");
      leg->AddEntry(projData0[ptbin],"HT0","lpe");
      leg->AddEntry(projData2[ptbin],"HT2","lpe");
      leg->Draw();

      // Do the actual fit
      TObjArray *mc = new TObjArray(2);   // MC histograms are put in this array
      mc->Add(projC[ptbin]);
      mc->Add(projB[ptbin]);

      fitResult0->cd(ptbin+1);
      TFractionFitter* fit = new TFractionFitter(projData0[ptbin], mc,"V"); // initialise
      fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
      fit->SetRangeX(46,56);                    // use only the first 15 bins in the fit
      Int_t status = fit->Fit();               // perform the fit
      std::cout << "fit status: " << status << std::endl;
      if (status == 0) {                       // check on fit status
	TH1F* result = (TH1F*) fit->GetPlot();
	projData0[ptbin]->Draw("Ep");
	result->Draw("same");
	}

      
      /* fitResult2->cd(ptbin+1);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        else out_efficiency<<leg_names[files]<<"&\t  "<<histos[4]->Integral()*1000. <<"&\t"<< histos[4]->Integral()/xsec[files] <<endl;
        if (files==nfiles-1) out_efficiency<<"Total BG"<<"&\t \t \t  "<<totalBG*1000.<< endl;
        files++;
    } while (files<nfiles);
    out_efficiency.close();

    for (int d=0; d<nhistos; d++) {
        discriminators[d] = Discr(discr_histos[d],signal_histos[d]);
    }

    bubblesort(discriminators, hist_names_sort,nhistos);

    out_discrimination.open(trigger+dir_name+"discrimination.txt");
    for (int d=0; d<nhistos; d++) {
        if (d==0) out_discrimination<<"Variable &\t d"<<endl;
        out_discrimination<<"$"<<hist_names_sort[d]<<"$"<<" & \t "<< std::setprecision(2)<< discriminators[d]<<endl;
    }
    out_discrimination.close();
    TLatex* tex = new TLatex(0.90,0.92,"13 TeV, PU = 20, bx = 25 ns, 1 pb^{-1}");
    tex->SetNDC();
    tex->SetTextAlign(35);
    tex->SetTextFont(42);
    tex->SetTextSize(0.04);
    tex->SetLineWidth(2);
    TLatex *tex1 = new TLatex(0.13,0.83,"CMS");
    tex1->SetNDC();
    tex1->SetTextAlign(20);
    tex1->SetTextFont(61);
    tex1->SetTextSize(0.06);
    tex1->SetLineWidth(2);
    TLatex* tex2 = new TLatex(0.22,0.77,"Work in progress");
    tex2->SetNDC();
    tex2->SetTextAlign(20);
    tex2->SetTextFont(52);
    tex2->SetTextSize(0.04);
    tex2->SetLineWidth(2);
    TLatex* tex_file = new TLatex(0.35,0.92,"Spring15, DoubleBtag");
    //	TLatex* tex_file = new TLatex(0.35,0.92,"Spring15, SingleBtag");
    tex_file->SetNDC();
    tex_file->SetTextAlign(35);
    tex_file->SetTextFont(42);
    tex_file->SetTextSize(0.04);
    tex_file->SetLineWidth(2);
    for (int i=0; i<nhistos; i++) {
        TString temp_str;
        temp_str.Form("%2.2f",Discr(discr_histos[i],signal_histos[i]));
        TString disc_value = temp_str.Prepend(", d = ");
        TLatex *disc_value_text = new TLatex(0.86,0.853,disc_value);
        disc_value_text->SetNDC();
        disc_value_text->SetTextAlign(35);
        disc_value_text->SetTextFont(42);
        disc_value_text->SetTextSize(0.04);
        disc_value_text->SetLineWidth(2);

        TCanvas *c1 = new TCanvas();
        c1->SetBottomMargin(.12);
        c1->cd();
        c1->SetLogy();
        Double_t xmin = signal_histos[i]->GetBinCenter(0);
        Double_t xmax = signal_histos[i]->GetBinCenter(signal_histos[i]->GetNbinsX());
        TH1F *frame = new TH1F("frame","",1,xmin,xmax);
        frame->SetMinimum(1e-4);
        frame->SetMaximum(1e10);
        frame->GetYaxis()->SetTitleOffset(0.9);
        frame->GetXaxis()->SetTitleOffset(0.91);
        frame->SetStats(0);
        frame->SetTitleFont(42,"x");
        frame->SetTitleFont(42,"y");
        frame->SetTitleSize(0.05, "XYZ");
        frame->SetXTitle(signal_histos[i]->GetXaxis()->GetTitle());

        if ((i<4))frame->SetYTitle("Events  /  20 GeV");
        else if (i==37)frame->SetYTitle("Events  /  20 GeV");
        else if (i==41)frame->SetYTitle("Events  /  20 GeV");
        else if ((i>=4) && (i<8)) frame->SetYTitle("Events  /  0.5");
        else if (i==15) frame->SetYTitle("Events  /  0.5");
        else if ((i>=8) && (i<12)) frame->SetYTitle("Events  /  0.2");
        else if (i==40) frame->SetYTitle("Events / 100 GeV");
        else if ((i==10)||(i==12)||(i==16)||(i==17)||(i==18) || (i==20) || (i==25) || (i==26) || (i==27) || (i==35) || (i==36) ||(i==38) || (i==39)) frame->SetYTitle("Events / 10 GeV");
        else if ((i==13)||(i==14)) frame->SetYTitle("Events / 0.1");
        else if ((i==21)||(i==22)) frame->SetYTitle("Events / 0.05");
        else if ((i==23)||(i==24)) frame->SetYTitle("Events / 0.01");
        else if (i==28) frame->SetYTitle("Events / 0.02");
        else if ((i==29) || (i==30) || (i==31) || (i==32) )frame->SetYTitle("Events / 0.1");
        else if ((i==33) || (i==33))frame->SetYTitle("Events / 0.01");
        else frame->SetYTitle("Events");
        frame->Draw();
        tex->Draw();
        tex1->Draw();
        tex2->Draw();
        tex_file->Draw();
        stacks[i]->Draw("nostacksame");
        leg->Draw("same");
        signal_histos[i]->Draw("same");
        disc_value_text->Draw();
        c1->Print(output_names[i]);
        c1->Delete();
    }

}

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



  /**********************
      ADD BRIANS PLOT
    ********************/

  NUM_TOT_POSITRONS = 1e+07; 
  BIAS = 1e+04; 
  BINNING_WEIGHT = POSITRONS_PER_SEC/(BIAS*NUM_TOT_POSITRONS);
  B = pow(1-pow(GAMMA_PLUS, -2.), .5);

  TFile* fileinelastic = new TFile("brian_1e7_pos_1e4_bias.root");
  TTree* Hits_Infoinelastic = (TTree *)fileinelastic->Get("Signal");

  Hits_Infoinelastic->SetBranchAddress("numHits", &numHit);
  Hits_Infoinelastic->SetBranchAddress("energyTot", &energy);
  Hits_Infoinelastic->SetBranchAddress("XPosition", &xPos);
  Hits_Infoinelastic->SetBranchAddress("YPosition", &yPos);
  Hits_Infoinelastic->SetBranchAddress("Particle_ID", &pID);
  Hits_Infoinelastic->SetBranchAddress("Theta", &theta);

  TH1D* hminelastic = new TH1D("M_{A'}^{2}" ,             // plot label
                       "Inelastic (only registering gammas)",    // title
                       nM2Bins,                   // x number of bins
                       nM2Min,                    // x lower bound
                       nM2Max);                   // x upper bound
  // go through all entries and fill the histograms
  nentries = Hits_Infoinelastic->GetEntries();
  for (int i=0; i<nentries; i++) {
    Hits_Infoinelastic->GetEntry(i);
    if (pID == 22) { // gammas only
      theta*= TMath::Pi()/180; //radians
      cout << "M^2 is: " << mSquared(energy, theta) << endl;
      cout << "Energy is: " << energy << endl;
           
      hminelastic->Fill(mSquared(energy, theta), BINNING_WEIGHT / dM2BinSize);
    }
  }
  cout << "DM2binsize is:" << dM2BinSize << endl;


  hminelastic->SetFillColor(kGreen);
  hminelastic->SetFillStyle(3001);
  hminelastic->GetXaxis()->SetTitle("M_{A'}^{2} (MeV^{2})");
  hminelastic->GetYaxis()->SetTitle("Photons per MeV^{2} per Second (MeV^{-2} s^{-1})");
  hminelastic->GetXaxis()->CenterTitle();
  hminelastic->GetYaxis()->CenterTitle();



 hs->Add(hminelastic);
  hs->Add(hYIMIN);
      hs->Add(hmyyy);
  hs->Add(hmepluseminusy);
    hs->Add(hmyy);



  // create canvas and draw histogram
  TCanvas* canvas = new TCanvas("canvas", "canvas", 700, 700);
  canvas->Divide(3,2);
  TPad* p;

  p = (TPad*)canvas->cd(1);
  //p->SetGrid();
  p->SetLogy();

 


  hs->Draw();

    hs->GetXaxis()->SetTitle("M_{A'}^{2} (MeV^{2})");
  hs->GetYaxis()->SetTitle("Photons per MeV^{2} per Second (MeV^{-2} s^{-1})");
  hs->GetXaxis()->CenterTitle();
  hs->GetYaxis()->CenterTitle();
  
  p = (TPad*)canvas->cd(4);
  hmepluseminusy->Draw();
  
  p = (TPad*)canvas->cd(3);
  hmyyy->Draw();
  p = (TPad*)canvas->cd(2);
  hmyy->Draw();

  p = (TPad*)canvas->cd(5);
  hYIMIN->Draw();

  p = (TPad*)canvas->cd(6);
  hminelastic->Draw();

/*
  p = (TPad*)canvas->cd(4);
  p->SetLogy();
  p->SetGrid();
  hm2->Draw();*/


}

//.........这里部分代码省略.........
      ptHatCorr -> Add(projpthatall,wt);

      // pt bin independent
      delPhi2535 -> Add(temp,wt); // still need weight from pthat
      //trigCount -> Scale(wt);     // "      "    "     "     "
      norm2535 += trigCount->GetBinContent(1);
      //delPhi2535 -> Scale(1./norm2535);
   
      // Analyze each ptH bin individually, adding to the overall hists
      for(Int_t ptbin=0; ptbin<numPtBins; ptbin++)
	{
	  // DEBUGcout << "pthbin: " << pthBin << " ptbin: " << ptbin << endl;
	  projDelPhi[ptbin] = mh3delPhi->ProjectionZ(Form("projDelPhi_%i",ptbin),mh3delPhi->GetXaxis()->FindBin(lowpt[ptbin]),mh3delPhi->GetXaxis()->FindBin(highpt[ptbin])-1,mh3delPhi->GetYaxis()->FindBin(hptCut),mh3delPhi->GetYaxis()->FindBin(hptMax));
	  projNpeY[ptbin]   = mh2npePt->ProjectionY(Form("projNpeY_%i",ptbin),mh2npePt->GetXaxis()->FindBin(lowpt[ptbin]),mh2npePt->GetXaxis()->FindBin(highpt[ptbin])-1);
	  projptHat[ptbin]  = mh2ptHatPt->ProjectionY(Form("projPtHat_%i",ptbin),mh2ptHatPt->GetXaxis()->FindBin(lowpt[ptbin]),mh2ptHatPt->GetXaxis()->FindBin(highpt[ptbin])-1);
	
	  delPhi[ptbin] -> Add(projDelPhi[ptbin],wt);
	  NpeY[ptbin] -> Add(projNpeY[ptbin],wt);
	  
	  // Calculate scaling Factor
	  Double_t Norm = NpeY[ptbin]->Integral();
	  ptNorm->SetBinContent(ptNorm->GetBin(ptbin+1),Norm);
	  totalNorm[ptbin] += Norm;
	 
	}
    }

  // For making plots

  ptHatC->cd(1);
  gPad-> SetLogy();
  ptHat->GetXaxis()->SetTitle("pT-Hat (GeV/c)");
  ptHat->SetTitle("Raw pT Hat");
  ptHat->Draw();
  ptHatC->cd(2);
  gPad-> SetLogy();
  ptHatCorr->GetXaxis()->SetTitle("pT-Hat (GeV/c)");
  ptHatCorr->SetTitle("Weighted pT Hat");
  ptHatCorr->Draw();
  
  for(Int_t ptbin=0; ptbin<numPtBins; ptbin++)
    {
      // Init necessary plotting tools
      lbl[ptbin] = new TPaveText(.2,.8,.5,.85,Form("NB NDC%i",ptbin));
      sprintf(textLabel,"%.1f < P_{T,e} < %.1f",lowpt[ptbin],highpt[ptbin]);
      lbl[ptbin]->AddText(textLabel);
      lbl[ptbin]->SetFillColor(kWhite);

      deltaPhi->cd(ptbin+1);
      delPhi[ptbin]->GetXaxis()->SetTitle("#Delta#phi_{eh}");
      // delPhi[ptbin]->Sumw2();
      //cout << totalNorm[ptbin] << endl;
      //delPhi[ptbin]->Scale(wt);
      delPhi[ptbin]->GetYaxis()->SetTitle("1/N_{NPE} #upoint dN/d(#Delta)#phi");
      delPhi[ptbin]->GetXaxis()->SetRangeUser(-3.5,3.5);
      if(ptbin == 0)
	{
	  if(mode == 1)
	    delPhi[ptbin]->SetTitle("Pythia NPE-had #Delta#phi - c/#bar{c}");
	  if(mode == 2)
	    delPhi[ptbin]->SetTitle("Pythia NPE-had #Delta#phi - b/#bar{b}");
	}
      else
	delPhi[ptbin]->SetTitle("");
      if(ptbin < 13){
	delPhi[ptbin]->Draw("E");

void plotTreeNorms(TTree *tree_, std::string selectString, bool do7TeV){

	// Create a map for plotting the pullsummaries:
	std::map < const char*, std::pair <double,double> > pullSummaryMap;
	int nPulls=0;

	TObjArray *l_branches = tree_->GetListOfBranches();
	int nBranches = l_branches->GetEntries();

	gStyle->SetPadTopMargin(0.01);

	TCanvas *c = new TCanvas("c","",960,800);

	std::string treename = tree_->GetName();
	c->SaveAs(Form("%s_normresiduals.pdf[",treename.c_str()));
	// File to store plots in 
	TFile *fOut = new TFile(Form("%s_normresiduals.root",treename.c_str()),"RECREATE");

        TH1F *bHd = new TH1F("bHd","",50,-1.0,1.0);
        TH1F *bHfd = new TH1F("bHfd","",50,-1.0,1.0);

	for (int iobj=0;iobj<nBranches;iobj++){

		TBranch *br =(TBranch*) l_branches->At(iobj);

		// Draw the normal histogram
		const char* name = br->GetName();

                // select only the normalizations
                string namestr(name);
                if(namestr.find("n_exp")==string::npos) continue;

                bool fitPull=true;
                bool fitPullf=true;

		double p_mean =0;

		int nToysInTree = tree_->GetEntries();
		// Find out if paramter is fitted value or constraint term.
                bool isFitted = true;
			
                p_mean = prenorms_[name].first;	// toy initial parameters from the datacards
                std::cout << "******* "<< name << " *******"<<std::endl;
                std::cout << p_mean << std::endl;
                std::cout << "******************************" <<std::endl;

                TH1F* bH = (TH1F*)bHd->Clone(Form("%s",name));
                TH1F* bHf = (TH1F*)bHfd->Clone(Form("%s_fail",name));
                
                const char* drawInput = Form("(%s-%f)/%f",name,p_mean,p_mean);
                tree_->Draw(Form("%s>>%s",drawInput,name),"");
                tree_->Draw(Form("%s>>%s_fail",drawInput,name),selectString.c_str(),"same");
                fitPull  = true;
                fitPullf = true;
                  
		bHf->SetLineColor(2);
		bH->GetXaxis()->SetTitle(bH->GetTitle());
		bH->GetYaxis()->SetTitle(Form("no toys (%d total)",nToysInTree));
		bH->GetYaxis()->SetTitleOffset(1.05);
		bH->GetXaxis()->SetTitleOffset(0.9);
		bH->GetYaxis()->SetTitleSize(0.05);
		bH->GetXaxis()->SetTitleSize(0.05);
		bH->GetXaxis()->SetTitle(Form("%s",name));
                
		
		bH->SetTitle("");	

		if ( bH->Integral() <=0 )  fitPull = false;
		if (fitPull) {bH->Fit("gaus"); bH->GetFunction("gaus")->SetLineColor(4);}
		
		if ( bHf->Integral() <=0 )  fitPullf = false;
		if (fitPullf) {bHf->Fit("gaus"); bHf->GetFunction("gaus")->SetLineColor(2);}

		c->Clear();
		
 		TPad pad1("t1","",0.01,0.01,0.66,0.95);
 		TPad pad2("t2","",0.70,0.20,0.98,0.80);

		pad1.SetNumber(1); pad2.SetNumber(2);

                if ( isFitted ) {pad2.Draw();}

		pad1.Draw();
		pad1.SetGrid(true);


		TLatex *titletext = new TLatex();titletext->SetNDC();

		
		c->cd(1); bH->Draw(); bHf->Draw("same");
		TLegend *legend = new TLegend(0.6,0.8,0.9,0.89);
		legend->SetFillColor(0);
		legend->AddEntry(bH,"All Toys","L");
		legend->AddEntry(bHf,selectString.c_str(),"L");
		legend->Draw();

		if (fitPull){
			c->cd(2);
			double gap;
			TLatex *tlatex = new TLatex(); tlatex->SetNDC(); 
//.........这里部分代码省略.........

void met_cuts(double pt1min, double pt2min, double METmin, int jets, int DataFlag){


  TLegend* leg = new TLegend(0.5,0.72,0.89,0.89);
  leg->SetNColumns(2);


  TCanvas *canvas = new TCanvas("c1n","",500,500);
 
  gPad->SetLogy();
  TCut mggmax = "mgg<200";
  TCut mggmin = "mgg>100";
  TCut pt1Cut = Form("pt1>%lf",pt1min);
  TCut pt2Cut = Form("pt2>%lf",pt2min);
  TCut METCut = Form("t1pfmet>%lf",METmin);
  
  TCut eveto1 = "eleveto1 == 1";
  TCut eveto2 = "eleveto2 == 1";
  TCut eveto = eveto1 && eveto2;
  TCut genmatch = "((genmatch1==1 && genmatch2==0)||(genmatch1==0 && genmatch2==1)||(genmatch1==0 && genmatch2==0))";  
  TCut metF = "((metF_GV==1) && (metF_HBHENoise==1) && (metF_HBHENoiseIso==1) && (metF_CSC==1) && (metF_eeBadSC==1))";  
  TCut Cut_Jets = Form("nJets<%d",jets);
  /* 
     TFile *M1 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M1.root","READ");  
     TFile *M10 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M10.root","READ");
     TFile *M100 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M100.root","READ");  
     TFile *M1000 =  TFile::Open("./50ns_betaV4/NewWeightDMHtoGG_M1000.root","READ");
  */

  TFile *data = TFile::Open("./25ns_2246inv/DoubleEG.root","READ");
  
  TFile *sig1 = TFile::Open("./25ns_2246inv/2HDM_mZP600.root","READ");
  TFile *sig2 = TFile::Open("./25ns_2246inv/2HDM_mZP800.root","READ");
  TFile *sig3 = TFile::Open("./25ns_2246inv/2HDM_mZP1000.root","READ");
  TFile *sig4 = TFile::Open("./25ns_2246inv/2HDM_mZP1200.root","READ");
  TFile *sig5 = TFile::Open("./25ns_2246inv/2HDM_mZP1400.root","READ");  
  
  /*
    TFile *bkg1 =  TFile::Open("./50ns_betaV4/DiPhoton.root","READ");  
    TFile *bkg2 =  TFile::Open("./50ns_betaV4/DYJetsToLL.root","READ");  
    TFile *bkg3 =  TFile::Open("./50ns_betaV4/GJets.root","READ");  
    TFile *bkg4 =  TFile::Open("./50ns_betaV4/GluGluHToGG.root","READ");  
    TFile *bkg5 =  TFile::Open("./50ns_betaV4/QCD.root","READ");  
    TFile *bkg6 =  TFile::Open("./50ns_betaV4/VH.root","READ");  
  */

  TFile *bkg1 =  TFile::Open("./25ns_2246inv/DiPhoton.root","READ");  
  TFile *bkg2 =  TFile::Open("./25ns_2246inv/DYJetsToLL.root","READ");  
  TFile *bkg3 =  TFile::Open("./25ns_2246inv/GJets.root","READ");  
  TFile *bkg4 =  TFile::Open("./25ns_2246inv/GluGluHToGG.root","READ");  
  TFile *bkg5 =  TFile::Open("./25ns_2246inv/QCD.root","READ");  
  TFile *bkg6 =  TFile::Open("./25ns_2246inv/VH.root","READ");  
  
  /*
    TTree *tree_M1 = (TTree*) M1->Get("DiPhotonTree");
    TTree *tree_M10 = (TTree*) M10->Get("DiPhotonTree");
    TTree *tree_M100 = (TTree*) M100->Get("DiPhotonTree");
    TTree *tree_M1000 = (TTree*) M1000->Get("DiPhotonTree");
  */

  TTree *tree_data = (TTree*) data->Get("DiPhotonTree");
  
  TTree *tree_sig1 = (TTree*) sig1->Get("DiPhotonTree");
  TTree *tree_sig2 = (TTree*) sig2->Get("DiPhotonTree");
  TTree *tree_sig3 = (TTree*) sig3->Get("DiPhotonTree");
  TTree *tree_sig4 = (TTree*) sig4->Get("DiPhotonTree");
  TTree *tree_sig5 = (TTree*) sig5->Get("DiPhotonTree");
  
  
  TTree *tree_bkg1 = (TTree*) bkg1->Get("DiPhotonTree");
  TTree *tree_bkg2 = (TTree*) bkg2->Get("DiPhotonTree");
  TTree *tree_bkg3 = (TTree*) bkg3->Get("DiPhotonTree");
  TTree *tree_bkg4 = (TTree*) bkg4->Get("DiPhotonTree");
  TTree *tree_bkg5 = (TTree*) bkg5->Get("DiPhotonTree");
  TTree *tree_bkg6 = (TTree*) bkg6->Get("DiPhotonTree");
  
  
  
  /* 
     tree_M1->Draw("(t1pfmet)>>h1(30,100,200)","weight*10"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
     TH1F *h1 =(TH1F*)gPad->GetPrimitive("h1");
     tree_M10->Draw("(t1pfmet)>>h2(30,100,200)","weight*10"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
     TH1F *h2 =(TH1F*)gPad->GetPrimitive("h2");
     tree_M100->Draw("(t1pfmet)>>h3(30,100,200)","weight*10"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
     TH1F *h3 =(TH1F*)gPad->GetPrimitive("h3");
     tree_M1000->Draw("(t1pfmet)>>h4(30,100,200)","weight*10"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
     TH1F *h4 =(TH1F*)gPad->GetPrimitive("h4");  
  */
  
  tree_data->Draw("(t1pfmet)>>hdata(60,0,900)",(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets && metF));
  TH1F *hdata =(TH1F*)gPad->GetPrimitive("hdata");

  tree_sig1->Draw("(t1pfmet)>>h1(60,0,900)","weight"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
  TH1F *h1 =(TH1F*)gPad->GetPrimitive("h1");
  tree_sig2->Draw("(t1pfmet)>>h2(60,0,900)","weight"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
  TH1F *h2 =(TH1F*)gPad->GetPrimitive("h2");
  tree_sig3->Draw("(t1pfmet)>>h3(60,0,900)","weight"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
  TH1F *h3 =(TH1F*)gPad->GetPrimitive("h3");
  tree_sig4->Draw("(t1pfmet)>>h4(60,0,900)","weight"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto && Cut_Jets));
  TH1F *h4 =(TH1F*)gPad->GetPrimitive("h4");
//.........这里部分代码省略.........