TUnfold: testUnfold5a.C Source File

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 example for multi-dimensional unfolding
00013 //
00014 
00015 #include <iostream>
00016 #include <map>
00017 #include <cmath>
00018 #include <TMath.h>
00019 #include <TRandom3.h>
00020 #include <TFile.h>
00021 #include <TTree.h>
00022 
00023 using namespace std;
00024 
00025 /*
00026   This file is part of TUnfold.
00027 
00028   TUnfold is free software: you can redistribute it and/or modify
00029   it under the terms of the GNU General Public License as published by
00030   the Free Software Foundation, either version 3 of the License, or
00031   (at your option) any later version.
00032 
00033   TUnfold is distributed in the hope that it will be useful,
00034   but WITHOUT ANY WARRANTY; without even the implied warranty of
00035   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00036   GNU General Public License for more details.
00037 
00038   You should have received a copy of the GNU General Public License
00039   along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.
00040 */
00041 
00042 ///////////////////////////////////////////////////////////////////////
00043 // 
00044 // Test program for the classes TUnfoldDensity and TUnfoldBinning
00045 //
00046 // A toy test of the TUnfold package
00047 //
00048 // This is an example of unfolding a two-dimensional distribution
00049 // also using an auxillary measurement to constrain some background
00050 //
00051 // The example comprizes several macros
00052 //   testUnfold5a.C   create root files with TTree objects for
00053 //                      signal, background and data
00054 //            -> write files  testUnfold5_signal.root
00055 //                            testUnfold5_background.root
00056 //                            testUnfold5_data.root
00057 //
00058 //   testUnfold5b.C   create a root file with the TUnfoldBinning objects
00059 //            -> write file  testUnfold5_binning.root
00060 //
00061 //   testUnfold5c.C   loop over trees and fill histograms based on the
00062 //                      TUnfoldBinning objects
00063 //            -> read  testUnfold5_binning.root
00064 //                     testUnfold5_signal.root
00065 //                     testUnfold5_background.root
00066 //                     testUnfold5_data.root
00067 //
00068 //            -> write testUnfold5_histograms.root
00069 //
00070 //   testUnfold5d.C   run the unfolding
00071 //            -> read  testUnfold5_histograms.root
00072 //            -> write testUnfold5_result.root
00073 //                     testUnfold5_result.ps
00074 // 
00075 ///////////////////////////////////////////////////////////////////////
00076 
00077 TRandom *g_rnd=0;
00078 
00079 class ToyEvent {
00080 public:
00081    void GenerateDataEvent(TRandom *rnd);
00082    void GenerateSignalEvent(TRandom *rnd);
00083    void GenerateBgrEvent(TRandom *rnd);
00084    // reconstructed quantities
00085    inline Double_t GetPtRec(void) const { return fPtRec; }
00086    inline Double_t GetEtaRec(void) const { return fEtaRec; }
00087    inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
00088    inline Bool_t IsTriggered(void) const { return fIsTriggered; }
00089 
00090    // generator level quantities
00091    inline Double_t GetPtGen(void) const { 
00092       if(IsSignal()) return fPtGen;
00093       else return -1.0;
00094    }
00095    inline Double_t GetEtaGen(void) const {
00096        if(IsSignal()) return fEtaGen;
00097        else return 999.0;
00098    }
00099    inline Bool_t IsSignal(void) const { return fIsSignal; }
00100 protected:
00101 
00102    void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
00103    void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
00104    void GenerateReco(TRandom *rnd);
00105 
00106    // reconstructed quantities
00107    Double_t fPtRec;
00108    Double_t fEtaRec;
00109    Double_t fDiscriminator;
00110    Bool_t fIsTriggered;
00111    // generated quantities
00112    Double_t fPtGen;
00113    Double_t fEtaGen;
00114    Bool_t fIsSignal;
00115 
00116    static Double_t kDataSignalFraction;
00117 
00118 };
00119 
00120 void testUnfold5a()
00121 {
00122   // random generator
00123   g_rnd=new TRandom3();
00124 
00125   // data and MC number of events
00126   const Int_t neventData         =  20000;
00127   Double_t const neventSignalMC  =2000000;
00128   Double_t const neventBgrMC     =2000000;
00129 
00130   Float_t etaRec,ptRec,discr,etaGen,ptGen;
00131   Int_t istriggered,issignal;
00132 
00133   //==================================================================  
00134   // Step 1: generate data TTree
00135 
00136   TFile *dataFile=new TFile("testUnfold5_data.root","recreate");
00137   TTree *dataTree=new TTree("data","event");
00138 
00139   dataTree->Branch("etarec",&etaRec,"etarec/F");
00140   dataTree->Branch("ptrec",&ptRec,"ptrec/F");
00141   dataTree->Branch("discr",&discr,"discr/F");
00142 
00143   // for real data, only the triggered events are available
00144   dataTree->Branch("istriggered",&istriggered,"istriggered/I");
00145   // data truth parameters
00146   dataTree->Branch("etagen",&etaGen,"etagen/F");
00147   dataTree->Branch("ptgen",&ptGen,"ptgen/F");
00148   dataTree->Branch("issignal",&issignal,"issignal/I");
00149 
00150   cout<<"fill data tree\n";
00151 
00152   Int_t nEvent=0,nTriggered=0;
00153   while(nTriggered<neventData) {
00154      ToyEvent event;
00155      event.GenerateDataEvent(g_rnd);
00156 
00157      etaRec=event.GetEtaRec();
00158      ptRec=event.GetPtRec();
00159      discr=event.GetDiscriminator();
00160      istriggered=event.IsTriggered() ? 1 : 0;
00161      etaGen=event.GetEtaGen();
00162      ptGen=event.GetPtGen();
00163      issignal=event.IsSignal() ? 1 : 0;
00164 
00165      dataTree->Fill();
00166 
00167      if(!(nEvent%100000)) cout<<"   data event "<<nEvent<<"\n";
00168 
00169      if(istriggered) nTriggered++;
00170      nEvent++;
00171 
00172   }
00173 
00174   dataTree->Write();
00175   delete dataTree;
00176   delete dataFile;
00177 
00178   //==================================================================  
00179   // Step 2: generate signal TTree
00180 
00181   TFile *signalFile=new TFile("testUnfold5_signal.root","recreate");
00182   TTree *signalTree=new TTree("signal","event");
00183 
00184   signalTree->Branch("etarec",&etaRec,"etarec/F");
00185   signalTree->Branch("ptrec",&ptRec,"ptrec/F");
00186   signalTree->Branch("discr",&discr,"discr/F");
00187   signalTree->Branch("istriggered",&istriggered,"istriggered/I");
00188   signalTree->Branch("etagen",&etaGen,"etagen/F");
00189   signalTree->Branch("ptgen",&ptGen,"ptgen/F");
00190 
00191   cout<<"fill signal tree\n";
00192   
00193   for(int ievent=0;ievent<neventSignalMC;ievent++) {
00194      ToyEvent event;
00195      event.GenerateSignalEvent(g_rnd);
00196 
00197      etaRec=event.GetEtaRec();
00198      ptRec=event.GetPtRec();
00199      discr=event.GetDiscriminator();
00200      istriggered=event.IsTriggered() ? 1 : 0;
00201      etaGen=event.GetEtaGen();
00202      ptGen=event.GetPtGen();
00203 
00204      if(!(ievent%100000)) cout<<"   signal event "<<ievent<<"\n";
00205 
00206      signalTree->Fill();
00207   }
00208 
00209   signalTree->Write();
00210   delete signalTree;
00211   delete signalFile;
00212 
00213   // ==============================================================
00214   // Step 3: generate background MC TTree
00215 
00216   TFile *bgrFile=new TFile("testUnfold5_background.root","recreate");
00217   TTree *bgrTree=new TTree("background","event");
00218 
00219   bgrTree->Branch("etarec",&etaRec,"etarec/F");
00220   bgrTree->Branch("ptrec",&ptRec,"ptrec/F");
00221   bgrTree->Branch("discr",&discr,"discr/F");
00222   bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
00223 
00224   cout<<"fill background tree\n";
00225 
00226   for(int ievent=0;ievent<neventBgrMC;ievent++) {
00227      ToyEvent event;
00228      event.GenerateBgrEvent(g_rnd);
00229      etaRec=event.GetEtaRec();
00230      ptRec=event.GetPtRec();
00231      discr=event.GetDiscriminator();
00232      istriggered=event.IsTriggered() ? 1 : 0;
00233 
00234      if(!(ievent%100000)) cout<<"   background event "<<ievent<<"\n";
00235 
00236      bgrTree->Fill();
00237   }
00238 
00239   bgrTree->Write();
00240   delete bgrTree;
00241   delete bgrFile;
00242 }
00243 
00244 Double_t ToyEvent::kDataSignalFraction=0.8;
00245 
00246 void ToyEvent::GenerateDataEvent(TRandom *rnd) {
00247    fIsSignal=rnd->Uniform()<kDataSignalFraction;
00248    if(IsSignal()) {
00249       GenerateSignalKinematics(rnd,kTRUE);
00250    } else {
00251       GenerateBgrKinematics(rnd,kTRUE);
00252    }
00253    GenerateReco(rnd);
00254 }
00255 
00256 void ToyEvent::GenerateSignalEvent(TRandom *rnd) {
00257    fIsSignal=1;
00258    GenerateSignalKinematics(rnd,kFALSE);
00259    GenerateReco(rnd);
00260 }
00261 
00262 void ToyEvent::GenerateBgrEvent(TRandom *rnd) {
00263    fIsSignal=0;
00264    GenerateBgrKinematics(rnd,kFALSE);
00265    GenerateReco(rnd);
00266 }
00267 
00268 void ToyEvent::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
00269    Double_t e_T0=0.5;
00270    Double_t e_T0_eta=0.0;
00271    Double_t e_n=2.0;
00272    Double_t e_n_eta=0.0;
00273    Double_t eta_p2=0.0;
00274    Double_t etaMax=4.0;
00275    if(isData) {
00276       e_T0=0.6;
00277       e_n=2.5;
00278       e_T0_eta=0.05;
00279       e_n_eta=-0.05;
00280       eta_p2=1.5;
00281    }
00282    if(eta_p2>0.0) {
00283       fEtaGen=TMath::Power(rnd->Uniform(),eta_p2)*etaMax;
00284       if(rnd->Uniform()>=0.5) fEtaGen= -fEtaGen;
00285    } else {
00286       fEtaGen=rnd->Uniform(-etaMax,etaMax);
00287    }
00288    Double_t n=e_n   + e_n_eta*fEtaGen;
00289    Double_t T0=e_T0 + e_T0_eta*fEtaGen;
00290    fPtGen=(TMath::Power(rnd->Uniform(),1./(1.-n))-1.)*T0;
00291    /*   static int print=100;
00292       if(print) {
00293          cout<<fEtaGen
00294             <<" "<<fPtGen
00295              <<"\n";
00296          print--;
00297          } */
00298 }
00299 
00300 void ToyEvent::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
00301    fPtGen=0.0;
00302    fEtaGen=0.0;
00303    fPtRec=rnd->Exp(2.5);
00304    fEtaRec=rnd->Uniform(-3.,3.);
00305 }
00306 
00307 void ToyEvent::GenerateReco(TRandom *rnd) {
00308    if(fIsSignal) {
00309       Double_t expEta=TMath::Exp(fEtaGen);
00310       Double_t eGen=fPtGen*(expEta+1./expEta);
00311       Double_t sigmaE=0.1*TMath::Sqrt(eGen)+(0.01+0.002*TMath::Abs(fEtaGen))
00312          *eGen;
00313       Double_t eRec;
00314       do {
00315          eRec=rnd->Gaus(eGen,sigmaE);
00316       } while(eRec<=0.0);
00317       Double_t sigmaEta=0.1+0.02*fEtaGen;
00318       fEtaRec=rnd->Gaus(fEtaGen,sigmaEta);
00319       fPtRec=eRec/(expEta+1./expEta);
00320       do {
00321          Double_t tauDiscr=0.08-0.04/(1.+fPtRec/10.0);
00322          Double_t sigmaDiscr=0.01;
00323          fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
00324       } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
00325       /* static int print=100;
00326          if(print) {
00327          cout<<fEtaGen<<" "<<fPtGen
00328              <<" -> "<<fEtaRec<<" "<<fPtRec
00329              <<"\n";
00330          print--;
00331          } */
00332    } else {
00333       do {
00334          Double_t tauDiscr=0.15-0.05/(1.+fPtRec/5.0)+0.1*fEtaRec;
00335          Double_t sigmaDiscr=0.02+0.01*fEtaRec;
00336          fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
00337       } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
00338    }
00339    fIsTriggered=(rnd->Uniform()<1./(TMath::Exp(-fPtRec+3.5)+1.));
00340 }