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

/*============================================================================*/
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 QAcentrality(const Char_t *fdata)
{
  style();

  TFile *fin = TFile::Open(fdata);
  TList *lin = (TList *)fin->Get("clist");
  lin->ls();
  TH1 *hin = (TH1 *)lin->FindObject("EvCentrDist");
  Float_t sum = 1.2 * hin->Integral(hin->FindBin(0.1), hin->FindBin(79.9));
  hin->Scale(1. / sum);
  SetHistoStyle(hin, 20, kRed+1);
  TCanvas *c = new TCanvas("cQAcentrality", "cQAcentrality", 800, 800);
  TH1 * hfr = c->DrawFrame(0., 0.005, 100., 0.015);
  hfr->SetTitle(";centrality percentile;events");
  hin->Draw("same");
  c->SaveAs(canvasPrefix+"centrality.pdf");

  TH2 *hinv0 = (TH2 *)lin->FindObject("V0");
  TCanvas *cv0 = new TCanvas("cQAcentralityV0", "cQAcentralityV0", 800, 800);
  cv0->SetLogx();
  cv0->SetLogz();
  //  TH1 * hfrv0 = cv0->DrawFrame(100., -0.5, 50000., 10.5);
  // DrawBinLabelsY(hfrv0, kTRUE);
  // hfrv0->SetTitle(";V0 signal;");
  //hinv0->Draw("same,col");
  hinv0->Draw("col");
  cv0->SaveAs(canvasPrefix+"centralityV0.pdf");
}

void fill_data(double e_min, double e_max, int mask)
{
	const char name_pattern[] = "danss_report_v4n-sect%d-calc.root";
	const char pos[3][20] = {"hUp_%d", "hMid_%d", "hDown_%d"};
	char fname[1024];
	char cpos[32];
	int i, j;
	TFile *f;
	TH1* h;
	double val, err;
	
	for (i=0; i<3; i++) {
		sprintf(fname, name_pattern, 3-i);
		f = new TFile(fname);
		if (!f->IsOpen()) return;
		for (j=0; j<3; j++) {
			sprintf(cpos, pos[j], mask);
			h = (TH1*) f->Get(cpos);
			if (!h) {
				printf("Something is wrong: %d %d %s.\n", i, j, cpos);
				return;
			}
			val = h->IntegralAndError(h->FindBin(e_min), h->FindBin(e_max), err);
			DataArray.cnt[3*i+j] = val;
			DataArray.ecnt[3*i+j] = err;
			printf("%d %d %f +- %f\n", i, j, val, err);
		}
		f->Close();
	}
}

