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

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;
}

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;
}

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);
}

//.........这里部分代码省略.........
    if (int1 <= 0.) continue;
    Float_t int2 = hcorr->Integral(i + 1, i + 1, 2, -1);
    hnum->SetBinContent(i + 1, int2);
    hnum->SetBinError(i + 1, TMath::Sqrt(int2));
    hden->SetBinContent(i + 1, int1);
    hden->SetBinError(i + 1, TMath::Sqrt(int1));
  }
  TCanvas *cEfficiency = new TCanvas("cEfficiency", "cEfficiency");
  cEfficiency->SetLogx();
  cEfficiency->SetLogy();

  heff->Divide(hnum, hden, 1., 1., "B");
  heff->Draw();
#if 0
  for (Int_t ii = 0; ii < heff->GetNbinsX(); ii++) {
    heff->SetBinContent(ii + 1, 1.);
    heff->SetBinError(ii + 1, 0.);
  }
#endif

  for (Int_t i = 0; i < hcorr->GetNbinsX(); i++) {
    hcorr->SetBinContent(i + 1, 1, 0.);
    hcorr->SetBinError(i + 1, 1, 0.);
  }
  for (Int_t i = 0; i < hcorr->GetNbinsY(); i++) {
    hcorr->SetBinContent(1, i + 1, 0.);
    hcorr->SetBinError(1, i + 1, 0.);
  }
  TH2 *hcorrfit = ReturnCorrFromFit(hcorr);
  // Docs from AliUnfolding
  //Int_t AliUnfolding::Unfold(TH2* correlation, TH1* efficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
  // unfolds with unfolding method fgMethodType
  //
  // parameters:
  //  correlation: response matrix as measured vs. generated
  //  efficiency:  (optional) efficiency that is applied on the unfolded spectrum, i.e. it has to be in unfolded variables. If 0 no efficiency is applied.
  //  measured:    the measured spectrum
  //  initialConditions: (optional) initial conditions for the unfolding. if 0 the measured spectrum is used as initial conditions.
  //  result:      target for the unfolded result
  //  check:       depends on the unfolding method, see comments in specific functions

  for (Int_t iloop = 0; iloop < nloop; iloop++) {
    if (bayesian) {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
      AliUnfolding::SetBayesianParameters(smooth, iter);
    } else {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
      AliUnfolding::SetChi2Regularization(AliUnfolding::RegularizationType(regul), weight);
    }
    AliUnfolding::SetSkip0BinInChi2(kTRUE);
    AliUnfolding::SetSkipBinsBegin(1);
    AliUnfolding::SetNbins(150, 150);
    AliUnfolding::Unfold(usecorrfit ? hcorrfit : hcorr, heff, hdata, hinit, hresu);
    hinit = (TH1 *)hresu->Clone(Form("hinit_%d", iloop));
  }

  printf("hdata->Integral(2, -1) = %f\n", hdata->Integral(2, -1));
  printf("hresu->Integral(2, -1) = %f\n", hresu->Integral(2, -1));


  TCanvas *cUnfolded = new TCanvas ("cUnfolded", "cUnfolded", 400, 800);
  cUnfolded->Divide(1, 2);
  cUnfolded->cd(1)->SetLogx();
  cUnfolded->cd(1)->SetLogy();
  hdata->Draw();
  hresu->Draw("same");
  htrue->Draw("same");
  cUnfolded->cd(2)->SetLogx();
  cUnfolded->cd(2)->DrawFrame(1., 0, 300., 10);
  TH1 *hrat = (TH1 *)hresu->Clone("hrat");
  hrat->Divide(htrue);
  hrat->Draw("same");

  TH1 *htrig = (TH1 *)hresu->Clone("htrig");
  htrig->Multiply(heff); 
  

  Float_t dndeta_resu = 0.;
  Float_t integr_resu = 0.;
  Float_t dndeta_trig = 0.;
  Float_t integr_trig = 0.;
  for (Int_t i = 1; i < hresu->GetNbinsX(); i++) {
    dndeta_resu += hresu->GetBinContent(i + 1) * hresu->GetBinLowEdge(i + 1);
    integr_resu += hresu->GetBinContent(i + 1);
    dndeta_trig += htrig->GetBinContent(i + 1) * htrig->GetBinLowEdge(i + 1);
    integr_trig += htrig->GetBinContent(i + 1);
  }

  cUnfolded->SaveAs("unfold_efficiency.pdf");

  integr_eff = integr_trig / integr_resu;
  integr_eff_err = TMath::Sqrt(integr_eff * (1. - integr_eff) / integr_resu);
  dndeta_eff = dndeta_trig / dndeta_resu;
  dndeta_eff_err = TMath::Sqrt(dndeta_eff * (1. - dndeta_eff) / dndeta_resu);

  printf("INEL > 0 efficiency: %.3f +- %.3f\n", integr_eff, integr_eff_err);
  printf("dN/dEta correction:  %.3f +- %.3f\n", dndeta_eff, dndeta_eff_err);

  return hresu;
}

