ROOT logo
// @(#)root/roostats:$Id: ProposalHelper.cxx 34109 2010-06-24 15:00:16Z moneta $
// Authors: Kevin Belasco        7/22/2009
// Authors: Kyle Cranmer         7/22/2009
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//_________________________________________________
/*
BEGIN_HTML
END_HTML
*/
//_________________________________________________

#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif
#ifndef RooStats_ProposalHelper
#include "RooStats/ProposalHelper.h"
#endif
#ifndef ROOSTATS_PdfProposal
#include "RooStats/PdfProposal.h"
#endif
#ifndef RooStats_RooStatsUtils
#include "RooStats/RooStatsUtils.h"
#endif
#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif
#ifndef ROO_DATA_SET
#include "RooDataSet.h"
#endif
#ifndef ROO_ABS_PDF
#include "RooAbsPdf.h"
#endif
#ifndef ROO_ADD_PDF
#include "RooAddPdf.h"
#endif
#ifndef ROO_KEYS_PDF
#include "RooNDKeysPdf.h"
#endif
#ifndef ROO_UNIFORM
#include "RooUniform.h"
#endif
#ifndef ROO_MSG_SERVICE
#include "RooMsgService.h"
#endif
#ifndef ROO_REAL_VAR
#include "RooRealVar.h"
#endif
#ifndef ROOT_TIterator
#include "TIterator.h"
#endif
#ifndef ROO_MULTI_VAR_GAUSSIAN
#include "RooMultiVarGaussian.h"
#endif
#ifndef ROO_CONST_VAR
#include "RooConstVar.h"
#endif
#ifndef ROOT_TString
#include "TString.h"
#endif

#include <map>

ClassImp(RooStats::ProposalHelper);

using namespace RooFit;
using namespace RooStats;

static const Double_t DEFAULT_UNI_FRAC = 0.10;
static const Double_t DEFAULT_CLUES_FRAC = 0.20;
//static const Double_t SIGMA_RANGE_DIVISOR = 6;
static const Double_t SIGMA_RANGE_DIVISOR = 5;
static const Int_t DEFAULT_CACHE_SIZE = 100;
//static const Option_t* CLUES_OPTIONS = "a";

ProposalHelper::ProposalHelper()
{
   fPdfProp = new PdfProposal();
   fVars = NULL;
   fOwnsPdfProp = kTRUE;
   fOwnsPdf = kFALSE;
   fOwnsCluesPdf = kFALSE;
   fOwnsVars = kFALSE;
   fUseUpdates = kFALSE;
   fPdf = NULL;
   fSigmaRangeDivisor = SIGMA_RANGE_DIVISOR;
   fCluesPdf = NULL;
   fUniformPdf = NULL;
   fClues = NULL;
   fCovMatrix = NULL;
   fCluesFrac = -1;
   fUniFrac = -1;
   fCacheSize = -1;
   fCluesOptions = NULL;
}

ProposalFunction* ProposalHelper::GetProposalFunction()
{
   if (fPdf == NULL)
      CreatePdf();
   // kbelasco: check here for memory leaks: does RooAddPdf make copies or
   // take ownership of components, coeffs
   RooArgList* components = new RooArgList();
   RooArgList* coeffs = new RooArgList();
   if (fCluesPdf == NULL)
      CreateCluesPdf();
   if (fCluesPdf != NULL) {
      if (fCluesFrac < 0)
         fCluesFrac = DEFAULT_CLUES_FRAC;
      printf("added clues from dataset %s with fraction %g\n",
            fClues->GetName(), fCluesFrac);
      components->add(*fCluesPdf);
      coeffs->add(RooConst(fCluesFrac));
   }
   if (fUniFrac > 0.) {
      CreateUniformPdf();
      components->add(*fUniformPdf);
      coeffs->add(RooConst(fUniFrac));
   }
   components->add(*fPdf);
   RooAddPdf* addPdf = new RooAddPdf("proposalFunction", "Proposal Density",
         *components, *coeffs);
   fPdfProp->SetPdf(*addPdf);
   fPdfProp->SetOwnsPdf(kTRUE);
   if (fCacheSize > 0)
      fPdfProp->SetCacheSize(fCacheSize);
   fOwnsPdfProp = kFALSE;
   return fPdfProp;
}

