ROOT Scripts

Merging histograms from signal and background

Here is an example script of how you might merge histograms and make a simple calculation from signal and background. Clearly and can be expanded and made more general but it's just something quick I put together last Friday. Will add better script once written...

{

TFile f1("out.root");
TFile f2("out_background.root");

TH1F *pass;
TH1F *total;
TH1F *pass2;
TH1F *total2;

pass = (TH1F*)f1.Get("hltAnalysis/EtaPass");
total = (TH1F*)f1.Get("hltAnalysis/EtaTotal");
pass2 = (TH1F*)f2.Get("hltAnalysis/EtaPass");
total2 = (TH1F*)f2.Get("hltAnalysis/EtaTotal");

double ratio = pass->Integral() / total->Integral();
double effTotal = ((2.0 * ratio) - (ratio * ratio));
double effTotalErr =  sqrt((effTotal * (1.0 - effTotal)) / total->Integral());

std::cout << "no background" << std::endl;
std::cout << "eff is : " << effTotal << "+/-" << effTotalErr << std::endl;


pass->Add(pass2);
total->Add(total2);

TH1F *eff = new TH1F("eff", "eff", pass->GetNbinsX(), 0, pass->GetXaxis()->GetXmax());

  int numBins = pass->GetNbinsX();
  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);
      eff->SetBinContent(i, effVal);
      eff->SetBinError(i, errVal);
    }

eff->Draw();

ratio = pass->Integral() / total->Integral();
effTotal = ((2.0 * ratio) - (ratio * ratio));
effTotalErr =  sqrt((effTotal * (1.0 - effTotal)) / total->Integral());

std::cout << "background" << std::endl;
std::cout << "eff is : " << effTotal << "+/-" << effTotalErr << std::endl;

}

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

Edit | Attach | Watch | Print version | History: r3 < r2 < r1 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r2 - 2007-07-02 - DaveEvans
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2021 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback