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#include "TDirectory.h"
31
33
34using std::ofstream;
35
37 fName( ChanName ), fInputFile( ChanInputFile )
38{
39 // create channel with given name and input file
40}
41
42namespace RooStats{
43 namespace HistFactory{
44 //BadChannel = Channel();
46 // BadChannel.Name = "BadChannel"; // = Channel(); //.Name = "BadChannel";
47 }
48}
49
50
52{
53 // add fully configured sample to channel
54
55 sample.SetChannelName( GetName() );
56 fSamples.push_back( sample );
57}
58
59void RooStats::HistFactory::Channel::Print( std::ostream& stream ) {
60 // print information of channel to given stream
61
62 stream << "\t Channel Name: " << fName
63 << "\t InputFile: " << fInputFile
64 << std::endl;
65
66 stream << "\t Data:" << std::endl;
67 fData.Print( stream );
68
69
70 stream << "\t statErrorConfig:" << std::endl;
71 fStatErrorConfig.Print( stream );
72
73
74 if( !fSamples.empty() ) {
75
76 stream << "\t Samples: " << std::endl;
77 for( unsigned int i = 0; i < fSamples.size(); ++i ) {
78 fSamples.at(i).Print( stream );
79 }
80 }
81
82
83 stream << "\t End of Channel " << fName << std::endl;
84
85
86}
87
88
89void RooStats::HistFactory::Channel::PrintXML( std::string const& directory, std::string const& prefix ) const {
90
91 // Create an XML file for this channel
92 cxcoutPHF << "Printing XML Files for channel: " << GetName() << std::endl;
93
94 std::string XMLName = prefix + fName + ".xml";
95 if(!directory.empty()) XMLName = directory + "/" + XMLName;
96
97 ofstream xml( XMLName.c_str() );
98
99 // Add the time
100 xml << "<!--" << std::endl;
101 xml << "This xml file created automatically on: " << std::endl;
102 // LM: use TTimeStamp since time_t does not work on Windows
103 TTimeStamp t;
104 UInt_t year = 0;
105 UInt_t month = 0;
106 UInt_t day = 0;
107 t.GetDate(true, 0, &year, &month, &day);
108 xml << year << '-'
109 << month << '-'
110 << day
111 << std::endl;
112 xml << "-->" << std::endl;
113
114 // Add the DOCTYPE
115 xml << "<!DOCTYPE Channel SYSTEM 'HistFactorySchema.dtd'> " << std::endl << std::endl;
116
117 // Add the Channel
118 xml << " <Channel Name=\"" << fName << "\" InputFile=\"" << fInputFile << "\" >" << std::endl << std::endl;
119
120 fData.PrintXML( xml );
121 for(auto const& data : fAdditionalData) {
122 data.PrintXML(xml);
123 }
124
125 fStatErrorConfig.PrintXML( xml );
126 /*
127 xml << " <StatErrorConfig RelErrorThreshold=\"" << fStatErrorConfig.GetRelErrorThreshold() << "\" "
128 << "ConstraintType=\"" << Constraint::Name( fStatErrorConfig.GetConstraintType() ) << "\" "
129 << "/> " << std::endl << std::endl;
130 */
131
132 for(auto const& sample : fSamples) {
133 sample.PrintXML( xml );
134 xml << std::endl << std::endl;
135 }
136
137 xml << std::endl;
138 xml << " </Channel> " << std::endl;
139 xml.close();
140
141 cxcoutPHF << "Finished printing XML files" << std::endl;
142
143}
144
145
147 // set data for this channel by specifying the name of the histogram,
148 // the external ROOT file and the path to the histogram inside the ROOT file
149
150 fData.SetHistoName( DataHistoName );
151 fData.SetInputFile( DataInputFile );
152 fData.SetHistoPath( DataHistoPath );
153
154}
155
156
157
159 // set data directly to some histogram
160 fData.SetHisto( hData );
161}
162
164
165 // For a NumberCounting measurement only
166 // Set the value of data in a particular channel
167 //
168 // Internally, this simply creates a 1-bin TH1F for you
169
170 std::string DataHistName = fName + "_data";
171
172 // Histogram has 1-bin (hard-coded)
173 TH1F* hData = new TH1F( DataHistName.c_str(), DataHistName.c_str(), 1, 0, 1 );
174 hData->SetBinContent( 1, val );
175
176 // Set the histogram of the internally held data
177 // node of this channel to this newly created histogram
178 SetData( hData );
179
180}
181
182
184
185 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
186 fStatErrorConfig.SetConstraintType( StatConstraintType );
187
188}
189
191
192 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
193 fStatErrorConfig.SetConstraintType( Constraint::GetType(StatConstraintType) );
194
195}
196
197
198
200
201 // Loop through all Samples and Systematics
202 // and collect all necessary histograms
203
204 // Handles to open files for collecting histograms
205 std::map<std::string,std::unique_ptr<TFile>> fileHandles;
206
207 // Get the Data Histogram:
208
209 if( !fData.GetInputFile().empty() ) {
210 fData.SetHisto( GetHistogram(fData.GetInputFile(),
211 fData.GetHistoPath(),
212 fData.GetHistoName(),
213 fileHandles) );
214 }
215
216 // Collect any histograms for additional Datasets
217 for(auto& data : fAdditionalData) {
218 if( !data.GetInputFile().empty() ) {
219 data.SetHisto( GetHistogram(data.GetInputFile(), data.GetHistoPath(), data.GetHistoName(), fileHandles) );
220 }
221 }
222
223 // Get the histograms for the samples:
224 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
225
227
228
229 // Get the nominal histogram:
230 cxcoutDHF << "Collecting Nominal Histogram" << std::endl;
231 TH1* Nominal = GetHistogram(sample.GetInputFile(),
232 sample.GetHistoPath(),
233 sample.GetHistoName(),
235
236 sample.SetHisto( Nominal );
237
238
239 // Get the StatError Histogram (if necessary)
240 if( sample.GetStatError().GetUseHisto() ) {
241 sample.GetStatError().SetErrorHist( GetHistogram(sample.GetStatError().GetInputFile(),
242 sample.GetStatError().GetHistoPath(),
243 sample.GetStatError().GetHistoName(),
244 fileHandles) );
245 }
246
247
248 // Get the HistoSys Variations:
249 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
250
252
253 histoSys.SetHistoLow( GetHistogram(histoSys.GetInputFileLow(),
254 histoSys.GetHistoPathLow(),
255 histoSys.GetHistoNameLow(),
256 fileHandles) );
257
258 histoSys.SetHistoHigh( GetHistogram(histoSys.GetInputFileHigh(),
259 histoSys.GetHistoPathHigh(),
260 histoSys.GetHistoNameHigh(),
261 fileHandles) );
262 } // End Loop over HistoSys
263
264
265 // Get the HistoFactor Variations:
266 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
267
269
270 histoFactor.SetHistoLow( GetHistogram(histoFactor.GetInputFileLow(),
271 histoFactor.GetHistoPathLow(),
272 histoFactor.GetHistoNameLow(),
273 fileHandles) );
274
275 histoFactor.SetHistoHigh( GetHistogram(histoFactor.GetInputFileHigh(),
276 histoFactor.GetHistoPathHigh(),
277 histoFactor.GetHistoNameHigh(),
278 fileHandles) );
279 } // End Loop over HistoFactor
280
281
282 // Get the ShapeSys Variations:
283 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
284
286
287 shapeSys.SetErrorHist( GetHistogram(shapeSys.GetInputFile(),
288 shapeSys.GetHistoPath(),
289 shapeSys.GetHistoName(),
290 fileHandles) );
291 } // End Loop over ShapeSys
292
293
294 // Get any initial shape for a ShapeFactor
295 for( unsigned int shapeFactorItr = 0; shapeFactorItr < sample.GetShapeFactorList().size(); ++shapeFactorItr ) {
296
298
299 // Check if we need an InitialShape
300 if( shapeFactor.HasInitialShape() ) {
301 TH1* hist = GetHistogram( shapeFactor.GetInputFile(), shapeFactor.GetHistoPath(),
302 shapeFactor.GetHistoName(), fileHandles );
303 shapeFactor.SetInitialShape( hist );
304 }
305
306 } // End Loop over ShapeFactor
307
308 } // End Loop over Samples
309}
310
311
313
314 // Check that all internal histogram pointers
315 // are properly configured (ie that they're not nullptr)
316
317 if( fData.GetHisto() == nullptr && !fData.GetInputFile().empty() ) {
318 cxcoutEHF << "Error: Data Histogram for channel " << GetName() << " is nullptr." << std::endl;
319 return false;
320 }
321
322 // Get the histograms for the samples:
323 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
324
325 const RooStats::HistFactory::Sample& sample = fSamples.at( sampItr );
326
327 // Get the nominal histogram:
328 if( sample.GetHisto() == nullptr ) {
329 cxcoutEHF << "Error: Nominal Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
330 return false;
331 }
332 else {
333
334 // Check if any bins are negative
335 std::vector<int> NegativeBinNumber;
336 std::vector<double> NegativeBinContent;
337 const TH1* histNominal = sample.GetHisto();
338 for(int ibin=1; ibin<=histNominal->GetNbinsX(); ++ibin) {
339 if(histNominal->GetBinContent(ibin) < 0) {
340 NegativeBinNumber.push_back(ibin);
341 NegativeBinContent.push_back(histNominal->GetBinContent(ibin));
342 }
343 }
344 if(!NegativeBinNumber.empty()) {
345 cxcoutWHF << "WARNING: Nominal Histogram " << histNominal->GetName() << " for Sample = " << sample.GetName()
346 << " in Channel = " << GetName() << " has negative entries in bin numbers = ";
347
348 for(unsigned int ibin=0; ibin<NegativeBinNumber.size(); ++ibin) {
349 if(ibin>0) std::cout << " , " ;
350 std::cout << NegativeBinNumber[ibin] << " : " << NegativeBinContent[ibin] ;
351 }
352 std::cout << std::endl;
353 }
354
355 }
356
357 // Get the StatError Histogram (if necessary)
358 if( sample.GetStatError().GetUseHisto() ) {
359 if( sample.GetStatError().GetErrorHist() == nullptr ) {
360 cxcoutEHF << "Error: Statistical Error Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
361 return false;
362 }
363 }
364
365
366 // Get the HistoSys Variations:
367 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
368
369 const RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoSysItr );
370
371 if( histoSys.GetHistoLow() == nullptr ) {
372 cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoSys.GetName()
373 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
374 return false;
375 }
376 if( histoSys.GetHistoHigh() == nullptr ) {
377 cxcoutEHF << "Error: HistoSyst High for Systematic " << histoSys.GetName()
378 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
379 return false;
380 }
381
382 } // End Loop over HistoSys
383
384
385 // Get the HistoFactor Variations:
386 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
387
388 const RooStats::HistFactory::HistoFactor& histoFactor = sample.GetHistoFactorList().at( histoFactorItr );
389
390 if( histoFactor.GetHistoLow() == nullptr ) {
391 cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoFactor.GetName()
392 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
393 return false;
394 }
395 if( histoFactor.GetHistoHigh() == nullptr ) {
396 cxcoutEHF << "Error: HistoSyst High for Systematic " << histoFactor.GetName()
397 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
398 return false;
399 }
400
401 } // End Loop over HistoFactor
402
403
404 // Get the ShapeSys Variations:
405 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
406
407 const RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeSysItr );
408
409 if( shapeSys.GetErrorHist() == nullptr ) {
410 cxcoutEHF << "Error: HistoSyst High for Systematic " << shapeSys.GetName()
411 << " in sample " << sample.GetName() << " is nullptr." << std::endl;
412 return false;
413 }
414
415 } // End Loop over ShapeSys
416
417 } // End Loop over Samples
418
419 return true;
420}
421
422
423
424/// Open a file and copy a histogram
425/// \param InputFile File where the histogram resides.
426/// \param HistoPath Path of the histogram in the file.
427/// \param HistoName Name of the histogram to retrieve.
428/// \param lsof List of open files. Helps to prevent opening and closing a file hundreds of times.
429TH1* RooStats::HistFactory::Channel::GetHistogram(std::string InputFile, std::string HistoPath, std::string HistoName, std::map<std::string,std::unique_ptr<TFile>>& lsof) {
430
431 cxcoutPHF << "Getting histogram " << InputFile << ":" << HistoPath << "/" << HistoName << std::endl;
432
433 auto& inFile = lsof[InputFile];
434 if (!inFile || !inFile->IsOpen()) {
435 inFile.reset( TFile::Open(InputFile.c_str()) );
436 if ( !inFile || !inFile->IsOpen() ) {
437 cxcoutEHF << "Error: Unable to open input file: " << InputFile << std::endl;
438 throw hf_exc();
439 }
440 cxcoutIHF << "Opened input file: " << InputFile << ": " << std::endl;
441 }
442
443 TDirectory* dir = inFile->GetDirectory(HistoPath.c_str());
444 if (dir == nullptr) {
445 cxcoutEHF << "Histogram path '" << HistoPath
446 << "' wasn't found in file '" << InputFile << "'." << std::endl;
447 throw hf_exc();
448 }
449
450 // Have to read histograms via keys, to ensure that the latest-greatest
451 // name cycle is read from file. Otherwise, they might come from memory.
452 auto key = dir->GetKey(HistoName.c_str());
453 if (key == nullptr) {
454 cxcoutEHF << "Histogram '" << HistoName
455 << "' wasn't found in file '" << InputFile
456 << "' in directory '" << HistoPath << "'." << std::endl;
457 throw hf_exc();
458 }
459
460 std::unique_ptr<TH1> hist(key->ReadObject<TH1>());
461 if( !hist ) {
462 cxcoutEHF << "Histogram '" << HistoName
463 << "' wasn't found in file '" << InputFile
464 << "' in directory '" << HistoPath << "'." << std::endl;
465 throw hf_exc();
466 }
467
468 TDirectory::TContext ctx{nullptr};
469 TH1 * ptr = static_cast<TH1 *>(hist->Clone());
470
471 if(!ptr){
472 std::cerr << "Not all necessary info are set to access the input file. Check your config" << std::endl;
473 std::cerr << "filename: " << InputFile
474 << "path: " << HistoPath
475 << "obj: " << HistoName << std::endl;
476 throw hf_exc();
477 }
478
479#ifdef DEBUG
480 std::cout << "Found Histogram: " << HistoName " at address: " << ptr
481 << " with integral " << ptr->Integral() << " and mean " << ptr->GetMean()
482 << std::endl;
483#endif
484
485 // Done
486 return ptr;
487
488}
#define cxcoutPHF
#define cxcoutDHF
#define cxcoutIHF
#define cxcoutWHF
#define cxcoutEHF
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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:59
void SetData(const RooStats::HistFactory::Data &data)
set data object
Definition Channel.h:53
void AddSample(RooStats::HistFactory::Sample sample)
Definition Channel.cxx:51
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:429
void SetStatErrorConfig(double RelErrorThreshold, Constraint::Type ConstraintType)
Definition Channel.cxx:183
void PrintXML(std::string const &directory, std::string const &prefix="") const
Definition Channel.cxx:89
Configuration for an *un*constrained, coherent shape variation of affected samples.
Configuration for a constrained, coherent shape variation of affected samples.
*Un*constrained bin-by-bin variation of affected histogram.
Constrained bin-by-bin variation of affected histogram.
TDirectory::TContext keeps track and restore the current directory.
Definition TDirectory.h:89
Describe directory structure in memory.
Definition TDirectory.h:45
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:4121
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:878
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:108
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:7559
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition TH1.cxx:7970
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 CodegenImpl.h:59