//.........这里部分代码省略.........
  TH1 *hbias = (TH1 *)hdata->Clone("hbias");
  hresu->SetMarkerStyle(20);
  hresu->SetMarkerColor(4);
  hresu->Reset();

  TH1 *hnum = hcorr->ProjectionY("hnum");
  TH1 *hden = hcorr->ProjectionY("hden");
  TH1 *heff = hcorr->ProjectionY("heff");
  hnum->Reset();
  hnum->Sumw2();
  hden->Reset();
  hden->Sumw2();
  heff->Reset();
  for (Int_t i = 0; i < heff->GetNbinsX(); i++) {
    Float_t int1 = hcorr->Integral(i + 1, i + 1, 0, -1);
    if (int1 <= 0.) continue;
    Float_t int2 = hcorr->Integral(i + 1, i + 1, 2, -1);
    hnum->SetBinContent(i + 1, int2);
    hnum->SetBinError(i + 1, TMath::Sqrt(int2));
    hden->SetBinContent(i + 1, int1);
    hden->SetBinError(i + 1, TMath::Sqrt(int1));
  }
  new TCanvas("cEfficiency");
  heff->Divide(hnum, hden, 1., 1., "B");
  heff->Draw();
#if 0
  for (Int_t ii = 0; ii < heff->GetNbinsX(); ii++) {
    heff->SetBinContent(ii + 1, 1.);
    heff->SetBinError(ii + 1, 0.);
  }
#endif
  
  for (Int_t i = 0; i < hcorr->GetNbinsX(); i++) {
    hcorr->SetBinContent(i + 1, 1, 0.);
    hcorr->SetBinError(i + 1, 1, 0.);
  }
  for (Int_t i = 0; i < hcorr->GetNbinsY(); i++) {
    hcorr->SetBinContent(1, i + 1, 0.);
    hcorr->SetBinError(1, i + 1, 0.);
  }
  TH2 *hcorrfit = ReturnCorrFromFit(hcorr);

  for (Int_t iloop = 0; iloop < nloop; iloop++) {
    if (bayesian) {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
      AliUnfolding::SetBayesianParameters(smooth, iter);
    } else {
      AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
      AliUnfolding::SetChi2Regularization(regul, weight);
    }
    AliUnfolding::SetSkip0BinInChi2(kTRUE);
    AliUnfolding::SetSkipBinsBegin(1);
    AliUnfolding::SetNbins(150, 150);
    AliUnfolding::Unfold(usecorrfit ? hcorrfit : hcorr, heff, hdata, hinit, hresu);
    hinit = (TH1 *)hresu->Clone(Form("hinit_%d", iloop));
  }

  printf("hdata->Integral(2, -1) = %f\n", hdata->Integral(2, -1));
  printf("hresu->Integral(2, -1) = %f\n", hresu->Integral(2, -1));
  
  
  TCanvas *cUnfolded = new TCanvas ("cUnfolded", "", 400, 800);
  cUnfolded->Divide(1, 2);
  cUnfolded->cd(1)->SetLogx();
  cUnfolded->cd(1)->SetLogy();
  hdata->Draw();
  hresu->Draw("same");
  htrue->Draw("same");
  cUnfolded->cd(2)->SetLogx();
  cUnfolded->cd(2)->DrawFrame(1., 0.75, 300., 1.25);
  TH1 *hrat = (TH1 *)hresu->Clone("hrat");
  hrat->Divide(htrue);
  hrat->Draw("same");

  TH1 *htrig = (TH1 *)hresu->Clone("htrig");
  htrig->Multiply(heff);

  Float_t dndeta_resu = 0.;
  Float_t integr_resu = 0.;
  Float_t dndeta_trig = 0.;
  Float_t integr_trig = 0.;
  for (Int_t i = 1; i < hresu->GetNbinsX(); i++) {
    dndeta_resu += hresu->GetBinContent(i + 1) * hresu->GetBinLowEdge(i + 1);
    integr_resu += hresu->GetBinContent(i + 1);
    dndeta_trig += htrig->GetBinContent(i + 1) * htrig->GetBinLowEdge(i + 1);
    integr_trig += htrig->GetBinContent(i + 1);
  }
  //  dndeta_resu /= integr_resu;
  //  dndeta_trig /= integr_trig;

  integr_eff = integr_trig / integr_resu;
  integr_eff_err = TMath::Sqrt(integr_eff * (1. - integr_eff) / integr_resu);
  dndeta_eff = dndeta_trig / dndeta_resu;
  dndeta_eff_err = TMath::Sqrt(dndeta_eff * (1. - dndeta_eff) / dndeta_resu);
  
  printf("INEL > 0 efficiency: %.3f +- %.3f\n", integr_eff, integr_eff_err);
  printf("dN/dEta correction:  %.3f +- %.3f\n", dndeta_eff, dndeta_eff_err);

  return hresu;
}