Double_t fitlang( char* hs, double klow = 18, double khigh = 40 ) {

  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 = h->FindBin(klow);
  int ib9 = h->FindBin(khigh);

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

TH1 *computeEffVsBEff(TH1 *effVsCut, TH1 *effVsCutB)
{
	TH1 *h = new TH1F(Form("%s_effVsBEff", effVsCut->GetName()), "effVsBEff",
	                  100, 0, 1);
	h->SetMaximum(1.5);
	h->SetMinimum(1e-3);

	unsigned int n = effVsCut->GetNbinsX();
	for(unsigned int bin = 1; bin <= n; bin++) {
		double eff = effVsCut->GetBinContent(bin);
		double error = effVsCut->GetBinError(bin);
		double effB = effVsCutB->GetBinContent(bin);

		h->SetBinContent(h->FindBin(effB), eff);
		h->SetBinError(h->FindBin(effB), error);
 		// FIXME: The error in effB is not propagated
	}

	return h;
}

TEST_F(EvalHistMethods, CreateHistogram1D) {
  evaluator->SetNormalizationBuffer(norm);
  evaluator->SetParameterBuffer(params);

  TH1* hist = evaluator->CreateHistogram();

  EXPECT_EQ(2, hist->GetNbinsX());
  ASSERT_FLOAT_EQ(1.0, hist->Integral("width"));
  ASSERT_FLOAT_EQ(1.6, hist->GetBinContent(hist->FindBin(0.25)));
  ASSERT_FLOAT_EQ(0.4, hist->GetBinContent(hist->FindBin(0.75)));

  delete hist;
}

void DrawReport(const char* psname, TObjArray* harr) 
{
  gStyle->SetOptFit(1);
  if (!harr) harr = &histoArr;
  TCanvas* cnv = new TCanvas("cl","cl",900,600);
  //
  TString psnm1 = psname;
  if (psnm1.IsNull()) psnm1 = "clusters.ps";
  TString psnm0 = psnm1.Data(); 
  psnm0 += "[";
  TString psnm2 = psnm1.Data(); 
  psnm2 += "]";
  cnv->Print(psnm0.Data());
  //
  TH1* clall = GetHistoClSize(0,kNPixAll,harr);
  clall->SetLineColor(kRed);
  clall->Draw();
  TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
  clSpl->SetLineColor(kBlue);
  clSpl->Draw("sames");
  gPad->Modified();
  gPad->Update();
  SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
  SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
  gPad->Modified();
  gPad->Update();
  gPad->SetLogy(1);
  //
  cnv->cd();
  cnv->Print(psnm1.Data());
  //
  // plot cluster sized from 1 to 10
  for (int i=1;i<=10;i++) {
    if (clall->GetBinContent(clall->FindBin(i))<100) continue;
    DrawNP(i,harr,cnv);
    cnv->Print(psnm1.Data());
  }
  cnv->Print(psnm2.Data());
}

//.........这里部分代码省略.........
      if(fMCTauMultiplicity == kMCTauTwoTaus && ntaus != 1) return true;
      if(fMCTauMultiplicity == kMCTauMoreThanTwoTaus && ntaus < 2) return true;
    }
    else {
      if(fMCTauMultiplicity == kMCTauZeroTaus && ntaus != 0) return true;
      if(fMCTauMultiplicity == kMCTauOneTau && ntaus != 1) return true;
      if(fMCTauMultiplicity == kMCTauTwoTaus && ntaus != 2) return true;
      if(fMCTauMultiplicity == kMCTauMoreThanTwoTaus && ntaus < 3) return true;
    }
   
    if(genTaus.size() > 0)
      theGenTau = genTaus[0];
  }
  cTauMCSelection.increment();

  bool originalMuonIsWMu = false;
  EmbeddingMuonCollection::Muon embeddingMuon;
  if(fIsEmbedded) {
    if(fMuons.size() != 1)
      throw std::runtime_error("Embedding muon collection size is not 1");
    embeddingMuon = fMuons.get(0);
    if(embeddingMuon.p4().Pt() < 41) return true;
    //std::cout << "Muon pt " << muon.p4().Pt() << std::endl;

    originalMuonIsWMu = std::abs(embeddingMuon.pdgId()) == 13 && std::abs(embeddingMuon.motherPdgId()) == 24 && std::abs(embeddingMuon.grandMotherPdgId()) == 6;
    if(!originalMuonIsWMu) return true;
  }

  cOnlyWMu.increment();

  // Per-event correction for W->tau->mu
  if(fIsEmbedded && fEmbeddingWTauMuWeights) {
    embeddingMuon.ensureValidity();
    int bin = fEmbeddingWTauMuWeights->FindBin(embeddingMuon.p4().Pt());
    //std::cout << embeddingMuon.p4().Pt() << " " << bin << std::endl;
    weight *= fEmbeddingWTauMuWeights->GetBinContent(bin);
    fEventCounter.setWeight(weight);
  }
  cWTauMuWeight.increment();

  // Muon trigger efficiency weighting
  if(fIsEmbedded && (fEmbeddingNormalizationMode == kIDTriggerEfficiency || fEmbeddingNormalizationMode == kIDTriggerEfficiencyTauBR)) {
    weight *= 1/embeddingMuon.triggerEfficiency();
    fEventCounter.setWeight(weight);
  }
  cTriggerEffWeight.increment();

  // Muon ID efficiency weighting
  if(fIsEmbedded) {
    weight *= 1/embeddingMuon.idEfficiency();
    fEventCounter.setWeight(weight);
  }
  cIdEffWeight.increment();

  // Weighting by arbitrary histogram, given in constructor argument
  if(fIsEmbedded && fEmbeddingMuonWeights) {
    embeddingMuon.ensureValidity();
    int bin = fEmbeddingMuonWeights->FindFixBin(embeddingMuon.p4().Pt());
    weight *= fEmbeddingMuonWeights->GetBinContent(bin);
    fEventCounter.setWeight(weight);
  }
  cMuonWeight.increment();

  // Jet selection
  /*
  size_t njets = fJets.size();

//.........这里部分代码省略.........
  cout << " Find cuts for all the centralities " << endl << endl;

  const int nval=9;//number of variables to fill
  char* centname[3] = {"5pStep","10pStep","20pStep"};
  char* varname[nval] = {"npart","ncoll","b","standard_ecc","reactionplane_ecc","participant_ecc","R_Ollitraut","R_Geo", "R_arith"};
  double vup[nval] = {499.5,2999.5,19.995,1,1,1,4,4,4};
  double vlo[nval] = {-0.5,-0.5,-0.005,-1,-1,-1,0,0,0};
  int vNb[nval] = {500,3000,2000,200,200,200,400,400,400};

  //initialize the histograms which are used to fill the distribution of variables for each centrality
  TH1* hvar[3][nval][20];
  for(int i=0; i<3; i++){
    int NC = 0;
    if(i==0) NC = 20;
    else if(i==1) NC = 10;
    else if(i==2) NC = 5;
    for(int ivar=0; ivar<nval; ivar++){
      for(int icen=0; icen<NC; icen++){
	sprintf(name,"hvar_%s_%s_cent%d",centname[i],varname[ivar],icen);
	hvar[i][ivar][icen] = new TH1D(name,name,vNb[ivar],vlo[ivar],vup[ivar]);
      }
    }
  }

  double qbbcsum;
  for(int i=0; i<n; i++){
    if(i%1000000==0) cout << i << endl;
    ch->GetEntry(i);

    if(!(nHitBbc_n>=2&&nHitBbc_s>=2)) continue;//BBC trigger condition

    qbbcsum = (qBBC_n+qBBC_s)/100.;

    int centbin5  = FindBin(20,bbcq5,qbbcsum);
    int centbin10 = FindBin(10,bbcq10,qbbcsum);
    int centbin20 = FindBin(5,bbcq20,qbbcsum);

    if(centbin5==-1) continue;
    if(centbin10==-1) continue;
    if(centbin20==-1) continue;

    //find the weight according to the corresponding average efficiency.
    double weight = hbbcqeff->GetBinContent(hbbcqeff->FindBin(qbbcsum));
    //5 percent step
    //
    hvar[0][0][centbin5]->Fill(npart,weight);
    hvar[0][1][centbin5]->Fill(ncoll,weight);
    hvar[0][2][centbin5]->Fill(b,weight);
    hvar[0][3][centbin5]->Fill(ecc_std,weight);
    hvar[0][4][centbin5]->Fill(ecc_rp,weight);
    hvar[0][5][centbin5]->Fill(ecc_part,weight);
    hvar[0][6][centbin5]->Fill(r_oll,weight);
    hvar[0][7][centbin5]->Fill(r_geo,weight);
    hvar[0][8][centbin5]->Fill(r_arith,weight);

    //10 percent step
    //
    hvar[1][0][centbin10]->Fill(npart,weight);
    hvar[1][1][centbin10]->Fill(ncoll,weight);
    hvar[1][2][centbin10]->Fill(b,weight);
    hvar[1][3][centbin10]->Fill(ecc_std,weight);
    hvar[1][4][centbin10]->Fill(ecc_rp,weight);
    hvar[1][5][centbin10]->Fill(ecc_part,weight);
    hvar[1][6][centbin10]->Fill(r_oll,weight);
    hvar[1][7][centbin10]->Fill(r_geo,weight);
    hvar[1][8][centbin10]->Fill(r_arith,weight);

/*============================================================================*/
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);
}

