testUnfold7b.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_unfold
/// \notebook
/// Test program for the classes TUnfoldDensity and TUnfoldBinning
///
/// A toy test of the TUnfold package
///
/// This example is documented in conference proceedings:
///
///   arXiv:1611.01927
/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
///
/// This is an example of unfolding a two-dimensional distribution
/// also using an auxiliary measurement to constrain some background
///
/// The example comprises several macros
///  - testUnfold7a.C   create root files with TTree objects for
///                      signal, background and data
///            - write files  testUnfold7_signal.root
///                           testUnfold7_background.root
///                           testUnfold7_data.root
///
///  - testUnfold7b.C   loop over trees and fill histograms based on the
///                      TUnfoldBinning objects
///            - read  testUnfold7binning.xml
///                    testUnfold7_signal.root
///                    testUnfold7_background.root
///                    testUnfold7_data.root
///
///            - write testUnfold7_histograms.root
///
///  - testUnfold7c.C   run the unfolding
///            - read  testUnfold7_histograms.root
///            - write testUnfold7_result.root
///                    testUnfold7_result.ps
///
/// \macro_output
/// \macro_code
///
///  **Version 17.6, in parallel to changes in TUnfold**
///
///  This file is part of TUnfold.
///
///  TUnfold is free software: you can redistribute it and/or modify
///  it under the terms of the GNU General Public License as published by
///  the Free Software Foundation, either version 3 of the License, or
///  (at your option) any later version.
///
///  TUnfold is distributed in the hope that it will be useful,
///  but WITHOUT ANY WARRANTY; without even the implied warranty of
///  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
///  GNU General Public License for more details.
///
///  You should have received a copy of the GNU General Public License
///  along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.
///
/// \author Stefan Schmitt DESY, 14.10.2008
 
 
/* below is the content of the file testUnfold7binning.xml,
   which is required as input to run this example.
 
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
<TUnfoldBinning>
<BinningNode name="fine" firstbin="1" factor="1">
  <Axis name="pt" lowEdge="0.">
   <Bin repeat="20" width="1."/>
   <Bin repeat="12" width="2.5"/>
   <Bin location="overflow" width="10"/>
  </Axis>
</BinningNode>
<BinningNode name="coarse" firstbin="1" factor="1">
  <Axis name="pt" lowEdge="0.">
   <Bin repeat="10" width="2"/>
   <Bin repeat="6" width="5"/>
   <Bin location="overflow" width="10"/>
  </Axis>
</BinningNode>
</TUnfoldBinning>
 
 */
 
#include <iostream>
#include <cmath>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TDOMParser.h>
#include <TXMLDocument.h>
#include "TUnfoldBinningXML.h"
 
