/// \file /// \ingroup tutorial_unfold /// \notebook /// Test program for the class TUnfoldSys. /// /// Simple toy tests of the TUnfold package /// /// Pseudo data (5000 events) are unfolded into three components /// The unfolding is performed once without and once with area constraint /// /// Ideally, the pulls may show that the result is biased if no constraint /// is applied. This is expected because the true data errors are not known, /// and instead the sqrt(data) errors are used. /// /// \macro_output /// \macro_code /// /// **Version 17.6, in parallel to changes in TUnfold** /// /// #### History: /// - Version 17.5, in parallel to changes in TUnfold /// - Version 17.4, in parallel to changes in TUnfold /// - Version 17.3, in parallel to changes in TUnfold /// - Version 17.2, in parallel to changes in TUnfold /// - Version 17.1, in parallel to changes in TUnfold /// - Version 16.1, parallel to changes in TUnfold /// - Version 16.0, parallel to changes in TUnfold /// - Version 15, use L-curve scan to scan the average correlation /// /// 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 . /// /// \author Stefan Schmitt DESY, 14.10.2008 #include #include #include #include #include #include #include #include #include #include "TUnfoldDensity.h" using namespace std; TRandom *rnd=nullptr; Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) { // choose an integer random number in the range [0,nmax] // (the generator level bin) // depending on the probabilities // probability[0],probability[1],...probability[nmax-1] Double_t f=rnd->Rndm(); Int_t r=0; while((r=probability[r])) { f -= probability[r]; r++; } return r; } Double_t GenerateRecEvent(const Double_t *shapeParm) { // return a coordinate (the reconstructed variable) // depending on shapeParm[] // shapeParm[0]: fraction of events with Gaussian distribution // shapeParm[1]: mean of Gaussian // shapeParm[2]: width of Gaussian // (1-shapeParm[0]): fraction of events with flat distribution // shapeParm[3]: minimum of flat component // shapeParm[4]: maximum of flat component Double_t f=rnd->Rndm(); Double_t r; if(fGaus(shapeParm[1],shapeParm[2]); } else { r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3]; } return r; } void testUnfold4(bool printInfo = false) { // switch off printing Info messages if (!printInfo) gErrorIgnoreLevel = kWarning; // switch on histogram errors TH1::SetDefaultSumw2(); // random generator rnd=new TRandom3(); // data and MC number of events Double_t const nData0= 500.0; Double_t const nMC0 = 50000.0; // Binning // reconstructed variable (0-10) Int_t const nDet=15; Double_t const xminDet=0.0; Double_t const xmaxDet=15.0; // signal binning (three shapes: 0,1,2) Int_t const nGen=3; Double_t const xminGen=-0.5; Double_t const xmaxGen= 2.5; // parameters // fraction of events per signal shape static const Double_t genFrac[]={0.3,0.6,0.1}; // signal shapes static const Double_t genShape[][5]= {{1.0,2.0,1.5,0.,15.}, {1.0,7.0,2.5,0.,15.}, {0.0,0.0,0.0,0.,15.}}; // define DATA histograms // observed data distribution TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet); // define MC histograms // matrix of migrations TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)", nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet); TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen); TH1D **histPullNC=new TH1D* [nGen]; TH1D **histPullArea=new TH1D* [nGen]; for(int i=0;iReset(); histGenDetMC->Reset(); Int_t nData=rnd->Poisson(nData0); for(Int_t i=0;iFill(yObs); } Int_t nMC=rnd->Poisson(nMC0); for(Int_t i=0;iFill(iGen,yObs); } /* for(Int_t ix=0;ix<=histGenDetMC->GetNbinsX()+1;ix++) { for(Int_t iy=0;iy<=histGenDetMC->GetNbinsY()+1;iy++) { cout<GetBinContent(ix,iy)<<"\n"; } } */ //======================== // unfolding TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz, TUnfold::kRegModeSize,TUnfold::kEConstraintNone); // define the input vector (the measured data distribution) unfold.SetInput(histDetDATA,0.0,1.0); // run the unfolding unfold.ScanLcurve(50,0.,0.,nullptr,nullptr,nullptr); // fill pull distributions without constraint unfold.GetOutput(histUnfold); for(int i=0;iFill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/ histUnfold->GetBinError(i+1)); } // repeat unfolding on the same data, now with Area constraint unfold.SetConstraint(TUnfold::kEConstraintArea); // run the unfolding unfold.ScanLcurve(50,0.,0.,nullptr,nullptr,nullptr); // fill pull distributions with constraint unfold.GetOutput(histUnfold); for(int i=0;iFill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/ histUnfold->GetBinError(i+1)); } } TCanvas output; output.Divide(3,2); gStyle->SetOptFit(1111); for(int i=0;iFit("gaus"); histPullNC[i]->Draw(); } for(int i=0;iFit("gaus"); histPullArea[i]->Draw(); } output.SaveAs("testUnfold4.ps"); }