ROOT Scripts

Merging histograms from signal and background

Here is an example script of how you might merge histograms of pass and total for both signal (and correctly weighted background) to make efficiency histograms combining both signal and background. Will update to work from simple flat ntuple rather than histograms once coded.

//
// As script to make efficiency histograms
// with signal (and background if in
// two separate files, given that background
// histograms already have correct weighting
// applied from analyzer).
//
// probably best to use root -b so it saves the
// results quickly if on a slow connection!
//

void merge(bool doBackground)
{

  // set up the tdr style
  gROOT->ProcessLine(".L ~/tdrStyle.c");
  setTDRStyle();

  // define a vector to hold in input 
  // ROOT files
  std::vector<TFile*> theRootFiles;

  // add the input root file(s)
  // note background only added if doBackground == true
  theRootFiles.push_back( new TFile("out.root") );
  if (doBackground == true) {
    theRootFiles.push_back( new TFile("out_background.root") );
  }

  // produce the efficiency histograms
  doHist("Eta", theRootFiles);
  doHist("Phi", theRootFiles);
  doHist("Pt", theRootFiles);
  doHist("EtaEB", theRootFiles);
  doHist("PhiEB", theRootFiles);
  doHist("PtEB", theRootFiles);
  doHist("EtaEE", theRootFiles);
  doHist("PhiEE", theRootFiles);
  doHist("PtEE", theRootFiles);

}

void doHist(TString &histName, std::vector<TFile*> &files)
{

  std::cout << "[doHist] doing " << histName << std::endl;

  // define placeholders for the signal
  // and background pass/fail hists
  TH1F *pass = new TH1F();
  TH1F *total = new TH1F();
  TH1F *pass2 = new TH1F();
  TH1F *total2 = new TH1F();

  // get pass/fail hists for signal 
  pass = (TH1F)files.at(0)->Get("hltAnalysis/" + histName + "Pass");
  total = (TH1F)files.at(0)->Get("hltAnalysis/" + histName + "Total");

  // if background file is present
  // do for background + add to signal
  if (files.size() == 2)
  {
    pass2 = (TH1F)files.at(1)->Get("hltAnalysis/" + histName + "Pass");
    total2 = (TH1F)files.at(1)->Get("hltAnalysis/" + histName + "Total");
    pass->Add(pass2);
    total->Add(total2);
  }

  // calculate the efficiency of signal (+background)
  double ratio = pass->Integral() / total->Integral();
  double effTotal = ((2.0 * ratio) - (ratio * ratio));
  double effTotalErr =  sqrt((effTotal * (1.0 - effTotal)) / total->Integral());
  std::cout << "eff is : " << effTotal << "+/-" << effTotalErr << std::endl;

  // plot the efficiency histogram for signal (+background)
  int numBins = pass->GetNbinsX();
  double xMax = pass->GetXaxis()->GetXmax();

  TH1F *eff = new TH1F("eff", histName, numBins, 0, xMax);

  for(int i = 1; i <= numBins; ++i)
  {
    double passCount = pass->GetBinContent(i); 
    double totalCount = total->GetBinContent(i);
    totalCount = ((totalCount == 0) ? 1 : totalCount);
    double ratio = passCount / totalCount;
    double effVal = ((2.0 * ratio) - (ratio * ratio));
    double errVal = sqrt((effVal * (1.0 - effVal)) / totalCount);

    if (!(passCount == totalCount && totalCount == 1))
    {
      eff->SetBinContent(i, effVal);
      eff->SetBinError(i, errVal);
    }
  }

  // save the efficiency histogram
  TCanvas *c1 = new TCanvas();
  c1->cd();
  eff->Draw();
  if (files.size() == 1) {
    c1->SaveAs(histName + "Eff.png");
  } else {
    c1->SaveAs(histName + "Eff_wbg.png"); 
  }

  // tidy up
  delete eff;
  delete c1;
  delete pass;
  delete total;
  delete pass2;
  delete total2;
}