void ProposalHelper::CreatePdf()
{
   // kbelasco: check here for memory leaks:
   // does RooMultiVarGaussian make copies of xVec and muVec?
   // or should we delete them?
   if (fVars == NULL) {
      coutE(InputArguments) << "ProposalHelper::CreatePdf(): " <<
         "Variables to create proposal function for are not set." << endl;
      return;
   }
   RooArgList* xVec = new RooArgList();
   RooArgList* muVec = new RooArgList();
   TIterator* it = fVars->createIterator();
   RooRealVar* r;
   RooRealVar* clone;
   while ((r = (RooRealVar*)it->Next()) != NULL) {
      xVec->add(*r);
      TString cloneName = TString::Format("%s%s", "mu__", r->GetName());
      clone = (RooRealVar*)r->clone(cloneName.Data());
      muVec->add(*clone);
      if (fUseUpdates)
         fPdfProp->AddMapping(*clone, *r);
   }
   if (fCovMatrix == NULL)
      CreateCovMatrix(*xVec);
   fPdf = new RooMultiVarGaussian("mvg", "MVG Proposal", *xVec, *muVec,
                                  *fCovMatrix);
   delete xVec;
   delete muVec;
   delete it;
}

void ProposalHelper::CreateCovMatrix(RooArgList& xVec)
{
   Int_t size = xVec.getSize();
   fCovMatrix = new TMatrixDSym(size);
   RooRealVar* r;
   for (Int_t i = 0; i < size; i++) {
      r = (RooRealVar*)xVec.at(i);
      Double_t range = r->getMax() - r->getMin();
      (*fCovMatrix)(i,i) = range / fSigmaRangeDivisor;
   }
}

void ProposalHelper::CreateCluesPdf()
{
   if (fClues != NULL) {
      if (fCluesOptions == NULL)
         fCluesPdf = new RooNDKeysPdf("cluesPdf", "Clues PDF", *fVars, *fClues);
      else
         fCluesPdf = new RooNDKeysPdf("cluesPdf", "Clues PDF", *fVars, *fClues,
               fCluesOptions);
   }
}