//.........这里部分代码省略.........
    // !!!!!!!!!!!!!!!!!!
    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; // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      else            grxE[count] = 0.50; // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    }
    else
      //      grxE[count] = (1.30-1.10)/2.0; // !!!!!!!!!!!!!!!!!!!!!!!!!!
      grxE[count] = 0.05;
    gry[count] = value;

    Double_t tmp1 = h1->Integral(h1->FindBin(fmini), h1->FindBin(fmaxi));
    Double_t tmp2 = h3_p->Integral(h3_p->FindBin(fmini), h3_p->FindBin(fmaxi));
    Double_t tmp_sg = tmp1 - tmp2;
    Double_t tmp_bg = tmp2;

    // if ((tmp_sg <= -tmp_bg) || (tmp_bg < 33.0)) {
    //   gry3[count] = 0.0;
    //   gry4[count] = 0.0;
    // }
    // else {
    gry3[count] = tmp_sg/tmp_bg;
    gry4[count] = tmp_sg/TMath::Sqrt(tmp_sg + tmp_bg);
    // }

void Fit::nll(int& ndim, double* gout, double& result, double* par, int flags) {
  size_t n = Fit::signals.size();
  result = 0;

  // make summed fit histogram
  TH1* hfit = (TH1*) signals[0].histogram->Clone("hfit");
  hfit->Reset();
  for (size_t i=0; i<n; i++) {
    TH1* h = Fit::signals.at(i).histogram;
    h->Scale(1.0/h->Integral());
    hfit->Add(h, Fit::norms[i] * par[i]);
  }

  // loop over bins
  if (hfit->IsA() == TH2F::Class()) {
    TH2F* h2 = dynamic_cast<TH2F*>(hfit);
    for (int i=1; i<h2->GetNbinsX(); i++) {
      if (i < h2->GetXaxis()->FindBin(Fit::r_range.min) || i > h2->GetXaxis()->FindBin(Fit::r_range.max)) {
        continue;
      }
      for (int j=1; j<h2->GetNbinsY(); j++) {
        if (j < h2->GetYaxis()->FindBin(Fit::e_range.min) || j > h2->GetYaxis()->FindBin(Fit::e_range.max)) {
          continue;
        }
        double nexp = h2->GetBinContent(i, j);
        double nobs = dynamic_cast<TH2F*>(Fit::data)->GetBinContent(i, j);
        result += (nexp - nobs * TMath::Log(TMath::Max(1e-12, nexp)));
        //std::cout << "- " << nexp << " " << nobs << " " << TMath::Log(nexp) << " // " << result << std::endl;
      }
    }
  }
  else {
    for (int i=1; i<hfit->GetNbinsX(); i++) {
      if (i < hfit->FindBin(Fit::e_range.min) || i > hfit->FindBin(Fit::e_range.max)) {
        continue;
      }
      double nexp = hfit->GetBinContent(i);
      double nobs = Fit::data->GetBinContent(i);
      result += (nexp - nobs * TMath::Log(TMath::Max(1e-12, nexp)));
      //std::cout << "- " << nexp << " " << nobs << " " << TMath::Log(nexp) << " // " << result << std::endl;
    }
  }

  // constraints
  for (size_t i=0; i<n; i++) {
    if (Fit::signals.at(i).constraint > 0) {
      result += 0.5 * TMath::Power((par[i]-1.0)
                / Fit::signals.at(i).constraint, 2);
    }
  }
  delete hfit;

#ifdef DEBUG
  // print parameters at each iteration
  std::cout << "+ ";
  for (size_t i=0; i<n; i++) {
    std::cout << par[i] * Fit::signals.at(i).rate << " (" << par[i] << ") \t";
  }
  std::cout << result << std::endl;
#endif
}

