00001 // Author: Stefan Schmitt 00002 // DESY, 14.10.2008 00003 00004 // Version 17.6, in parallel to changes in TUnfold 00005 // 00006 // History: 00007 // Version 17.5, in parallel to changes in TUnfold 00008 // Version 17.4, in parallel to changes in TUnfold 00009 // Version 17.3, in parallel to changes in TUnfold 00010 // Version 17.2, in parallel to changes in TUnfold 00011 // Version 17.1, in parallel to changes in TUnfold 00012 // Version 17.0, updated for changed methods in TUnfold 00013 // Version 16.1, parallel to changes in TUnfold 00014 // Version 16.0, parallel to changes in TUnfold 00015 // Version 15, with automatic L-curve scan, simplified example 00016 // Version 14, with changes in TUnfoldSys.cxx 00017 // Version 13, with changes to TUnfold.C 00018 // Version 12, with improvements to TUnfold.cxx 00019 // Version 11, print chi**2 and number of degrees of freedom 00020 // Version 10, with bug-fix in TUnfold.cxx 00021 // Version 9, with bug-fix in TUnfold.cxx, TUnfold.h 00022 // Version 8, with bug-fix in TUnfold.cxx, TUnfold.h 00023 // Version 7, with bug-fix in TUnfold.cxx, TUnfold.h 00024 // Version 6a, fix problem with dynamic array allocation under windows 00025 // Version 6, re-include class MyUnfold in the example 00026 // Version 5, move class MyUnfold to seperate files 00027 // Version 4, with bug-fix in TUnfold.C 00028 // Version 3, with bug-fix in TUnfold.C 00029 // Version 2, with changed ScanLcurve() arguments 00030 // Version 1, remove L curve analysis, use ScanLcurve() method instead 00031 // Version 0, L curve analysis included here 00032 00033 /* 00034 This file is part of TUnfold. 00035 00036 TUnfold is free software: you can redistribute it and/or modify 00037 it under the terms of the GNU General Public License as published by 00038 the Free Software Foundation, either version 3 of the License, or 00039 (at your option) any later version. 00040 00041 TUnfold is distributed in the hope that it will be useful, 00042 but WITHOUT ANY WARRANTY; without even the implied warranty of 00043 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00044 GNU General Public License for more details. 00045 00046 You should have received a copy of the GNU General Public License 00047 along with TUnfold. If not, see <http://www.gnu.org/licenses/>. 00048 */ 00049 00050 #include <TMath.h> 00051 #include <TCanvas.h> 00052 #include <TRandom3.h> 00053 #include <TFitter.h> 00054 #include <TF1.h> 00055 #include <TStyle.h> 00056 #include <TVector.h> 00057 #include <TGraph.h> 00058 #include "TUnfold.h" 00059 00060 using namespace std; 00061 00062 /////////////////////////////////////////////////////////////////////// 00063 // 00064 // Test program as an example for a user specific regularisation scheme 00065 // 00066 // (1) Generate Monte Carlo and Data events 00067 // The events consist of 00068 // signal 00069 // background 00070 // 00071 // The signal is a resonance. It is generated with a Breit-Wigner, 00072 // smeared by a Gaussian 00073 // 00074 // (2) Unfold the data. The result is: 00075 // The background level 00076 // The shape of the resonance, corrected for detector effects 00077 // 00078 // The regularisation is done on the curvature, excluding the bins 00079 // near the peak. 00080 // 00081 // (3) produce some plots 00082 // 00083 /////////////////////////////////////////////////////////////////////// 00084 00085 TRandom *rnd=0; 00086 00087 // generate an event 00088 // output: 00089 // negative mass: background event 00090 // positive mass: signal event 00091 Double_t GenerateEvent(Double_t bgr, // relative fraction of background 00092 Double_t mass, // peak position 00093 Double_t gamma) // peak width 00094 { 00095 Double_t t; 00096 if(rnd->Rndm()>bgr) { 00097 // generate signal event 00098 // with positive mass 00099 do { 00100 do { 00101 t=rnd->Rndm(); 00102 } while(t>=1.0); 00103 t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass; 00104 } while(t<=0.0); 00105 return t; 00106 } else { 00107 // generate background event 00108 // generate events following a power-law distribution 00109 // f(E) = K * TMath::power((E0+E),N0) 00110 static Double_t const E0=2.4; 00111 static Double_t const N0=2.9; 00112 do { 00113 do { 00114 t=rnd->Rndm(); 00115 } while(t>=1.0); 00116 // the mass is returned negative 00117 // In our example a convenient way to indicate it is a background event. 00118 t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0; 00119 } while(t>=0.0); 00120 return t; 00121 } 00122 } 00123 00124 // smear the event to detector level 00125 // input: 00126 // mass on generator level (mTrue>0 !) 00127 // output: 00128 // mass on detector level 00129 Double_t DetectorEvent(Double_t mTrue) { 00130 // smear by double-gaussian 00131 static Double_t frac=0.1; 00132 static Double_t wideBias=0.03; 00133 static Double_t wideSigma=0.5; 00134 static Double_t smallBias=0.0; 00135 static Double_t smallSigma=0.1; 00136 if(rnd->Rndm()>frac) { 00137 return rnd->Gaus(mTrue+smallBias,smallSigma); 00138 } else { 00139 return rnd->Gaus(mTrue+wideBias,wideSigma); 00140 } 00141 } 00142 00143 int testUnfold2() 00144 { 00145 // switch on histogram errors 00146 TH1::SetDefaultSumw2(); 00147 00148 // random generator 00149 rnd=new TRandom3(); 00150 00151 // data and MC luminosity, cross-section 00152 Double_t const luminosityData=100000; 00153 Double_t const luminosityMC=1000000; 00154 Double_t const crossSection=1.0; 00155 00156 Int_t const nDet=250; 00157 Int_t const nGen=100; 00158 Double_t const xminDet=0.0; 00159 Double_t const xmaxDet=10.0; 00160 Double_t const xminGen=0.0; 00161 Double_t const xmaxGen=10.0; 00162 00163 //============================================ 00164 // generate MC distribution 00165 // 00166 TH1D *histMgenMC=new TH1D("MgenMC",";mass(gen)",nGen,xminGen,xmaxGen); 00167 TH1D *histMdetMC=new TH1D("MdetMC",";mass(det)",nDet,xminDet,xmaxDet); 00168 TH2D *histMdetGenMC=new TH2D("MdetgenMC",";mass(det);mass(gen)",nDet,xminDet,xmaxDet, 00169 nGen,xminGen,xmaxGen); 00170 Int_t neventMC=rnd->Poisson(luminosityMC*crossSection); 00171 for(Int_t i=0;i<neventMC;i++) { 00172 Double_t mGen=GenerateEvent(0.3, // relative fraction of background 00173 4.0, // peak position in MC 00174 0.2); // peak width in MC 00175 Double_t mDet=DetectorEvent(TMath::Abs(mGen)); 00176 // the generated mass is negative for background 00177 // and positive for signal 00178 // so it will be filled in the underflow bin 00179 // this is very convenient for the unfolding: 00180 // the unfolded result will contain the number of background 00181 // events in the underflow bin 00182 00183 // generated MC distribution (for comparison only) 00184 histMgenMC->Fill(mGen,luminosityData/luminosityMC); 00185 // reconstructed MC distribution (for comparison only) 00186 histMdetMC->Fill(mDet,luminosityData/luminosityMC); 00187 00188 // matrix describing how the generator input migrates to the 00189 // reconstructed level. Unfolding input. 00190 // NOTE on underflow/overflow bins: 00191 // (1) the detector level under/overflow bins are used for 00192 // normalisation ("efficiency" correction) 00193 // in our toy example, these bins are populated from tails 00194 // of the initial MC distribution. 00195 // (2) the generator level underflow/overflow bins are 00196 // unfolded. In this example: 00197 // underflow bin: background events reconstructed in the detector 00198 // overflow bin: signal events generated at masses > xmaxDet 00199 // for the unfolded result these bins will be filled 00200 // -> the background normalisation will be contained in the underflow bin 00201 histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC); 00202 } 00203 00204 //============================================ 00205 // generate data distribution 00206 // 00207 TH1D *histMgenData=new TH1D("MgenData",";mass(gen)",nGen,xminGen,xmaxGen); 00208 TH1D *histMdetData=new TH1D("MdetData",";mass(det)",nDet,xminDet,xmaxDet); 00209 Int_t neventData=rnd->Poisson(luminosityData*crossSection); 00210 for(Int_t i=0;i<neventData;i++) { 00211 Double_t mGen=GenerateEvent(0.4, // relative fraction of background 00212 3.8, // peak position 00213 0.15); // peak width 00214 Double_t mDet=DetectorEvent(TMath::Abs(mGen)); 00215 // generated data mass for comparison plots 00216 // for real data, we do not have this histogram 00217 histMgenData->Fill(mGen); 00218 00219 // reconstructed mass, unfolding input 00220 histMdetData->Fill(mDet); 00221 } 00222 00223 //========================================================================= 00224 // set up the unfolding 00225 TUnfold unfold(histMdetGenMC,TUnfold::kHistMapOutputVert, 00226 TUnfold::kRegModeNone); 00227 // regularisation 00228 //---------------- 00229 // the regularisation is done on the curvature (2nd derivative) of 00230 // the output distribution 00231 // 00232 // One has to exclude the bins near the peak of the Breit-Wigner, 00233 // because there the curvature is high 00234 // (and the regularisation eventually could enforce a small 00235 // curvature, thus biasing result) 00236 // 00237 // in real life, the parameters below would have to be optimized, 00238 // depending on the data peak position and width 00239 // Or maybe one finds a different regularisation scheme... this is 00240 // just an example... 00241 Double_t estimatedPeakPosition=3.8; 00242 Int_t nPeek=3; 00243 TUnfold::ERegMode regMode=TUnfold::kRegModeCurvature; 00244 // calculate bin number correspoinding to estimated peak position 00245 Int_t iPeek=(Int_t)(nGen*(estimatedPeakPosition-xminGen)/(xmaxGen-xminGen) 00246 // offset 1.5 00247 // accounts for start bin 1 00248 // and rounding errors +0.5 00249 +1.5); 00250 // regularize output bins 1..iPeek-nPeek 00251 unfold.RegularizeBins(1,1,iPeek-nPeek,regMode); 00252 // regularize output bins iPeek+nPeek..nGen 00253 unfold.RegularizeBins(iPeek+nPeek,1,nGen-(iPeek+nPeek),regMode); 00254 00255 // unfolding 00256 //----------- 00257 00258 // set input distribution and bias scale (=0) 00259 if(unfold.SetInput(histMdetData,0.0)>=10000) { 00260 std::cout<<"Unfolding result may be wrong\n"; 00261 } 00262 00263 // do the unfolding here 00264 Double_t tauMin=0.0; 00265 Double_t tauMax=0.0; 00266 Int_t nScan=30; 00267 Int_t iBest; 00268 TSpline *logTauX,*logTauY; 00269 TGraph *lCurve; 00270 // this method scans the parameter tau and finds the kink in the L curve 00271 // finally, the unfolding is done for the "best" choice of tau 00272 iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY); 00273 std::cout<<"tau="<<unfold.GetTau()<<"\n"; 00274 std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L() 00275 <<" / "<<unfold.GetNdf()<<"\n"; 00276 00277 // save point corresponding to the kink in the L curve as TGraph 00278 Double_t t[1],x[1],y[1]; 00279 logTauX->GetKnot(iBest,t[0],x[0]); 00280 logTauY->GetKnot(iBest,t[0],y[0]); 00281 TGraph *bestLcurve=new TGraph(1,x,y); 00282 TGraph *bestLogTauX=new TGraph(1,t,x); 00283 00284 //============================================================ 00285 // extract unfolding results into histograms 00286 00287 // set up a bin map, excluding underflow and overflow bins 00288 // the binMap relates the the output of the unfolding to the final 00289 // histogram bins 00290 Int_t *binMap=new Int_t[nGen+2]; 00291 for(Int_t i=1;i<=nGen;i++) binMap[i]=i; 00292 binMap[0]=-1; 00293 binMap[nGen+1]=-1; 00294 00295 TH1D *histMunfold=new TH1D("Unfolded",";mass(gen)",nGen,xminGen,xmaxGen); 00296 unfold.GetOutput(histMunfold,binMap); 00297 TH1D *histMdetFold=new TH1D("FoldedBack","mass(det)",nDet,xminDet,xmaxDet); 00298 unfold.GetFoldedOutput(histMdetFold); 00299 00300 // store global correlation coefficients 00301 TH1D *histRhoi=new TH1D("rho_I","mass",nGen,xminGen,xmaxGen); 00302 unfold.GetRhoI(histRhoi,binMap); 00303 00304 delete[] binMap; 00305 binMap=0; 00306 00307 //===================================================================== 00308 // plot some histograms 00309 TCanvas output; 00310 00311 // produce some plots 00312 output.Divide(3,2); 00313 00314 // Show the matrix which connects input and output 00315 // There are overflow bins at the bottom, not shown in the plot 00316 // These contain the background shape. 00317 // The overflow bins to the left and right contain 00318 // events which are not reconstructed. These are necessary for proper MC 00319 // normalisation 00320 output.cd(1); 00321 histMdetGenMC->Draw("BOX"); 00322 00323 // draw generator-level distribution: 00324 // data (red) [for real data this is not available] 00325 // MC input (black) [with completely wrong peak position and shape] 00326 // unfolded data (blue) 00327 output.cd(2); 00328 histMunfold->SetLineColor(kBlue); 00329 histMunfold->Draw(); 00330 histMgenData->SetLineColor(kRed); 00331 histMgenData->Draw("SAME"); 00332 histMgenMC->Draw("SAME HIST"); 00333 00334 // show detector level distributions 00335 // data (red) 00336 // MC (black) 00337 // unfolded data (blue) 00338 output.cd(3); 00339 histMdetFold->SetLineColor(kBlue); 00340 histMdetFold->Draw(); 00341 histMdetData->SetLineColor(kRed); 00342 histMdetData->Draw("SAME"); 00343 histMdetMC->Draw("SAME HIST"); 00344 00345 // show correlation coefficients 00346 // all bins outside the peak are found to be highly correlated 00347 // But they are compatible with zero anyway 00348 // If the peak shape is fitted, 00349 // these correlations have to be taken into account, see example 00350 output.cd(4); 00351 histRhoi->Draw(); 00352 00353 // show rhoi_max(tau) distribution 00354 output.cd(5); 00355 logTauX->Draw(); 00356 bestLogTauX->SetMarkerColor(kRed); 00357 bestLogTauX->Draw("*"); 00358 00359 output.cd(6); 00360 lCurve->Draw("AL"); 00361 bestLcurve->SetMarkerColor(kRed); 00362 bestLcurve->Draw("*"); 00363 00364 output.SaveAs("testUnfold2.ps"); 00365 return 0; 00366 }