// Test program for the class TUnfoldSys // Author: Stefan Schmitt // DESY, 14.10.2008 // Version 16, parallel to changes in TUnfold // // History: // Version 15, simple example including background subtraction #include #include #include #include #include #include #include #include #include using namespace std; /////////////////////////////////////////////////////////////////////// // // Test program for the class TUnfoldSys // // the goal is to unfold the underlying "true" distribution of a variable Pt // // the reconstructed Pt is measured in 24 bins from 4 to 28 // the generator-level Pt is unfolded into 10 bins from 6 to 26 // plus underflow bin from 0 to 6 // plus overflow bin above 26 // there are two background sources // bgr1 and bgr2 // the signal has a finite trigger efficiency at a threshold of 8 GeV // // two types of systematic error are studied: // SYS1: variation of the input shape // SYS2: variation of the trigger efficiency at threshold // // Finally, the unfolding is compared to a "bin-by-bin" correction method // /////////////////////////////////////////////////////////////////////// TRandom *rnd=0; Double_t GenerateEvent(const Double_t *parm, const Double_t *triggerParm, Double_t *intLumi, Bool_t *triggerFlag, Double_t *ptGen,Int_t *iType) { // generate an event // input: // parameters for the event generator // return value: // reconstructed Pt // output to pointers: // integrated luminosity // several variables only accessible on generator level // // the parm array defines the physical parameters // parm[0]: background source 1 fraction // parm[1]: background source 2 fraction // parm[2]: lower limit of generated Pt distribution // parm[3]: upper limit of generated Pt distribution // parm[4]: exponent for Pt distribution signal // parm[5]: exponent for Pt distribution background 1 // parm[6]: exponent for Pt distribution background 2 // parm[7]: resolution parameter a goes with sqrt(Pt) // parm[8]: resolution parameter b goes with Pt // triggerParm[0]: trigger threshold turn-on position // triggerParm[1]: trigger threshold turn-on width // triggerParm[2]: trigger efficiency for high pt // // intLumi is advanced bu 1 for each *generated* event // for data, several events may be generated, until one passes the trigger // // some generator-level quantities are also returned: // triggerFlag: whether the event passed the trigger threshold // ptGen: the generated pt // iType: which type of process was simulated // // the "triggered" also has another meaning: // if(triggerFlag==0) only triggered events are returned // // Usage to generate data events // ptObs=GenerateEvent(parm,triggerParm,0,0,0) // // Usage to generate MC events // ptGen=GenerateEvent(parm,triggerParm,&triggerFlag,&ptGen,&iType); // Double_t ptObs; Bool_t isTriggered=kFALSE; do { Int_t itype; Double_t ptgen; Double_t f=rnd->Rndm(); // decide whether this is background or signal itype=0; if(fRndm(); if(a1 == 0.0) { Double_t x0=TMath::Log(parm[2]); ptgen=TMath::Exp(t*(TMath::Log(parm[3])-x0)+x0); } else { Double_t x0=pow(parm[2],a1); ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1); } if(iType) *iType=itype; if(ptGen) *ptGen=ptgen; // smearing in Pt with large asymmetric tail Double_t sigma= TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen); ptObs=rnd->BreitWigner(ptgen,sigma); // decide whether event was triggered Double_t triggerProb = triggerParm[2]/(1.+TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1])); isTriggered= rnd->Rndm()SetOptFit(1111); // random generator rnd=new TRandom3(); // data and MC luminosities Double_t const lumiData= 30000; Double_t const lumiMC =1000000; // Binning // reconstructed mass Int_t const nDet=24; Double_t const xminDet=4.0; Double_t const xmaxDet=28.0; // signal contributions Int_t const nGen=10; Double_t const xminGen= 6.0; Double_t const xmaxGen=26.0; //======================== // Step 1: fill histograms // =================== data true parameters ============================ // // in real life these are unknown // static Double_t parm_data[]={ 0.05, // fraction of background 1 (on generator level) 0.05, // fraction of background 2 (on generator level) 3.5, // lower Pt cut (generator level) 100.,// upper Pt cut (generator level) -3.0,// signal exponent 0.1, // background 1 exponent -0.5, // background 2 exponent 0.2, // energy resolution a term 0.01, // energy resolution b term }; static Double_t triggerParm_data[]={8.0, // threshold is 8 GeV 4.0, // width is 4 GeV 0.95 // high Pt efficiency os 95% }; //============================================ // generate data distribution TH1D *histDetDATA=new TH1D("PT(data)",";Pt(obs)",nDet,xminDet,xmaxDet); // second data distribution, with generator-level bins // for bin-by-bin correction study TH1D *histDetDATAbbb=new TH1D("PT(data,bbb)",";Pt(obs,bbb)", nGen,xminGen,xmaxGen); // in real life we do not have this distribution... TH1D *histGenDATA=new TH1D("PT(MC,signal,gen)",";Pt(gen)", nGen,xminGen,xmaxGen); Double_t intLumi=0.0; while(intLumiFill(ptObs); histDetDATAbbb->Fill(ptObs); } // fill generator-level signal histogram for data // for real data, this does not exist if(iTypeGen==0) histGenDATA->Fill(ptGen); } // =================== MC parameters ============================ // default settings static Double_t parm_MC[]={ 0.05, // fraction of background 1 (on generator level) 0.05, // fraction of background 2 (on generator level) 3.5, // lower Pt cut (generator level) 100.,// upper Pt cut (generator level) -4.0,// signal exponent !!! UNKNOWN !!! steeper than in data // to illustrate bin-by-bin bias 0.1, // background 1 exponent -0.5, // background 2 exponent 0.2, // energy resolution a term 0.01, // energy resolution b term }; static Double_t triggerParm_MC[]={8.0, // threshold is 8 GeV 4.0, // width is 4 GeV 0.95 // high Pt efficiency os 95% }; // study trigger systematics: change parameters for trigger threshold static Double_t triggerParm_MCSYS1[]={7.7, // threshold is 7 GeV 3.7, // width is 3 GeV 0.95 // high Pt efficiency is 9% }; // study dependence on initial signal shape static Double_t parm_MC_SYS2[]={ 0.0, // fraction of background: not needed 0.0, // fraction of background: not needed 3.5, // lower Pt cut (generator level) 100.,// upper Pt cut (generator level) -2.0, // signal exponent changed 0.1, // background 1 exponent -0.5, // background 2 exponent 0.2, // energy resolution a term 0.01, // energy resolution b term }; //============================================ // generate MC distributions // detector-level histograms TH1D *histDetMC[3]; histDetMC[0]=new TH1D("PT(MC,signal)",";Pt(obs)",nDet,xminDet,xmaxDet); histDetMC[1]=new TH1D("PT(MC,bgr1)",";Pt(obs)",nDet,xminDet,xmaxDet); histDetMC[2]=new TH1D("PT(MC,bgr2)",";Pt(obs)",nDet,xminDet,xmaxDet); TH1D *histDetMCall=new TH1D("PT(MC,all)",";Pt(obs)",nDet,xminDet,xmaxDet); TH1D *histDetMCbgr=new TH1D("PT(MC,bgr)",";Pt(obs)",nDet,xminDet,xmaxDet); // another set of detector level histograms, this time with // generator-level binning, to try the bin-by-bin correction TH1D *histDetMCbbb[3]; histDetMCbbb[0]=new TH1D("PT(MC,sig)",";Pt(obs,bbb)",nGen,xminGen,xmaxGen); histDetMCbbb[1]=new TH1D("PT(MC,bgr1)bbb",";Pt(obs,bbb)",nGen,xminGen,xmaxGen); histDetMCbbb[2]=new TH1D("PT(MC,bgr2)bbb",";Pt(obs,bbb)",nGen,xminGen,xmaxGen); // true signal distribution, two different binnings TH1D *histGenMC=new TH1D("PT(MC,signal,gen)bbb",";Pt(gen)", nGen,xminGen,xmaxGen); TH2D *histGenDet=new TH2D("PT(MC,signal,gen,det)",";Pt(gen);Pt(det)", nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet); // weighting factor to accomodate for the different luminosity in data and MC Double_t lumiWeight = lumiData/lumiMC; intLumi=0.0; while(intLumiFill(ptObs,lumiWeight); histDetMCbbb[iTypeGen]->Fill(ptObs,lumiWeight); // fill total MC prediction (signal+bgr1+bgr2) histDetMCall->Fill(ptObs,lumiWeight); // fill bgr only MC prediction (bgr1+bgr2) if(iTypeGen>0) { histDetMCbgr->Fill(ptObs,lumiWeight); } } // generator level distributions (signal only) if(iTypeGen==0) { histGenMC->Fill(ptGen,lumiWeight); } // matrix dectribing the migrations (signal only) if(iTypeGen==0) { // if event was triggered, fill histogram with ptObs if(isTriggered) { histGenDet->Fill(ptGen,ptObs,lumiWeight); } else { // not triggered: fill detector-level underflow bin histGenDet->Fill(ptGen,-999.,lumiWeight); } } } // generate another MC, this time with changed trigger settings // fill 2D histogram and histograms for bin-by-bin study TH2D *histGenDetSYS1=new TH2D("PT(MC,signal,gen,det,sys1)", ";Pt(gen);Pt(obs)", nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet); TH1D *histGenSYS1=new TH1D("PT(MC,signal,gen,sys1)", ";Pt(obs)", nGen,xminGen,xmaxGen); TH1D *histDetSYS1bbb=new TH1D("PT(MC,signal,det,sys1,bbb)", ";Pt(obs)", nGen,xminGen,xmaxGen); intLumi=0.; while(intLumiFill(ptGen,ptObs,lumiWeight); } else { // not triggered: fill detector-level underflow bin histGenDetSYS1->Fill(ptGen,-999.,lumiWeight); } } // generator level distribution for bin-by-bin study if(iTypeGen==0) { histGenSYS1->Fill(ptGen,lumiWeight); // detector level for bin-by-bin study if(isTriggered) { histDetSYS1bbb->Fill(ptObs,lumiWeight); } } } // generate a third MC, this time with changed initial distribution // only the 2D histogram is filled TH2D *histGenDetSYS2=new TH2D("PT(MC,signal,gen,det,sys2)", ";Pt(gen);Pt(det)", nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet); TH1D *histGenSYS2=new TH1D("PT(MC,signal,gen,sys2)", ";Pt(obs)", nGen,xminGen,xmaxGen); TH1D *histDetSYS2bbb=new TH1D("PT(MC,signal,det,sys2,bbb)", ";Pt(obs)", nGen,xminGen,xmaxGen); intLumi=0.0; while(intLumiFill(ptGen,ptObs,lumiWeight); } else { // not triggered: fill detector-level underflow bin histGenDetSYS2->Fill(ptGen,-999.,lumiWeight); } } if(iTypeGen==0) { // generator level distribution for bin-by-bin study histGenSYS2->Fill(ptGen,lumiWeight); // detector level for bin-by-bin study if(isTriggered) { histDetSYS2bbb->Fill(ptObs,lumiWeight); } } } //======================== // Step 2: unfolding // // this is based on a matrix (TH2D) which describes the connection // between // nDet Detector bins // nGen+2 Generator level bins // // why are there nGen+2 Generator level bins? // the two extra bins are the underflow and overflow bin // These are unfolded as well. Such bins are needed to // accomodate for migrations from outside the phase-space into the // observed phase-space. // // In addition there are overflow/underflow bins // for the reconstructed variable. These are used to count events // which are -NOT- reconstructed. Either because they migrate // out of the reconstructed phase-space. Or because there are detector // inefficiencies. // TUnfoldSys unfold(histGenDet,TUnfold::kHistMapOutputHoriz, TUnfold::kRegModeSize,TUnfold::kEConstraintNone); // define the input vector (the measured data distribution) unfold.SetInput(histDetDATA); // subtract background, normalized to data luminosity // with 10% scale error each Double_t scale_bgr=1.0; Double_t dscale_bgr=0.2; unfold.SubtractBackground(histDetMC[1],"background1",scale_bgr,dscale_bgr); unfold.SubtractBackground(histDetMC[2],"background2",scale_bgr,dscale_bgr); // add systematic errors // trigger threshold unfold.AddSysError(histGenDetSYS1,"trigger_SYS1", TUnfold::kHistMapOutputHoriz, TUnfoldSys::kSysErrModeMatrix); // calorimeter response unfold.AddSysError(histGenDetSYS2,"signalshape_SYS2", TUnfold::kHistMapOutputHoriz, TUnfoldSys::kSysErrModeMatrix); // run the unfolding Int_t nScan=30; TSpline *logTauX,*logTauY; TGraph *lCurve; // this method scans the parameter tau and finds the kink in the L curve // finally, the unfolding is done for the best choice of tau Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY); // save graphs with one point to visualize best choice of tau Double_t t[1],x[1],y[1]; logTauX->GetKnot(iBest,t[0],x[0]); logTauY->GetKnot(iBest,t[0],y[0]); TGraph *bestLcurve=new TGraph(1,x,y); TGraph *bestLogTauLogChi2=new TGraph(1,t,x); // get unfolding output // Reminder: there are // nGen+2 Generator level bins // But we are not interested in the additional bins. // The mapping is required to get the proper errors and correlations // Here a mapping is defined from the nGen+2 bins to nGen bins. Int_t *binMap=new Int_t[nGen+2]; binMap[0] = -1; // underflow bin should be ignored binMap[nGen+1]=-1; // overflow bin should be ignored for(Int_t i=1;i<=nGen;i++) { binMap[i]=i; // map inner bins to output histogram bins } // this returns the unfolding output // The errors included at this point are: // * statistical errros // * background subtraction errors TH1D *histUnfoldStatBgr=new TH1D("PT(unfold,signal,stat+bgr)",";Pt(gen)", nGen,xminGen,xmaxGen); unfold.GetOutput(histUnfoldStatBgr,binMap); // retreive error matrix of statistical errors only TH2D *histEmatStat=new TH2D("ErrMat(stat)",";Pt(gen);Pt(gen)", nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen); unfold.GetEmatrixInput(histEmatStat,binMap); // retreive full error matrix // This includes also all systematic errors defined above TH2D *histEmatTotal=new TH2D("ErrMat(total)",";Pt(gen);Pt(gen)", nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen); unfold.GetEmatrixTotal(histEmatTotal,binMap); // create two copies of the unfolded data, one with statistical errors // the other with total errors TH1D *histUnfoldStat=new TH1D("PT(unfold,signal,stat)",";Pt(gen)", nGen,xminGen,xmaxGen); TH1D *histUnfoldTotal=new TH1D("PT(unfold,signal,total)",";Pt(gen)", nGen,xminGen,xmaxGen); for(Int_t i=0;iGetBinContent(i); // histogram with unfolded data and stat errors histUnfoldStat->SetBinContent(i,c); histUnfoldStat->SetBinError (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i))); // histogram with unfolded data and total errors histUnfoldTotal->SetBinContent(i,c); histUnfoldTotal->SetBinError (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i))); } // create histogram with correlation matrix TH2D *histCorr=new TH2D("Corr(total)",";Pt(gen);Pt(gen)", nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen); for(Int_t i=0;iGetBinContent(i,i)); if(ei<=0.0) continue; for(Int_t j=0;jGetBinContent(j,j)); if(ej<=0.0) continue; histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej); } } delete [] binMap; //======================== // Step 3: plots TCanvas output; output.Divide(3,2); output.cd(1); histDetDATA->SetMinimum(0.0); histDetDATA->Draw("E"); histDetMCall->SetMinimum(0.0); histDetMCall->SetLineColor(kBlue); histDetMCbgr->SetLineColor(kRed); histDetMC[1]->SetLineColor(kCyan); histDetMCall->Draw("SAME HIST"); histDetMCbgr->Draw("SAME HIST"); histDetMC[1]->Draw("SAME HIST"); output.cd(2); histUnfoldTotal->SetMinimum(0.0); histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5); // outer error: total error histUnfoldTotal->Draw("E"); // middle error: stat+bgr histUnfoldStatBgr->Draw("SAME E1"); // inner error: stat only histUnfoldStat->Draw("SAME E1"); histGenDATA->Draw("SAME HIST"); histGenMC->SetLineColor(kBlue); histGenMC->Draw("SAME HIST"); output.cd(3); histGenDet->SetLineColor(kBlue); histGenDet->Draw("BOX"); // show tau as a function of chi**2 output.cd(4); logTauX->Draw(); bestLogTauLogChi2->SetMarkerColor(kRed); bestLogTauLogChi2->Draw("*"); // show the L curve output.cd(5); lCurve->Draw("AL"); bestLcurve->SetMarkerColor(kRed); bestLcurve->Draw("*"); // show correlation matrix output.cd(6); histCorr->Draw("BOX"); output.SaveAs("testUnfold3.ps"); // step 4: compare results to // the so-called bin-by-bin "correction" for(Int_t i=1;i<=nGen;i++) { // data contribution in this bin Double_t data=histDetDATAbbb->GetBinContent(i); Double_t errData=histDetDATAbbb->GetBinError(i); // subtract background contribution Double_t data_bgr=data; Double_t errData_bgr=errData; for(Int_t j=1;j<=2;j++) { data_bgr -= scale_bgr*histDetMCbbb[j]->GetBinContent(i); errData_bgr = TMath::Sqrt(errData_bgr*errData_bgr+ scale_bgr*histDetMCbbb[j]->GetBinError(i)* scale_bgr*histDetMCbbb[j]->GetBinError(i)+ dscale_bgr*histDetMCbbb[j]->GetBinContent(i)* dscale_bgr*histDetMCbbb[j]->GetBinContent(i) ); } // "correct" the data, using the Monte Carlo and neglecting off-diagonals Double_t fCorr=(histGenMC->GetBinContent(i)/ histDetMCbbb[0]->GetBinContent(i)); Double_t data_bbb= data_bgr *fCorr; // stat only error Double_t errData_stat_bbb = errData*fCorr; // stat plus background subtraction Double_t errData_statbgr_bbb = errData_bgr*fCorr; // estimate systematic error by repeating the exercise // using the MC with systematic shifts applied Double_t fCorr_1=(histGenSYS1->GetBinContent(i)/ histDetSYS1bbb->GetBinContent(i)); Double_t shift_sys1= data_bgr*fCorr_1 - data_bbb; Double_t fCorr_2=(histGenSYS2->GetBinContent(i)/ histDetSYS2bbb->GetBinContent(i)); Double_t shift_sys2= data_bgr*fCorr_2 - data_bbb; cout<GetBinContent(i); Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i); Double_t errData_statbgr_unfold=histUnfoldStatBgr->GetBinError(i); Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i); // compare std::cout<<"Bin "<GetBinContent(i) <<" unfold: "<GetBinContent(i) <<" binbybin: "<