00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <TMath.h>
00017 #include <TCanvas.h>
00018 #include <TRandom3.h>
00019 #include <TFitter.h>
00020 #include <TF1.h>
00021 #include <TStyle.h>
00022 #include <TVector.h>
00023 #include <TGraph.h>
00024
00025 #include "TUnfoldDensity.h"
00026
00027 using namespace std;
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 TRandom *rnd=0;
00062
00063 Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) {
00064
00065
00066
00067
00068 Double_t f=rnd->Rndm();
00069 Int_t r=0;
00070 while((r<nmax)&&(f>=probability[r])) {
00071 f -= probability[r];
00072 r++;
00073 }
00074 return r;
00075 }
00076
00077 Double_t GenerateRecEvent(const Double_t *shapeParm) {
00078
00079
00080
00081
00082
00083
00084
00085
00086 Double_t f=rnd->Rndm();
00087 Double_t r;
00088 if(f<shapeParm[0]) {
00089 r=rnd->Gaus(shapeParm[1],shapeParm[2]);
00090 } else {
00091 r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3];
00092 }
00093 return r;
00094 }
00095
00096 void testUnfold4()
00097 {
00098
00099 TH1::SetDefaultSumw2();
00100
00101
00102 rnd=new TRandom3();
00103
00104
00105 Double_t const nData0= 500.0;
00106 Double_t const nMC0 = 50000.0;
00107
00108
00109
00110 Int_t const nDet=15;
00111 Double_t const xminDet=0.0;
00112 Double_t const xmaxDet=15.0;
00113
00114
00115 Int_t const nGen=3;
00116 Double_t const xminGen=-0.5;
00117 Double_t const xmaxGen= 2.5;
00118
00119
00120
00121 static const Double_t genFrac[]={0.3,0.6,0.1};
00122
00123
00124 static const Double_t genShape[][5]=
00125 {{1.0,2.0,1.5,0.,15.},
00126 {1.0,7.0,2.5,0.,15.},
00127 {0.0,0.0,0.0,0.,15.}};
00128
00129
00130
00131 TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet);
00132
00133
00134
00135 TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)",
00136 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
00137
00138 TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen);
00139
00140 TH1D **histPullNC=new TH1D* [nGen];
00141 TH1D **histPullArea=new TH1D* [nGen];
00142 for(int i=0;i<nGen;i++) {
00143 histPullNC[i]=new TH1D(TString::Format("PullNC%d",i),"pull",15,-3.,3.);
00144 histPullArea[i]=new TH1D(TString::Format("PullArea%d",i),"pull",15,-3.,3.);
00145 }
00146
00147
00148 cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
00149
00150 for(int itoy=0;itoy<1000;itoy++) {
00151 if(!(itoy %10)) cout<<"toy iteration: "<<itoy<<"\n";
00152 histDetDATA->Reset();
00153 histGenDetMC->Reset();
00154
00155 Int_t nData=rnd->Poisson(nData0);
00156 for(Int_t i=0;i<nData;i++) {
00157 Int_t iGen=GenerateGenEvent(nGen,genFrac);
00158 Double_t yObs=GenerateRecEvent(genShape[iGen]);
00159 histDetDATA->Fill(yObs);
00160 }
00161
00162 Int_t nMC=rnd->Poisson(nMC0);
00163 for(Int_t i=0;i<nMC;i++) {
00164 Int_t iGen=GenerateGenEvent(nGen,genFrac);
00165 Double_t yObs=GenerateRecEvent(genShape[iGen]);
00166 histGenDetMC->Fill(iGen,yObs);
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176 TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz,
00177 TUnfold::kRegModeSize,TUnfold::kEConstraintNone);
00178
00179 unfold.SetInput(histDetDATA,0.0,1.0);
00180
00181
00182 unfold.ScanLcurve(50,0.,0.,0,0,0);
00183
00184
00185 unfold.GetOutput(histUnfold);
00186
00187 for(int i=0;i<nGen;i++) {
00188 histPullNC[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
00189 histUnfold->GetBinError(i+1));
00190 }
00191
00192
00193 unfold.SetConstraint(TUnfold::kEConstraintArea);
00194
00195
00196 unfold.ScanLcurve(50,0.,0.,0,0,0);
00197
00198
00199 unfold.GetOutput(histUnfold);
00200
00201 for(int i=0;i<nGen;i++) {
00202 histPullArea[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
00203 histUnfold->GetBinError(i+1));
00204 }
00205
00206 }
00207 TCanvas output;
00208 output.Divide(3,2);
00209
00210 gStyle->SetOptFit(1111);
00211
00212 for(int i=0;i<nGen;i++) {
00213 output.cd(i+1);
00214 histPullNC[i]->Fit("gaus");
00215 histPullNC[i]->Draw();
00216 }
00217 for(int i=0;i<nGen;i++) {
00218 output.cd(i+4);
00219 histPullArea[i]->Fit("gaus");
00220 histPullArea[i]->Draw();
00221 }
00222 output.SaveAs("testUnfold4.ps");
00223 }