void ProposalHelper::CreateUniformPdf()
{
   fUniformPdf = new RooUniform("uniform", "Uniform Proposal PDF",
         RooArgSet(*fVars));
}
 ProposalHelper.cxx:1
 ProposalHelper.cxx:2
 ProposalHelper.cxx:3
 ProposalHelper.cxx:4
 ProposalHelper.cxx:5
 ProposalHelper.cxx:6
 ProposalHelper.cxx:7
 ProposalHelper.cxx:8
 ProposalHelper.cxx:9
 ProposalHelper.cxx:10
 ProposalHelper.cxx:11
 ProposalHelper.cxx:12
 ProposalHelper.cxx:13
 ProposalHelper.cxx:14
 ProposalHelper.cxx:15
 ProposalHelper.cxx:16
 ProposalHelper.cxx:17
 ProposalHelper.cxx:18
 ProposalHelper.cxx:19
 ProposalHelper.cxx:20
 ProposalHelper.cxx:21
 ProposalHelper.cxx:22
 ProposalHelper.cxx:23
 ProposalHelper.cxx:24
 ProposalHelper.cxx:25
 ProposalHelper.cxx:26
 ProposalHelper.cxx:27
 ProposalHelper.cxx:28
 ProposalHelper.cxx:29
 ProposalHelper.cxx:30
 ProposalHelper.cxx:31
 ProposalHelper.cxx:32
 ProposalHelper.cxx:33
 ProposalHelper.cxx:34
 ProposalHelper.cxx:35
 ProposalHelper.cxx:36
 ProposalHelper.cxx:37
 ProposalHelper.cxx:38
 ProposalHelper.cxx:39
 ProposalHelper.cxx:40
 ProposalHelper.cxx:41
 ProposalHelper.cxx:42
 ProposalHelper.cxx:43
 ProposalHelper.cxx:44
 ProposalHelper.cxx:45
 ProposalHelper.cxx:46
 ProposalHelper.cxx:47
 ProposalHelper.cxx:48
 ProposalHelper.cxx:49
 ProposalHelper.cxx:50
 ProposalHelper.cxx:51
 ProposalHelper.cxx:52
 ProposalHelper.cxx:53
 ProposalHelper.cxx:54
 ProposalHelper.cxx:55
 ProposalHelper.cxx:56
 ProposalHelper.cxx:57
 ProposalHelper.cxx:58
 ProposalHelper.cxx:59
 ProposalHelper.cxx:60
 ProposalHelper.cxx:61
 ProposalHelper.cxx:62
 ProposalHelper.cxx:63
 ProposalHelper.cxx:64
 ProposalHelper.cxx:65
 ProposalHelper.cxx:66
 ProposalHelper.cxx:67
 ProposalHelper.cxx:68
 ProposalHelper.cxx:69
 ProposalHelper.cxx:70
 ProposalHelper.cxx:71
 ProposalHelper.cxx:72
 ProposalHelper.cxx:73
 ProposalHelper.cxx:74
 ProposalHelper.cxx:75
 ProposalHelper.cxx:76
 ProposalHelper.cxx:77
 ProposalHelper.cxx:78
 ProposalHelper.cxx:79
 ProposalHelper.cxx:80
 ProposalHelper.cxx:81
 ProposalHelper.cxx:82
 ProposalHelper.cxx:83
 ProposalHelper.cxx:84
 ProposalHelper.cxx:85
 ProposalHelper.cxx:86
 ProposalHelper.cxx:87
 ProposalHelper.cxx:88
 ProposalHelper.cxx:89
 ProposalHelper.cxx:90
 ProposalHelper.cxx:91
 ProposalHelper.cxx:92
 ProposalHelper.cxx:93
 ProposalHelper.cxx:94
 ProposalHelper.cxx:95
 ProposalHelper.cxx:96
 ProposalHelper.cxx:97
 ProposalHelper.cxx:98
 ProposalHelper.cxx:99
 ProposalHelper.cxx:100
 ProposalHelper.cxx:101
 ProposalHelper.cxx:102
 ProposalHelper.cxx:103
 ProposalHelper.cxx:104
 ProposalHelper.cxx:105
 ProposalHelper.cxx:106
 ProposalHelper.cxx:107
 ProposalHelper.cxx:108
 ProposalHelper.cxx:109
 ProposalHelper.cxx:110
 ProposalHelper.cxx:111
 ProposalHelper.cxx:112
 ProposalHelper.cxx:113
 ProposalHelper.cxx:114
 ProposalHelper.cxx:115
 ProposalHelper.cxx:116
 ProposalHelper.cxx:117
 ProposalHelper.cxx:118
 ProposalHelper.cxx:119
 ProposalHelper.cxx:120
 ProposalHelper.cxx:121
 ProposalHelper.cxx:122
 ProposalHelper.cxx:123
 ProposalHelper.cxx:124
 ProposalHelper.cxx:125
 ProposalHelper.cxx:126
 ProposalHelper.cxx:127
 ProposalHelper.cxx:128
 ProposalHelper.cxx:129
 ProposalHelper.cxx:130
 ProposalHelper.cxx:131
 ProposalHelper.cxx:132
 ProposalHelper.cxx:133
 ProposalHelper.cxx:134
 ProposalHelper.cxx:135
 ProposalHelper.cxx:136
 ProposalHelper.cxx:137
 ProposalHelper.cxx:138
 ProposalHelper.cxx:139
 ProposalHelper.cxx:140
 ProposalHelper.cxx:141
 ProposalHelper.cxx:142
 ProposalHelper.cxx:143
 ProposalHelper.cxx:144
 ProposalHelper.cxx:145
 ProposalHelper.cxx:146
 ProposalHelper.cxx:147
 ProposalHelper.cxx:148
 ProposalHelper.cxx:149
 ProposalHelper.cxx:150
 ProposalHelper.cxx:151
 ProposalHelper.cxx:152
 ProposalHelper.cxx:153
 ProposalHelper.cxx:154
 ProposalHelper.cxx:155
 ProposalHelper.cxx:156
 ProposalHelper.cxx:157
 ProposalHelper.cxx:158
 ProposalHelper.cxx:159
 ProposalHelper.cxx:160
 ProposalHelper.cxx:161
 ProposalHelper.cxx:162
 ProposalHelper.cxx:163
 ProposalHelper.cxx:164
 ProposalHelper.cxx:165
 ProposalHelper.cxx:166
 ProposalHelper.cxx:167
 ProposalHelper.cxx:168
 ProposalHelper.cxx:169
 ProposalHelper.cxx:170
 ProposalHelper.cxx:171
 ProposalHelper.cxx:172
 ProposalHelper.cxx:173
 ProposalHelper.cxx:174
 ProposalHelper.cxx:175
 ProposalHelper.cxx:176
 ProposalHelper.cxx:177
 ProposalHelper.cxx:178
 ProposalHelper.cxx:179
 ProposalHelper.cxx:180
 ProposalHelper.cxx:181
 ProposalHelper.cxx:182
 ProposalHelper.cxx:183
 ProposalHelper.cxx:184
 ProposalHelper.cxx:185
 ProposalHelper.cxx:186
 ProposalHelper.cxx:187
 ProposalHelper.cxx:188
 ProposalHelper.cxx:189
 ProposalHelper.cxx:190
 ProposalHelper.cxx:191
 ProposalHelper.cxx:192
 ProposalHelper.cxx:193
 ProposalHelper.cxx:194
 ProposalHelper.cxx:195
 ProposalHelper.cxx:196