00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #include <TError.h>
00051 #include <TMath.h>
00052 #include <TCanvas.h>
00053 #include <TRandom3.h>
00054 #include <TFitter.h>
00055 #include <TF1.h>
00056 #include <TStyle.h>
00057 #include <TVector.h>
00058 #include <TGraph.h>
00059
00060 #include "TUnfoldDensity.h"
00061
00062
00063
00064 using namespace std;
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 TRandom *rnd=0;
00108
00109 TH2 *gHistInvEMatrix;
00110
00111 TVirtualFitter *gFitter=0;
00112
00113 void chisquare_corr(Int_t &npar, Double_t * , Double_t &f, Double_t *u, Int_t ) {
00114
00115
00116
00117
00118 Double_t x;
00119 TH1 *hfit = (TH1*)gFitter->GetObjectFit();
00120 TF1 *f1 = (TF1*)gFitter->GetUserFunc();
00121
00122 f1->InitArgs(&x,u);
00123 npar = f1->GetNpar();
00124 f = 0;
00125
00126 Int_t npfit = 0;
00127 Int_t nPoints=hfit->GetNbinsX();
00128 Double_t *df=new Double_t[nPoints];
00129 for (Int_t i=0;i<nPoints;i++) {
00130 x = hfit->GetBinCenter(i+1);
00131 TF1::RejectPoint(kFALSE);
00132 df[i] = f1->EvalPar(&x,u)-hfit->GetBinContent(i+1);
00133 if (TF1::RejectedPoint()) df[i]=0.0;
00134 else npfit++;
00135 }
00136 for (Int_t i=0;i<nPoints;i++) {
00137 for (Int_t j=0;j<nPoints;j++) {
00138 f += df[i]*df[j]*gHistInvEMatrix->GetBinContent(i+1,j+1);
00139 }
00140 }
00141 delete[] df;
00142 f1->SetNumberFitPoints(npfit);
00143 }
00144
00145 Double_t bw_func(Double_t *x,Double_t *par) {
00146 Double_t dm=x[0]-par[1];
00147 return par[0]/(dm*dm+par[2]*par[2]);
00148 }
00149
00150
00151
00152
00153
00154
00155 Double_t GenerateEvent(Double_t bgr,
00156 Double_t mass,
00157 Double_t gamma)
00158 {
00159 Double_t t;
00160 if(rnd->Rndm()>bgr) {
00161
00162
00163 do {
00164 do {
00165 t=rnd->Rndm();
00166 } while(t>=1.0);
00167 t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass;
00168 } while(t<=0.0);
00169 return t;
00170 } else {
00171
00172
00173
00174 static Double_t const E0=2.4;
00175 static Double_t const N0=2.9;
00176 do {
00177 do {
00178 t=rnd->Rndm();
00179 } while(t>=1.0);
00180
00181
00182 t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0;
00183 } while(t>=0.0);
00184 return t;
00185 }
00186 }
00187
00188
00189
00190
00191
00192
00193 Double_t DetectorEvent(Double_t mTrue) {
00194
00195 static Double_t frac=0.1;
00196 static Double_t wideBias=0.03;
00197 static Double_t wideSigma=0.5;
00198 static Double_t smallBias=0.0;
00199 static Double_t smallSigma=0.1;
00200 if(rnd->Rndm()>frac) {
00201 return rnd->Gaus(mTrue+smallBias,smallSigma);
00202 } else {
00203 return rnd->Gaus(mTrue+wideBias,wideSigma);
00204 }
00205 }
00206
00207 int testUnfold1()
00208 {
00209
00210 TH1::SetDefaultSumw2();
00211
00212
00213 gStyle->SetOptFit(1111);
00214
00215
00216 rnd=new TRandom3();
00217
00218
00219 Double_t const luminosityData=100000;
00220 Double_t const luminosityMC=1000000;
00221 Double_t const crossSection=1.0;
00222
00223 Int_t const nDet=250;
00224 Int_t const nGen=100;
00225 Double_t const xminDet=0.0;
00226 Double_t const xmaxDet=10.0;
00227 Double_t const xminGen=0.0;
00228 Double_t const xmaxGen=10.0;
00229
00230
00231
00232
00233 TH1D *histMgenMC=new TH1D("MgenMC",";mass(gen)",nGen,xminGen,xmaxGen);
00234 TH1D *histMdetMC=new TH1D("MdetMC",";mass(det)",nDet,xminDet,xmaxDet);
00235 TH2D *histMdetGenMC=new TH2D("MdetgenMC",";mass(det);mass(gen)",
00236 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
00237 Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);
00238 for(Int_t i=0;i<neventMC;i++) {
00239 Double_t mGen=GenerateEvent(0.3,
00240 4.0,
00241 0.2);
00242 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
00243
00244
00245
00246
00247
00248
00249
00250
00251 histMgenMC->Fill(mGen,luminosityData/luminosityMC);
00252
00253 histMdetMC->Fill(mDet,luminosityData/luminosityMC);
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC);
00269 }
00270
00271
00272
00273
00274
00275 TH2D *histMdetGenSysMC=new TH2D("MdetgenSysMC",";mass(det);mass(gen)",
00276 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
00277 neventMC=rnd->Poisson(luminosityMC*crossSection);
00278 for(Int_t i=0;i<neventMC;i++) {
00279 Double_t mGen=GenerateEvent
00280 (0.5,
00281 3.6,
00282 0.15);
00283 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
00284 histMdetGenSysMC->Fill(mDet,mGen,luminosityData/luminosityMC);
00285 }
00286
00287
00288
00289
00290 TH1D *histMgenData=new TH1D("MgenData",";mass(gen)",nGen,xminGen,xmaxGen);
00291 TH1D *histMdetData=new TH1D("MdetData",";mass(det)",nDet,xminDet,xmaxDet);
00292 Int_t neventData=rnd->Poisson(luminosityData*crossSection);
00293 for(Int_t i=0;i<neventData;i++) {
00294 Double_t mGen=GenerateEvent(0.4,
00295 3.8,
00296 0.15);
00297 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
00298
00299
00300 histMgenData->Fill(mGen);
00301
00302
00303 histMdetData->Fill(mDet);
00304 }
00305
00306
00307
00308 TH1D *histDensityGenData=new TH1D("DensityGenData",";mass(gen)",
00309 nGen,xminGen,xmaxGen);
00310 TH1D *histDensityGenMC=new TH1D("DensityGenMC",";mass(gen)",
00311 nGen,xminGen,xmaxGen);
00312 for(Int_t i=1;i<=nGen;i++) {
00313 histDensityGenData->SetBinContent(i,histMgenData->GetBinContent(i)/
00314 histMgenData->GetBinWidth(i));
00315 histDensityGenMC->SetBinContent(i,histMgenMC->GetBinContent(i)/
00316 histMgenMC->GetBinWidth(i));
00317 }
00318
00319
00320
00321
00322 TUnfoldDensity unfold(histMdetGenMC,TUnfold::kHistMapOutputVert);
00323
00324
00325
00326
00327
00328
00329
00330 if(unfold.SetInput(histMdetData)>=10000) {
00331 std::cout<<"Unfolding result may be wrong\n";
00332 }
00333
00334
00335
00336
00337
00338 Int_t nScan=30;
00339
00340 Double_t tauMin=0.0;
00341 Double_t tauMax=0.0;
00342 Int_t iBest;
00343 TSpline *logTauX,*logTauY;
00344 TGraph *lCurve;
00345
00346
00347 #ifdef VERBOSE_LCURVE_SCAN
00348 Int_t oldinfo=gErrorIgnoreLevel;
00349 gErrorIgnoreLevel=kInfo;
00350 #endif
00351
00352
00353 iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
00354
00355
00356 #ifdef VERBOSE_LCURVE_SCAN
00357 gErrorIgnoreLevel=oldinfo;
00358 #endif
00359
00360
00361
00362
00363
00364 Double_t SYS_ERROR1_MSTART=6;
00365 Double_t SYS_ERROR1_SIZE=0.1;
00366 TH2D *histMdetGenSys1=new TH2D("Mdetgensys1",";mass(det);mass(gen)",
00367 nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
00368 for(Int_t i=0;i<=nDet+1;i++) {
00369 if(histMdetData->GetBinCenter(i)>=SYS_ERROR1_MSTART) {
00370 for(Int_t j=0;j<=nGen+1;j++) {
00371 histMdetGenSys1->SetBinContent(i,j,SYS_ERROR1_SIZE);
00372 }
00373 }
00374 }
00375 unfold.AddSysError(histMdetGenSysMC,"SYSERROR_MC",TUnfold::kHistMapOutputVert,
00376 TUnfoldSys::kSysErrModeMatrix);
00377 unfold.AddSysError(histMdetGenSys1,"SYSERROR1",TUnfold::kHistMapOutputVert,
00378 TUnfoldSys::kSysErrModeRelative);
00379
00380
00381
00382
00383 std::cout<<"tau="<<unfold.GetTau()<<"\n";
00384 std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
00385 <<" / "<<unfold.GetNdf()<<"\n";
00386 std::cout<<"chi**2(sys)="<<unfold.GetChi2Sys()<<"\n";
00387
00388
00389
00390
00391
00392 Double_t t[1],x[1],y[1];
00393 logTauX->GetKnot(iBest,t[0],x[0]);
00394 logTauY->GetKnot(iBest,t[0],y[0]);
00395 TGraph *bestLcurve=new TGraph(1,x,y);
00396 TGraph *bestLogTauLogChi2=new TGraph(1,t,x);
00397
00398
00399
00400
00401
00402 TH1 *histMunfold=unfold.GetOutput("Unfolded");
00403
00404
00405 TH1 *histMdetFold=unfold.GetFoldedOutput("FoldedBack");
00406
00407
00408
00409
00410
00411
00412
00413 TH2 *histEmatTotal=unfold.GetEmatrixTotal("EmatTotal");
00414
00415
00416 TH1D *histTotalError=
00417 new TH1D("TotalError",";mass(gen)",nGen,xminGen,xmaxGen);
00418 for(Int_t bin=1;bin<=nGen;bin++) {
00419 histTotalError->SetBinContent(bin,histMunfold->GetBinContent(bin));
00420 histTotalError->SetBinError
00421 (bin,TMath::Sqrt(histEmatTotal->GetBinContent(bin,bin)));
00422 }
00423
00424
00425
00426
00427
00428
00429 TH1 *histRhoi=unfold.GetRhoItotal("rho_I",
00430 0,
00431 0,
00432 "*[UO]",
00433 kTRUE,
00434 &gHistInvEMatrix
00435 );
00436
00437
00438
00439
00440
00441
00442 gFitter=TVirtualFitter::Fitter(histMunfold);
00443 gFitter->SetFCN(chisquare_corr);
00444
00445 TF1 *bw=new TF1("bw",bw_func,xminGen,xmaxGen,3);
00446 bw->SetParameter(0,1000.);
00447 bw->SetParameter(1,3.8);
00448 bw->SetParameter(2,0.2);
00449
00450
00451
00452 histMunfold->Fit(bw,"UE");
00453
00454
00455
00456 TCanvas output;
00457 output.Divide(3,2);
00458
00459
00460
00461
00462
00463
00464
00465 output.cd(1);
00466 histMdetGenMC->Draw("BOX");
00467
00468
00469
00470
00471
00472 output.cd(2);
00473 histTotalError->SetLineColor(kBlue);
00474 histTotalError->Draw("E");
00475 histMunfold->SetLineColor(kGreen);
00476 histMunfold->Draw("SAME E1");
00477 histDensityGenData->SetLineColor(kRed);
00478 histDensityGenData->Draw("SAME");
00479 histDensityGenMC->Draw("SAME HIST");
00480
00481
00482
00483
00484
00485 output.cd(3);
00486 histMdetFold->SetLineColor(kBlue);
00487 histMdetFold->Draw();
00488 histMdetMC->Draw("SAME HIST");
00489
00490 TH1 *histInput=unfold.GetInput("Minput",";mass(det)");
00491
00492 histInput->SetLineColor(kRed);
00493 histInput->Draw("SAME");
00494
00495
00496 output.cd(4);
00497 histRhoi->Draw();
00498
00499
00500 output.cd(5);
00501 logTauX->Draw();
00502 bestLogTauLogChi2->SetMarkerColor(kRed);
00503 bestLogTauLogChi2->Draw("*");
00504
00505
00506 output.cd(6);
00507 lCurve->Draw("AL");
00508 bestLcurve->SetMarkerColor(kRed);
00509 bestLcurve->Draw("*");
00510
00511 output.SaveAs("testUnfold1.ps");
00512
00513 return 0;
00514 }