00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <TMath.h>
00018 #include <TCanvas.h>
00019 #include <TRandom3.h>
00020 #include <TFitter.h>
00021 #include <TF1.h>
00022 #include <TStyle.h>
00023 #include <TVector.h>
00024 #include <TGraph.h>
00025
00026 #include "TUnfoldDensity.h"
00027
00028 using namespace std;
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
00062
00063
00064
00065
00066
00067
00068
00069
00070 TRandom *rnd=0;
00071
00072 Double_t GenerateEvent(const Double_t *parm,
00073 const Double_t *triggerParm,
00074 Double_t *intLumi,
00075 Bool_t *triggerFlag,
00076 Double_t *ptGen,Int_t *iType)
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
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 Double_t ptObs;
00119 Bool_t isTriggered=kFALSE;
00120 do {
00121 Int_t itype;
00122 Double_t ptgen;
00123 Double_t f=rnd->Rndm();
00124
00125 itype=0;
00126 if(f<parm[0]) itype=1;
00127 else if(f<parm[0]+parm[1]) itype=2;
00128
00129
00130 Double_t a=parm[4+itype];
00131 Double_t a1=a+1.0;
00132 Double_t t=rnd->Rndm();
00133 if(a1 == 0.0) {
00134 Double_t x0=TMath::Log(parm[2]);
00135 ptgen=TMath::Exp(t*(TMath::Log(parm[3])-x0)+x0);
00136 } else {
00137 Double_t x0=pow(parm[2],a1);
00138 ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1);
00139 }
00140 if(iType) *iType=itype;
00141 if(ptGen) *ptGen=ptgen;
00142
00143
00144 Double_t sigma=
00145 TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
00146 ptObs=rnd->BreitWigner(ptgen,sigma);
00147
00148
00149 Double_t triggerProb =
00150 triggerParm[2]/(1.+TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
00151 isTriggered= rnd->Rndm()<triggerProb;
00152 (*intLumi) ++;
00153 } while((!triggerFlag) && (!isTriggered));
00154
00155 if(triggerFlag) *triggerFlag=isTriggered;
00156 return ptObs;
00157 }
00158
00159
00160 void testUnfold3()
00161 {
00162
00163 TH1::SetDefaultSumw2();
00164
00165
00166 gStyle->SetOptFit(1111);
00167
00168
00169 rnd=new TRandom3();
00170
00171
00172 Double_t const lumiData= 30000;
00173 Double_t const lumiMC =1000000;
00174
00175
00176
00177
00178
00179 Int_t const nDet=24;
00180 Double_t const xminDet=4.0;
00181 Double_t const xmaxDet=28.0;
00182
00183
00184 Int_t const nGen=10;
00185 Double_t const xminGen= 6.0;
00186 Double_t const xmaxGen=26.0;
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 TH1D *histUnfoldInput=
00215 new TH1D("unfolding input rec",";ptrec",nDet,xminDet,xmaxDet);
00216 TH2D *histUnfoldMatrix=
00217 new TH2D("unfolding matrix",";ptgen;ptrec",
00218 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
00219 TH1D *histUnfoldBgr1=
00220 new TH1D("unfolding bgr1 rec",";ptrec",nDet,xminDet,xmaxDet);
00221 TH1D *histUnfoldBgr2=
00222 new TH1D("unfolding bgr2 rec",";ptrec",nDet,xminDet,xmaxDet);
00223 TH2D *histUnfoldMatrixSys=
00224 new TH2D("unfolding matrix sys",";ptgen;ptrec",
00225 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
00226
00227
00228 TH1D *histBbbInput=
00229 new TH1D("bbb input rec",";ptrec",nGen,xminGen,xmaxGen);
00230 TH1D *histBbbSignalRec=
00231 new TH1D("bbb signal rec",";ptrec",nGen,xminGen,xmaxGen);
00232 TH1D *histBbbBgr1=
00233 new TH1D("bbb bgr1 rec",";ptrec",nGen,xminGen,xmaxGen);
00234 TH1D *histBbbBgr2=
00235 new TH1D("bbb bgr2 rec",";ptrec",nGen,xminGen,xmaxGen);
00236 TH1D *histBbbBgrPt=
00237 new TH1D("bbb bgrPt rec",";ptrec",nGen,xminGen,xmaxGen);
00238 TH1D *histBbbSignalGen=
00239 new TH1D("bbb signal gen",";ptgen",nGen,xminGen,xmaxGen);
00240 TH1D *histBbbSignalRecSys=
00241 new TH1D("bbb signal sys rec",";ptrec",nGen,xminGen,xmaxGen);
00242 TH1D *histBbbBgrPtSys=
00243 new TH1D("bbb bgrPt sys rec",";ptrec",nGen,xminGen,xmaxGen);
00244 TH1D *histBbbSignalGenSys=
00245 new TH1D("bbb signal sys gen",";ptgen",nGen,xminGen,xmaxGen);
00246
00247
00248 TH1D *histDataTruth=
00249 new TH1D("DATA truth gen",";ptgen",nGen,xminGen,xmaxGen);
00250 TH1D *histDetMC=
00251 new TH1D("MC prediction rec",";ptrec",nDet,xminDet,xmaxDet);
00252
00253
00254
00255
00256 static Double_t parm_DATA[]={
00257 0.05,
00258 0.05,
00259 3.5,
00260 100.,
00261 -3.0,
00262 0.1,
00263 -0.5,
00264 0.2,
00265 0.01,
00266 };
00267 static Double_t triggerParm_DATA[]={8.0,
00268 4.0,
00269 0.95
00270 };
00271
00272 Double_t intLumi=0.0;
00273 while(intLumi<lumiData) {
00274 Int_t iTypeGen;
00275 Bool_t isTriggered;
00276 Double_t ptGen;
00277 Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
00278 &isTriggered,&ptGen,&iTypeGen);
00279 if(isTriggered) {
00280
00281 histUnfoldInput->Fill(ptObs);
00282
00283
00284 histBbbInput->Fill(ptObs);
00285 }
00286
00287 if(iTypeGen==0) histDataTruth->Fill(ptGen);
00288 }
00289
00290
00291
00292
00293
00294
00295 static Double_t parm_MC[]={
00296 0.05,
00297 0.05,
00298 3.5,
00299 100.,
00300 -4.0,
00301
00302 0.1,
00303 -0.5,
00304 0.2,
00305 0.01,
00306 };
00307 static Double_t triggerParm_MC[]={8.0,
00308 4.0,
00309 0.95
00310 };
00311
00312
00313 Double_t lumiWeight = lumiData/lumiMC;
00314 intLumi=0.0;
00315 while(intLumi<lumiMC) {
00316 Int_t iTypeGen;
00317 Bool_t isTriggered;
00318 Double_t ptGen;
00319 Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
00320 &ptGen,&iTypeGen);
00321 if(!isTriggered) ptObs=0.0;
00322
00323
00324
00325 if(iTypeGen==0) {
00326 histUnfoldMatrix->Fill(ptGen,ptObs,lumiWeight);
00327 } else if(iTypeGen==1) {
00328 histUnfoldBgr1->Fill(ptObs,lumiWeight);
00329 } else if(iTypeGen==2) {
00330 histUnfoldBgr2->Fill(ptObs,lumiWeight);
00331 }
00332
00333
00334 if(iTypeGen==0) {
00335 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
00336 histBbbSignalRec->Fill(ptObs,lumiWeight);
00337 } else {
00338 histBbbBgrPt->Fill(ptObs,lumiWeight);
00339 }
00340 histBbbSignalGen->Fill(ptGen,lumiWeight);
00341 } else if(iTypeGen==1) {
00342 histBbbBgr1->Fill(ptObs,lumiWeight);
00343 } else if(iTypeGen==2) {
00344 histBbbBgr2->Fill(ptObs,lumiWeight);
00345 }
00346
00347
00348 histDetMC ->Fill(ptObs,lumiWeight);
00349 }
00350
00351
00352
00353
00354
00355 static Double_t parm_MC_SYS[]={
00356 0.05,
00357 0.05,
00358 3.5,
00359 100.,
00360 -2.0,
00361 0.1,
00362 -0.5,
00363 0.2,
00364 0.01,
00365 };
00366
00367 intLumi=0.0;
00368 while(intLumi<lumiMC) {
00369 Int_t iTypeGen;
00370 Bool_t isTriggered;
00371 Double_t ptGen;
00372 Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
00373 &isTriggered,&ptGen,&iTypeGen);
00374 if(!isTriggered) ptObs=0.0;
00375
00376
00377 if(iTypeGen==0) {
00378 histUnfoldMatrixSys->Fill(ptGen,ptObs);
00379 }
00380
00381
00382 if(iTypeGen==0) {
00383 if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
00384 histBbbSignalRecSys->Fill(ptObs);
00385 } else {
00386 histBbbBgrPtSys->Fill(ptObs);
00387 }
00388 histBbbSignalGenSys->Fill(ptGen);
00389 }
00390 }
00391
00392
00393 cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
00394
00395
00396
00397
00398
00399
00400
00401
00402 TUnfold::ERegMode regMode =
00403 TUnfold::kRegModeCurvature;
00404
00405
00406 TUnfold::EConstraint constraintMode=
00407 TUnfold::kEConstraintArea;
00408
00409
00410 TUnfoldDensity::EDensityMode densityFlags=
00411 TUnfoldDensity::kDensityModeBinWidth;
00412
00413
00414 TUnfoldDensity unfold(histUnfoldMatrix,TUnfold::kHistMapOutputHoriz,
00415 regMode,constraintMode,densityFlags);
00416
00417
00418 unfold.SetInput(histUnfoldInput);
00419
00420
00421
00422 Double_t scale_bgr=1.0;
00423 Double_t dscale_bgr=0.1;
00424 unfold.SubtractBackground(histUnfoldBgr1,"background1",scale_bgr,dscale_bgr);
00425 unfold.SubtractBackground(histUnfoldBgr2,"background2",scale_bgr,dscale_bgr);
00426
00427
00428 unfold.AddSysError(histUnfoldMatrixSys,"signalshape_SYS",
00429 TUnfold::kHistMapOutputHoriz,
00430 TUnfoldSys::kSysErrModeMatrix);
00431
00432
00433 Int_t nScan=30;
00434 TSpline *logTauX,*logTauY;
00435 TGraph *lCurve;
00436
00437
00438 Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);
00439
00440 cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
00441 <<" / "<<unfold.GetNdf()<<"\n";
00442
00443
00444 Double_t t[1],x[1],y[1];
00445 logTauX->GetKnot(iBest,t[0],x[0]);
00446 logTauY->GetKnot(iBest,t[0],y[0]);
00447 TGraph *bestLcurve=new TGraph(1,x,y);
00448 TGraph *bestLogTauLogChi2=new TGraph(1,t,x);
00449
00450
00451
00452
00453
00454
00455
00456
00457 TH1 *histUnfoldOutput=unfold.GetOutput("PT(unfold,stat+bgrerr)");
00458
00459
00460 TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");
00461
00462
00463
00464 TH2 *histEmatTotal=unfold.GetEmatrixTotal("unfolding total error matrix");
00465
00466
00467
00468 TH1 *histUnfoldStat=new TH1D("PT(unfold,staterr)",";Pt(gen)",
00469 nGen,xminGen,xmaxGen);
00470 TH1 *histUnfoldTotal=new TH1D("PT(unfold,totalerr)",";Pt(gen)",
00471 nGen,xminGen,xmaxGen);
00472 for(Int_t i=0;i<nGen+2;i++) {
00473 Double_t c=histUnfoldOutput->GetBinContent(i);
00474
00475
00476 histUnfoldStat->SetBinContent(i,c);
00477 histUnfoldStat->SetBinError
00478 (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i)));
00479
00480
00481 histUnfoldTotal->SetBinContent(i,c);
00482 histUnfoldTotal->SetBinError
00483 (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i)));
00484 }
00485
00486
00487 TH2D *histCorr=new TH2D("Corr(total)",";Pt(gen);Pt(gen)",
00488 nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
00489 for(Int_t i=0;i<nGen+2;i++) {
00490 Double_t ei,ej;
00491 ei=TMath::Sqrt(histEmatTotal->GetBinContent(i,i));
00492 if(ei<=0.0) continue;
00493 for(Int_t j=0;j<nGen+2;j++) {
00494 ej=TMath::Sqrt(histEmatTotal->GetBinContent(j,j));
00495 if(ej<=0.0) continue;
00496 histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej);
00497 }
00498 }
00499
00500
00501 TH1 *histDetNormBgr1=unfold.GetBackground("bgr1 normalized",
00502 "background1");
00503 TH1 *histDetNormBgrTotal=unfold.GetBackground("bgr total normalized");
00504
00505
00506
00507
00508 TCanvas output;
00509 output.Divide(3,2);
00510 output.cd(1);
00511
00512 histUnfoldInput->SetMinimum(0.0);
00513 histUnfoldInput->Draw("E");
00514 histDetMC->SetMinimum(0.0);
00515 histDetMC->SetLineColor(kBlue);
00516 histDetNormBgrTotal->SetLineColor(kRed);
00517 histDetNormBgr1->SetLineColor(kCyan);
00518 histDetMC->Draw("SAME HIST");
00519 histDetNormBgr1->Draw("SAME HIST");
00520 histDetNormBgrTotal->Draw("SAME HIST");
00521
00522 output.cd(2);
00523
00524 histUnfoldTotal->SetMinimum(0.0);
00525 histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5);
00526
00527 histUnfoldTotal->Draw("E");
00528
00529 histUnfoldOutput->Draw("SAME E1");
00530
00531 histUnfoldStat->Draw("SAME E1");
00532 histDataTruth->Draw("SAME HIST");
00533 histBbbSignalGen->SetLineColor(kBlue);
00534 histBbbSignalGen->Draw("SAME HIST");
00535
00536 output.cd(3);
00537
00538 histUnfoldMatrix->SetLineColor(kBlue);
00539 histUnfoldMatrix->Draw("BOX");
00540
00541
00542 output.cd(4);
00543 logTauX->Draw();
00544 bestLogTauLogChi2->SetMarkerColor(kRed);
00545 bestLogTauLogChi2->Draw("*");
00546
00547
00548 output.cd(5);
00549 lCurve->Draw("AL");
00550 bestLcurve->SetMarkerColor(kRed);
00551 bestLcurve->Draw("*");
00552
00553
00554 output.cd(6);
00555 histCorr->Draw("BOX");
00556
00557 output.SaveAs("testUnfold3.ps");
00558
00559
00560
00561
00562 std::cout<<"bin truth result (stat) (bgr) (sys)\n";
00563 std::cout<<"===+=====+=========+========+========+=======\n";
00564 for(Int_t i=1;i<=nGen;i++) {
00565
00566 Double_t data=histBbbInput->GetBinContent(i);
00567 Double_t errData=histBbbInput->GetBinError(i);
00568
00569
00570 Double_t data_bgr=data
00571 - scale_bgr*(histBbbBgr1->GetBinContent(i)
00572 + histBbbBgr2->GetBinContent(i)
00573 + histBbbBgrPt->GetBinContent(i));
00574 Double_t errData_bgr=TMath::Sqrt
00575 (errData*errData+
00576 TMath::Power(dscale_bgr*histBbbBgr1->GetBinContent(i),2)+
00577 TMath::Power(scale_bgr*histBbbBgr1->GetBinError(i),2)+
00578 TMath::Power(dscale_bgr*histBbbBgr2->GetBinContent(i),2)+
00579 TMath::Power(scale_bgr*histBbbBgr2->GetBinError(i),2)+
00580 TMath::Power(dscale_bgr*histBbbBgrPt->GetBinContent(i),2)+
00581 TMath::Power(scale_bgr*histBbbBgrPt->GetBinError(i),2));
00582
00583 Double_t fCorr=(histBbbSignalGen->GetBinContent(i)/
00584 histBbbSignalRec->GetBinContent(i));
00585 Double_t data_bbb= data_bgr *fCorr;
00586
00587 Double_t errData_stat_bbb = errData*fCorr;
00588
00589 Double_t errData_statbgr_bbb = errData_bgr*fCorr;
00590
00591
00592
00593 Double_t fCorr_sys=(histBbbSignalGenSys->GetBinContent(i)/
00594 histBbbSignalRecSys->GetBinContent(i));
00595 Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;
00596
00597
00598 Double_t errData_total_bbb=
00599 TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
00600 +shift_sys*shift_sys);
00601
00602
00603 Double_t data_unfold= histUnfoldStat->GetBinContent(i);
00604 Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i);
00605 Double_t errData_statbgr_unfold=histUnfoldOutput->GetBinError(i);
00606 Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);
00607
00608
00609
00610 std::cout<<TString::Format
00611 ("%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
00612 i,histDataTruth->GetBinContent(i),data_unfold,
00613 errData_stat_unfold,TMath::Sqrt(errData_statbgr_unfold*
00614 errData_statbgr_unfold-
00615 errData_stat_unfold*
00616 errData_stat_unfold),
00617 TMath::Sqrt(errData_total_unfold*
00618 errData_total_unfold-
00619 errData_statbgr_unfold*
00620 errData_statbgr_unfold))<<"\n";
00621 std::cout<<TString::Format
00622 (" %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
00623 data_bbb,errData_stat_bbb,TMath::Sqrt(errData_statbgr_bbb*
00624 errData_statbgr_bbb-
00625 errData_stat_bbb*
00626 errData_stat_bbb),
00627 TMath::Sqrt(errData_total_bbb*
00628 errData_total_bbb-
00629 errData_statbgr_bbb*
00630 errData_statbgr_bbb))
00631 <<"\n\n";
00632 }
00633 }