C++ TH1::GetBinWidth方法代码示例

TH1 *
YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth)
{
  /* find lowest edge in histo */
  Int_t binlo;
  Double_t lo;
  for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
    if (h->GetBinContent(ibin) != 0.) {
      binlo = ibin;
      lo = h->GetBinLowEdge(ibin);
      break;
    }
  }
  
  Int_t nbins = (lo - min) / binwidth;
  if(nbins<1)
	return 0x0;		
  TH1 *hlo = new TH1F("hlo", "", nbins, min, lo);
  
  /* integrate function in histogram bins */
  Double_t cont, err, width;
  for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
    width = hlo->GetBinWidth(ibin + 1);
    cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
    err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
    hlo->SetBinContent(ibin + 1, cont / width);
    hlo->SetBinError(ibin + 1, err / width);
  }

  return hlo;
}

TH1 *
YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth)
{
  /* find highest edge in histo */
  Int_t binhi;
  Double_t hi;
  for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
    if (h->GetBinContent(ibin) != 0.) {
      binhi = ibin + 1;
      hi = h->GetBinLowEdge(ibin + 1);
      break;
    }
  }
  if(max<hi) {
	Printf("Warning! You should probably set a higher max value (Max = %f, hi = %f)", max, hi);
    return 0x0;
  }
  Int_t nbins = (max - hi) / binwidth;
 if(nbins<1)
	return 0x0;	 
  TH1 *hhi = new TH1F("hhi", "", nbins, hi, max);
  
  /* integrate function in histogram bins */
  Double_t cont, err, width;
  for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
    width = hhi->GetBinWidth(ibin + 1);
    cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
    err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
    hhi->SetBinContent(ibin + 1, cont / width);
    hhi->SetBinError(ibin + 1, err / width);
  }

  return hhi;
}

Double_t fitgp0( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){
    cout << hs << " does not exist\n";
    return 0;
  }

  h->SetMarkerStyle(21);
  h->SetMarkerSize(0.8);
  h->SetStats(1);
  gStyle->SetOptFit(101);

  gROOT->ForceStyle();

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1 = h->GetBinCenter(1);
  double x9 = h->GetBinCenter(nb);

  // create a TF1 with the range from x1 to x9 and 4 parameters
  TF1 *gp0Fcn = new TF1( "gp0Fcn", gp0Fit, x1, x9, 4 );

  gp0Fcn->SetParName( 0, "mean" );
  gp0Fcn->SetParName( 1, "sigma" );
  gp0Fcn->SetParName( 2, "area" );
  gp0Fcn->SetParName( 3, "BG" );

  gp0Fcn->SetNpx(500);
  gp0Fcn->SetLineWidth(4);
  gp0Fcn->SetLineColor(kMagenta);
  gp0Fcn->SetLineColor(kGreen);

  // set start values for some parameters:
  gp0Fcn->SetParameter( 0, xmax ); // peak position
  gp0Fcn->SetParameter( 1, 4*dx ); // width
  gp0Fcn->SetParameter( 2, nn ); // N
  gp0Fcn->SetParameter( 3, bg );

  // N: not drawing
  // Q: quiet
  // R: use specified range
  h->Fit( "gp0Fcn", "NQR", "ep" );

  return gp0Fcn->GetParameter(1);

}

RooHistN::RooHistN(const TH1 &data, Double_t nominalBinWidth, Double_t nSigma, RooAbsData::ErrorType etype, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the contents of the specified TH1 object
  // which may have fixed or variable bin widths. Error bars are
  // calculated using Poisson 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 addBin(), 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 input histogram's name and title
  SetName(data.GetName());
  SetTitle(data.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  // TH1::GetYaxis() is not const (why!?)
  setYAxisLabel(const_cast<TH1&>(data).GetYaxis()->GetTitle());
  
  // initialize our contents from the input histogram's contents
  Int_t nbin= data.GetNbinsX();
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data.GetBinCenter(bin);
    Stat_t y= data.GetBinContent(bin);
    Stat_t dy = data.GetBinError(bin) ;
    if (etype==RooAbsData::Poisson) {
      addBin(x,roundBin(y),data.GetBinWidth(bin),xErrorFrac);
    } else {
      addBinWithError(x,y,dy,dy,data.GetBinWidth(bin),xErrorFrac);
    }
  }
  // add over/underflow bins to our event count
  _entries+= data.GetBinContent(0) + data.GetBinContent(nbin+1);
}

Double_t fitfulllang( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){
    cout << hs << " does not exist\n";
    return 0;
  }

  double aa = h->GetEntries();//normalization

  // find peak:
  int ipk = h->GetMaximumBin();
  double xpk = h->GetBinCenter(ipk);
  double sm = xpk / 9; // sigma
  double ns = sm; // noise

  // fit range:
  int ib0 = ipk/2;
  int ib9 = h->GetNbinsX() - 1;

  double x0 = h->GetBinLowEdge(ib0);
  double x9 = h->GetBinLowEdge(ib9) + h->GetBinWidth(ib9);

  // create a TF1 with the range from x0 to x9 and 4 parameters
  TF1 *fitFcn = new TF1( "fitFcn", fitLandauGauss, x0, x9, 4 );

  fitFcn->SetParName( 0, "peak" );
  fitFcn->SetParName( 1, "sigma" );
  fitFcn->SetParName( 2, "area" );
  fitFcn->SetParName( 3, "smear" );

  fitFcn->SetNpx(500);
  fitFcn->SetLineWidth(4);
  fitFcn->SetLineColor(kMagenta);

  // set start values:
  fitFcn->SetParameter( 0, xpk ); // peak position, defined above
  fitFcn->SetParameter( 1, sm ); // width
  fitFcn->SetParameter( 2, aa ); // area
  fitFcn->SetParameter( 3, ns ); // noise

  h->Fit("fitFcn", "NQR", "ep" );// R = range from fitFcn
  return fitFcn->GetParameter(0);
}

