TUnfold: testUnfold7a.C Source File

00001 // Author: Stefan Schmitt
00002 // DESY, July 2016
00003 
00004 //  Version 17.6, example study of unfolding methods
00005 
00006 #include <iostream>
00007 #include <map>
00008 #include <cmath>
00009 #include <TMath.h>
00010 #include <TRandom3.h>
00011 #include <TFile.h>
00012 #include <TTree.h>
00013 #include <TLorentzVector.h>
00014 
00015 #define MASS1 0.511E-3
00016 
00017 using namespace std;
00018 
00019 /*
00020   This file is part of TUnfold.
00021 
00022   TUnfold is free software: you can redistribute it and/or modify
00023   it under the terms of the GNU General Public License as published by
00024   the Free Software Foundation, either version 3 of the License, or
00025   (at your option) any later version.
00026 
00027   TUnfold is distributed in the hope that it will be useful,
00028   but WITHOUT ANY WARRANTY; without even the implied warranty of
00029   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00030   GNU General Public License for more details.
00031 
00032   You should have received a copy of the GNU General Public License
00033   along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.
00034 */
00035 
00036 ///////////////////////////////////////////////////////////////////////
00037 // 
00038 // Test program for the classes TUnfoldDensity and TUnfoldBinning
00039 //
00040 // A toy test of the TUnfold package
00041 //
00042 // This example is documented in conference proceedings
00043 //   arXiv:1611.01927
00044 // 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
00045 //
00046 // This is an example of unfolding a two-dimensional distribution
00047 // also using an auxillary measurement to constrain some background
00048 //
00049 // The example comprizes several macros
00050 //   testUnfold7a.C   create root files with TTree objects for
00051 //                      signal, background and data
00052 //            -> write files  testUnfold7_signal.root
00053 //                            testUnfold7_background.root
00054 //                            testUnfold7_data.root
00055 //
00056 //   testUnfold7b.C   loop over trees and fill histograms based on the
00057 //                      TUnfoldBinning objects
00058 //            -> read  testUnfold7binning.xml
00059 //                     testUnfold7_signal.root
00060 //                     testUnfold7_background.root
00061 //                     testUnfold7_data.root
00062 //
00063 //            -> write testUnfold7_histograms.root
00064 //
00065 //   testUnfold7c.C   run the unfolding
00066 //            -> read  testUnfold7_histograms.root
00067 //            -> write testUnfold7_result.root
00068 //                     testUnfold7_result.ps
00069 // 
00070 ///////////////////////////////////////////////////////////////////////
00071 
00072 TRandom *g_rnd=0;
00073 
00074 class ToyEvent7 {
00075 public:
00076    void GenerateDataEvent(TRandom *rnd);
00077    void GenerateSignalEvent(TRandom *rnd);
00078    void GenerateBgrEvent(TRandom *rnd);
00079    // reconstructed quantities
00080    inline Double_t GetMRec(int i) const { return fMRec[i]; }
00081    inline Double_t GetPtRec(int i) const { return fPtRec[i]; }
00082    inline Double_t GetEtaRec(int i) const { return fEtaRec[i]; }
00083    inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
00084    inline Double_t GetPhiRec(int i) const { return fPhiRec[i]; }
00085    inline Bool_t IsTriggered(void) const { return fIsTriggered; }
00086 
00087    // generator level quantities
00088    inline Double_t GetMGen(int i) const { 
00089       if(IsSignal()) return fMGen[i];
00090       else return -1.0;
00091    }
00092    inline Double_t GetPtGen(int i) const { 
00093       if(IsSignal()) return fPtGen[i];
00094       else return -1.0;
00095    }
00096    inline Double_t GetEtaGen(int i) const {
00097        if(IsSignal()) return fEtaGen[i];
00098        else return 999.0;
00099    }
00100    inline Double_t GetPhiGen(int i) const {
00101        if(IsSignal()) return fPhiGen[i];
00102        else return 999.0;
00103    }
00104    inline Bool_t IsSignal(void) const { return fIsSignal; }
00105 protected:
00106 
00107    void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
00108    void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
00109    void GenerateReco(TRandom *rnd);
00110 
00111    // reconstructed quantities
00112    Double_t fMRec[3];
00113    Double_t fPtRec[3];
00114    Double_t fEtaRec[3];
00115    Double_t fPhiRec[3];
00116    Double_t fDiscriminator;
00117    Bool_t fIsTriggered;
00118    // generated quantities
00119    Double_t fMGen[3];
00120    Double_t fPtGen[3];
00121    Double_t fEtaGen[3];
00122    Double_t fPhiGen[3];
00123    Bool_t fIsSignal;
00124 public:
00125    static Double_t kDataSignalFraction;
00126    static Double_t kMCSignalFraction;
00127 
00128 };
00129 
00130 void testUnfold7a()
00131 {
00132   // random generator
00133   g_rnd=new TRandom3(4711);
00134 
00135   // data and MC number of events
00136   Double_t muData0=5000.;
00137   // luminosity error
00138   Double_t muData=muData0*g_rnd->Gaus(1.0,0.03);
00139   // stat error
00140   Int_t neventData      = g_rnd->Poisson( muData);
00141 
00142   // generated number of MC events
00143   Int_t neventSigmc     = 250000;
00144   Int_t neventBgrmc     = 100000;
00145 
00146   Float_t etaRec[3],ptRec[3],phiRec[3],mRec[3],discr;
00147   Float_t etaGen[3],ptGen[3],phiGen[3],mGen[3];
00148   Float_t weight;
00149   Int_t istriggered,issignal;
00150 
00151   //==================================================================  
00152   // Step 1: generate data TTree
00153 
00154   TFile *dataFile=new TFile("testUnfold7_data.root","recreate");
00155   TTree *dataTree=new TTree("data","event");
00156 
00157   dataTree->Branch("etarec",etaRec,"etarec[3]/F");
00158   dataTree->Branch("ptrec",ptRec,"ptrec[3]/F");
00159   dataTree->Branch("phirec",phiRec,"phirec[3]/F");
00160   dataTree->Branch("mrec",mRec,"mrec[3]/F");
00161   dataTree->Branch("discr",&discr,"discr/F");
00162 
00163   // for real data, only the triggered events are available
00164   dataTree->Branch("istriggered",&istriggered,"istriggered/I");
00165   // data truth parameters
00166   dataTree->Branch("etagen",etaGen,"etagen[3]/F");
00167   dataTree->Branch("ptgen",ptGen,"ptgen[3]/F");
00168   dataTree->Branch("phigen",phiGen,"phigen[3]/F");
00169   dataTree->Branch("mgen",mGen,"mgen[3]/F");
00170   dataTree->Branch("issignal",&issignal,"issignal/I");
00171 
00172   cout<<"fill data tree\n";
00173 
00174   //Int_t nEvent=0,nTriggered=0;
00175   for(int ievent=0;ievent<neventData;ievent++) {
00176      ToyEvent7 event;
00177      event.GenerateDataEvent(g_rnd);
00178      for(int i=0;i<3;i++) {
00179         etaRec[i]=event.GetEtaRec(i);
00180         ptRec[i]=event.GetPtRec(i);
00181         phiRec[i]=event.GetPhiRec(i);
00182         mRec[i]=event.GetMRec(i);
00183         etaGen[i]=event.GetEtaGen(i);
00184         ptGen[i]=event.GetPtGen(i);
00185         phiGen[i]=event.GetPhiGen(i);
00186         mGen[i]=event.GetMGen(i);
00187      }
00188      discr=event.GetDiscriminator();
00189      istriggered=event.IsTriggered() ? 1 : 0;
00190      issignal=event.IsSignal() ? 1 : 0;
00191 
00192      dataTree->Fill();
00193 
00194      if(!(ievent%100000)) cout<<"   data event "<<ievent<<"\n";
00195 
00196      //if(istriggered) nTriggered++;
00197      //nEvent++;
00198 
00199   }
00200 
00201   dataTree->Write();
00202   delete dataTree;
00203   delete dataFile;
00204 
00205   //==================================================================  
00206   // Step 2: generate signal TTree
00207 
00208   TFile *signalFile=new TFile("testUnfold7_signal.root","recreate");
00209   TTree *signalTree=new TTree("signal","event");
00210 
00211   signalTree->Branch("etarec",etaRec,"etarec[3]/F");
00212   signalTree->Branch("ptrec",ptRec,"ptrec[3]/F");
00213   signalTree->Branch("phirec",ptRec,"phirec[3]/F");
00214   signalTree->Branch("mrec",mRec,"mrec[3]/F");
00215   signalTree->Branch("discr",&discr,"discr/F");
00216 
00217   // for real data, only the triggered events are available
00218   signalTree->Branch("istriggered",&istriggered,"istriggered/I");
00219   // data truth parameters
00220   signalTree->Branch("etagen",etaGen,"etagen[3]/F");
00221   signalTree->Branch("ptgen",ptGen,"ptgen[3]/F");
00222   signalTree->Branch("phigen",phiGen,"phigen[3]/F");
00223   signalTree->Branch("weight",&weight,"weight/F");
00224   signalTree->Branch("mgen",mGen,"mgen[3]/F");
00225 
00226   cout<<"fill signal tree\n";
00227 
00228   weight=ToyEvent7::kMCSignalFraction*muData0/neventSigmc;
00229   
00230   for(int ievent=0;ievent<neventSigmc;ievent++) {
00231      ToyEvent7 event;
00232      event.GenerateSignalEvent(g_rnd);
00233 
00234      for(int i=0;i<3;i++) {
00235         etaRec[i]=event.GetEtaRec(i);
00236         ptRec[i]=event.GetPtRec(i);
00237         phiRec[i]=event.GetPhiRec(i);
00238         mRec[i]=event.GetMRec(i);
00239         etaGen[i]=event.GetEtaGen(i);
00240         ptGen[i]=event.GetPtGen(i);
00241         phiGen[i]=event.GetPhiGen(i);
00242         mGen[i]=event.GetMGen(i);
00243      }
00244      discr=event.GetDiscriminator();
00245      istriggered=event.IsTriggered() ? 1 : 0;
00246 
00247      if(!(ievent%100000)) cout<<"   signal event "<<ievent<<"\n";
00248 
00249      signalTree->Fill();
00250   }
00251 
00252   signalTree->Write();
00253   delete signalTree;
00254   delete signalFile;
00255 
00256   // ==============================================================
00257   // Step 3: generate background MC TTree
00258 
00259   TFile *bgrFile=new TFile("testUnfold7_background.root","recreate");
00260   TTree *bgrTree=new TTree("background","event");
00261 
00262   bgrTree->Branch("etarec",&etaRec,"etarec[3]/F");
00263   bgrTree->Branch("ptrec",&ptRec,"ptrec[3]/F");
00264   bgrTree->Branch("phirec",&phiRec,"phirec[3]/F");
00265   bgrTree->Branch("mrec",&mRec,"mrec[3]/F");
00266   bgrTree->Branch("discr",&discr,"discr/F");
00267   bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
00268   signalTree->Branch("weight",&weight,"weight/F");
00269 
00270   cout<<"fill background tree\n";
00271 
00272   weight=(1.-ToyEvent7::kMCSignalFraction)*muData0/neventBgrmc;
00273 
00274   for(int ievent=0;ievent<neventBgrmc;ievent++) {
00275      ToyEvent7 event;
00276      event.GenerateBgrEvent(g_rnd);
00277      for(int i=0;i<3;i++) {
00278         etaRec[i]=event.GetEtaRec(i);
00279         ptRec[i]=event.GetPtRec(i);
00280         phiRec[i]=event.GetPhiRec(i);
00281      }
00282      discr=event.GetDiscriminator();
00283      istriggered=event.IsTriggered() ? 1 : 0;
00284 
00285      if(!(ievent%100000)) cout<<"   background event "<<ievent<<"\n";
00286 
00287      bgrTree->Fill();
00288   }
00289 
00290   bgrTree->Write();
00291   delete bgrTree;
00292   delete bgrFile;
00293 }
00294 
00295 Double_t ToyEvent7::kDataSignalFraction=0.75;
00296 Double_t ToyEvent7::kMCSignalFraction=0.75;
00297 
00298 void ToyEvent7::GenerateDataEvent(TRandom *rnd) {
00299    fIsSignal=rnd->Uniform()<kDataSignalFraction;
00300    if(IsSignal()) {
00301       GenerateSignalKinematics(rnd,kTRUE);
00302    } else {
00303       GenerateBgrKinematics(rnd,kTRUE);
00304    }
00305    GenerateReco(rnd);
00306 }
00307 
00308 void ToyEvent7::GenerateSignalEvent(TRandom *rnd) {
00309    fIsSignal=1;
00310    GenerateSignalKinematics(rnd,kFALSE);
00311    GenerateReco(rnd);
00312 }
00313 
00314 void ToyEvent7::GenerateBgrEvent(TRandom *rnd) {
00315    fIsSignal=0;
00316    GenerateBgrKinematics(rnd,kFALSE);
00317    GenerateReco(rnd);
00318 }
00319 
00320 void ToyEvent7::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
00321 
00322    // fake decay of Z0 to two fermions
00323    double M0=91.1876;
00324    double Gamma=2.4952;
00325    // generated mass
00326    do {
00327       fMGen[2]=rnd->BreitWigner(M0,Gamma);
00328    } while(fMGen[2]<=0.0);
00329 
00330    double N_ETA=3.0;
00331    double MU_PT=5.;
00332    double SIGMA_PT=2.0;
00333    double DECAY_A=0.2;
00334    if(isData) {
00335       //N_ETA=2.5;
00336       MU_PT=6.;
00337       SIGMA_PT=1.8;
00338       //DECAY_A=0.5;
00339    }
00340    fEtaGen[2]=TMath::Power(rnd->Uniform(0,1.5),N_ETA);
00341    if(rnd->Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
00342    fPhiGen[2]=rnd->Uniform(-M_PI,M_PI);
00343    do {
00344       fPtGen[2]=rnd->Landau(MU_PT,SIGMA_PT);
00345    } while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
00346    //========================== decay
00347    TLorentzVector sum;
00348    sum.SetPtEtaPhiM(fPtGen[2],fEtaGen[2],fPhiGen[2],fMGen[2]);
00349    // boost into lab-frame
00350    TVector3 boost=sum.BoostVector();
00351    // decay in rest-frame
00352 
00353    TLorentzVector p[3];
00354    double m=MASS1;
00355    double costh;
00356    do {
00357       double r=rnd->Uniform(-1.,1.);
00358       costh=r*(1.+DECAY_A*r*r);
00359    } while(fabs(costh)>=1.0);
00360    double phi=rnd->Uniform(-M_PI,M_PI);
00361    double e=0.5*sum.M();
00362    double ptot=TMath::Sqrt(e+m)*TMath::Sqrt(e-m);
00363    double pz=ptot*costh;
00364    double pt=TMath::Sqrt(ptot+pz)*TMath::Sqrt(ptot-pz);
00365    double px=pt*cos(phi);
00366    double py=pt*sin(phi);
00367    p[0].SetXYZT(px,py,pz,e);
00368    p[1].SetXYZT(-px,-py,-pz,e);
00369    for(int i=0;i<2;i++) {
00370       p[i].Boost(boost);
00371    }
00372    p[2]=p[0]+p[1];
00373    for(int i=0;i<3;i++) {
00374       fPtGen[i]=p[i].Pt();
00375       fEtaGen[i]=p[i].Eta();
00376       fPhiGen[i]=p[i].Phi();
00377       fMGen[i]=p[i].M();
00378    }
00379 }
00380 
00381 void ToyEvent7::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
00382    for(int i=0;i<3;i++) {
00383       fPtGen[i]=0.0;
00384       fEtaGen[i]=0.0;
00385       fPhiGen[i]=0.0;
00386    }
00387    TLorentzVector p[3];
00388    for(int i=0;i<2;i++) {
00389       p[i].SetPtEtaPhiM(rnd->Exp(15.0),rnd->Uniform(-3.,3.),
00390                         rnd->Uniform(-M_PI,M_PI),MASS1);
00391    }
00392    p[2]=p[0]+p[1];
00393    for(int i=0;i<3;i++) {
00394       fPtRec[i]=p[i].Pt();
00395       fEtaRec[i]=p[i].Eta();
00396       fPhiRec[i]=p[i].Phi();
00397       fMRec[i]=p[i].M();
00398    }
00399 }
00400 
00401 void ToyEvent7::GenerateReco(TRandom *rnd) {
00402    if(fIsSignal) {
00403       TLorentzVector p[3];
00404       for(int i=0;i<2;i++) {
00405          Double_t expEta=TMath::Exp(fEtaGen[i]);
00406          Double_t coshEta=(expEta+1./expEta);
00407          Double_t eGen=fPtGen[i]*coshEta;
00408          Double_t sigmaE=
00409             0.1*TMath::Sqrt(eGen)
00410             +1.0*coshEta
00411             +0.01*eGen;
00412          Double_t eRec;
00413          do {
00414             eRec=rnd->Gaus(eGen,sigmaE);
00415          } while(eRec<=0.0);
00416          Double_t sigmaEta=0.1+0.02*TMath::Abs(fEtaGen[i]);
00417          p[i].SetPtEtaPhiM(eRec/(expEta+1./expEta),
00418                            rnd->Gaus(fEtaGen[i],sigmaEta),
00419                            remainder(rnd->Gaus(fPhiGen[i],0.03),2.*M_PI),
00420                            MASS1);
00421       }
00422       p[2]=p[0]+p[1];
00423       for(int i=0;i<3;i++) {
00424          fPtRec[i]=p[i].Pt();
00425          fEtaRec[i]=p[i].Eta();
00426          fPhiRec[i]=p[i].Phi();
00427          fMRec[i]=p[i].M();
00428       }
00429    }
00430    if(fIsSignal) {
00431       do {
00432          Double_t tauDiscr=0.08-0.04/(1.+fPtRec[2]/10.0);
00433          Double_t sigmaDiscr=0.01;
00434          fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
00435       } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
00436    } else {
00437       do {
00438          Double_t tauDiscr=0.15-0.05/(1.+fPtRec[2]/5.0)+0.1*fEtaRec[2];
00439          Double_t sigmaDiscr=0.02+0.01*fEtaRec[2];
00440          fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
00441       } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
00442    }
00443    fIsTriggered=false;
00444    for(int i=0;i<2;i++) {
00445       if(rnd->Uniform()<0.92/(TMath::Exp(-(fPtRec[i]-15.5)/2.5)+1.)) fIsTriggered=true;
00446    }
00447 }