void Fit(TString SIGNAL,TString HISTO,double scaleSGN)
{
  gROOT->ForceStyle();
  TFile *fDat = TFile::Open("Histo_flatTree_data_tmva"+SIGNAL+".root");
  TFile *fBkg = TFile::Open("Histo_flatTree_qcd_weights_tmva"+SIGNAL+".root");
  TFile *fSgn = TFile::Open("Histo_flatTree_"+SIGNAL+"_weights_tmva"+SIGNAL+".root");
  
  TH1 *hDat = (TH1*)fDat->Get(HISTO);
  TH1 *hBkgRaw = (TH1*)fBkg->Get(HISTO);
  TH1 *hSgn = (TH1*)fSgn->Get(HISTO);
  TH1 *hDat_JESlo = (TH1*)fDat->Get(HISTO+"_JESlo");
  TH1 *hBkgRaw_JESlo = (TH1*)fBkg->Get(HISTO+"_JESlo");
  TH1 *hSgn_JESlo = (TH1*)fSgn->Get(HISTO+"_JESlo");
  TH1 *hDat_JESup = (TH1*)fDat->Get(HISTO+"_JESup");
  TH1 *hBkgRaw_JESup = (TH1*)fBkg->Get(HISTO+"_JESup");
  TH1 *hSgn_JESup = (TH1*)fSgn->Get(HISTO+"_JESup");
  
  TH1F *hBkg = (TH1F*)hBkgRaw->Clone("Bkg");
  TH1F *hBkg_JESlo = (TH1F*)hBkgRaw_JESlo->Clone("Bkg_JESlo");
  TH1F *hBkg_JESup = (TH1F*)hBkgRaw_JESup->Clone("Bkg_JESup");
  
  hBkg->Smooth(2);
  hBkg_JESlo->Smooth(2);
  hBkg_JESup->Smooth(2);
  hSgn->Smooth(2);
  hSgn_JESlo->Smooth(2);
  hSgn_JESup->Smooth(2);
  
  double lumi = 4967;
  hBkg->Scale(lumi);
  hBkg_JESlo->Scale(lumi);
  hBkg_JESup->Scale(lumi);
  double k_factor = hDat->Integral()/hBkg->Integral();
  double k_factor_JESlo = hDat->Integral()/hBkg_JESlo->Integral();
  double k_factor_JESup = hDat->Integral()/hBkg_JESup->Integral();
  hBkg->Scale(k_factor);
  cout<<"Signal entries = "<<hSgn->GetEntries()<<endl;
  hSgn->Scale(lumi/scaleSGN);
  hBkg_JESlo->Scale(k_factor_JESlo);
  hSgn_JESlo->Scale(lumi/scaleSGN);
  hBkg_JESup->Scale(k_factor_JESup);
  hSgn_JESup->Scale(lumi/scaleSGN);
  hSgn_JESlo->Scale(hSgn->Integral()/hSgn_JESlo->Integral());
  hSgn_JESup->Scale(hSgn->Integral()/hSgn_JESup->Integral());
  
  TH1 *hBkg_STATlo = (TH1*)hBkg->Clone(HISTO+"_STATlo");
  TH1 *hSgn_STATlo = (TH1*)hSgn->Clone(HISTO+"_STATlo");
  TH1 *hBkg_STATup = (TH1*)hBkg->Clone(HISTO+"_STATup");
  TH1 *hSgn_STATup = (TH1*)hSgn->Clone(HISTO+"_STATup");
  
  float y1,e1;
  for(int i=0;i<hBkg->GetNbinsX();i++) {
    y1 = hBkg->GetBinContent(i+1);
    e1 = hBkg->GetBinError(i+1);
    hBkg_STATlo->SetBinContent(i+1,y1-e1);
    hBkg_STATup->SetBinContent(i+1,y1+e1);
    y1 = hSgn->GetBinContent(i+1);
    e1 = hSgn->GetBinError(i+1);
    hSgn_STATlo->SetBinContent(i+1,y1-e1);
    hSgn_STATup->SetBinContent(i+1,y1+e1);
  }
  hBkg_STATlo->Scale(hBkg->Integral()/hBkg_STATlo->Integral());
  hBkg_STATup->Scale(hBkg->Integral()/hBkg_STATup->Integral());
  hSgn_STATlo->Scale(hSgn->Integral()/hSgn_STATlo->Integral());
  hSgn_STATup->Scale(hSgn->Integral()/hSgn_STATup->Integral());
  
  double xMIN = hBkg->GetBinLowEdge(1);
  double xMAX = hBkg->GetBinLowEdge(hBkg->GetNbinsX()+1);
  double xMIN2 = hDat->GetBinLowEdge(hDat->FindFirstBinAbove(0.5));
  double xMAX2 = hDat->GetBinLowEdge(hDat->FindLastBinAbove(0.5)+1);
  RooRealVar x("x","x",xMIN2,xMAX2);
  
  RooDataHist data("data","dataset with x",x,hDat);
  RooDataHist bkg("qcd","bkg with x",x,hBkg);
  RooDataHist sgn("signal","sgn with x",x,hSgn);
  
  RooHistPdf bkgPDF("bkgPDF","bkgPDF",x,bkg,0);
  RooHistPdf sgnPDF("sgnPDF","sgnPDF",x,sgn,0);
  
  RooRealVar f("f","f",0,0.,1.);
  
  RooAddPdf model("model","model",RooArgList(sgnPDF,bkgPDF),RooArgList(f));
  
  RooFitResult* r = model.fitTo(data,Save());
  r->Print("v");
  double N = hDat->Integral();
  double B = hBkg->Integral();
  double S = hSgn->Integral();
  double m = f.getVal();
  double e = f.getError();
  cout<<"k-factor = "<<k_factor<<endl;
  cout<<N<<" "<<B<<" "<<S<<endl;
  cout<<"Total cross section =       "<<N/lumi<<" pb"<<endl;
  cout<<"Model cross section =       "<<S/lumi<<" pb"<<endl;
  cout<<"Fitted signal strength =    "<<m*N/S<<endl;
  cout<<"Fitted signal error =       "<<e*N/S<<endl;
  double p = 0.95;
  double xup = (N/S)*(m+sqrt(2.)*e*TMath::ErfInverse((1-p)*TMath::Erf(m/e)+p));
  cout<<"Bayesian Upper limit =      "<<xup<<endl;
  RooPlot* frame1 = x.frame();  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      TString selection =  "(" + cuts + "&&";
      if (isam == 18) selection+= "MuonEle && isMuEGTriggerPass";
      else if (isam == 19) selection += "TwoEle && isDoubleElectronTriggerPass";
      else if (isam == 20) selection += "TwoMuon && isDoubleMuTriggerPass";
      selection += ")";

      TString var = string(histnames[i])+">>+"+histName ;
      tmpTree->Draw(var, selection, "goff");


      tmpfile->Close();
      delete tmpfile;
    }

    std::cout << "tempHist name = " << tempHist->GetName() << std::endl;
    std::cout << "tempHist: Entries = " << tempHist->GetEntries() << ", Integral = " << tempHist->Integral() << std::endl;
    double threshold = (tempHist->Integral())/20;
 
   //Now, look through the bins of the histograms, combining when necessary
    int nBinsUsed = 0;
    double *xBins = new double[histbins[i]+1];
    xBins[0] = histmin[i];  //First bin always starts at minimum
    double binSum = 0;

    for (int iBin = 1; iBin <= histbins[i]; ++iBin) {

      binSum += tempHist->GetBinContent(iBin);
      std::cout << "binSum = " << binSum << ", iBin = " << iBin << std::endl;


      if (binSum >= threshold) {
        std::cout << "Setting bin edge!" << std::endl;
        ++nBinsUsed;
        xBins[nBinsUsed] = tempHist->GetBinLowEdge(iBin+1);
        binSum = 0.;
      }

    }

    //Now check: if binSum > 0, that means the last bin was too small to be a bin
    //by itself.  If so, we need to make sure it's combined with the bin on it's left by setting the bin edge to the histogram max
    if (binSum > 0) {
      xBins[nBinsUsed] = tempHist->GetBinLowEdge(histbins[i]+1);
    }

    std::cout << "  -->Proposed binning: " << histbins[i] << ", " << histmin[i] << ", " << histmax[i] << std::endl;
    std::cout << "  -->Using " << nBinsUsed << " bins: ";
    for (int iBin = 0; iBin < nBinsUsed; ++iBin) std::cout << xBins[iBin] << ", ";
    std::cout << xBins[nBinsUsed] << std::endl;


    for (int j = 0; j < numhists; j++)
    {
      currentDir->cd(); //This is because the "Draw" command only finds
      profiles_Data[numhists*i+j] = new TProfile(Form("profiles_Data_%s_%s",histnames[i],histnames[j]),"",nBinsUsed,xBins);
      profiles_Data[numhists*i+j]->Sumw2();
      profiles_Data[numhists*i+j]->SetLineWidth(2);
      profiles_MC[numhists*i+j] = new TProfile(Form("profiles_MC_%s_%s",histnames[i],histnames[j]),"",nBinsUsed,xBins);
      profiles_MC[numhists*i+j]->Sumw2();
      profiles_MC[numhists*i+j]->SetLineWidth(2);
    }

    delete[] xBins;

  }