double fittp0sigma( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){ cout << hs << " does not exist\n"; return 0; }

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1 = h->GetBinCenter(1);
  double x9 = h->GetBinCenter(nb);

  // create a TF1 with the range from x1 to x9 and 5 parameters
  TF1 *tp0Fcn = new TF1( "tp0Fcn", tp0Fit, x1, x9, 5 );

  tp0Fcn->SetParName( 0, "mean" );
  tp0Fcn->SetParName( 1, "sigma" );
  tp0Fcn->SetParName( 2, "nu" );
  tp0Fcn->SetParName( 3, "area" );
  tp0Fcn->SetParName( 4, "BG" );
   
  // set start values for some parameters:
  tp0Fcn->SetParameter( 0, xmax ); // peak position
  tp0Fcn->SetParameter( 1, 4*dx ); // width
  tp0Fcn->SetParameter( 2, 2.2 ); // nu
  tp0Fcn->SetParameter( 3, nn ); // N
  tp0Fcn->SetParameter( 4, bg );
    
  h->Fit( "tp0Fcn", "Q R", "ep" );
  // h->Fit("tp0Fcn","V+","ep");
  TF1 *fit = h->GetFunction("tp0Fcn");
  return fit->GetParameter(1);
}

//.........这里部分代码省略.........
//             sig->SetTitle( tit.ReplaceAll("Input variable", "Regression target" ) );
//          }
//       }
      sig->SetTitle( "" );            
      

      // finally plot and overlay
      Float_t sc = 1.1;
      if (countPad == 1) sc = 1.3;
      sig->SetMaximum( TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*sc );
      sig->Draw( "hist" );
      cPad->SetLeftMargin( 0.17 );

      sig->GetYaxis()->SetTitleOffset( 1.50 );
      if (!isRegression) {
         bgd->Draw("histsame");
         TString ytit = TString("(1/N) ") + sig->GetYaxis()->GetTitle();
         ytit = TString("Fraction of Events");
         sig->GetYaxis()->SetTitle( ytit ); // histograms are normalised
      }

      if (countPad == 1) sig->GetXaxis()->SetTitle("Leading Lepton p_{T} [GeV/c]");
      if (countPad == 2) sig->GetXaxis()->SetTitle("Trailing Lepton p_{T} [GeV/c]");
      if (countPad == 3) sig->GetXaxis()->SetTitle("#Delta#phi(l,l)");
      if (countPad == 4) sig->GetXaxis()->SetTitle("#Delta R(l,l)");
      if (countPad == 5) sig->GetXaxis()->SetTitle("Dilepton Mass [GeV/c^{2}]");
      if (countPad == 6) sig->GetXaxis()->SetTitle("Dilepton Flavor Final State");
      if (countPad == 7) sig->GetXaxis()->SetTitle("M_{T} (Higgs) [GeV/c^{2}]");
      if (countPad == 8) sig->GetXaxis()->SetTitle("#Delta#phi(Dilepton System, MET)");
      if (countPad == 9) sig->GetXaxis()->SetTitle("#Delta#phi(Dilepton System, Jet)");


      // Draw legend
//       if (countPad == 1 && !isRegression) {
         TLegend *legend= new TLegend( cPad->GetLeftMargin(), 
                                       1-cPad->GetTopMargin()-.15, 
                                       cPad->GetLeftMargin()+.4, 
                                       1-cPad->GetTopMargin() );

         if(countPad == 1 || countPad == 2 ||countPad == 3 ||countPad == 4 ||countPad == 5 ||countPad == 7  ) {
           legend= new TLegend( 0.50, 
                                1-cPad->GetTopMargin()-.15, 
                                0.90, 
                                1-cPad->GetTopMargin() );
         }

         legend->SetFillStyle(0);
         legend->AddEntry(sig,"Signal","F");
         legend->AddEntry(bgd,"Background","F");
         legend->SetBorderSize(0);
         legend->SetMargin( 0.3 );
         legend->SetTextSize( 0.03 );
         legend->Draw("same");
//       } 

      // redraw axes
      sig->Draw("sameaxis");

      // text for overflows
      Int_t    nbin = sig->GetNbinsX();
      Double_t dxu  = sig->GetBinWidth(0);
      Double_t dxo  = sig->GetBinWidth(nbin+1);
      TString uoflow = "";
      if (isRegression) {
         uoflow = Form( "U/O-flow: %.1f%% / %.1f%%", 
                        sig->GetBinContent(0)*dxu*100, sig->GetBinContent(nbin+1)*dxo*100 );
      }
      else {
         uoflow = Form( "U/O-flow (S,B): (%.1f, %.1f)%% / (%.1f, %.1f)%%", 
                        sig->GetBinContent(0)*dxu*100, bgd->GetBinContent(0)*dxu*100,
                        sig->GetBinContent(nbin+1)*dxo*100, bgd->GetBinContent(nbin+1)*dxo*100 );
      }
  
      TText* t = new TText( 0.98, 0.14, uoflow );
      t->SetNDC();
      t->SetTextSize( 0.040 );
      t->SetTextAngle( 90 );
