void binary( string inFile = "allSim.root", int cSpecies = 1, bool cutDedx = false ){

	TCanvas * c = new TCanvas( "c", "c", 800, 800 );
	string outName = "rpBinaryPid.pdf";
	c->Print( (outName+"[").c_str() );

	TFile * f = new TFile( inFile.c_str(), "READ" );

	string rOutName = "rootBinaryPid.root";
	TFile * fOut = new TFile( rOutName.c_str(), "RECREATE" );

	vector<double>effVsP;
	vector<double>pureVsP;

	c->Divide( 2, 2 );
	for ( int i = 0; i < 70; i++ ){

		stringstream sstr;
		sstr << "h_dedx_tof_p3_b" << i;
		TH2* sum = (TH2*)f->Get( sstr.str().c_str() );
		
		sstr.str("");
		sstr << "h_dedx_tof_p0_b" << i;

		TH2* p0 = (TH2*)f->Get( sstr.str().c_str() );
		sstr.str("");
		sstr << "h_dedx_tof_p1_b" << i;
		TH2* p1 = (TH2*)f->Get( sstr.str().c_str() );
		sstr.str("");
		sstr << "h_dedx_tof_p2_b" << i;
		TH2* p2 = (TH2*)f->Get( sstr.str().c_str() );

		pidBinary * pid = new pidBinary( sum, cSpecies, p0, p1, p2 );
		pid->cutDedx( cutDedx );

		c->cd( 3 );
		gPad->SetLogz(1);
		sum->Draw("colz");
		

		c->cd( 4 );
		gPad->SetLogx(1);
		TH1* pX = sum->ProjectionX();
		TH1* pY = sum->ProjectionY();
		pY->SetFillColor( kBlue );
		pY->SetLineColor( kBlue );
		pY->Draw("hbar");


		c->cd(1);
		gPad->SetLogy(1);
		pX->SetFillColor( kBlue );
		pX->SetLineColor( kBlue );
		pX->Draw("h");
		c->cd(2);
		
		sstr.str("");
		sstr << "eff_" << i;
		TH1D* eff = pid->efficiency( sstr.str(), 0.0, 5.0, 0.1 );
		sstr.str("");
		sstr << "pure_" << i;
		TH1D* pure = pid->purity( sstr.str(), 0.0, 5.0, 0.1 );

		gStyle->SetOptStat( 0 );
		eff->SetTitle( "Efficiecy (Blue), Purity (Red)" );
		eff->GetYaxis()->SetRangeUser(0, 1.05);
		eff->SetLineWidth( 2 );
		eff->Draw();
		pure->SetLineColor( kRed );
		pure->SetLineWidth( 2 );
		pure->Draw("same");

		effVsP.push_back( pid->efficiency() );
		pureVsP.push_back( pid->purity( ) );

		c->Print( outName.c_str());
	}

	int nBins = (3.7 - 0.2) / 0.05;
	TH1D * hEffVsP = new TH1D( "hEffVsP", "Efficiency Vs. P; P [GeV]", nBins, 0.2, 3.7 );
	for ( int i = 0; i < effVsP.size(); i++ ){
		hEffVsP->SetBinContent( i, effVsP[ i ] );
	}
	TH1D * hPureVsP = new TH1D( "hPureVsP", "Purity Vs. P; P [GeV]", nBins, 0.2, 3.7 );
	for ( int i = 0; i < pureVsP.size(); i++ ){
		hPureVsP->SetBinContent( i, pureVsP[ i ] );
	}

	c->Divide( 1 );
	c->cd( 0 );
	hEffVsP->GetYaxis()->SetRangeUser( 0.0, 1.05);
	hEffVsP->SetLineWidth( 2 );
	hEffVsP->SetTitle( "Efficiency (Blue), Purity (Red)" );
	hEffVsP->Draw( "");
	hPureVsP->SetLineColor( kRed );
	hPureVsP->SetLineWidth( 2 );
	hPureVsP->Draw( "same" );

	c->Print( outName.c_str());

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

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

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

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

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

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

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

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

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

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

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

      h->Draw();

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

      double peak = -1;

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

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

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

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

void plot_efficiencies( TFile* file, Int_t type = 2, TDirectory* BinDir=0)
{
   // input:   - Input file (result from TMVA),
   //          - type = 1 --> plot efficiency(B) versus eff(S)
   //                 = 2 --> plot rejection (B) versus efficiency (S)
   //                 = 3 --> plot 1/eff(B) versus efficiency (S)

   Bool_t __PLOT_LOGO__  = kTRUE;
   Bool_t __SAVE_IMAGE__ = kTRUE;

   // the coordinates
   Float_t x1 = 0;
   Float_t x2 = 1;
   Float_t y1 = 0;
   Float_t y2 = 0.8;

   // reverse order if "rejection"
   if (type == 2) {
      Float_t z = y1;
      y1 = 1 - y2;
      y2 = 1 - z;
      //      cout << "--- type==2: plot background rejection versus signal efficiency" << endl;
   } else if (type == 3) {
      y1 = 0;
      y2 = -1; // will be set to the max found in the histograms

   } else {
      //  cout << "--- type==1: plot background efficiency versus signal efficiency" << endl;
   }
   // create canvas
   TCanvas* c = new TCanvas( "c", "the canvas", 200, 0, 650, 500 );

   // global style settings
   c->SetGrid();
   c->SetTicks();

   // legend
   Float_t x0L = 0.107,     y0H = 0.899;
   Float_t dxL = 0.457-x0L, dyH = 0.22;
   if (type == 2) {
      x0L = 0.15;
      y0H = 1 - y0H + dyH + 0.07;
   }
   TLegend *legend = new TLegend( x0L, y0H-dyH, x0L+dxL, y0H );
   //legend->SetTextSize( 0.05 );
   legend->SetHeader( "MVA Method:" );
   legend->SetMargin( 0.4 );

   TString xtit = "Signal efficiency";
   TString ytit = "Background efficiency";
   if (type == 2) ytit = "Background rejection";
   if (type == 3) ytit = "1/(Background eff.)";
   TString ftit = ytit + " versus " + xtit;

   TString hNameRef = "effBvsS";
   if (type == 2) hNameRef = "rejBvsS";
   if (type == 3) hNameRef = "invBeffvsSeff";


   if (TString(BinDir->GetName()).Contains("multicut")){
      ftit += "  Bin: ";
      ftit += (BinDir->GetTitle());
   }

   TList xhists;
   TList xmethods;
   UInt_t xnm = TMVAGlob::GetListOfMethods( xmethods );
   TIter xnext(&xmethods);
   // loop over all methods
   TKey *xkey;
   while ((xkey = (TKey*)xnext())) {
      TDirectory * mDir = (TDirectory*)xkey->ReadObj();
      TList titles;
      UInt_t ninst = TMVAGlob::GetListOfTitles(mDir,titles);
      TIter nextTitle(&titles);
      TKey *titkey;
      TDirectory *titDir;
      while ((titkey = TMVAGlob::NextKey(nextTitle,"TDirectory"))) {
         titDir = (TDirectory *)titkey->ReadObj();
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,titDir);
         TIter nextKey( titDir->GetListOfKeys() );
         TKey *hkey2;
         while ((hkey2 = TMVAGlob::NextKey(nextKey,"TH1"))) {
            TH1 *h = (TH1*)hkey2->ReadObj();
            TString hname = h->GetName();
            if (hname.Contains( hNameRef ) && hname.BeginsWith( "MVA_" )) {
               if (type==3 && h->GetMaximum() > y2) y2 = h->GetMaximum()*1.1;
            }
         }
      }
   }


   // draw empty frame
   if(gROOT->FindObject("frame")!=0) gROOT->FindObject("frame")->Delete();
   TH2F* frame = new TH2F( "frame", ftit, 500, x1, x2, 500, y1, y2 );
   frame->GetXaxis()->SetTitle( xtit );
   frame->GetYaxis()->SetTitle( ytit );
   TMVAGlob::SetFrameStyle( frame, 1.0 );
//.........这里部分代码省略.........

RooHistN::RooHistN(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the asymmetry between the specified TH1 objects
  // which may have fixed or variable bin widths, but which must both have
  // the same binning. The asymmetry is calculated as (1-2)/(1+2). Error bars are
  // calculated using Binomial statistics. Prints a warning and rounds
  // any bins with non-integer contents. Use the optional parameter to
  // specify the confidence level in units of sigma to use for
  // calculating error bars. The nominal bin width specifies the
  // default used by addAsymmetryBin(), and is used to set the relative
  // normalization of bins with different widths. If not set, the
  // nominal bin width is calculated as range/nbins.

  initialize();
  // copy the first input histogram's name and title
  SetName(data1.GetName());
  SetTitle(data1.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data1).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
		     data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
  // initialize our contents from the input histogram contents
  Int_t nbin= data1.GetNbinsX();
  if(data2.GetNbinsX() != nbin) {
    coutE(InputArguments) << "RooHistN::RooHistN: histograms have different number of bins" << endl;
    return;
  }
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data1.GetBinCenter(bin);
    if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
      coutW(InputArguments) << "RooHistN::RooHistN: histograms have different centers for bin " << bin << endl;
    }
    Stat_t y1= data1.GetBinContent(bin);
    Stat_t y2= data2.GetBinContent(bin);
    addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac);
  }
  // we do not have a meaningful number of entries
  _entries= -1;
}

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

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

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

      if (log_scale)
      {
        p->SetLogz();
      }
      p->SetGridx(0);
      p->SetGridy(0);

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

      TH1 *h = NULL;

      if (log_scale)
        h = new TH2F(hname,
                     Form(";Sample ID;ADC"), 31, -.5,
                     30.5,
                     //                128+64, 0, 3096);
                     1 << 10, 0, 1 << 14);
      //          else
      //            h = new TH2F(hname,
      //                Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100,
      //                -.050, .5,128,0,2048);

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

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

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

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

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

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

      //          return;
    }
  }

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

