source: huonglan/ANA2011/getrejS.C

{
  TFile *f = TFile::Open("res.root");
  TH1F *hl = (TH1F*)f->Get("sv0L");
  TH1F *hb = (TH1F*)f->Get("sv0B");

  TH1F *hipl = (TH1F*)f->Get("ip3dsv1L");
  TH1F *hipb = (TH1F*)f->Get("ip3dsv1B");

  int nb   = hl->GetNbinsX();
  int bi50 = hl->FindBin(5.85);
  cout << "Cut Bin " << bi50 << endl;
  cout << "epsilon_b = " << hb->Integral(bi50,nb+1)/hb->Integral(0,nb+1) << endl;
  cout << "R_light   = " << hl->Integral(0,nb+1)/hl->Integral(bi50,nb+1) << endl;
  double effsv[nb],eeffsv[nb];
  double R_lightsv[nb], eRsv[nb];
  double nblabelsv = hb->Integral(0,nb+1);
  double nllabelsv = hl->Integral(0,nb+1);

  for (int i = 1; i <= nb; i++) {
    effsv[i-1] = hb->Integral(i,nb+1)/nblabelsv;
    eeffsv[i-1] = 0;
    double intl = hl->Integral(i,nb+1);
    if (intl != 0) {
      R_lightsv[i-1] = nllabelsv/intl;
      eRsv[i-1] = R_lightsv[i-1]*sqrt((R_lightsv[i-1]-1)/nllabelsv);
    }
    else {
      R_lightsv[i-1] = 60000.;
      eRsv[i-1] = 0;
    }
  }

  TGraphErrors *heffbsv = new TGraphErrors(nb, effsv, R_lightsv, eeffsv, eRsv);

  int nbip = hipl->GetNbinsX();
  double nblabel = hipb->Integral(0,nbip+1);
  double nllabel = hipl->Integral(0,nbip+1);

  double eff[nbip],eeff[nbip];
  double R_light[nbip], eR[nbip];
  int bi50 = 0, bi57 = 0, bi60 = 0, bi70 = 0, bi01 = 0;
  for (int i = 1; i <= nbip; i++) {
    eff[i-1] = hipb->Integral(i,nbip+1)/nblabel;
    eeff[i-1] = 0;
    double intl = hipl->Integral(i,nbip+1);
    if (intl != 0) {
      R_light[i-1] = nllabel/intl;
      eR[i-1] = R_light[i-1]*sqrt((R_light[i-1]-1)/nllabel);
    }
    else {
      R_light[i-1] = 60000.;
      eR[i-1] = 0;
    }
    if (i > 1) {
      if (eff[i-1] <= 0.5 && eff[i-2] >= 0.5)   {bi50 = i-1;}
      if (eff[i-1] <= 0.57 && eff[i-2] >= 0.57) {bi57 = i-1;}
      if (eff[i-1] <= 0.6 && eff[i-2] >= 0.6)   {bi60 = i-1;}
      if (eff[i-1] <= 0.7 && eff[i-2] >= 0.7)   {bi70 = i-1;}
      if (1/R_light[i-1] <= 0.01 && 1/R_light[i-2] >= 0.01) {bi01 = i-1;} 
    }
  }

  double ru50  = (R_light[bi50]*(eff[bi50-1]-0.5)+R_light[bi50-1]*(0.5-eff[bi50]))/(eff[bi50-1]-eff[bi50]);
  double ru57  = (R_light[bi57]*(eff[bi57-1]-0.57)+R_light[bi57-1]*(0.57-eff[bi57]))/(eff[bi57-1]-eff[bi57]);
  double ru60  = (R_light[bi60]*(eff[bi60-1]-0.6)+R_light[bi60-1]*(0.6-eff[bi60]))/(eff[bi60-1]-eff[bi60]);
  double ru70  = (R_light[bi70]*(eff[bi70-1]-0.7)+R_light[bi70-1]*(0.7-eff[bi70]))/(eff[bi70-1]-eff[bi70]);
  double eru50 = (eR[bi50]*(eff[bi50-1]-0.5)+eR[bi50-1]*(0.5-eff[bi50]))/(eff[bi50-1]-eff[bi50]);
  double eru57 = (eR[bi57]*(eff[bi57-1]-0.57)+eR[bi57-1]*(0.57-eff[bi57]))/(eff[bi57-1]-eff[bi57]);
  double eru60 = (eR[bi60]*(eff[bi60-1]-0.6)+eR[bi60-1]*(0.6-eff[bi60]))/(eff[bi60-1]-eff[bi60]);
  double eru70 = (eR[bi70]*(eff[bi70-1]-0.7)+eR[bi70-1]*(0.7-eff[bi70]))/(eff[bi70-1]-eff[bi70]);
  cout <<" Ru(50%) = "<<ru50<<" +- "<<eru50<<" cut at "<<hipb->GetXaxis()->GetBinCenter(bi50)<<" "<<bi50<<endl;
  cout <<" Ru(57%) = "<<ru57<<" +- "<<eru57<<" cut at "<<hipb->GetXaxis()->GetBinCenter(bi57)<<" "<<bi57<<endl;
  cout <<" Ru(60%) = "<<ru60<<" +- "<<eru60<<" cut at "<<hipb->GetXaxis()->GetBinCenter(bi60)<<" "<<bi60<<endl;
  cout <<" Ru(70%) = "<<ru70<<" +- "<<eru70<<" cut at "<<hipb->GetXaxis()->GetBinCenter(bi70)<<" "<<bi70<<endl;
  //
  double epsb01  = (eff[bi01]*(1./R_light[bi01-1]-0.01)+eff[bi01-1]*(0.01-1./R_light[bi01]))/(1./R_light[bi01-1]-1./R_light[bi01]);
  double eepsb01 = (eeff[bi01]*(1./R_light[bi01-1]-0.01)+eeff[bi01-1]*(0.01-1./R_light[bi01]))/(1./R_light[bi01-1]-1./R_light[bi01]);
  cout <<" epsilon_b(Ru=100) = "<<epsb01<<" +- "<<eepsb01<<" for a cut at (center) "<<hb->GetXaxis()->GetBinCenter(bi01)<<endl;


  TGraphErrors *heffb = new TGraphErrors(nbip, eff, R_light, eeff, eR);
  heffb->SetMarkerStyle(20);
  heffb->SetMarkerColor(2);
  heffb->SetMarkerSize(0.5);
  heffb->Draw("LAP");

  heffbsv->Draw("LP");
}