//       t->AppendPad();    

      // save canvas to file
      if (countPad%noPadPerCanv==0) {
         TString fname = Form( "plots/%s_c%i", outfname.Data(), countCanvas );
         TMVAGlob::plot_logo();
         TMVAGlob::imgconv( canv, fname );
         createNewFig = kFALSE;
      }
      else {
         createNewFig = kTRUE;
      }
   }
   
   if (createNewFig) {
      TString fname = Form( "plots/%s_c%i", outfname.Data(), countCanvas );
      TMVAGlob::plot_logo();
      TMVAGlob::imgconv( canv, fname );
      createNewFig = kFALSE;
   }

   return;
}

void likelihoodrefs( TDirectory *lhdir ) {
   Bool_t newCanvas = kTRUE;
   
   const UInt_t maxCanvas = 200;
   TCanvas** c = new TCanvas*[maxCanvas];
   Int_t width  = 670;
   Int_t height = 380;

   // avoid duplicated printing
   std::vector<std::string> hasBeenUsed;
   const TString titName = lhdir->GetName();
   UInt_t ic = -1;   

   TIter next(lhdir->GetListOfKeys());
   TKey *key;
   while ((key = TMVAGlob::NextKey(next,"TH1"))) { // loop over all TH1
      TH1 *h = (TH1*)key->ReadObj();
      TH1F *b( 0 );
      TString hname( h->GetName() );

      // avoid duplicated plotting
      Bool_t found = kFALSE;
      for (UInt_t j = 0; j < hasBeenUsed.size(); j++) {
         if (hasBeenUsed[j] == hname.Data()) found = kTRUE;
      }
      if (!found) {

         // draw original plots
         if (hname.EndsWith("_sig_nice")) {

            if (newCanvas) {
               char cn[20];
               sprintf( cn, "cv%d_%s", ic+1, titName.Data() );
               ++ic;
               TString n = hname;	  
               c[ic] = new TCanvas( cn, Form( "%s reference for variable: %s", 
                                              titName.Data(),(n.ReplaceAll("_sig","")).Data() ), 
                                    ic*50+50, ic*20, width, height ); 
               c[ic]->Divide(2,1);
               newCanvas = kFALSE;
            }      

            // signal
            Int_t color = 4; 
            TPad * cPad = (TPad*)c[ic]->cd(1);
            TString plotname = hname;

            h->SetMaximum(h->GetMaximum()*1.3);
            h->SetMinimum( 0 );
            h->SetMarkerColor(color);
            h->SetMarkerSize( 0.7 );
            h->SetMarkerStyle( 24 );
            h->SetLineWidth(1);
            h->SetLineColor(color);
            color++;
            h->Draw("e1");
            Double_t hSscale = 1.0/(h->GetSumOfWeights()*h->GetBinWidth(1));

            TLegend *legS= new TLegend( cPad->GetLeftMargin(), 
                                        1-cPad->GetTopMargin()-.14, 
                                        cPad->GetLeftMargin()+.77, 
                                        1-cPad->GetTopMargin() );
            legS->SetBorderSize(1);
            legS->AddEntry(h,"Input data (signal)","p");

            // background
            TString bname( hname );	
            b = (TH1F*)lhdir->Get( bname.ReplaceAll("_sig","_bgd") );
            cPad = (TPad*)c[ic]->cd(2);
            color = 2;
            b->SetMaximum(b->GetMaximum()*1.3);
            b->SetMinimum( 0 );
            b->SetLineWidth(1);
            b->SetLineColor(color);
            b->SetMarkerColor(color);
            b->SetMarkerSize( 0.7 );
            b->SetMarkerStyle( 24 );
            b->Draw("e1");       
            Double_t hBscale = 1.0/(b->GetSumOfWeights()*b->GetBinWidth(1));
            TLegend *legB= new TLegend( cPad->GetLeftMargin(), 
                                        1-cPad->GetTopMargin()-.14, 
                                        cPad->GetLeftMargin()+.77, 
                                        1-cPad->GetTopMargin() );
            legB->SetBorderSize(1);
            legB->AddEntry(b,"Input data (backgr.)","p");

            // register
            hasBeenUsed.push_back( bname.Data() );

            // the PDFs --------------

            // check for splines
            h = 0;
            b = 0;
            TString pname = hname; pname.ReplaceAll("_nice","");            
            for (int i=0; i<= 5; i++) {
               TString hspline = pname + Form( "_smoothed_hist_from_spline%i", i );
               h = (TH1F*)lhdir->Get( hspline );
               if (h) {
                  b = (TH1F*)lhdir->Get( hspline.ReplaceAll("_sig","_bgd") );
//.........这里部分代码省略.........

//
//----------------------------------------------------------------------
//
int fittp0( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){

    cout << hs << " does not exist\n";

  }

  else{
   
    h->SetMarkerStyle(21);
    h->SetMarkerSize(0.8);
    h->SetStats(1);
    gStyle->SetOptFit(101);

    gROOT->ForceStyle();

    double dx = h->GetBinWidth(1);
    double nmax = h->GetBinContent(h->GetMaximumBin());
    double xmax = h->GetBinCenter(h->GetMaximumBin());
    double nn = 7*nmax;

    int nb = h->GetNbinsX();
    double n1 = h->GetBinContent(1);
    double n9 = h->GetBinContent(nb);
    double bg = 0.5*(n1+n9);

    double x1 = h->GetBinCenter(1);
    double x9 = h->GetBinCenter(nb);
    cout << hs << ": " << x1 << " - " << x9 << endl;

    // create a TF1 with the range from x1 to x9 and 5 parameters

    TF1 *tp0Fcn = new TF1( "tp0Fcn", tp0Fit, x1, x9, 5 );

    tp0Fcn->SetParName( 0, "mean" );
    tp0Fcn->SetParName( 1, "sigma" );
    tp0Fcn->SetParName( 2, "nu" );
    tp0Fcn->SetParName( 3, "area" );
    tp0Fcn->SetParName( 4, "BG" );

    tp0Fcn->SetNpx(500);
    tp0Fcn->SetLineWidth(4);
    tp0Fcn->SetLineColor(kMagenta);
    tp0Fcn->SetLineColor(kGreen);
   
    // set start values for some parameters:

    cout << hs << " " << dx << ", " << nn << ", " << xmax << endl;

    tp0Fcn->SetParameter( 0, xmax ); // peak position
    tp0Fcn->SetParameter( 1, 4*dx ); // width
    tp0Fcn->SetParameter( 2, 2.2 ); // nu
    tp0Fcn->SetParameter( 3, nn ); // N
    tp0Fcn->SetParameter( 4, bg );
    
    h->Fit( "tp0Fcn", "R", "ep" );
    // h->Fit("tp0Fcn","V+","ep");

    h->Draw("histepsame");  // data again on top
  }
}

//.........这里部分代码省略.........
        else {
          if(emptyBins > maxEmptyBins) maxEmptyBins = emptyBins;
          emptyBins = 0;
        }
      }
      rebinFactor[pdf.first] = std::max(maxEmptyBins, emptyBins);
    }

    for(auto& pdf : pdfs) rebinFactor[pdf.first] = std::max(rebinFactor[pdf.first], rebinFactor[switchQG(pdf.first)]);	// Use same rebin factor in quark and gluon pdf (otherwise bias if second derivate of pdf is non-zero)
    
    for(auto& pdf : pdfs){
      if(rebinFactor[pdf.first] > 19)     std::cout << "This pdf has a lot of emtpy bins and fluctuations:" << std::endl;
      if(rebinFactor[pdf.first] > 9)      rebin(pdf.second, 20);
      else if(rebinFactor[pdf.first] > 4) rebin(pdf.second, 10);
      else if(rebinFactor[pdf.first] > 3) rebin(pdf.second, 5);
      else if(rebinFactor[pdf.first] > 1) rebin(pdf.second, 4);
      else if(rebinFactor[pdf.first] > 0) rebin(pdf.second, 2);
    }


    // Normalization of the pdf's
    for(auto& pdf : pdfs) pdf.second->Scale(1./pdf.second->Integral(0, pdf.second->GetNbinsX() + 1));			// Scale to integral=1 (also include underflow/overflow)


    // Try to average out leftover fluctuations and empty bins
    for(auto& pdf : pdfs){
      if(pdf.first.Contains("gluon")) continue;

      TH1* tempQ = (TH1*) pdf.second->Clone();
      TH1* tempG = (TH1*) pdfs[switchQG(pdf.first)]->Clone();
      for(int i = 1; i < pdf.second->GetNbinsX() + 1; ++i){								// Do not consider underflow/overflow
        float contentQ = tempQ->GetBinContent(i);
        float contentG = tempG->GetBinContent(i);
        float width    = tempQ->GetBinWidth(i);

        if((1.5*contentQ < tempQ->GetBinContent(i-1) && 1.5*contentQ < tempQ->GetBinContent(i+1)) ||			// Try to average out some extreme fluctuations (i.e. only allow a difference of max 50% between two neighbouring bins)
           (1.5*contentG < tempG->GetBinContent(i-1) && 1.5*contentG < tempG->GetBinContent(i+1)) ||
           (contentQ < 1.5*tempQ->GetBinContent(i-1) && contentQ < 1.5*tempQ->GetBinContent(i+1)) ||
           (contentG < 1.5*tempG->GetBinContent(i-1) && contentG < 1.5*tempG->GetBinContent(i+1))){
          if(i-1 > 0 && i+1 < pdf.second->GetNbinsX() + 1){
            contentQ += tempQ->GetBinContent(i-1) + tempQ->GetBinContent(i+1);
            contentG += tempG->GetBinContent(i-1) + tempG->GetBinContent(i+1);
            width    += tempQ->GetBinWidth(i-1)   + tempQ->GetBinWidth(i+1);
          }
        }

        int j = 1;
        while(contentQ <= 0 || contentG <= 0){										// Average empty bins
          if(tempQ->Integral(0, i-j) <= 0) break;									// but not when surpassing the extreme edges of the pdf (see next part)
          if(tempG->Integral(0, i-j) <= 0) break;
          if(tempQ->Integral(i+j, tempQ->GetNbinsX()+1) <= 0) break;
          if(tempG->Integral(i+j, tempG->GetNbinsX()+1) <= 0) break;
          if(i-j == 0) break;
          if(i+j == pdf.second->GetNbinsX()) break;
          contentQ += tempQ->GetBinContent(i-j) + tempQ->GetBinContent(i+j);
          contentG += tempG->GetBinContent(i-j) + tempG->GetBinContent(i+j);
          width    += tempQ->GetBinWidth(i-j)   + tempQ->GetBinWidth(i+j);
          ++j;
        }

        pdf.second->SetBinContent(i, contentQ/width);
        pdfs[switchQG(pdf.first)]->SetBinContent(i, contentG/width);
      }
      delete tempQ;
      delete tempG;
    }

void ptBestFit(float BIN_SIZE=5.0,bool BLIND=false,TString MASS,TString NAME)
{
  gROOT->ProcessLine(".x ../../common/styleCMSTDR.C");
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
  gROOT->ForceStyle();
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gROOT->SetBatch(1);
  gStyle->SetPadRightMargin(0.04);
  gStyle->SetPadLeftMargin(0.16);
  gStyle->SetPadTopMargin(0.06);
  gStyle->SetPadBottomMargin(0.10);
  gStyle->SetTitleFont(42,"XY");
  gStyle->SetTitleSize(0.0475,"XY");
  gStyle->SetTitleOffset(0.9,"X");
  gStyle->SetTitleOffset(1.5,"Y");
  gStyle->SetLabelSize(0.0375,"XY");

  RooMsgService::instance().setSilentMode(kTRUE);
  for(int i=0;i<2;i++) {
    RooMsgService::instance().setStreamStatus(i,kFALSE);
  }
  float XMIN = 80;
  float XMAX = 200; 

  TFile *f1 = TFile::Open("datacards/datacard_m"+MASS+"_"+NAME+".root");
  TFile *f2 = TFile::Open("combine/mlfit.vbfHbb_"+NAME+"_mH"+MASS+".root");
  TFile *f3 = TFile::Open("root/sig_shapes_workspace_B80-200.root");
  TFile *f4 = TFile::Open("root/data_shapes_workspace_"+NAME+".root");

  RooWorkspace *w = (RooWorkspace*)f1->Get("w");
  //w->Print();
  RooAbsPdf *bkg_model = (RooAbsPdf*)w->pdf("model_s");
  RooFitResult *res_s  = (RooFitResult*)f2->Get("fit_s"); 
  RooFitResult *res_b  = (RooFitResult*)f2->Get("fit_b");
  RooRealVar *rFit     = dynamic_cast<RooRealVar *>(res_s->floatParsFinal()).find("r");
  RooDataSet *data     = (RooDataSet*)w->data("data_obs");
  
  int nparS=0,nparB=0;
  cout << res_s->floatParsFinal().getSize() << endl;
  cout << res_b->floatParsFinal().getSize() << endl;
  nparS = res_s->floatParsFinal().getSize();
  nparB = res_b->floatParsFinal().getSize();  
  float chi2sumS = 0.;
  float chi2sumB = 0.;
  int nparsum = 0;
//  if (BLIND) {
//    res_b->Print();
//  }
//  else {
//    res_s->Print();
//  }
  
  w->allVars().assignValueOnly(res_s->floatParsFinal());
//  w->Print();
//  w->allVars()->Print();

  RooWorkspace *wSig = (RooWorkspace*)f3->Get("w"); 
  RooWorkspace *wDat = (RooWorkspace*)f4->Get("w"); 

  const RooSimultaneous *sim = dynamic_cast<const RooSimultaneous *> (bkg_model);
  const RooAbsCategoryLValue &cat = (RooAbsCategoryLValue &) sim->indexCat();
  TList *datasets = data->split(cat,true);
  TIter next(datasets);
  //int count = 0; 
  for(RooAbsData *ds = (RooAbsData*)next();ds != 0; ds = (RooAbsData*)next()) {
	 //if (count > 0) return 0;
	 //count++;
    RooAbsPdf *pdfi = sim->getPdf(ds->GetName());
    RooArgSet *obs = (RooArgSet*)pdfi->getObservables(ds);
    RooRealVar *x = dynamic_cast<RooRealVar *>(obs->first());

    RooRealVar *yield_vbf = (RooRealVar*)wSig->var("yield_signalVBF_mass"+MASS+"_"+TString(ds->GetName()));
    RooRealVar *yield_gf  = (RooRealVar*)wSig->var("yield_signalGF_mass"+MASS+"_"+TString(ds->GetName()));
    TString ds_name(ds->GetName());
    //----- get the QCD normalization -----------
    RooRealVar *qcd_norm_final = dynamic_cast<RooRealVar *>(res_s->floatParsFinal()).find("CMS_vbfbb_qcd_norm_"+ds_name);
    RooRealVar *qcd_yield      = (RooRealVar*)wDat->var("yield_data_"+ds_name);

    float Nqcd  = exp(log(1.5)*qcd_norm_final->getVal())*qcd_yield->getVal();
    float eNqcd = log(1.5)*qcd_norm_final->getError()*Nqcd;
    cout<<"QCD normalization = "<<Nqcd<<" +/- "<<eNqcd<<endl;
    
    TH1 *hCoarse = (TH1*)ds->createHistogram("coarseHisto_"+ds_name,*x);
    float norm = hCoarse->Integral();
  
	 int rebin = BIN_SIZE/hCoarse->GetBinWidth(1);
    hCoarse->Rebin(rebin);

    float MIN_VAL = TMath::Max(0.9*hCoarse->GetBinContent(hCoarse->GetMinimumBin()),1.0);
    float MAX_VAL = 1.3*hCoarse->GetBinContent(hCoarse->GetMaximumBin());
    RooDataHist ds_coarse("ds_coarse_"+ds_name,"ds_coarse_"+ds_name,*x,hCoarse);

    TH1F *hBlind = (TH1F*)hCoarse->Clone("blindHisto_"+ds_name);
    for(int i=0;i<hBlind->GetNbinsX();i++) {
      double x0 = hBlind->GetBinCenter(i+1);
      if (x0 > 100 && x0 < 150) {
        hBlind->SetBinContent(i+1,0);
        hBlind->SetBinError(i+1,0);
      }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      c3->cd(count + 1)->SetGrid();
      Printf("%s", s1namegen.Data());
      hg = (TH1D *)PullHisto(l, s1namegen.Data(), bf[i], bl[i], ptmean);
      hg->SetLineColor(kMagenta);
      hg->GetXaxis()->SetTitle("inv. mass, GeV/c^2");
      hg->Draw("hist");
      c4->cd(count + 1)->SetGrid();
      ht = (TH1D *)PullHisto(l, s1nametrue.Data(), bf[i], bl[i], ptmean);
      ht->SetLineColor(kMagenta-5);
      ht->GetXaxis()->SetTitle("inv. mass, GeV/c^2");
      ht->Draw("hist");
    }
    c2->cd(count + 1)->SetGrid();
    TH1 *hh = (TH1 *)h1->Clone("hh");
    hh->SetLineColor(kRed+1);
    hh->Add(h3_p, -1);
    /// !!!!!!!!!!!!!!!!!!!!!!
    //////////    if ((ilist == 3) && (count < 2)) hh->Reset();
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    hh->Draw("hist");

    // !!!!!!!!!!!!!!!!!!
    ff->SetParameters(0.1, 1.02, 0.004, -25000., 0., 0., 0.);
    ff->SetLineColor(hh->GetLineColor());
    ff->SetLineWidth(1);
    //    ff->SetLineStyle(kDashed);
    // where fit
    Double_t fmin = 1.02-2*0.004;
    Double_t fmax = 1.02+2*0.004;
    //    Double_t fmin = 0.995;
    //    Double_t fmax = 1.185;
    // !!!!!!!!!!!!!!!!!!
    Bool_t hisfun = kFALSE; // kFALSE = integral from function
    Double_t   hisfun_k = 1.0/hh->GetBinWidth(10);
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if (binhaluska)
      if (i > 9) hisfun_k = 0.5/hh->GetBinWidth(10);
    Printf("======= %f", hisfun_k);
    // !!!!!!!!!!!!!!!!!!
    // wehere integral (his or fun)
    Double_t fmini = 1.02-2*0.004;
    Double_t fmaxi = 1.02+2*0.004;
    hh->Fit(ff, "Q", "", fmin, fmax);
    hh->Fit(ff, "Q", "", fmin, fmax);
    fitStatus = hh->Fit(ff, "Q", "", fmin, fmax);
    TF1 *pp3 = new TF1("pp3", "[0]+x*[1]+x*x*[2]+x*x*x*[3]", fmin, fmax);
    pp3->SetParameters(ff->GetParameter(3), ff->GetParameter(4),
                       ff->GetParameter(5), ff->GetParameter(6));
    pp3->SetLineWidth(1);
    pp3->SetLineColor(h3_p->GetLineColor());
    pp3->Draw("same");
    // ff->SetRange(fmin, fmax);
    // ff->DrawCopy("same");

    value              = hh->Integral(hh->FindBin(fmini), hh->FindBin(fmaxi));
    if (!hisfun) value = ff->Integral(fmini, fmaxi)*hisfun_k -
                   pp3->Integral(fmini, fmaxi)*hisfun_k;
    if (value < 0) value = 0;

    if ((fitStatus != 0) || (ff->GetParameter(2) > 0.1)) {
      printf(" SKIP Data");
      value = 0;
    }
    grx[count] = ptmean;
    if (binhaluska) {
      if (count < 10) grxE[count] = 0.25; // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

void MakePi0Analysis(){

  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadRightMargin(0.01);
  gStyle->SetPadTopMargin(0.09);
  gStyle->SetPadBottomMargin(0.11);
  //gStyle->SetOptStat(0);
  //gStyle->SetOptTitle(1);
  //gStyle->SetPadTickX(1);
  //gStyle->SetPadTickY(1);
  TGaxis::SetMaxDigits(3);
  
  Double_t bins[] = {1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0, 3.2, 3.4, 3.6,3.8,4.0,4.5, 5.0, 5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,10.0,11,12}; //Nicolas LHC12d PHI
  const Int_t binnum = sizeof(bins)/sizeof(Double_t) - 1;
  
  TString cwd = gSystem->WorkingDirectory();

  TH1F* fHistPeakMean  = new TH1F("fHistPeakMean","",binnum,bins);
  TH1F* fHistPeakWidth = new TH1F("fHistPeakWidth","",binnum,bins);
  TH1F* fHistRawYield  = new TH1F("fHistRawYield","",binnum,bins);
  TH1F* fHistRawYieldPerEvent  = new TH1F("fHistRawYieldPerEvent","",binnum,bins);
  
  TCanvas *c1[binnum];
  TCanvas *c2[binnum];
  TCanvas *c3;

  TFile* ef     = TFile::Open("../AnalysisResults.root");
  TList* list   = (TList*)ef->Get("list_kINT7_Pi0");
  TList *listEv = (TList*)ef->Get("list_kINT7_Event");
  
  TH1F* fHistEvents = (TH1F*)listEv->FindObject("fHistAnalyzedEvents");
  Double_t nEvt = fHistEvents->GetEntries();
  
  TString tofName[] = {"","_TOFcut1","_TOFcut2"};
  TString fitName[] = {"CrystalBall","AsymmGauss"};
  TString modName[] = {"","_M1","_M3"};
  
  Int_t ntof = 3;
  
  Double_t signal_range_min=0.11;
  Double_t signal_range_max=0.15;
  Double_t bg_range_min=0.050;
  Double_t bg_range_max=0.250;
  
  for(Int_t itof=0; itof<ntof; ++itof){
    for(Int_t iMod=0; iMod<1; iMod++){
      for(Int_t iFit=0; iFit<2; iFit++){
	for(Int_t iPol=0; iPol<3; iPol++){
	  
	  Int_t pol = iPol;

	  TString fNameSame = "fHistMassTwoGammas"+tofName[itof]+modName[iMod];
	  TString fNameMix  = "fHistMixMassTwoGammas"+tofName[itof]+modName[iMod];
	  
	  THnSparse *fTHnSparseRecPi0;
	  THnSparse *fTHnSparseRecMixPi0;
	  THnSparse *fTHnSparseGeneratePi0;
	  
	  TH2* hPi0A08    = 0x0;
	  TH2* hMixPi0A08 = 0x0;
	  hPi0A08    = (TH2*)list->FindObject(fNameSame)->Clone();
	  hMixPi0A08 = (TH2*)list->FindObject(fNameMix)->Clone();

	  Int_t sbin = 0;
	  Int_t lbin = 0;
      
	  TH1F* fHistFinalRatio[binnum];
	  TH1F* fHistFinalSame[binnum];
	  TH1F* fHistFinalBG[binnum];
	  TH1F* fHistFinalSignal[binnum];
	  TH1F* fHistOnlyFitSignal[binnum];
	  
	  TF1* fFitFinalRatio[binnum];
	  TF1* fFitFinalSignal[binnum];
	  TF1* fFitOnlyFitSignal[binnum];

	  for(Int_t i=0; i<binnum; i++){    
	    
	    Double_t bin_width = (bins[i+1]-bins[i]);
	    Double_t bin_mean  = (bins[i+1]+bins[i])/2;
	    
	    Int_t rebin = 1;
	    
	    if(bin_mean>5.0){
	      rebin=10;
	    }
	    else{
	      rebin=5;
	    }
	    
	    Int_t sproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i]);
	    Int_t lproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i+1])-1;
	    
	    TH1* fHistBasicSame = (TH1*)hPi0A08   ->ProjectionX(Form("fHistBasicSame_No%d",i+1),sproj_bin,lproj_bin,""); 
	    TH1* fHistBasicMix  = (TH1*)hMixPi0A08->ProjectionX(Form("fHistBasicMix_No%d",i+1),sproj_bin,lproj_bin,""); 
	    fHistBasicSame->Rebin(rebin);
	    fHistBasicMix->Rebin(rebin);
	    fHistBasicSame->Sumw2();
	    fHistBasicMix->Sumw2();
	    fHistBasicSame->GetYaxis()->SetTitle(Form("dN/dM per %.0f MeV/c^{2}",fHistBasicSame->GetBinWidth(1)*1000));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            legend->AddEntry(bgd,TString("Background") + ((htype == CompareType) ? " (test sample)" : ""), "F");
            legend->SetBorderSize(1);
            legend->SetMargin( (htype == CompareType ? 0.2 : 0.3) );
            legend->Draw("same");

            // overlay signal and background histograms
            sig->Draw("samehist");
            bgd->Draw("samehist");
   
            if (htype == CompareType) {
               // if overtraining check, load additional histograms
               TH1* sigOv = 0;
               TH1* bgdOv = 0;

               TString ovname = hname += "_Train";
               sigOv = dynamic_cast<TH1*>(titDir->Get( ovname + "_S" ));
               bgdOv = dynamic_cast<TH1*>(titDir->Get( ovname + "_B" ));
      
               if (sigOv == 0 || bgdOv == 0) {
                  cout << "+++ Problem in \"mvas.C\": overtraining check histograms do not exist" << endl;
               }
               else {
                  cout << "--- Found comparison histograms for overtraining check" << endl;

                  TLegend *legend2= new TLegend( 1 - c->GetRightMargin() - 0.42, 1 - c->GetTopMargin() - 0.12,
                                                 1 - c->GetRightMargin(), 1 - c->GetTopMargin() );
                  legend2->SetFillStyle( 1 );
                  legend2->SetBorderSize(1);
                  legend2->AddEntry(sigOv,"Signal (training sample)","P");
                  legend2->AddEntry(bgdOv,"Background (training sample)","P");
                  legend2->SetMargin( 0.1 );
                  legend2->Draw("same");
               }
               Int_t col = sig->GetLineColor();
               sigOv->SetMarkerColor( col );
               sigOv->SetMarkerSize( 0.7 );
               sigOv->SetMarkerStyle( 20 );
               sigOv->SetLineWidth( 1 );
               sigOv->SetLineColor( col );
               sigOv->Draw("e1same");
      
               col = bgd->GetLineColor();
               bgdOv->SetMarkerColor( col );
               bgdOv->SetMarkerSize( 0.7 );
               bgdOv->SetMarkerStyle( 20 );
               bgdOv->SetLineWidth( 1 );
               bgdOv->SetLineColor( col );
               bgdOv->Draw("e1same");

               ymax = TMath::Max( ymax, TMath::Max( sigOv->GetMaximum(), bgdOv->GetMaximum() )*maxMult );
               frame->GetYaxis()->SetLimits( 0, ymax );
      
               // for better visibility, plot thinner lines
               sig->SetLineWidth( 1 );
               bgd->SetLineWidth( 1 );

               // perform K-S test
               cout << "--- Perform Kolmogorov-Smirnov tests" << endl;
               Double_t kolS = sig->KolmogorovTest( sigOv );
               Double_t kolB = bgd->KolmogorovTest( bgdOv );
               cout << "--- Goodness of signal (background) consistency: " << kolS << " (" << kolB << ")" << endl;

               TString probatext = Form( "Kolmogorov-Smirnov test: signal (background) probability = %5.3g (%5.3g)", kolS, kolB );
               TText* tt = new TText( 0.12, 0.74, probatext );
               tt->SetNDC(); tt->SetTextSize( 0.032 ); tt->AppendPad(); 
            }

            // redraw axes
            frame->Draw("sameaxis");

            // text for overflows
            Int_t    nbin = sig->GetNbinsX();
            Double_t dxu  = sig->GetBinWidth(0);
            Double_t dxo  = sig->GetBinWidth(nbin+1);
            TString uoflow = Form( "U/O-flow (S,B): (%.1f, %.1f)%% / (%.1f, %.1f)%%", 
                                   sig->GetBinContent(0)*dxu*100, bgd->GetBinContent(0)*dxu*100,
                                   sig->GetBinContent(nbin+1)*dxo*100, bgd->GetBinContent(nbin+1)*dxo*100 );
            TText* t = new TText( 0.975, 0.115, uoflow );
            t->SetNDC();
            t->SetTextSize( 0.030 );
            t->SetTextAngle( 90 );
            t->AppendPad();    
   
            // update canvas
            c->Update();

            // save canvas to file

            TMVAGlob::plot_logo(1.058);
            if (Save_Images) {
               if      (htype == MVAType)     TMVAGlob::imgconv( c, Form("plots/mva_%s",     methodTitle.Data()) );
               else if (htype == ProbaType)   TMVAGlob::imgconv( c, Form("plots/proba_%s",   methodTitle.Data()) ); 
               else if (htype == CompareType) TMVAGlob::imgconv( c, Form("plots/overtrain_%s", methodTitle.Data()) ); 
               else                           TMVAGlob::imgconv( c, Form("plots/rarity_%s",  methodTitle.Data()) ); 
            }
            countCanvas++;
         }
      }
   }
}

Bool_t fitsHere(TLegend *l,Double_t x1, Double_t y1, Double_t x2, Double_t y2)
{
    Bool_t fits = true;
    TList *list = l->GetListOfPrimitives();
    for (Int_t k = 0; list->At(k) != 0 && fits; k++)
    {
        TObject *obj = ((TLegendEntry*)(list->At(k)))->GetObject();
        if (obj == 0) continue;
        TClass *cl = obj->IsA();

        //Histogram, drawn as a histogram
        if (cl->InheritsFrom("TH1") && !cl->InheritsFrom("TH2") && !cl->InheritsFrom("TH3")
         && cl != TProfile::Class() && ((TH1*)obj)->GetMarkerColor() == kWhite)
        {
            Int_t where = 0;
            TH1 *h = (TH1*)obj;
            for (Int_t i = 1; i <= h->GetNbinsX() && fits; i++)
            {
                if (h->GetBinLowEdge(i) + h->GetBinWidth(i) < x1) continue;   //to the left of the legend
                if (h->GetBinLowEdge(i)                     > x2) continue;   //to the right of the legend
                if (h->GetBinContent(i) > y1 && h->GetBinContent(i) < y2) fits = false;   //inside the legend
                if (h->GetBinContent(i) < y1)
                {
                    if (where == 0) where = -1;             //below the legend
                    if (where == 1) fits = false;           //a previous bin was above it so there's a vertical line through it
                }
                if (h->GetBinContent(i) > y2)
                {
                    if (where == 0) where = 1;              //above the legend
                    if (where == -1) fits = false;          //a previous bin was below it so there's a vertical line through it
                }
            }
            continue;
        }
        //Histogram, drawn with Draw("P")
        else if (cl->InheritsFrom("TH1") && !cl->InheritsFrom("TH2") && !cl->InheritsFrom("TH3")
              && cl != TProfile::Class())
        //Probably TProfile would be the same but I haven't tested it
        {
            TH1 *h = (TH1*)obj;
            for (Int_t i = 1; i <= h->GetNbinsX() && fits; i++)
            {
                if (h->GetBinLowEdge(i) + h->GetBinWidth(i)/2 < x1) continue;
                if (h->GetBinLowEdge(i)                       > x2) continue;
                if (h->GetBinContent(i) > y1 && h->GetBinContent(i) < y2) fits = false;
                if (h->GetBinContent(i) + h->GetBinError(i) > y2 && h->GetBinContent(i) - h->GetBinError(i) < y2) fits = false;
                if (h->GetBinContent(i) + h->GetBinError(i) > y1 && h->GetBinContent(i) - h->GetBinError(i) < y1) fits = false;
            }
        }
        else if (cl->InheritsFrom("TF1") && !cl->InheritsFrom("TF2"))
        {
            //TF1 *f = (TF1*)obj;
            //Double_t max = f->GetMaximum(x1,x2);
            //Double_t min = f->GetMinimum(x1,x2);
            //if (min < y2 && max > y1) fits = false;
        }
        // else if (cl->InheritsFrom(...... add more objects here
        else
        {
            cout << "Don't know how to place the legend around objects of type " << obj->ClassName() << "." << endl
                 << "Add this class into placeLegend.C if you want it to work properly." << endl
                 << "The legend will still be placed around any other objects." << endl;
        }
    }
    return fits;
}