void SetStyle(TH1& h, double size, int color, int style, int fillstyle=0, int linestyle=1){
	h.SetMarkerSize(size);
	h.SetMarkerColor(color);
	h.SetLineColor(color);
	h.SetMarkerStyle(style);
	h.SetFillStyle(fillstyle);
	h.SetLineStyle(linestyle);
	h.GetXaxis()->SetTitleFont(42);
	h.GetYaxis()->SetTitleFont(42);
	h.GetXaxis()->SetTitleSize(0.048);
	h.GetYaxis()->SetTitleSize(0.048);
	h.GetXaxis()->CenterTitle();
	h.GetYaxis()->CenterTitle();
}

void SetXRange(TH1 &h, double xmin, double xmax){
	h.GetXaxis()->SetRangeUser(xmin,xmax);
}

/*============================================================================*/
void gaus2peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[8],epar[8],x[2],y[2];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x4);
  y[1] = hist->GetBinContent(bin2);
  TGraph *g = new TGraph(2,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[6]=fpol1->GetParameter(0);
  par[7]=fpol1->GetParameter(1);

  //--**-- Two Gaussian Peak estimation without background --**--//
  fgaus1 = new TF1("m1","gaus",x1,x2);
  fgaus2 = new TF1("m2","gaus",x3,x4);
  hist->Fit(fgaus1,"R+QN");
  hist->Fit(fgaus2,"R+QN");
  fgaus1->GetParameters(&par[0]);
  fgaus2->GetParameters(&par[3]);

  //--**-- Final Peak Fit with Background --**--//
  func = new TF1("m","gaus(0)+gaus(3)+pol1(6)",x1,x4);
  func->SetParameters(par);

  
  hist->Fit(func,"R+QN");
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  epar[3]=func->GetParError(3);
  epar[4]=func->GetParError(4);
  epar[5]=func->GetParError(5);
 
  Double_t fwhm1 = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm1 = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N10 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r10 = epar[0]/par[0];
  Double_t r12 = epar[2]/par[2];
  Double_t eN10= N10*TMath::Sqrt(r10*r10+r12*r12);

  Double_t fwhm2 = par[5]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm2 = epar[5]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N20 = par[3]*(TMath::Sqrt(TMath::TwoPi())*par[5]);
  Double_t r20 = epar[3]/par[3];
  Double_t r22 = epar[5]/par[5];
  Double_t eN20= N20*TMath::Sqrt(r20*r20+r22*r22);

  //printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
  //        par[1],epar[1],fwhm1,efwhm1,N10,eN10);
  //printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
  //        par[4],epar[4],fwhm2,efwhm2,N20,eN20);
  printf("%11.4f %11.4f %11.0f %11.0f\n",
          par[1],epar[1],N10,eN10);
  printf("%11.4f %11.4f %11.0f %11.0f\n",
          par[4],epar[4],N20,eN20);
}

