#ifndef ROO_MSG_SERVICE
#include "RooMsgService.h"
#endif
#include "RooStats/SamplingDistribution.h"
#include "RooNumber.h"
#include "TMath.h"
#include <algorithm>
#include <iostream>
#include <cmath>
#include <limits>
using namespace std ;
ClassImp(RooStats::SamplingDistribution)
using namespace RooStats;
SamplingDistribution::SamplingDistribution( const char *name, const char *title,
std::vector<Double_t>& samplingDist, const char * varName) :
TNamed(name,title)
{
fSamplingDist = samplingDist;
fSampleWeights.resize(fSamplingDist.size(),1.0) ;
fVarName = varName;
}
SamplingDistribution::SamplingDistribution( const char *name, const char *title,
std::vector<Double_t>& samplingDist, std::vector<Double_t>& sampleWeights, const char * varName) :
TNamed(name,title)
{
fSamplingDist = samplingDist;
fSampleWeights = sampleWeights;
fVarName = varName;
}
SamplingDistribution::SamplingDistribution( const char *name, const char *title, const char * varName) :
TNamed(name,title)
{
fVarName = varName;
}
SamplingDistribution::SamplingDistribution(
const char *name,
const char *title,
RooDataSet& dataSet,
const char * _columnName,
const char * varName
) : TNamed(name, title) {
if( dataSet.numEntries() == 0 || !dataSet.get()->first() ) {
if( varName ) fVarName = varName;
return;
}
TString columnName( _columnName );
if( !columnName.Length() ) {
columnName.Form( "%s_TS0", name );
if( !dataSet.get()->find(columnName) ) {
columnName = dataSet.get()->first()->GetName();
}
}
if( !varName ) {
fVarName = (*dataSet.get())[columnName].GetTitle();
}else{
fVarName = varName;
}
for(Int_t i=0; i < dataSet.numEntries(); i++) {
fSamplingDist.push_back(dataSet.get(i)->getRealValue(columnName));
fSampleWeights.push_back(dataSet.weight());
}
}
SamplingDistribution::SamplingDistribution( ) :
TNamed("SamplingDistribution_DefaultName","SamplingDistribution")
{
}
SamplingDistribution::~SamplingDistribution()
{
fSamplingDist.clear();
fSampleWeights.clear();
}
void SamplingDistribution::Add(const SamplingDistribution* other)
{
if(!other) return;
std::vector<double> newSamplingDist = other->fSamplingDist;
std::vector<double> newSampleWeights = other->fSampleWeights;
fSamplingDist.reserve(fSamplingDist.size()+newSamplingDist.size());
fSampleWeights.reserve(fSampleWeights.size()+newSampleWeights.size());
for(unsigned int i=0; i<newSamplingDist.size(); ++i){
fSamplingDist.push_back(newSamplingDist[i]);
fSampleWeights.push_back(newSampleWeights[i]);
}
if(GetVarName().Length() == 0 && other->GetVarName().Length() > 0)
fVarName = other->GetVarName();
if(strlen(GetName()) == 0 && strlen(other->GetName()) > 0)
SetName(other->GetName());
if(strlen(GetTitle()) == 0 && strlen(other->GetTitle()) > 0)
SetTitle(other->GetTitle());
}
Double_t SamplingDistribution::Integral(Double_t low, Double_t high, Bool_t normalize, Bool_t lowClosed, Bool_t
highClosed) const
{
double error = 0;
return IntegralAndError(error, low,high, normalize, lowClosed, highClosed);
}
void SamplingDistribution::SortValues() const {
unsigned int n = fSamplingDist.size();
std::vector<unsigned int> index(n);
TMath::SortItr(fSamplingDist.begin(), fSamplingDist.end(), index.begin(), false );
fSumW = std::vector<double>( n );
fSumW2 = std::vector<double>( n );
std::vector<double> sortedDist( n);
std::vector<double> sortedWeights( n);
for(unsigned int i=0; i <n; i++) {
unsigned int j = index[i];
if (i > 0) {
fSumW[i] += fSumW[i-1];
fSumW2[i] += fSumW2[i-1];
}
fSumW[i] += fSampleWeights[j];
fSumW2[i] += fSampleWeights[j]*fSampleWeights[j];
sortedDist[i] = fSamplingDist[ j] ;
sortedWeights[i] = fSampleWeights[ j] ;
}
fSamplingDist = sortedDist;
fSampleWeights = sortedWeights;
}
Double_t SamplingDistribution::IntegralAndError(Double_t & error, Double_t low, Double_t high, Bool_t normalize, Bool_t lowClosed, Bool_t
highClosed) const
{
int n = fSamplingDist.size();
if( n == 0 ) {
error = numeric_limits<Double_t>::infinity();
return 0.0;
}
if (int(fSumW.size()) != n)
SortValues();
int indexLow = -1;
int indexHigh = -1;
if (lowClosed) {
indexLow = std::lower_bound( fSamplingDist.begin(), fSamplingDist.end() , low) - fSamplingDist.begin() -1;
}
else {
indexLow = std::upper_bound( fSamplingDist.begin(), fSamplingDist.end() , low) - fSamplingDist.begin() - 1;
}
if (highClosed) {
indexHigh = std::upper_bound( fSamplingDist.begin(), fSamplingDist.end() , high) - fSamplingDist.begin() -1;
}
else {
indexHigh = std::lower_bound( fSamplingDist.begin(), fSamplingDist.end() , high) - fSamplingDist.begin() -1;
}
assert(indexLow < n && indexHigh < n);
double sum = 0;
double sum2 = 0;
if (indexHigh >= 0) {
sum = fSumW[indexHigh];
sum2 = fSumW2[indexHigh];
if (indexLow >= 0) {
sum -= fSumW[indexLow];
sum2 -= fSumW2[indexLow];
}
}
if(normalize) {
double norm = fSumW.back();
double norm2 = fSumW2.back();
sum /= norm;
error = std::sqrt( sum2 * (1. - 2. * sum) + norm2 * sum * sum ) / norm;
}
else {
error = std::sqrt(sum2);
}
return sum;
}
Double_t SamplingDistribution::CDF(Double_t x) const {
return Integral(-RooNumber::infinity(), x, kTRUE, kTRUE, kTRUE);
}
Double_t SamplingDistribution::InverseCDF(Double_t pvalue)
{
Double_t dummy=0;
return InverseCDF(pvalue,0,dummy);
}
Double_t SamplingDistribution::InverseCDF(Double_t pvalue,
Double_t sigmaVariation,
Double_t& inverseWithVariation)
{
if (fSumW.size() != fSamplingDist.size())
SortValues();
if (!TMath::AreEqualRel(fSumW.back(), fSumW2.back(), 1.E-6) )
Warning("InverseCDF","Estimation of Quantiles (InverseCDF) for weighted events is not yet supported");
int nominal = (unsigned int) (pvalue*fSamplingDist.size());
if(nominal <= 0) {
inverseWithVariation = -1.*RooNumber::infinity();
return -1.*RooNumber::infinity();
}
else if(nominal >= (Int_t)fSamplingDist.size()-1 ) {
inverseWithVariation = RooNumber::infinity();
return RooNumber::infinity();
}
else if(pvalue < 0.5){
int delta = (int)(sigmaVariation*sqrt(1.0*nominal));
int variation = nominal+delta;
if(variation>=(Int_t)fSamplingDist.size()-1)
inverseWithVariation = RooNumber::infinity();
else if(variation<=0)
inverseWithVariation = -1.*RooNumber::infinity();
else
inverseWithVariation = fSamplingDist[ variation ];
return fSamplingDist[nominal];
}
else if(pvalue >= 0.5){
int delta = (int)(sigmaVariation*sqrt(1.0*fSamplingDist.size()- nominal));
int variation = nominal+delta;
if(variation>=(Int_t)fSamplingDist.size()-1)
inverseWithVariation = RooNumber::infinity();
else if(variation<=0)
inverseWithVariation = -1.*RooNumber::infinity();
else
inverseWithVariation = fSamplingDist[ variation+1 ];
return fSamplingDist[nominal+1];
}
else{
std::cout << "problem in SamplingDistribution::InverseCDF" << std::endl;
}
inverseWithVariation = RooNumber::infinity();
return RooNumber::infinity();
}
Double_t SamplingDistribution::InverseCDFInterpolate(Double_t pvalue)
{
if (fSumW.size() != fSamplingDist.size())
SortValues();
if (!TMath::AreEqualRel(fSumW.back(), fSumW2.back(), 1.E-6) )
Warning("InverseCDFInterpolate","Estimation of Quantiles (InverseCDF) for weighted events is not yet supported.");
int nominal = (unsigned int) (pvalue*fSamplingDist.size());
if(nominal <= 0) {
return -1.*RooNumber::infinity();
}
if(nominal >= (Int_t)fSamplingDist.size()-1 ) {
return RooNumber::infinity();
}
Double_t upperX = fSamplingDist[nominal+1];
Double_t upperY = ((Double_t) (nominal+1))/fSamplingDist.size();
Double_t lowerX = fSamplingDist[nominal];
Double_t lowerY = ((Double_t) nominal)/fSamplingDist.size();
return (upperX-lowerX)/(upperY-lowerY)*(pvalue-lowerY)+lowerX;
}
SamplingDistribution.cxx:1 SamplingDistribution.cxx:2 SamplingDistribution.cxx:3 SamplingDistribution.cxx:4 SamplingDistribution.cxx:5 SamplingDistribution.cxx:6 SamplingDistribution.cxx:7 SamplingDistribution.cxx:8 SamplingDistribution.cxx:9 SamplingDistribution.cxx:10 SamplingDistribution.cxx:11 SamplingDistribution.cxx:12 SamplingDistribution.cxx:13 SamplingDistribution.cxx:14 SamplingDistribution.cxx:15 SamplingDistribution.cxx:16 SamplingDistribution.cxx:17 SamplingDistribution.cxx:18 SamplingDistribution.cxx:19 SamplingDistribution.cxx:20 SamplingDistribution.cxx:21 SamplingDistribution.cxx:22 SamplingDistribution.cxx:23 SamplingDistribution.cxx:24 SamplingDistribution.cxx:25 SamplingDistribution.cxx:26 SamplingDistribution.cxx:27 SamplingDistribution.cxx:28 SamplingDistribution.cxx:29 SamplingDistribution.cxx:30 SamplingDistribution.cxx:31 SamplingDistribution.cxx:32 SamplingDistribution.cxx:33 SamplingDistribution.cxx:34 SamplingDistribution.cxx:35 SamplingDistribution.cxx:36 SamplingDistribution.cxx:37 SamplingDistribution.cxx:38 SamplingDistribution.cxx:39 SamplingDistribution.cxx:40 SamplingDistribution.cxx:41 SamplingDistribution.cxx:42 SamplingDistribution.cxx:43 SamplingDistribution.cxx:44 SamplingDistribution.cxx:45 SamplingDistribution.cxx:46 SamplingDistribution.cxx:47 SamplingDistribution.cxx:48 SamplingDistribution.cxx:49 SamplingDistribution.cxx:50 SamplingDistribution.cxx:51 SamplingDistribution.cxx:52 SamplingDistribution.cxx:53 SamplingDistribution.cxx:54 SamplingDistribution.cxx:55 SamplingDistribution.cxx:56 SamplingDistribution.cxx:57 SamplingDistribution.cxx:58 SamplingDistribution.cxx:59 SamplingDistribution.cxx:60 SamplingDistribution.cxx:61 SamplingDistribution.cxx:62 SamplingDistribution.cxx:63 SamplingDistribution.cxx:64 SamplingDistribution.cxx:65 SamplingDistribution.cxx:66 SamplingDistribution.cxx:67 SamplingDistribution.cxx:68 SamplingDistribution.cxx:69 SamplingDistribution.cxx:70 SamplingDistribution.cxx:71 SamplingDistribution.cxx:72 SamplingDistribution.cxx:73 SamplingDistribution.cxx:74 SamplingDistribution.cxx:75 SamplingDistribution.cxx:76 SamplingDistribution.cxx:77 SamplingDistribution.cxx:78 SamplingDistribution.cxx:79 SamplingDistribution.cxx:80 SamplingDistribution.cxx:81 SamplingDistribution.cxx:82 SamplingDistribution.cxx:83 SamplingDistribution.cxx:84 SamplingDistribution.cxx:85 SamplingDistribution.cxx:86 SamplingDistribution.cxx:87 SamplingDistribution.cxx:88 SamplingDistribution.cxx:89 SamplingDistribution.cxx:90 SamplingDistribution.cxx:91 SamplingDistribution.cxx:92 SamplingDistribution.cxx:93 SamplingDistribution.cxx:94 SamplingDistribution.cxx:95 SamplingDistribution.cxx:96 SamplingDistribution.cxx:97 SamplingDistribution.cxx:98 SamplingDistribution.cxx:99 SamplingDistribution.cxx:100 SamplingDistribution.cxx:101 SamplingDistribution.cxx:102 SamplingDistribution.cxx:103 SamplingDistribution.cxx:104 SamplingDistribution.cxx:105 SamplingDistribution.cxx:106 SamplingDistribution.cxx:107 SamplingDistribution.cxx:108 SamplingDistribution.cxx:109 SamplingDistribution.cxx:110 SamplingDistribution.cxx:111 SamplingDistribution.cxx:112 SamplingDistribution.cxx:113 SamplingDistribution.cxx:114 SamplingDistribution.cxx:115 SamplingDistribution.cxx:116 SamplingDistribution.cxx:117 SamplingDistribution.cxx:118 SamplingDistribution.cxx:119 SamplingDistribution.cxx:120 SamplingDistribution.cxx:121 SamplingDistribution.cxx:122 SamplingDistribution.cxx:123 SamplingDistribution.cxx:124 SamplingDistribution.cxx:125 SamplingDistribution.cxx:126 SamplingDistribution.cxx:127 SamplingDistribution.cxx:128 SamplingDistribution.cxx:129 SamplingDistribution.cxx:130 SamplingDistribution.cxx:131 SamplingDistribution.cxx:132 SamplingDistribution.cxx:133 SamplingDistribution.cxx:134 SamplingDistribution.cxx:135 SamplingDistribution.cxx:136 SamplingDistribution.cxx:137 SamplingDistribution.cxx:138 SamplingDistribution.cxx:139 SamplingDistribution.cxx:140 SamplingDistribution.cxx:141 SamplingDistribution.cxx:142 SamplingDistribution.cxx:143 SamplingDistribution.cxx:144 SamplingDistribution.cxx:145 SamplingDistribution.cxx:146 SamplingDistribution.cxx:147 SamplingDistribution.cxx:148 SamplingDistribution.cxx:149 SamplingDistribution.cxx:150 SamplingDistribution.cxx:151 SamplingDistribution.cxx:152 SamplingDistribution.cxx:153 SamplingDistribution.cxx:154 SamplingDistribution.cxx:155 SamplingDistribution.cxx:156 SamplingDistribution.cxx:157 SamplingDistribution.cxx:158 SamplingDistribution.cxx:159 SamplingDistribution.cxx:160 SamplingDistribution.cxx:161 SamplingDistribution.cxx:162 SamplingDistribution.cxx:163 SamplingDistribution.cxx:164 SamplingDistribution.cxx:165 SamplingDistribution.cxx:166 SamplingDistribution.cxx:167 SamplingDistribution.cxx:168 SamplingDistribution.cxx:169 SamplingDistribution.cxx:170 SamplingDistribution.cxx:171 SamplingDistribution.cxx:172 SamplingDistribution.cxx:173 SamplingDistribution.cxx:174 SamplingDistribution.cxx:175 SamplingDistribution.cxx:176 SamplingDistribution.cxx:177 SamplingDistribution.cxx:178 SamplingDistribution.cxx:179 SamplingDistribution.cxx:180 SamplingDistribution.cxx:181 SamplingDistribution.cxx:182 SamplingDistribution.cxx:183 SamplingDistribution.cxx:184 SamplingDistribution.cxx:185 SamplingDistribution.cxx:186 SamplingDistribution.cxx:187 SamplingDistribution.cxx:188 SamplingDistribution.cxx:189 SamplingDistribution.cxx:190 SamplingDistribution.cxx:191 SamplingDistribution.cxx:192 SamplingDistribution.cxx:193 SamplingDistribution.cxx:194 SamplingDistribution.cxx:195 SamplingDistribution.cxx:196 SamplingDistribution.cxx:197 SamplingDistribution.cxx:198 SamplingDistribution.cxx:199 SamplingDistribution.cxx:200 SamplingDistribution.cxx:201 SamplingDistribution.cxx:202 SamplingDistribution.cxx:203 SamplingDistribution.cxx:204 SamplingDistribution.cxx:205 SamplingDistribution.cxx:206 SamplingDistribution.cxx:207 SamplingDistribution.cxx:208 SamplingDistribution.cxx:209 SamplingDistribution.cxx:210 SamplingDistribution.cxx:211 SamplingDistribution.cxx:212 SamplingDistribution.cxx:213 SamplingDistribution.cxx:214 SamplingDistribution.cxx:215 SamplingDistribution.cxx:216 SamplingDistribution.cxx:217 SamplingDistribution.cxx:218 SamplingDistribution.cxx:219 SamplingDistribution.cxx:220 SamplingDistribution.cxx:221 SamplingDistribution.cxx:222 SamplingDistribution.cxx:223 SamplingDistribution.cxx:224 SamplingDistribution.cxx:225 SamplingDistribution.cxx:226 SamplingDistribution.cxx:227 SamplingDistribution.cxx:228 SamplingDistribution.cxx:229 SamplingDistribution.cxx:230 SamplingDistribution.cxx:231 SamplingDistribution.cxx:232 SamplingDistribution.cxx:233 SamplingDistribution.cxx:234 SamplingDistribution.cxx:235 SamplingDistribution.cxx:236 SamplingDistribution.cxx:237 SamplingDistribution.cxx:238 SamplingDistribution.cxx:239 SamplingDistribution.cxx:240 SamplingDistribution.cxx:241 SamplingDistribution.cxx:242 SamplingDistribution.cxx:243 SamplingDistribution.cxx:244 SamplingDistribution.cxx:245 SamplingDistribution.cxx:246 SamplingDistribution.cxx:247 SamplingDistribution.cxx:248 SamplingDistribution.cxx:249 SamplingDistribution.cxx:250 SamplingDistribution.cxx:251 SamplingDistribution.cxx:252 SamplingDistribution.cxx:253 SamplingDistribution.cxx:254 SamplingDistribution.cxx:255 SamplingDistribution.cxx:256 SamplingDistribution.cxx:257 SamplingDistribution.cxx:258 SamplingDistribution.cxx:259 SamplingDistribution.cxx:260 SamplingDistribution.cxx:261 SamplingDistribution.cxx:262 SamplingDistribution.cxx:263 SamplingDistribution.cxx:264 SamplingDistribution.cxx:265 SamplingDistribution.cxx:266 SamplingDistribution.cxx:267 SamplingDistribution.cxx:268 SamplingDistribution.cxx:269 SamplingDistribution.cxx:270 SamplingDistribution.cxx:271 SamplingDistribution.cxx:272 SamplingDistribution.cxx:273 SamplingDistribution.cxx:274 SamplingDistribution.cxx:275 SamplingDistribution.cxx:276 SamplingDistribution.cxx:277 SamplingDistribution.cxx:278 SamplingDistribution.cxx:279 SamplingDistribution.cxx:280 SamplingDistribution.cxx:281 SamplingDistribution.cxx:282 SamplingDistribution.cxx:283 SamplingDistribution.cxx:284 SamplingDistribution.cxx:285 SamplingDistribution.cxx:286 SamplingDistribution.cxx:287 SamplingDistribution.cxx:288 SamplingDistribution.cxx:289 SamplingDistribution.cxx:290 SamplingDistribution.cxx:291 SamplingDistribution.cxx:292 SamplingDistribution.cxx:293 SamplingDistribution.cxx:294 SamplingDistribution.cxx:295 SamplingDistribution.cxx:296 SamplingDistribution.cxx:297 SamplingDistribution.cxx:298 SamplingDistribution.cxx:299 SamplingDistribution.cxx:300 SamplingDistribution.cxx:301 SamplingDistribution.cxx:302 SamplingDistribution.cxx:303 SamplingDistribution.cxx:304 SamplingDistribution.cxx:305 SamplingDistribution.cxx:306 SamplingDistribution.cxx:307 SamplingDistribution.cxx:308 SamplingDistribution.cxx:309 SamplingDistribution.cxx:310 SamplingDistribution.cxx:311 SamplingDistribution.cxx:312 SamplingDistribution.cxx:313 SamplingDistribution.cxx:314 SamplingDistribution.cxx:315 SamplingDistribution.cxx:316 SamplingDistribution.cxx:317 SamplingDistribution.cxx:318 SamplingDistribution.cxx:319 SamplingDistribution.cxx:320 SamplingDistribution.cxx:321 SamplingDistribution.cxx:322 SamplingDistribution.cxx:323 SamplingDistribution.cxx:324 SamplingDistribution.cxx:325 SamplingDistribution.cxx:326 SamplingDistribution.cxx:327 SamplingDistribution.cxx:328 SamplingDistribution.cxx:329 SamplingDistribution.cxx:330 SamplingDistribution.cxx:331 SamplingDistribution.cxx:332 SamplingDistribution.cxx:333 SamplingDistribution.cxx:334 SamplingDistribution.cxx:335 SamplingDistribution.cxx:336 SamplingDistribution.cxx:337 SamplingDistribution.cxx:338 SamplingDistribution.cxx:339 SamplingDistribution.cxx:340 SamplingDistribution.cxx:341 SamplingDistribution.cxx:342 SamplingDistribution.cxx:343 SamplingDistribution.cxx:344 SamplingDistribution.cxx:345 SamplingDistribution.cxx:346 SamplingDistribution.cxx:347 SamplingDistribution.cxx:348 SamplingDistribution.cxx:349 SamplingDistribution.cxx:350 SamplingDistribution.cxx:351 SamplingDistribution.cxx:352 SamplingDistribution.cxx:353 SamplingDistribution.cxx:354 SamplingDistribution.cxx:355 SamplingDistribution.cxx:356 SamplingDistribution.cxx:357 SamplingDistribution.cxx:358 SamplingDistribution.cxx:359 SamplingDistribution.cxx:360 SamplingDistribution.cxx:361 SamplingDistribution.cxx:362 SamplingDistribution.cxx:363 SamplingDistribution.cxx:364 SamplingDistribution.cxx:365 SamplingDistribution.cxx:366 SamplingDistribution.cxx:367 SamplingDistribution.cxx:368 SamplingDistribution.cxx:369 SamplingDistribution.cxx:370 SamplingDistribution.cxx:371 SamplingDistribution.cxx:372 SamplingDistribution.cxx:373 SamplingDistribution.cxx:374 SamplingDistribution.cxx:375 SamplingDistribution.cxx:376 SamplingDistribution.cxx:377 SamplingDistribution.cxx:378 SamplingDistribution.cxx:379 SamplingDistribution.cxx:380 SamplingDistribution.cxx:381 SamplingDistribution.cxx:382 SamplingDistribution.cxx:383 SamplingDistribution.cxx:384 SamplingDistribution.cxx:385 SamplingDistribution.cxx:386 SamplingDistribution.cxx:387 SamplingDistribution.cxx:388 SamplingDistribution.cxx:389 SamplingDistribution.cxx:390 SamplingDistribution.cxx:391 SamplingDistribution.cxx:392 SamplingDistribution.cxx:393 SamplingDistribution.cxx:394 SamplingDistribution.cxx:395 SamplingDistribution.cxx:396 SamplingDistribution.cxx:397 SamplingDistribution.cxx:398 SamplingDistribution.cxx:399 SamplingDistribution.cxx:400 SamplingDistribution.cxx:401 SamplingDistribution.cxx:402 SamplingDistribution.cxx:403 SamplingDistribution.cxx:404 SamplingDistribution.cxx:405 SamplingDistribution.cxx:406 SamplingDistribution.cxx:407 SamplingDistribution.cxx:408 SamplingDistribution.cxx:409 SamplingDistribution.cxx:410 SamplingDistribution.cxx:411 SamplingDistribution.cxx:412 SamplingDistribution.cxx:413 SamplingDistribution.cxx:414 SamplingDistribution.cxx:415 SamplingDistribution.cxx:416 SamplingDistribution.cxx:417 SamplingDistribution.cxx:418 SamplingDistribution.cxx:419 SamplingDistribution.cxx:420 SamplingDistribution.cxx:421 SamplingDistribution.cxx:422 SamplingDistribution.cxx:423 SamplingDistribution.cxx:424 SamplingDistribution.cxx:425 SamplingDistribution.cxx:426 SamplingDistribution.cxx:427 SamplingDistribution.cxx:428 SamplingDistribution.cxx:429 SamplingDistribution.cxx:430 SamplingDistribution.cxx:431 SamplingDistribution.cxx:432