// The probability to find a value of the test statistic equal or
// worse than the one obtained with the observed yields
void ToyExperiments::printGlobalPValueOfObservation(unsigned int nExperiments) const {
  std::cout << "Performing " << nExperiments << " toy experiments to compute global p-value of observation ...  " << std::flush;

  // The test statistic
  TH1* hTestStat = new TH1D("hTestStat","",10000*Parameters::nBins(),0,1000);
  std::vector<unsigned int> obs(Parameters::nBins(),0);

  // Minimal value (in standard deviations) allowed for
  //correlated fluctuation
  const double minCorr = findMinValidRandomNumberForCorrelatedUncertainties();

  // Throw mean values
  for(unsigned int p = 0; p < nExperiments; ++p) {
    // Throw one (normalized) random number for correlated
    // uncertainties that is valid in all bins
    double rCorr = rand_->Gaus(0.,1.);
    while( rCorr <= minCorr ) {
      rCorr = rand_->Gaus(0.,1.);
    }

    // Loop over all bins and get individual predictions
    for(unsigned int bin = 0; bin < Parameters::nBins(); ++bin) {
      double prediction = -1.;
      bool negativePrediction = true;
      while( negativePrediction ) {
	// Throw one (normalized) random number for uncorrelated
	// uncertainties that is valid in this bin only
	double rUncorr = rand_->Gaus(0.,1.);
	
	// Scale the normalized random numbers by the uncertainties' size
	// to obtain variation of yield
	double uncorrVar = rUncorr * uncorrelatedUncerts_.at(bin);
	double corrVar   = rCorr   * correlatedUncerts_.at(bin);
	
	// Add variations to yield
	prediction = meanPredictions_.at(bin) + uncorrVar + corrVar;
	
	// Check if prediction is positive
	if( prediction >= 0. ) {
	  negativePrediction = false;
	}
      }

      // Throw predicted yields from Poisson with
      // this mean
      obs.at(bin) = rand_->Poisson(prediction);
    } // End of loop over bins
    
    hTestStat->Fill(testStatistic(obs));
  }
  std::cout << "ok" << std::endl;
  
  std::cout << "\n\n----- Global p-value of observation -----" << std::endl;
  double tSObs = testStatistic(observedYields_);
  int binTSObs = hTestStat->FindBin(tSObs);
  int binTSTot = hTestStat->GetNbinsX()+1; // Include overflows
  double gpVal = hTestStat->Integral(binTSObs,binTSTot)/hTestStat->Integral(1,binTSTot);
  if( binTSObs == binTSTot ) {
    std::cerr << "  ERROR: observed value of test statistic out of histogram range" << std::endl;
  } else if( hTestStat->GetBinContent(binTSTot) > 0 ) {
    std::cerr << "  WARNING: distribution of test statistic has overflows (N = " << hTestStat->GetBinContent(binTSTot) << ")" << std::endl;
  }
  std::cout << "         observed value of test statistic : " << tSObs << std::endl;
  std::cout << "  global p-value (including correlations) : " << gpVal  << " (" << TMath::NormQuantile(1.-gpVal) << "sig)" << std::endl;

  hTestStat->Draw();
  gPad->SetLogy();
  gPad->SaveAs("TestStatistic.pdf");

  delete hTestStat;
}