using std::cout;
 
 
void testUnfold7b()
{
  // switch on histogram errors
  TH1::SetDefaultSumw2();
 
  //=======================================================
  // Step 1: open file to save histograms and binning schemes
 
  TFile *outputFile=new TFile("testUnfold7_histograms.root","recreate");
 
  //=======================================================
  // Step 2: read binning from XML
  //         and save them to output file
 
  TUnfoldBinning *fineBinningRoot,*coarseBinningRoot;
 
  // read binning schemes in XML format
 
  TDOMParser parser;
  TString dir = gSystem->UnixPathName(gSystem->GetDirName(__FILE__));
  Int_t error = parser.ParseFile(dir+"/testUnfold7binning.xml");
  if(error) {
     cout<<"error="<<error<<" from TDOMParser\n";
     cout<<"==============================================================\n";
     cout<<"Maybe the file testUnfold7binning.xml is missing?\n";
     cout<<"The content of the file is included in the comments section\n";
     cout<<"of this macro \"testUnfold7b.C\"\n";
     cout<<"==============================================================\n";
  }
  TXMLDocument const *XMLdocument=parser.GetXMLDocument();
  fineBinningRoot=
     TUnfoldBinningXML::ImportXML(XMLdocument,"fine");
  coarseBinningRoot=
     TUnfoldBinningXML::ImportXML(XMLdocument,"coarse");
 
  // write binning schemes to output file
  fineBinningRoot->Write();
  coarseBinningRoot->Write();
 
  if(fineBinningRoot) {
     fineBinningRoot->PrintStream(cout);
  } else {
     cout<<"could not read 'detector' binning\n";
  }
  if(coarseBinningRoot) {
     coarseBinningRoot->PrintStream(cout);
  } else {
     cout<<"could not read 'generator' binning\n";
  }
 
  TUnfoldBinning const *fineBinning=fineBinningRoot;//->FindNode("ptfine");
  TUnfoldBinning const *coarseBinning=coarseBinningRoot;//->FindNode("ptcoarse");
 
 
  //=======================================================
  // Step 3: book and fill data histograms
 
  Float_t ptRec[3],ptGen[3],weight;
  Int_t isTriggered,isSignal;
 
  outputFile->cd();
 
  TH1 *histDataRecF=fineBinning->CreateHistogram("histDataRecF");
  TH1 *histDataRecC=coarseBinning->CreateHistogram("histDataRecC");
  TH1 *histDataBgrF=fineBinning->CreateHistogram("histDataBgrF");
  TH1 *histDataBgrC=coarseBinning->CreateHistogram("histDataBgrC");
  TH1 *histDataGen=coarseBinning->CreateHistogram("histDataGen");
 
  for (auto histo : {histDataRecF, histDataRecC, histDataBgrF, histDataBgrC, histDataGen}) {
     histo->SetDirectory(outputFile);
  }
 
  TFile *dataFile=new TFile("testUnfold7_data.root");
  TTree *dataTree=(TTree *) dataFile->Get("data");
 
  if(!dataTree) {
     cout<<"could not read 'data' tree\n";
  }
 
  dataTree->ResetBranchAddresses();
  dataTree->SetBranchAddress("ptrec",ptRec);
  //dataTree->SetBranchAddress("discr",&discr);
  // for real data, only the triggered events are available
  dataTree->SetBranchAddress("istriggered",&isTriggered);
  // data truth parameters
  dataTree->SetBranchAddress("ptgen",ptGen);
  dataTree->SetBranchAddress("issignal",&isSignal);
  dataTree->SetBranchStatus("*",1);
 
  cout<<"loop over data events\n";
 
#define VAR_REC (ptRec[2])
#define VAR_GEN (ptGen[2])
 
  for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
     if(dataTree->GetEntry(ievent)<=0) break;
     // fill histogram with reconstructed quantities
     if(isTriggered) {
        int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
        int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
        histDataRecF->Fill(binF);
        histDataRecC->Fill(binC);
        if(!isSignal) {
           histDataBgrF->Fill(binF);
           histDataBgrC->Fill(binC);
        }
     }
     // fill histogram with data truth parameters
     if(isSignal) {
        int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
        histDataGen->Fill(binGen);
     }
  }
 
  delete dataTree;
  delete dataFile;
 
  //=======================================================
  // Step 4: book and fill histogram of migrations
  //         it receives events from both signal MC and background MC
 
  TH2 *histMcsigGenRecF=TUnfoldBinning::CreateHistogramOfMigrations
     (coarseBinning,fineBinning,"histMcsigGenRecF");
  TH2 *histMcsigGenRecC=TUnfoldBinning::CreateHistogramOfMigrations
     (coarseBinning,coarseBinning,"histMcsigGenRecC");
  TH1 *histMcsigRecF=fineBinning->CreateHistogram("histMcsigRecF");
  TH1 *histMcsigRecC=coarseBinning->CreateHistogram("histMcsigRecC");
  TH1 *histMcsigGen=coarseBinning->CreateHistogram("histMcsigGen");
  for (auto histo :
       std::initializer_list<TH1 *>{histMcsigGenRecF, histMcsigGenRecC, histMcsigRecF, histMcsigRecC, histMcsigGen}) {
     histo->SetDirectory(outputFile);
  }
 
  TFile *signalFile=new TFile("testUnfold7_signal.root");
  TTree *signalTree=(TTree *) signalFile->Get("signal");
 
  if(!signalTree) {
     cout<<"could not read 'signal' tree\n";
  }
 
  signalTree->ResetBranchAddresses();
  signalTree->SetBranchAddress("ptrec",&ptRec);
  //signalTree->SetBranchAddress("discr",&discr);
  signalTree->SetBranchAddress("istriggered",&isTriggered);
  signalTree->SetBranchAddress("ptgen",&ptGen);
  signalTree->SetBranchAddress("weight",&weight);
  signalTree->SetBranchStatus("*",1);
 
  cout<<"loop over MC signal events\n";
 
  for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
     if(signalTree->GetEntry(ievent)<=0) break;
 
     int binC=0,binF=0;
     if(isTriggered) {
        binF=fineBinning->GetGlobalBinNumber(VAR_REC);
        binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
     }
     int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
     histMcsigGenRecF->Fill(binGen,binF,weight);
     histMcsigGenRecC->Fill(binGen,binC,weight);
     histMcsigRecF->Fill(binF,weight);
     histMcsigRecC->Fill(binC,weight);
     histMcsigGen->Fill(binGen,weight);
  }
 
  delete signalTree;
  delete signalFile;
 
  TH1 *histMcbgrRecF=fineBinning->CreateHistogram("histMcbgrRecF");
  TH1 *histMcbgrRecC=coarseBinning->CreateHistogram("histMcbgrRecC");
  histMcbgrRecF->SetDirectory(outputFile);
  histMcbgrRecC->SetDirectory(outputFile);
 
  TFile *bgrFile=new TFile("testUnfold7_background.root");
  TTree *bgrTree=(TTree *) bgrFile->Get("background");
 
  if(!bgrTree) {
     cout<<"could not read 'background' tree\n";
  }
 
  bgrTree->ResetBranchAddresses();
  bgrTree->SetBranchAddress("ptrec",&ptRec);
  //bgrTree->SetBranchAddress("discr",&discr);
  bgrTree->SetBranchAddress("istriggered",&isTriggered);
  bgrTree->SetBranchAddress("weight",&weight);
  bgrTree->SetBranchStatus("*",1);
 
  cout<<"loop over MC background events\n";
 
  for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
     if(bgrTree->GetEntry(ievent)<=0) break;
 
     // here, for background only reconstructed quantities are known
     // and only the reconstructed events are relevant
     if(isTriggered) {
        int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
        int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
        histMcbgrRecF->Fill(binF,weight);
        histMcbgrRecC->Fill(binC,weight);
     }
  }
 
  delete bgrTree;
  delete bgrFile;
 
  outputFile->Write();
  delete outputFile;
 
}