void CombinePreScale()
{

    double jetptbin[] = {27, 33, 39, 47, 55, 64,74, 84, 97, 114, 133, 153, 174, 196, 220, 245, 272, 300, 429, 692, 1000};   
    int nbins = sizeof(jetptbin)/sizeof(double)-1; 

    TString kDir="/scratch/xuq7/RpA/TreeAna";
    TString algo ="akPu3PF"; //"AkPu3PF" ;
    TString residual = "NoResidual"; //"NoResidual" ;
    TString coll = "PPbNoGplus" ; // or "PbP" ;
    TString class = "HFsumEta4Bin1" ; // "" for inclusive ;
    TString effTab = "TrkEffHIN12017v5TrkCorr2DCut" ;  //"HistIterTrkCorrtestFilterCut"; //  "Trk" ; "HistIterTrkCorrtest" ;
    const int Nfile = 5 ;
    TFile * f ;
    TString name[]={"jetptEta","jetptphi","jetEtaphi","jetptchMax","jetptchSum","jetptneuMax","jetptneuSum","jetptphoMax","jetptphoSum","jetptchMaxpt","jetptchSumpt","jetptneuMaxpt","jetptneuSumpt","jetptphoMaxpt","jetptphoSumpt","jetptSumSumpt","jetptSumSumrawpt","jetptneuMaxr","jetptchN","jetptneuN","jetptphoN"};
    TString etaname[]={"jetpt","rawptJES","jetptchMaxpt","jetptneuMaxr","jetptSumSumpt"};
    int nname=sizeof(name)/sizeof(TString);
    int netaname=sizeof(etaname)/sizeof(TString);
    f = TFile::Open(Form("%s/AllTrigLumiDATAPPb%sDiJetMass.root",kDir.Data(),algo.Data()),"readonly");
    TString outname ;
    outname ="Datacombined.root";
    TFile * outf = new TFile(Form("%s/%s", kDir.Data(), outname.Data()), "RECREATE");
//    f =  TFile::Open("AllTrigLumiDATAPPbakPu3PFDiJetMass.root", "readonly");
for(int m=0;m<nname;m++){
    TH1 *hCombined;
    TH1 *hTrkPt[5];
  
                hTrkPt[0] = (TH1*)f->Get(Form("Jet20%s_0-100%%",name[m].Data()));
                hTrkPt[1] = (TH1*)f->Get(Form("Jet40%s_0-100%%",name[m].Data()));
                hTrkPt[2] = (TH1*)f->Get(Form("Jet60%s_0-100%%",name[m].Data()));
                hTrkPt[3] = (TH1*)f->Get(Form("Jet80%s_0-100%%",name[m].Data()));
                hTrkPt[4] = (TH1*)f->Get(Form("Jet100%s_0-100%%",name[m].Data()));
    for( int i = 0; i<5; i++)
    {
        if(i==0){ 
	hCombined = (TH1*)hTrkPt[i]->Clone(Form("%sCombinedSpectra",name[m].Data()));
	hCombined->SetTitle("");
	}
        else hCombined->Add(hTrkPt[i]);
    }
    hCombined->Write();
}
     double etaMin[8] = {-1.,-2.2,-1.2,-0.7,-0.3,0.3,0.7,1.2};
    double etaMax[8] = { 1.,-1.2,-0.7,-0.3, 0.3,0.7,1.2,2.2};
TH1* hTrk[8][5];
TH1* hCombinedEtaBin[8];
for(int m=0;m<netaname;m++){
	    for(int ieta=0; ieta<8; ieta++){
            if(ieta==0){
                hTrk[ieta][0] = (TH1*)f->Get(Form("Jet20%s_0-100%%",etaname[m].Data()));
                hTrk[ieta][1] = (TH1*)f->Get(Form("Jet40%s_0-100%%",etaname[m].Data()));
                hTrk[ieta][2] = (TH1*)f->Get(Form("Jet60%s_0-100%%",etaname[m].Data()));
                hTrk[ieta][3] = (TH1*)f->Get(Form("Jet80%s_0-100%%",etaname[m].Data()));
                hTrk[ieta][4] = (TH1*)f->Get(Form("Jet100%s_0-100%%",etaname[m].Data()));
            }
            else {
                hTrk[ieta][0] = (TH1*)f->Get(Form("Jet20%sEtaBin%.f_%.f_Cen0-100%%",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
                hTrk[ieta][1] = (TH1*)f->Get(Form("Jet40%sEtaBin%.f_%.f_Cen0-100%%",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
                hTrk[ieta][2] = (TH1*)f->Get(Form("Jet60%sEtaBin%.f_%.f_Cen0-100%%",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
                hTrk[ieta][3] = (TH1*)f->Get(Form("Jet80%sEtaBin%.f_%.f_Cen0-100%%",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
                hTrk[ieta][4] = (TH1*)f->Get(Form("Jet100%sEtaBin%.f_%.f_Cen0-100%%",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
                }

 for(int i = 0; i<5; i++)
    {
        if(i==0) {
	hCombinedEtaBin[ieta] = (TH1*)hTrk[ieta][i]->Clone(Form("%sCombinedSpectraInEtaBin%.f_%.f",etaname[m].Data(), etaMin[ieta]*10, etaMax[ieta]*10));
	hCombinedEtaBin[ieta]->SetTitle("");
	}
        else hCombinedEtaBin[ieta]->Add(hTrk[ieta][i]);
    }
    hCombinedEtaBin[ieta]->Write();
}
}

    outf->Close();//write into a root file
    
}

void plot_fitter_validation()
{
  TFile* fin = new TFile("fitter_validation_cafana.root");

  const int dcpCols[] = {kBlack, kRed, kGreen+2, kBlue};
  const std::string dcpStrs[] = {"0", "#pi/2", "#pi", "3#pi/2"};

  for(std::string hcStr: {"fhc", "rhc"}){
    const std::string HCStr = (hcStr == "fhc") ? "FHC" : "RHC";

    for(std::string chanStr: {"numu", "nue"}){
      const std::string CHANStr = (chanStr == "numu") ? "#nu_{#mu}" : "#nu_{e}";

      for(std::string hieStr: {"nh", "ih"}){
        const std::string HIEStr = (hieStr == "nh") ? "NH" : "IH";

        new TCanvas;

        for(int deltaIdx2 = 0; deltaIdx2 < 4; ++deltaIdx2){
          // For neutrinos 3pi/2 is the tallest histogram, draw it first, for
          // antineutrinos we need pi/2 first.
          const int deltaIdx = (hcStr == "fhc") ? 3-deltaIdx2 : (deltaIdx2+1)%4;

          const std::string dcpStr = TString::Format("%gpi", deltaIdx/2.).Data();

          TH1* h = (TH1*)fin->Get((chanStr+"_"+hcStr+"_"+hieStr+"_"+dcpStr).c_str());
          h->SetLineColor(dcpCols[deltaIdx]);
          h->Draw("same");

          h->SetTitle(("5 yrs "+HCStr+" "+CHANStr+" "+HIEStr).c_str());
        } // end for deltaIdx

        TLegend* leg = new TLegend(.6, .6, .85, .85);
        leg->SetFillStyle(0);
        for(int deltaIdx = 0; deltaIdx < 4; ++deltaIdx){
          TH1* dummy = new TH1F("", "", 1, 0, 1);
          dummy->SetLineColor(dcpCols[deltaIdx]);
          leg->AddEntry(dummy, ("#delta_{CP}="+dcpStrs[deltaIdx]).c_str(), "l");
        }
        leg->Draw("same");

        gPad->Print((hcStr+"_"+chanStr+"_"+hieStr+".pdf").c_str());
      } // end for hieStr
    } // end for chanStr
  } // end for hcStr


  TGraph* gNH = (TGraph*)fin->Get("sens_nh");
  TGraph* gIH = (TGraph*)fin->Get("sens_ih");
  TGraph* gNHOscErr = (TGraph*)fin->Get("sens_nh_oscerr");
  TGraph* gIHOscErr = (TGraph*)fin->Get("sens_ih_oscerr");

  TGraph* gNHFlux[10];
  TGraph* gIHFlux[10];
  for(int i = 0; i < 10; ++i){
    gNHFlux[i] = (TGraph*)fin->Get(TString::Format("sens_nh_flux%d", i).Data());
    gIHFlux[i] = (TGraph*)fin->Get(TString::Format("sens_ih_flux%d", i).Data());
  }

  TGraph* gNHXSec[10];
  TGraph* gIHXSec[10];
  for(int i = 0; i < 10; ++i){
    gNHXSec[i] = (TGraph*)fin->Get(TString::Format("sens_nh_xsec%d", i).Data());
    gIHXSec[i] = (TGraph*)fin->Get(TString::Format("sens_ih_xsec%d", i).Data());
  }

  TH2* axes = new TH2F("", ";#delta_{CP} / #pi;#sigma = #sqrt{#Delta#chi^{2}}", 100, 0, 2, 100, 0, 8);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->SetTitleOffset(.75);
  axes->GetYaxis()->CenterTitle();
  axes->Draw();

  gNH->Draw("l same");
  gIH->Draw("l same");
  gNHOscErr->Draw("l same");
  gIHOscErr->Draw("l same");

  TLegend* leg = new TLegend(.4, .65, .6, .875);
  leg->SetFillStyle(0);
  leg->AddEntry(gNH, "NH", "l");
  leg->AddEntry(gIH, "IH", "l");
  leg->AddEntry(gNHOscErr, "NH osc err", "l");
  leg->AddEntry(gIHOscErr, "IH osc err", "l");
  leg->Draw();

  gPad->Print("mcd.pdf");


  new TCanvas;
  axes->Draw();
  for(int i = 0; i < 10; ++i){
    gNHFlux[i]->Draw("l same");
    gIHFlux[i]->Draw("l same");
  }
  gNH->Draw("l same");
  gIH->Draw("l same");

  leg = new TLegend(.4, .65, .6, .875);
  leg->SetFillStyle(0);
  leg->AddEntry(gNH, "NH", "l");
//.........这里部分代码省略.........

/*============================================================================*/
void gaus1peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[5],epar[5],x[4],y[4];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  setcanvas(1);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x2;
  x[2] = x3;
  x[3] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x2);
  y[1] = hist->GetBinContent(bin2);
  Int_t bin3 = hist->FindBin(x3);
  y[2] = hist->GetBinContent(bin3);
  Int_t bin4 = hist->FindBin(x4);
  y[3] = hist->GetBinContent(bin4);
  TGraph *g = new TGraph(4,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[3]=fpol1->GetParameter(0);
  par[4]=fpol1->GetParameter(1);

  //--**-- Gaussian Peak estimation without background --**--//
  TF1 *fgaus = new TF1("GAUS","gaus",x2,x3);
  hist->Fit(fgaus,"RQN");
  fgaus->GetParameters(&par[0]);

  //--**-- Final Peak Fit with Background --**--//
  TF1 *func = new TF1("FGAUS","gaus(0)+pol1(3)",x1,x4);
  func->SetParameters(par);
  hist->Fit(func,"R+QN");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  Double_t fwhm = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N0 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r0 = epar[0]/par[0];
  Double_t r2 = epar[2]/par[2];
  Double_t eN0= N0*TMath::Sqrt(r0*r0+r2*r2);
  printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
          par[1],epar[1],fwhm,efwhm,N0,eN0);
  //printf("%11.4f %11.4f %11.0f %11.0f\n",
  //        par[1],epar[1],N0,eN0);
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
}

void gluinoMass(double lumi=-1., double maxInstLumi=-1.) {
	
  if (lumi<0)
    lumi=877.;
  if (maxInstLumi<0)
    maxInstLumi=1300.;
  LimitPlots plots(lumi);
  
  //mchamp index 0 is used, corresponds to 0th mass point = 100 GeV
  plots.calculateCrossSections(7,4,0,39,9);
	
  // expected limit (1 and 2 sigma bands)
  TGraph* g_exp = plots.getExpMassLimitGluino();
  TGraphAsymmErrors* g_exp1 = plots.getExpMassLimitGluino1Sig();
  TGraphAsymmErrors* g_exp2 = plots.getExpMassLimitGluino2Sig();
  
  // three points on counting expt curve
  TGraph* g_gluino = plots.getMassLimitGluino();
  TGraph* g_stop = plots.getMassLimitStop();
  
  // one point from lifetime fit
  TGraph* g_tp = plots.getMassLimitGluinoTP();
  
  // theory prediction
  TGraph* g_thGluino = plots.getGluinoTheory();
  TGraph* g_thStop = plots.getStopTheory();
  
  TCanvas* canvas = new TCanvas("canvas");
  
  //canvas->SetGrid();
  canvas->SetLogy();
  
  TH1 * h;
  h = canvas->DrawFrame(300., .02, 1000., 1e2);
  //h->SetTitle("Beamgap Expt;m_{#tilde{g}} [GeV/c^{2}]; Stopped HSCP Cross Section #times BR [pb]");
  h->SetTitle("Beamgap Expt;m_{#tilde{g}} [GeV/c^{2}]; #sigma(pp #rightarrow #tilde{g}#tilde{g}) #times BR(#tilde{g} #rightarrow g#tilde{#chi}^{0}) [pb]");
  
  // not covered region
  TBox* nc = new TBox(100., .1, 150., 5e2);
  nc->SetFillStyle(3354);
  nc->SetFillColor(kRed-4);
  //nc->Draw();
  
  // details
  TPaveText* blurb = new TPaveText(300., 2, 550., 1e2);
  blurb->AddText("CMS Preliminary 2012");

  std::stringstream label;
  label<<"#int L dt = "<<lumi<<" fb^{-1}";
  blurb->AddText(label.str().c_str());
  label.str("");
  double peakInstLumi=maxInstLumi;
  int exponent=30;
  while (peakInstLumi>10) {
    peakInstLumi/=10.;
    ++exponent;
  }
  label<<"L^{max}_{inst} = "<<peakInstLumi<<" x 10^{"<<exponent<<"} cm^{-2}s^{-1}";
  blurb->AddText(label.str().c_str());
  label.str("");
  label << "#sqrt{s} = " << ENERGY << " TeV";
  blurb->AddText(label.str().c_str());
  blurb->AddText("m_{#tilde{g}} - m_{#tilde{#chi}^{0}} = 100 GeV/c^{2}");
  //blurb->AddText("m_{#tilde{t}} - m_{#tilde{#chi}^{0}} = 200 GeV/c^{2}");	
  blurb->SetTextFont(42);
  blurb->SetBorderSize(0);
  blurb->SetFillColor(0);
  blurb->SetShadowColor(0);
  blurb->SetTextAlign(12);
  blurb->SetTextSize(0.032);

  
  // legend
  TBox *legbg = new TBox(600., 2., 900., 1e2);
  legbg->Draw();
  TLegend *leg = new TLegend(600., 2., 900., 1e2,"95% C.L. Limits","");
  leg->SetTextSize(0.028);
  leg->SetBorderSize(0);
  leg->SetTextFont(42);
  leg->SetFillColor(0);
  leg->AddEntry(g_exp, "Expected: 10 #mus - 1000 s Counting Exp. ", "l");
  leg->AddEntry(g_exp1, "Expected #pm1#sigma: 10 #mus - 1000 s Counting Exp. ", "f");
  leg->AddEntry(g_exp2, "Expected #pm2#sigma: 10 #mus - 1000 s Counting Exp. ", "f");
  //	  leg->AddEntry(graph3, "Obs.: 10^{6} s Counting Exp.", "l");
  leg->AddEntry(g_gluino, "Obs.: 10 #mus - 1000 s Counting Exp. ", "l");
  leg->AddEntry(g_tp, "Obs.: 10 #mus Timing Profile ", "l");
  //leg->AddEntry(g_stop, "Obs.: 10 #mus - 1000 s Counting Exp. (#tilde{t})", "l");
  //leg->AddEntry(graph_em, "Obs.: 10 #mus - 1000 s Counting Exp. (EM only)", "l");
  //  leg->AddEntry(graph1, "Obs.: 570 ns Counting Exp.", "l");
  leg->Draw();
  
  
  
  // 2 sigma expected band
  g_exp2->SetLineColor(0);
  g_exp2->SetLineStyle(0);
  g_exp2->SetLineWidth(0);
  g_exp2->SetFillColor(5);
  g_exp2->SetFillStyle(1001);
  g_exp2->Draw("3");
//.........这里部分代码省略.........

void makeSplitQCDhist_PythiaBinned(vector<string> folders, const string histname, 
											const string htrange, const string htbinlabel,
											const hist_t histinfo)
{
	const float scaleTo = fDATA_LUMI; // pb

	TLegend *leg  = new TLegend(0.6,0.65,0.9,0.9);
	leg->SetTextFont(42);
	vector<TH1*> hists;
	new TCanvas();
	gPad->SetLogy();
	gPad->SetTickx();
	gPad->SetTicky();

		stringstream title;
		title << htrange << " [" << histinfo.name << "];" << histinfo.title;
		cout << title.str() << endl;
		TH1* Hist = GetHist(histname);

		Hist->SetTitle(title.str().c_str());
		Hist->SetMarkerStyle(20);
		Hist->SetStats(0);
		Hist->Rebin(histinfo.rebin);
		Hist->Draw("P");


/*
	for (unsigned i = 0; i < folders.size(); ++i)
	{	
		string njet("");
		if (i==0) njet += "[2-3]";
		else if (i==1) njet += "[4-5]";
		else if (i==2) njet += "[6-7]";
		else if (i==3) njet += "#geq 8";
*/		/*if (i==0) njet += "3";
		else if (i==1) njet += "4";
		else if (i==2) njet += "5";
		else if (i==3) njet += "6";
		else if (i==4) njet += "7";
		else if (i==5) njet += "8";*/

/*		stringstream title, histName;
		title << htrange << ";" << histinfo.title;
		cout << title.str() << endl;
		histName << folders.at(i) << "/" << histname;;
		hists.push_back(GetHist(histName.str()));
		hists.at(i)->SetTitle(title.str().c_str());
		hists.at(i)->SetMarkerStyle(20+i);
		hists.at(i)->SetMarkerColor(1+i*2);
		hists.at(i)->SetStats(0);
		hists.at(i)->Rebin(histinfo.rebin);
		if (histinfo.normalizeByBinWidth) NormByBinWidth(hists.at(i));

		stringstream legname;
		legname << "Njets " << njet;
		leg->AddEntry(hists.at(i), legname.str().c_str());
		if (i==0) hists.at(i)->Draw("P");
		else hists.at(i)->Draw("same P");
	}
	leg->Draw();
*/
}

void BasePlotter::writeStacked(string filename, const HistogramContainer& histContainer, string extension) const
{
    // Check histContainer consistency
    if( ! histContainer.check() ) 
        throw 1;

    if( extension[0] == '.' )
        extension.erase(0,1);
    
    if( filename.find(".root") == string::npos )
    filename.append( ".root" );
    
    system(("mkdir -p " + configContainer.outputDir).c_str()); 
    TFile* f = new TFile((configContainer.outputDir + filename).c_str(), "UPDATE");
    
    unsigned int nSamples = histContainer.reducedNames.size();
    TCanvas *c = new TCanvas(("Canv_stacked_" + histContainer.containerName).c_str(), "", 600, 600 + (120 * configContainer.plotRatio));
    TPad *histPad=nullptr, *ratioPad=nullptr;
    if( configContainer.plotRatio || configContainer.plotSignificance )
    {
        histPad = new TPad("pad1", "pad1", 0, 0.25, 1, 1);
        histPad->Draw();     
        ratioPad = new TPad("pad2", "pad2", 0, 0, 1, 0.3); 
        ratioPad->Draw();
        histPad->cd();
    }    
    
    TLegend* leg = nullptr;
    vector<TLatex*> latexVector;
    setCanvasAttributes( nSamples, leg, latexVector);
    
    vector<TH1F*> hStack;
    TH1* hData = nullptr, *hSignal = nullptr;
    // Add backgrounds to stack
    for( int iSample = nSamples-1; iSample > -1; --iSample )
    {
        if( histContainer.sampleType[iSample] == SampleType::MC || histContainer.sampleType[iSample] == SampleType::FAKELEPTON || histContainer.sampleType[iSample] == SampleType::MCFAKELEPTON )
        {
            TH1F *temp = (TH1F*) histContainer.histograms[iSample]->Clone(("stack_hist_"+histContainer.reducedNames[iSample]).c_str());
            if( hStack.size() > 0 ) temp->Add(hStack[hStack.size()-1]);
            temp->SetLineColor(histContainer.color[iSample]);
            temp->SetFillColor(histContainer.color[iSample]);
            temp->SetFillStyle(1001);
            hStack.push_back(temp);
            string legendEntry = histContainer.reducedNames[iSample];
            replace(legendEntry.begin(), legendEntry.end(), '_', ' ');
            leg->AddEntry(temp,legendEntry.c_str(),"f");
        }
    }
    
    // Add signal to stack
    for( unsigned int iSample = 0; iSample < nSamples; ++iSample )
    {
        if( histContainer.sampleType[iSample] == SampleType::SIGNAL )
        {
            hSignal = histContainer.histograms[iSample];
            hSignal->SetLineColor(histContainer.color[iSample]);
            string legendEntry = histContainer.reducedNames[iSample];
            replace(legendEntry.begin(), legendEntry.end(), '_', ' ');
            if( configContainer.signalStacked )
            {
                TH1F *temp = (TH1F*) histContainer.histograms[iSample]->Clone(("stack_hist_"+histContainer.reducedNames[iSample]).c_str());
                if( hStack.size() > 0 ) temp->Add(hStack[hStack.size()-1]);
                temp->SetFillColor(histContainer.color[iSample]);
                temp->SetFillStyle(1001);
                hStack.push_back(temp);
                leg->AddEntry(temp,legendEntry.c_str(),"f");
            }
            else
            {
                hSignal->SetLineWidth(2);
                leg->AddEntry(hSignal,legendEntry.c_str(),"l");
            }
        }
    }
    
    // Add data
    for( unsigned int iSample = 0; iSample < nSamples; ++iSample )
    {
        if( histContainer.sampleType[iSample] == SampleType::DATA && configContainer.unblind )
        {
            histContainer.histograms[iSample]->SetLineColor(histContainer.color[iSample]);
            histContainer.histograms[iSample]->SetMarkerColor(histContainer.color[iSample]);
            histContainer.histograms[iSample]->SetMarkerSize(1);
            histContainer.histograms[iSample]->SetMarkerStyle(20);
            histContainer.histograms[iSample]->SetLineWidth(2);
            string legendEntry = histContainer.reducedNames[iSample];
            replace(legendEntry.begin(), legendEntry.end(), '_', ' ');
            leg->AddEntry(histContainer.histograms[iSample],legendEntry.c_str(),"lp");
            hData = histContainer.histograms[iSample];
        }
    }

    // Set y-range
    float hMax = 0.;
    float hMin = 0.;
    bool maxMakesSense = true;
    if( configContainer.drawNormalized )
        hMax = 1;
    else
//.........这里部分代码省略.........

// input: - Input file (result from TMVA)
//        - use of TMVA plotting TStyle
void deviations( TString fin = "TMVAReg.root", 
                 HistType htype = MVAType, Bool_t showTarget, Bool_t useTMVAStyle = kTRUE )
{
   // set style and remove existing canvas'
   TMVAGlob::Initialize( useTMVAStyle );
   gStyle->SetNumberContours(999);

   // switches
   const Bool_t Save_Images     = kTRUE;

   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );  

   // define Canvas layout here!
   Int_t xPad = 1; // no of plots in x
   Int_t yPad = 1; // no of plots in y
   Int_t noPad = xPad * yPad ; 
   const Int_t width = 650;   // size of canvas

   // this defines how many canvases we need
   TCanvas* c[100];

   // counter variables
   Int_t countCanvas = 0;

   // search for the right histograms in full list of keys
   //    TList* methods = new TMap();

   TIter next(file->GetListOfKeys());
   TKey *key(0);   
   while ((key = (TKey*)next())) {

      if (!TString(key->GetName()).BeginsWith("Method_")) continue;
      if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) continue;

      TString methodName;
      TMVAGlob::GetMethodName(methodName,key);
      cout << "--- Plotting deviation for method: " << methodName << endl;

      TObjString *mN  = new TObjString( methodName );
      TDirectory* mDir = (TDirectory*)key->ReadObj();

      TList* jobNames = new TList();

      TIter keyIt(mDir->GetListOfKeys());
      TKey *titkey;
      while ((titkey = (TKey*)keyIt())) {

         if (!gROOT->GetClass(titkey->GetClassName())->InheritsFrom("TDirectory")) continue;

         TDirectory *titDir = (TDirectory *)titkey->ReadObj();

         TObjString *jN = new TObjString( titDir->GetName() );
         if (!jobNames->Contains( jN )) jobNames->Add( jN );
         else delete jN;
	    
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,titDir);

         TString hname = "MVA_" + methodTitle;
         TIter   dirKeyIt( titDir->GetListOfKeys() );
         TKey*   dirKey;

         Int_t countPlots = 0;
         while ((dirKey = (TKey*)dirKeyIt())){
            if (dirKey->ReadObj()->InheritsFrom("TH2F")) {
               TString s(dirKey->ReadObj()->GetName());
               if (s.Contains("_reg_") && 
                   ( (showTarget && s.Contains("_tgt")) || (!showTarget && !s.Contains("_tgt")) ) && 
                   s.Contains( (htype == CompareType ? "train" : "test" ))) {
                  c[countCanvas] = new TCanvas( Form("canvas%d", countCanvas+1), 
                                                Form( "Regression output deviation versus %s for method: %s",
                                                      (showTarget ? "target" : "input variables"), methodName.Data() ),
                                                countCanvas*50+100, (countCanvas+1)*20, width, (Int_t)width*0.72 ); 
                  c[countCanvas]->SetRightMargin(0.10); // leave space for border
                  TH1* h = (TH1*)dirKey->ReadObj();
                  h->SetTitle( Form("Output deviation for method: %s (%s sample)", 
                                    hname.Data(), (htype == CompareType ? &