Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Channel.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, George Lewis
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11////////////////////////////////////////////////////////////////////////////////
12/** \class RooStats::HistFactory::Channel
13 * \ingroup HistFactory
14 This class encapsulates all information for the statistical interpretation of one experiment.
15 It can be combined with other channels (e.g. for the combination of multiple experiments, or
16 to constrain nuisance parameters with information obtained in a control region).
17 A channel contains one or more samples which describe the contribution from different processes
18 to this measurement.
19*/
20
21
22
24#include "HFMsgService.h"
25#include <cstdlib>
26
27#include "TFile.h"
28#include "TKey.h"
29#include "TTimeStamp.h"
30
32
33using namespace std;
34
35RooStats::HistFactory::Channel::Channel(std::string ChanName, std::string ChanInputFile) :
36 fName( ChanName ), fInputFile( ChanInputFile )
37{
38 // create channel with given name and input file
39}
40
41namespace RooStats{
42 namespace HistFactory{
43 //BadChannel = Channel();
45 // BadChannel.Name = "BadChannel"; // = Channel(); //.Name = "BadChannel";
46 }
47}
48
49
51{
52 // add fully configured sample to channel
53
54 sample.SetChannelName( GetName() );
55 fSamples.push_back( sample );
56}
57
58void RooStats::HistFactory::Channel::Print( std::ostream& stream ) {
59 // print information of channel to given stream
60
61 stream << "\t Channel Name: " << fName
62 << "\t InputFile: " << fInputFile
63 << std::endl;
64
65 stream << "\t Data:" << std::endl;
66 fData.Print( stream );
67
68
69 stream << "\t statErrorConfig:" << std::endl;
70 fStatErrorConfig.Print( stream );
71
72
73 if( !fSamples.empty() ) {
74
75 stream << "\t Samples: " << std::endl;
76 for( unsigned int i = 0; i < fSamples.size(); ++i ) {
77 fSamples.at(i).Print( stream );
78 }
79 }
80
81
82 stream << "\t End of Channel " << fName << std::endl;
83
84
85}
86
87
88void RooStats::HistFactory::Channel::PrintXML( std::string const& directory, std::string const& prefix ) const {
89
90 // Create an XML file for this channel
91 cxcoutPHF << "Printing XML Files for channel: " << GetName() << std::endl;
92
93 std::string XMLName = prefix + fName + ".xml";
94 if(!directory.empty()) XMLName = directory + "/" + XMLName;
95
96 ofstream xml( XMLName.c_str() );
97
98 // Add the time
99 xml << "<!--" << std::endl;
100 xml << "This xml file created automatically on: " << std::endl;
101 // LM: use TTimeStamp since time_t does not work on Windows
102 TTimeStamp t;
103 UInt_t year = 0;
104 UInt_t month = 0;
105 UInt_t day = 0;
106 t.GetDate(true, 0, &year, &month, &day);
107 xml << year << '-'
108 << month << '-'
109 << day
110 << std::endl;
111 xml << "-->" << std::endl;
112
113 // Add the DOCTYPE
114 xml << "<!DOCTYPE Channel SYSTEM 'HistFactorySchema.dtd'> " << std::endl << std::endl;
115
116 // Add the Channel
117 xml << " <Channel Name=\"" << fName << "\" InputFile=\"" << fInputFile << "\" >" << std::endl << std::endl;
118
119 fData.PrintXML( xml );
120 for(auto const& data : fAdditionalData) {
121 data.PrintXML(xml);
122 }
123
124 fStatErrorConfig.PrintXML( xml );
125 /*
126 xml << " <StatErrorConfig RelErrorThreshold=\"" << fStatErrorConfig.GetRelErrorThreshold() << "\" "
127 << "ConstraintType=\"" << Constraint::Name( fStatErrorConfig.GetConstraintType() ) << "\" "
128 << "/> " << std::endl << std::endl;
129 */
130
131 for(auto const& sample : fSamples) {
132 sample.PrintXML( xml );
133 xml << std::endl << std::endl;
134 }
135
136 xml << std::endl;
137 xml << " </Channel> " << std::endl;
138 xml.close();
139
140 cxcoutPHF << "Finished printing XML files" << std::endl;
141
142}
143
144
145void RooStats::HistFactory::Channel::SetData( std::string DataHistoName, std::string DataInputFile, std::string DataHistoPath ) {
146 // set data for this channel by specifying the name of the histogram,
147 // the external ROOT file and the path to the histogram inside the ROOT file
148
149 fData.SetHistoName( DataHistoName );
150 fData.SetInputFile( DataInputFile );
151 fData.SetHistoPath( DataHistoPath );
152
153}
154
155
156
158 // set data directly to some histogram
159 fData.SetHisto( hData );
160}
161
163
164 // For a NumberCounting measurement only
165 // Set the value of data in a particular channel
166 //
167 // Internally, this simply creates a 1-bin TH1F for you
168
169 std::string DataHistName = fName + "_data";
170
171 // Histogram has 1-bin (hard-coded)
172 TH1F* hData = new TH1F( DataHistName.c_str(), DataHistName.c_str(), 1, 0, 1 );
173 hData->SetBinContent( 1, val );
174
175 // Set the histogram of the internally held data
176 // node of this channel to this newly created histogram
177 SetData( hData );
178
179}
180
181
182void RooStats::HistFactory::Channel::SetStatErrorConfig( double StatRelErrorThreshold, Constraint::Type StatConstraintType ) {
183
184 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
185 fStatErrorConfig.SetConstraintType( StatConstraintType );
186
187}
188
189void RooStats::HistFactory::Channel::SetStatErrorConfig( double StatRelErrorThreshold, std::string StatConstraintType ) {
190
191 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
192 fStatErrorConfig.SetConstraintType( Constraint::GetType(StatConstraintType) );
193
194}
195
196
197
199
200 // Loop through all Samples and Systematics
201 // and collect all necessary histograms
202
203 // Handles to open files for collecting histograms
204 std::map<std::string,std::unique_ptr<TFile>> fileHandles;
205
206 // Get the Data Histogram:
207
208 if( !fData.GetInputFile().empty() ) {
209 fData.SetHisto( GetHistogram(fData.GetInputFile(),
210 fData.GetHistoPath(),
211 fData.GetHistoName(),
212 fileHandles) );
213 }
214
215 // Collect any histograms for additional Datasets
216 for(auto& data : fAdditionalData) {
217 if( !data.GetInputFile().empty() ) {
218 data.SetHisto( GetHistogram(data.GetInputFile(), data.GetHistoPath(), data.GetHistoName(), fileHandles) );
219 }
220 }
221
222 // Get the histograms for the samples:
223 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
224
225 RooStats::HistFactory::Sample& sample = fSamples.at( sampItr );
226
227
228 // Get the nominal histogram:
229 cxcoutDHF << "Collecting Nominal Histogram" << std::endl;
230 TH1* Nominal = GetHistogram(sample.GetInputFile(),
231 sample.GetHistoPath(),
232 sample.GetHistoName(),
233 fileHandles);
234
235 sample.SetHisto( Nominal );
236
237
238 // Get the StatError Histogram (if necessary)
239 if( sample.GetStatError().GetUseHisto() ) {
240 sample.GetStatError().SetErrorHist( GetHistogram(sample.GetStatError().GetInputFile(),
241 sample.GetStatError().GetHistoPath(),
242 sample.GetStatError().GetHistoName(),
243 fileHandles) );
244 }
245
246
247 // Get the HistoSys Variations:
248 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
249
250 RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoSysItr );
251
252 histoSys.SetHistoLow( GetHistogram(histoSys.GetInputFileLow(),
253 histoSys.GetHistoPathLow(),
254 histoSys.GetHistoNameLow(),
255 fileHandles) );
256
257 histoSys.SetHistoHigh( GetHistogram(histoSys.GetInputFileHigh(),
258 histoSys.GetHistoPathHigh(),
259 histoSys.GetHistoNameHigh(),
260 fileHandles) );
261 } // End Loop over HistoSys
262
263
264 // Get the HistoFactor Variations:
265 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
266
267 RooStats::HistFactory::HistoFactor& histoFactor = sample.GetHistoFactorList().at( histoFactorItr );
268
269 histoFactor.SetHistoLow( GetHistogram(histoFactor.GetInputFileLow(),
270 histoFactor.GetHistoPathLow(),
271 histoFactor.GetHistoNameLow(),
272 fileHandles) );
273
274 histoFactor.SetHistoHigh( GetHistogram(histoFactor.GetInputFileHigh(),
275 histoFactor.GetHistoPathHigh(),
276 histoFactor.GetHistoNameHigh(),
277 fileHandles) );
278 } // End Loop over HistoFactor
279
280
281 // Get the ShapeSys Variations:
282 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
283
284 RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeSysItr );
285
286 shapeSys.SetErrorHist( GetHistogram(shapeSys.GetInputFile(),
287 shapeSys.GetHistoPath(),
288 shapeSys.GetHistoName(),
289 fileHandles) );
290 } // End Loop over ShapeSys
291
292
293 // Get any initial shape for a ShapeFactor
294 for( unsigned int shapeFactorItr = 0; shapeFactorItr < sample.GetShapeFactorList().size(); ++shapeFactorItr ) {
295
296 RooStats::HistFactory::ShapeFactor& shapeFactor = sample.GetShapeFactorList().at( shapeFactorItr );
297
298 // Check if we need an InitialShape
299 if( shapeFactor.HasInitialShape() ) {
300 TH1* hist = GetHistogram( shapeFactor.GetInputFile(), shapeFactor.GetHistoPath(),
301 shapeFactor.GetHistoName(), fileHandles );
302 shapeFactor.SetInitialShape( hist );
303 }
304
305 } // End Loop over ShapeFactor
306
307 } // End Loop over Samples
308}
309
310
312
313 // Check that all internal histogram pointers
314 // are properly configured (ie that they're not nullptr)
315
316 if( fData.GetHisto() == nullptr && !fData.GetInputFile().empty() ) {
317 cxcoutEHF << "Error: Data Histogram for channel " << GetName() << " is nullptr." << std::endl;
318 return false;
319 }
320
321 // Get the histograms for the samples:
322 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
323
324 const RooStats::HistFactory::Sample& sample = fSamples.at( sampItr );
325
326 // Get the nominal histogram:
327 if( sample.GetHisto() == nullptr ) {
328 cxcoutEHF << "Error: Nominal Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
329 return false;
330 }
331 else {
332
333 // Check if any bins are negative
334 std::vector<int> NegativeBinNumber;
335 std::vector<double> NegativeBinContent;
336 const TH1* histNominal = sample.GetHisto();
337 for(int ibin=1; ibin<=histNominal->GetNbinsX(); ++ibin) {
338 if(histNominal->GetBinContent(ibin) < 0) {
339 NegativeBinNumber.push_back(ibin);
340 NegativeBinContent.push_back(histNominal->GetBinContent(ibin));
341 }
342 }
343 if(!NegativeBinNumber.empty()) {
344 cxcoutWHF << "WARNING: Nominal Histogram " << histNominal->GetName() << " for Sample = " << sample.GetName()
345 << " in Channel = " << GetName() << " has negative entries in bin numbers = ";
346
347 for(unsigned int ibin=0; ibin<NegativeBinNumber.size(); ++ibin) {
348 if(ibin>0) std::cout << " , " ;
349 std::cout << NegativeBinNumber[ibin] << " : " << NegativeBinContent[ibin] ;
350 }
351 std::cout << std::endl;
352 }
353
354 }
355
356 // Get the StatError Histogram (if necessary)
357 if( sample.GetStatError().GetUseHisto() ) {
358 if( sample.GetStatError().GetErrorHist() == nullptr ) {
359 cxcoutEHF << "Error: Statistical Error Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
360 return false;
361 }
362 }
363
364
365 // Get the HistoSys Variations:
366 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
367
368 const RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoSysItr );
369
370 if( histoSys.GetHistoLow() == nullptr ) {
371 cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoSys.GetName()
372 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
373 return false;
374 }
375 if( histoSys.GetHistoHigh() == nullptr ) {
376 cxcoutEHF << "Error: HistoSyst High for Systematic " << histoSys.GetName()
377 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
378 return false;
379 }
380
381 } // End Loop over HistoSys
382
383
384 // Get the HistoFactor Variations:
385 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
386
387 const RooStats::HistFactory::HistoFactor& histoFactor = sample.GetHistoFactorList().at( histoFactorItr );
388
389 if( histoFactor.GetHistoLow() == nullptr ) {
390 cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoFactor.GetName()
391 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
392 return false;
393 }
394 if( histoFactor.GetHistoHigh() == nullptr ) {
395 cxcoutEHF << "Error: HistoSyst High for Systematic " << histoFactor.GetName()
396 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
397 return false;
398 }
399
400 } // End Loop over HistoFactor
401
402
403 // Get the ShapeSys Variations:
404 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
405
406 const RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeSysItr );
407
408 if( shapeSys.GetErrorHist() == nullptr ) {
409 cxcoutEHF << "Error: HistoSyst High for Systematic " << shapeSys.GetName()
410 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
411 return false;
412 }
413
414 } // End Loop over ShapeSys
415
416 } // End Loop over Samples
417
418 return true;
419}
420
421
422
423/// Open a file and copy a histogram
424/// \param InputFile File where the histogram resides.
425/// \param HistoPath Path of the histogram in the file.
426/// \param HistoName Name of the histogram to retrieve.
427/// \param lsof List of open files. Helps to prevent opening and closing a file hundreds of times.
428TH1* RooStats::HistFactory::Channel::GetHistogram(std::string InputFile, std::string HistoPath, std::string HistoName, std::map<std::string,std::unique_ptr<TFile>>& lsof) {
429
430 cxcoutPHF << "Getting histogram " << InputFile << ":" << HistoPath << "/" << HistoName << std::endl;
431
432 auto& inFile = lsof[InputFile];
433 if (!inFile || !inFile->IsOpen()) {
434 inFile.reset( TFile::Open(InputFile.c_str()) );
435 if ( !inFile || !inFile->IsOpen() ) {
436 cxcoutEHF << "Error: Unable to open input file: " << InputFile << std::endl;
437 throw hf_exc();
438 }
439 cxcoutIHF << "Opened input file: " << InputFile << ": " << std::endl;
440 }
441
442 TDirectory* dir = inFile->GetDirectory(HistoPath.c_str());
443 if (dir == nullptr) {
444 cxcoutEHF << "Histogram path '" << HistoPath
445 << "' wasn't found in file '" << InputFile << "'." << std::endl;
446 throw hf_exc();
447 }
448
449 // Have to read histograms via keys, to ensure that the latest-greatest
450 // name cycle is read from file. Otherwise, they might come from memory.
451 auto key = dir->GetKey(HistoName.c_str());
452 if (key == nullptr) {
453 cxcoutEHF << "Histogram '" << HistoName
454 << "' wasn't found in file '" << InputFile
455 << "' in directory '" << HistoPath << "'." << std::endl;
456 throw hf_exc();
457 }
458
459 auto hist = key->ReadObject<TH1>();
460 if( !hist ) {
461 cxcoutEHF << "Histogram '" << HistoName
462 << "' wasn't found in file '" << InputFile
463 << "' in directory '" << HistoPath << "'." << std::endl;
464 throw hf_exc();
465 }
466
467
468 TH1 * ptr = static_cast<TH1 *>(hist->Clone());
469
470 if(!ptr){
471 std::cerr << "Not all necessary info are set to access the input file. Check your config" << std::endl;
472 std::cerr << "filename: " << InputFile
473 << "path: " << HistoPath
474 << "obj: " << HistoName << std::endl;
475 throw hf_exc();
476 }
477
478 ptr->SetDirectory(nullptr);
479
480
481#ifdef DEBUG
482 std::cout << "Found Histogram: " << HistoName " at address: " << ptr
483 << " with integral " << ptr->Integral() << " and mean " << ptr->GetMean()
484 << std::endl;
485#endif
486
487 // Done
488 return ptr;
489
490}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
This class encapsulates all information for the statistical interpretation of one experiment.
Definition Channel.h:30
void Print(std::ostream &=std::cout)
Definition Channel.cxx:58
void SetData(const RooStats::HistFactory::Data &data)
set data object
Definition Channel.h:53
void AddSample(RooStats::HistFactory::Sample sample)
Definition Channel.cxx:50
TH1 * GetHistogram(std::string InputFile, std::string HistoPath, std::string HistoName, std::map< std::string, std::unique_ptr< TFile > > &lsof)
Open a file and copy a histogram.
Definition Channel.cxx:428
void SetStatErrorConfig(double RelErrorThreshold, Constraint::Type ConstraintType)
Definition Channel.cxx:182
void PrintXML(std::string const &directory, std::string const &prefix="") const
Definition Channel.cxx:88
Configuration for an *un*constrained, coherent shape variation of affected samples.
Configuration for a constrained, coherent shape variation of affected samples.
const std::string & GetHistoNameHigh() const
const std::string & GetHistoNameLow() const
const std::string & GetHistoPathLow() const
const std::string & GetInputFileHigh() const
const std::string & GetInputFileLow() const
const std::string & GetHistoPathHigh() const
std::string GetHistoName() const
get histogram name
Definition Sample.h:92
std::string GetName() const
get name of sample
Definition Sample.h:82
const TH1 * GetHisto() const
Definition Sample.cxx:89
void SetChannelName(const std::string &ChannelName)
set name of associated channel
Definition Sample.h:104
void SetHisto(TH1 *histo)
Definition Sample.h:45
std::string GetHistoPath() const
get histogram path
Definition Sample.h:97
RooStats::HistFactory::StatError & GetStatError()
Definition Sample.h:124
std::vector< RooStats::HistFactory::ShapeFactor > & GetShapeFactorList()
Definition Sample.h:113
std::vector< RooStats::HistFactory::HistoFactor > & GetHistoFactorList()
Definition Sample.h:111
std::string GetInputFile() const
get input ROOT file
Definition Sample.h:87
std::vector< RooStats::HistFactory::HistoSys > & GetHistoSysList()
Definition Sample.h:110
std::vector< RooStats::HistFactory::ShapeSys > & GetShapeSysList()
Definition Sample.h:112
*Un*constrained bin-by-bin variation of affected histogram.
const std::string & GetHistoPath() const
const std::string & GetInputFile() const
const std::string & GetHistoName() const
Constrained bin-by-bin variation of affected histogram.
std::string GetHistoPath() const
const TH1 * GetErrorHist() const
std::string GetHistoName() const
std::string GetInputFile() const
const std::string & GetHistoPath() const
const TH1 * GetErrorHist() const
const std::string & GetInputFile() const
const std::string & GetHistoName() const
Describe directory structure in memory.
Definition TDirectory.h:45
virtual TDirectory * GetDirectory(const char *namecycle, Bool_t printError=false, const char *funcname="GetDirectory")
Find a directory using apath.
virtual TKey * GetKey(const char *, Short_t=9999) const
Definition TDirectory.h:221
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:4070
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual void SetDirectory(TDirectory *dir)
By default, when a histogram is created, it is added to the list of histogram objects in the current ...
Definition TH1.cxx:8902
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition TH1.cxx:7500
virtual Int_t GetNbinsX() const
Definition TH1.h:297
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition TH1.cxx:7906
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9187
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5025
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
The TTimeStamp encapsulates seconds and ns since EPOCH.
Definition TTimeStamp.h:45
UInt_t GetDate(Bool_t inUTC=kTRUE, Int_t secOffset=0, UInt_t *year=nullptr, UInt_t *month=nullptr, UInt_t *day=nullptr) const
Return date in form of 19971224 (i.e.
Type GetType(const std::string &Name)
Namespace for the RooStats classes.
Definition Asimov.h:19