Merging histograms with different weighting

This example is quite specific... Will add something more general when time allows.

struct qcdBin {
  TFile * file_;
  double sigma_;
  long eventsPass_;
  long eventsTotal_;
  double weight_;
};

void loop(){

// set up the tdr style
gROOT->ProcessLine(".L ~/tdrStyle.c");
setTDRStyle();

// set up the qcd bins
qcdBin pt30_50;
        pt30_50.file_ =  new TFile("30-50.root");
        pt30_50.sigma_ =  0.163 * pow(10, -3);
        pt30_50.eventsPass_ = 25;
        pt30_50.eventsTotal_ = 174876;
        pt30_50.weight_ = -1;  // weight is not initialized

qcdBin pt50_80;
        pt50_80.file_ =  new TFile("50-80.root");
        pt50_80.sigma_ =  0.0216 * pow(10, -3);
        pt50_80.eventsPass_ = 518;
        pt50_80.eventsTotal_ = 1885806;
        pt50_80.weight_ = -1;  // weight is not initialized

qcdBin pt80_120;
        pt80_120.file_ =  new TFile("80-120.root");
        pt80_120.sigma_ =  0.00308 * pow(10, -3);
        pt80_120.eventsPass_ = 186;
        pt80_120.eventsTotal_ = 691098;
        pt80_120.weight_ = -1;  // weight is not initialized

qcdBin pt120_170;
        pt120_170.file_ =  new TFile("120-170.root");
        pt120_170.sigma_ =  0.000494 * pow(10, -3);
        pt120_170.eventsPass_ = 101;
        pt120_170.eventsTotal_ = 357434;
        pt120_170.weight_ = -1;  // weight is not initialized

qcdBin signal;
        pt120_170.file_ =  new TFile("signal.root");
        pt120_170.sigma_ =  -1;          // signal events 
        pt120_170.eventsPass_ = -1;     // don't need weighting
        pt120_170.eventsTotal_ = -1;
        pt120_170.weight_ = 1;           // weight is 1              



TH1F *totalWeighted2 = new TH1F("totalWeighted2", "totalWeighted method2", 50, 0., 100.);
        totalWeighted2->GetYaxis()->SetRangeUser(0, 3000);
        totalWeighted2->Sumw2();

doWeightMethod2(pt30_50, *totalWeighted2);
doWeightMethod2(pt50_80, *totalWeighted2);
doWeightMethod2(pt80_120, *totalWeighted2);
doWeightMethod2(pt120_170, *totalWeighted2);

TCanvas *c2 = new TCanvas();
//c2->SetLogy();
c2->cd();
totalWeighted2->Draw();

}

// use Sumw2
void doWeightMethod2(qcdBin &thisBin, TH1F &theHist)
{

        double eff;
        double weight;
        double intSigma = 10 * pow(10, 12);

        TH1F *tmp;
        tmp = (TH1F*)thisBin.file_->Get("ProbeBuilder__SimpleEcalProbe/_tagProbeInvMass");

        if (thisBin.weight_ == -1) {
          eff = (static_cast<double>(thisBin.eventsPass_) / static_cast<double>(thisBin.eventsTotal_));
          weight = (eff * thisBin.sigma_ * intSigma) / static_cast<double>(thisBin.eventsPass_);
        } else weight = thisBin.weight_;

        //std::cout << "bin with sigma " << thisBin.sigma_ << std::endl;
        //std::cout << "   has weight " << weight << std::endl;

        for (Int_t bin = 1; bin < tmp->GetNbinsX(); ++bin)
        {
                double tmpContent =  tmp->GetBinContent(bin);
                for (unsigned int i = 0; i < tmpContent; ++i)
                {
                        theHist.Fill(tmp->GetBinCenter(bin), weight);
                }
        }

}

-- DaveEvans - 02 Jul 2007