float PurityFit::fit_specific( const TH1* h_, const TH1* sig_, const TH1* bkg_, 
        string name, // unique name of the result
        string outname , // output file name, where to save results
        map<string,float> *pars	 // to gather additional params
        )
{

    // 1) perform preliminary checks
    if ( h_ == NULL ) { cout<<"[PurityFit]::[fit_specific]::[ERROR] no target histogram"<<endl; return -1;}
    if ( sig_ == NULL ) { cout<<"[PurityFit]::[fit_specific]::[ERROR] no sig histogram"<<endl; return -1;}
    if ( bkg_ == NULL ) { cout<<"[PurityFit]::[fit_specific]::[ERROR] no bkg histogram"<<endl; return -1;}

    if (sig_ -> Integral() == 0 ) { cout<<"[PurityFit]::[fit_specific]::[ERROR] sig integrall is NULL"<<endl; return -2;}
    if (bkg_ -> Integral() == 0 ) { cout<<"[PurityFit]::[fit_specific]::[ERROR] bkg integrall is NULL"<<endl; return -2;}

    if (verbose_ >0) cout <<"[PurityFit]::[fit_specific]::[INFO] fitting "<<h_->GetName() << " " << sig_->GetName()<<" "<<bkg_->GetName()<<endl;
    // 1.5) Clone
    TH1 * sig = (TH1*)sig_ -> Clone(Form("%s_fitspecific_clone",sig_->GetName()));
    TH1 * bkg = (TH1*)bkg_ -> Clone(Form("%s_fitspecific_clone",bkg_->GetName()));
    TH1 * h = (TH1*)h_ -> Clone(Form("%s_fitspecific_clone",h_->GetName()));
    // 1.6) check no negative entries otherwise 0
    for(int i=1;i<=sig->GetNbinsX();++i)
        if(sig->GetBinContent(i) <0) sig->SetBinContent(i,0);
    for(int i=1;i<=bkg->GetNbinsX();++i)
        if(bkg->GetBinContent(i) <0) bkg->SetBinContent(i,0);

    // 2) scale templates: template normalization is meaningless
    sig->Sumw2();	
    bkg->Sumw2();

    if (bkg->Integral() <=0)
    {
    for(int i=1;i<=bkg->GetNbinsX();++i)
        bkg->SetBinContent(i,1.0);
    }

    // 2.5) fit background with exponential
    int   nbins = h->GetNbinsX();
    float xmin = h->GetBinLowEdge(1);
    float xmax = h->GetBinLowEdge(nbins+1);

    double prevChi2=-1;

    vector<TCanvas*> cbkgs;
    vector<TF1> expos;

    double cached_pars[50]; //cache
    for(int i=0;i< sizeof(cached_pars)/sizeof(double); ++i) cached_pars[i]=0;

    int poln=0;
    for( poln=0;poln<5; ++poln)
    {
        string formula = "TMath::Exp(-[0]*x) * ( ";
        for(int i=0;i<=poln;++i)
                {
                    if (i>0 ) formula += " + ";
                    formula += Form("[%d]",i+1);
                    if (i>0) formula += Form("*TMath::Power(x,%d)",i);
                }
        formula += ")";
        cout <<"Considering Formula"<<formula<<endl;
        TF1 expo(Form("expo%d",poln),formula.c_str(),xmin,xmax);
        expo.SetParameter(0,0.01);
        expo.SetParameter(1,0.05);

        expo.SetParLimits(0,TMath::Min(bkg->Integral()*.1,0.001), bkg->Integral()*10);
        expo.SetParLimits(1,1e-9, 1.);

        if (poln>0)for(int k=0;k < poln+1;++k) expo.SetParameter(k, cached_pars[k]);

        bkg->Fit( &expo ,"QN") ;
        bkg->Fit( &expo ,"QNM") ;
        expos.push_back(expo);

        //cache parameters
        for(int k=0;k < poln+1;++k) cached_pars[k]=expo.GetParameter(k);

        double chi2=expo.GetChisquare();
        double prob = 0 ; 
        if (poln > 0)
        {
            int n = 0;
            for(int i=1;i<=bkg->GetNbinsX();++i) if ( bkg->GetBinContent(i)> 0 ) ++n;
            int dof = n - (poln+2) + 1;
            double f = (prevChi2 - chi2) / ( chi2 /dof) ;
            prob= 1.- TMath::FDistI (f , 1, dof  ) ;
        }
        prevChi2 = chi2;

        cout<<"----------- BKG PARAMETERS ARE -------"<<endl;
        cout << "Prob = "<<prob<<endl;
        cout << "chi2 = "<<chi2<<endl;
        cout<<" 0 : "<< expo.GetParameter(0) <<endl;
        for(int i=0; i<=poln;++i)
            cout<<" "<< i+1 <<" : "<< expo.GetParameter(i+1) <<endl;
        cout<<"--------------------------------------"<<endl;
        TCanvas *cbkg=new TCanvas( (string(name)+"_bkgfit"+ Form("pol%d",poln)).c_str(),"Canvas"); cbkgs.push_back(cbkg);
        bkg->Clone("bkgForDraw")->Draw("P E"); // FIXME, memory leak
        expo.DrawClone("L SAME");

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

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;
}