Logo ROOT   6.18/05
Reference Guide
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 <stdlib.h>
25
26#include "TFile.h"
27#include "TTimeStamp.h"
28
30
31using namespace std;
32
34 fName( "" )
35{
36 // standard constructor
37}
38
40 fName( other.fName ),
41 fInputFile( other.fInputFile ),
42 fHistoPath( other.fHistoPath ),
43 fData( other.fData ),
44 fAdditionalData( other.fAdditionalData ),
45 fStatErrorConfig( other.fStatErrorConfig ),
46 fSamples( other.fSamples )
47{ ; }
48
49
50RooStats::HistFactory::Channel::Channel(std::string ChanName, std::string ChanInputFile) :
51 fName( ChanName ), fInputFile( ChanInputFile )
52{
53 // create channel with given name and input file
54}
55
56namespace RooStats{
57 namespace HistFactory{
58 //BadChannel = Channel();
60 // BadChannel.Name = "BadChannel"; // = Channel(); //.Name = "BadChannel";
61 }
62}
63
64
66{
67 // add fully configured sample to channel
68
69 sample.SetChannelName( GetName() );
70 fSamples.push_back( sample );
71}
72
73void RooStats::HistFactory::Channel::Print( std::ostream& stream ) {
74 // print information of channel to given stream
75
76 stream << "\t Channel Name: " << fName
77 << "\t InputFile: " << fInputFile
78 << std::endl;
79
80 stream << "\t Data:" << std::endl;
81 fData.Print( stream );
82
83
84 stream << "\t statErrorConfig:" << std::endl;
85 fStatErrorConfig.Print( stream );
86
87
88 if( fSamples.size() != 0 ) {
89
90 stream << "\t Samples: " << std::endl;
91 for( unsigned int i = 0; i < fSamples.size(); ++i ) {
92 fSamples.at(i).Print( stream );
93 }
94 }
95
96
97 stream << "\t End of Channel " << fName << std::endl;
98
99
100}
101
102
103void RooStats::HistFactory::Channel::PrintXML( std::string directory, std::string prefix ) {
104
105 // Create an XML file for this channel
106 std::cout << "Printing XML Files for channel: " << GetName() << std::endl;
107
108 std::string XMLName = prefix + fName + ".xml";
109 if( directory != "" ) XMLName = directory + "/" + XMLName;
110
111 ofstream xml( XMLName.c_str() );
112
113 // Add the time
114 xml << "<!--" << std::endl;
115 xml << "This xml file created automatically on: " << std::endl;
116 // LM: use TTimeStamp since time_t does not work on Windows
117 TTimeStamp t;
118 UInt_t year = 0;
119 UInt_t month = 0;
120 UInt_t day = 0;
121 t.GetDate(true, 0, &year, &month, &day);
122 xml << year << '-'
123 << month << '-'
124 << day
125 << std::endl;
126 xml << "-->" << std::endl;
127
128 // Add the DOCTYPE
129 xml << "<!DOCTYPE Channel SYSTEM 'HistFactorySchema.dtd'> " << std::endl << std::endl;
130
131 // Add the Channel
132 xml << " <Channel Name=\"" << fName << "\" InputFile=\"" << fInputFile << "\" >" << std::endl << std::endl;
133
134 fData.PrintXML( xml );
135 /*
136 xml << " <Data HistoName=\"" << fData.GetHistoName() << "\" "
137 << "InputFile=\"" << fData.GetInputFile() << "\" "
138 << "HistoPath=\"" << fData.GetHistoPath() << "\" "
139 << " /> " << std::endl << std::endl;
140 */
141
142 fStatErrorConfig.PrintXML( xml );
143 /*
144 xml << " <StatErrorConfig RelErrorThreshold=\"" << fStatErrorConfig.GetRelErrorThreshold() << "\" "
145 << "ConstraintType=\"" << Constraint::Name( fStatErrorConfig.GetConstraintType() ) << "\" "
146 << "/> " << std::endl << std::endl;
147 */
148
149 for( unsigned int i = 0; i < fSamples.size(); ++i ) {
150 fSamples.at(i).PrintXML( xml );
151 xml << std::endl << std::endl;
152 }
153
154 xml << std::endl;
155 xml << " </Channel> " << std::endl;
156 xml.close();
157
158 std::cout << "Finished printing XML files" << std::endl;
159
160}
161
162
163
164void RooStats::HistFactory::Channel::SetData( std::string DataHistoName, std::string DataInputFile, std::string DataHistoPath ) {
165 // set data for this channel by specifying the name of the histogram,
166 // the external ROOT file and the path to the histogram inside the ROOT file
167
168 fData.SetHistoName( DataHistoName );
169 fData.SetInputFile( DataInputFile );
170 fData.SetHistoPath( DataHistoPath );
171
172}
173
174
175
177 // set data directly to some histogram
178 fData.SetHisto( hData );
179}
180
182
183 // For a NumberCounting measurement only
184 // Set the value of data in a particular channel
185 //
186 // Internally, this simply creates a 1-bin TH1F for you
187
188 std::string DataHistName = fName + "_data";
189
190 // Histogram has 1-bin (hard-coded)
191 TH1F* hData = new TH1F( DataHistName.c_str(), DataHistName.c_str(), 1, 0, 1 );
192 hData->SetBinContent( 1, val );
193
194 // Set the histogram of the internally held data
195 // node of this channel to this newly created histogram
196 SetData( hData );
197
198}
199
200
201void RooStats::HistFactory::Channel::SetStatErrorConfig( double StatRelErrorThreshold, Constraint::Type StatConstraintType ) {
202
203 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
204 fStatErrorConfig.SetConstraintType( StatConstraintType );
205
206}
207
208void RooStats::HistFactory::Channel::SetStatErrorConfig( double StatRelErrorThreshold, std::string StatConstraintType ) {
209
210 fStatErrorConfig.SetRelErrorThreshold( StatRelErrorThreshold );
211 fStatErrorConfig.SetConstraintType( Constraint::GetType(StatConstraintType) );
212
213}
214
215
216
218
219 // Loop through all Samples and Systematics
220 // and collect all necessary histograms
221
222 // Get the Data Histogram:
223
224 if( fData.GetInputFile() != "" ) {
225 fData.SetHisto( GetHistogram(fData.GetInputFile(),
226 fData.GetHistoPath(),
227 fData.GetHistoName()) );
228 }
229
230 // Collect any histograms for additional Datasets
231 for( unsigned int i=0; i < fAdditionalData.size(); ++i) {
232 RooStats::HistFactory::Data& data = fAdditionalData.at(i);
233 if( data.GetInputFile() != "" ) {
234 data.SetHisto( GetHistogram(data.GetInputFile(), data.GetHistoPath(),data.GetHistoName()) );
235 }
236 }
237
238 // Get the histograms for the samples:
239 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
240
241 RooStats::HistFactory::Sample& sample = fSamples.at( sampItr );
242
243
244 // Get the nominal histogram:
245 std::cout << "Collecting Nominal Histogram" << std::endl;
246 TH1* Nominal = GetHistogram(sample.GetInputFile(),
247 sample.GetHistoPath(),
248 sample.GetHistoName());
249
250 sample.SetHisto( Nominal );
251
252
253 // Get the StatError Histogram (if necessary)
254 if( sample.GetStatError().GetUseHisto() ) {
255 sample.GetStatError().SetErrorHist( GetHistogram(sample.GetStatError().GetInputFile(),
256 sample.GetStatError().GetHistoPath(),
257 sample.GetStatError().GetHistoName()) );
258 }
259
260
261 // Get the HistoSys Variations:
262 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
263
264 RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoSysItr );
265
266 histoSys.SetHistoLow( GetHistogram(histoSys.GetInputFileLow(),
267 histoSys.GetHistoPathLow(),
268 histoSys.GetHistoNameLow()) );
269
270 histoSys.SetHistoHigh( GetHistogram(histoSys.GetInputFileHigh(),
271 histoSys.GetHistoPathHigh(),
272 histoSys.GetHistoNameHigh()) );
273 } // End Loop over HistoSys
274
275
276 // Get the HistoFactor Variations:
277 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
278
279 RooStats::HistFactory::HistoFactor& histoFactor = sample.GetHistoFactorList().at( histoFactorItr );
280
281 histoFactor.SetHistoLow( GetHistogram(histoFactor.GetInputFileLow(),
282 histoFactor.GetHistoPathLow(),
283 histoFactor.GetHistoNameLow()) );
284
285 histoFactor.SetHistoHigh( GetHistogram(histoFactor.GetInputFileHigh(),
286 histoFactor.GetHistoPathHigh(),
287 histoFactor.GetHistoNameHigh()) );
288 } // End Loop over HistoFactor
289
290
291 // Get the ShapeSys Variations:
292 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
293
294 RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeSysItr );
295
296 shapeSys.SetErrorHist( GetHistogram(shapeSys.GetInputFile(),
297 shapeSys.GetHistoPath(),
298 shapeSys.GetHistoName()) );
299 } // End Loop over ShapeSys
300
301
302 // Get any initial shape for a ShapeFactor
303 for( unsigned int shapeFactorItr = 0; shapeFactorItr < sample.GetShapeFactorList().size(); ++shapeFactorItr ) {
304
305 RooStats::HistFactory::ShapeFactor& shapeFactor = sample.GetShapeFactorList().at( shapeFactorItr );
306
307 // Check if we need an InitialShape
308 if( shapeFactor.HasInitialShape() ) {
309 TH1* hist = GetHistogram( shapeFactor.GetInputFile(), shapeFactor.GetHistoPath(),
310 shapeFactor.GetHistoName() );
311 shapeFactor.SetInitialShape( hist );
312 }
313
314 } // End Loop over ShapeFactor
315
316 } // End Loop over Samples
317
318 return;
319
320}
321
322
324
325 // Check that all internal histogram pointers
326 // are properly configured (ie that they're not NULL)
327
328 try {
329
330 if( fData.GetHisto() == NULL && fData.GetInputFile() != "" ) {
331 std::cout << "Error: Data Histogram for channel " << GetName() << " is NULL." << std::endl;
332 throw hf_exc();
333 }
334
335 // Get the histograms for the samples:
336 for( unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr ) {
337
338 RooStats::HistFactory::Sample& sample = fSamples.at( sampItr );
339
340 // Get the nominal histogram:
341 if( sample.GetHisto() == NULL ) {
342 std::cout << "Error: Nominal Histogram for sample " << sample.GetName() << " is NULL." << std::endl;
343 throw hf_exc();
344 }
345 else {
346
347 // Check if any bins are negative
348 std::vector<int> NegativeBinNumber;
349 std::vector<double> NegativeBinContent;
350 const TH1* histNominal = sample.GetHisto();
351 for(int ibin=1; ibin<=histNominal->GetNbinsX(); ++ibin) {
352 if(histNominal->GetBinContent(ibin) < 0) {
353 NegativeBinNumber.push_back(ibin);
354 NegativeBinContent.push_back(histNominal->GetBinContent(ibin));
355 }
356 }
357 if(NegativeBinNumber.size()>0) {
358 std::cout << "WARNING: Nominal Histogram " << histNominal->GetName() << " for Sample = " << sample.GetName()
359 << " in Channel = " << GetName() << " has negative entries in bin numbers = ";
360
361 for(unsigned int ibin=0; ibin<NegativeBinNumber.size(); ++ibin) {
362 if(ibin>0) std::cout << " , " ;
363 std::cout << NegativeBinNumber[ibin] << " : " << NegativeBinContent[ibin] ;
364 }
365 std::cout << std::endl;
366 }
367
368 }
369
370 // Get the StatError Histogram (if necessary)
371 if( sample.GetStatError().GetUseHisto() ) {
372 if( sample.GetStatError().GetErrorHist() == NULL ) {
373 std::cout << "Error: Statistical Error Histogram for sample " << sample.GetName() << " is NULL." << std::endl;
374 throw hf_exc();
375 }
376 }
377
378
379 // Get the HistoSys Variations:
380 for( unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr ) {
381
382 RooStats::HistFactory::HistoSys& histoSys = sample.GetHistoSysList().at( histoSysItr );
383
384 if( histoSys.GetHistoLow() == NULL ) {
385 std::cout << "Error: HistoSyst Low for Systematic " << histoSys.GetName()
386 << " in sample " << sample.GetName() << " is NULL." << std::endl;
387 throw hf_exc();
388 }
389 if( histoSys.GetHistoHigh() == NULL ) {
390 std::cout << "Error: HistoSyst High for Systematic " << histoSys.GetName()
391 << " in sample " << sample.GetName() << " is NULL." << std::endl;
392 throw hf_exc();
393 }
394
395 } // End Loop over HistoSys
396
397
398 // Get the HistoFactor Variations:
399 for( unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr ) {
400
401 RooStats::HistFactory::HistoFactor& histoFactor = sample.GetHistoFactorList().at( histoFactorItr );
402
403 if( histoFactor.GetHistoLow() == NULL ) {
404 std::cout << "Error: HistoSyst Low for Systematic " << histoFactor.GetName()
405 << " in sample " << sample.GetName() << " is NULL." << std::endl;
406 throw hf_exc();
407 }
408 if( histoFactor.GetHistoHigh() == NULL ) {
409 std::cout << "Error: HistoSyst High for Systematic " << histoFactor.GetName()
410 << " in sample " << sample.GetName() << " is NULL." << std::endl;
411 throw hf_exc();
412 }
413
414 } // End Loop over HistoFactor
415
416
417 // Get the ShapeSys Variations:
418 for( unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr ) {
419
420 RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at( shapeSysItr );
421
422 if( shapeSys.GetErrorHist() == NULL ) {
423 std::cout << "Error: HistoSyst High for Systematic " << shapeSys.GetName()
424 << " in sample " << sample.GetName() << " is NULL." << std::endl;
425 throw hf_exc();
426 }
427
428 } // End Loop over ShapeSys
429
430 } // End Loop over Samples
431
432 }
433 catch(std::exception& e)
434 {
435 std::cout << e.what() << std::endl;
436 return false;
437 }
438
439 return true;
440
441
442
443
444}
445
446
447
448
449TH1* RooStats::HistFactory::Channel::GetHistogram(std::string InputFile, std::string HistoPath, std::string HistoName) {
450
451 std::cout << "Getting histogram. "
452 << " InputFile " << InputFile
453 << " HistoPath " << HistoPath
454 << " HistoName " << HistoName
455 << std::endl;
456
457 // TFile* file = TFile::Open( InputFile.c_str() );
458
459 TFile* inFile = TFile::Open( InputFile.c_str() );
460 if( !inFile ) {
461 std::cout << "Error: Unable to open input file: " << InputFile << std::endl;
462 throw hf_exc();
463 }
464
465 std::cout << "Opened input file: " << InputFile << ": " << inFile << std::endl;
466
467 std::string HistNameFull = HistoPath + HistoName;
468
469 if( HistoPath != std::string("") ) {
470 if( HistoPath[ HistoPath.length()-1 ] != std::string("/") ) {
471 std::cout << "WARNING: Histogram path is set to: " << HistoPath
472 << " but it should end with a '/' " << std::endl;
473 std::cout << "Total histogram path is now: " << HistNameFull << std::endl;
474 }
475 }
476
477 TH1* hist = NULL;
478 try{
479 hist = dynamic_cast<TH1*>( inFile->Get( HistNameFull.c_str() ) );
480 }
481 catch(std::exception& e)
482 {
483 std::cout << "Failed to cast object to TH1*" << std::endl;
484 std::cout << e.what() << std::endl;
485 throw hf_exc();
486 }
487 if( !hist ) {
488 std::cout << "Failed to get histogram: " << HistNameFull
489 << " in file: " << InputFile << std::endl;
490 throw hf_exc();
491 }
492
493
494 TH1 * ptr = (TH1 *) hist->Clone();
495
496 if(!ptr){
497 std::cerr << "Not all necessary info are set to access the input file. Check your config" << std::endl;
498 std::cerr << "filename: " << InputFile
499 << "path: " << HistoPath
500 << "obj: " << HistoName << std::endl;
501 throw hf_exc();
502 }
503 else {
504 ptr->SetDirectory(0); // for the current histogram h
505 }
506
507
508#ifdef DEBUG
509 std::cout << "Found Histogram: " << HistoName " at address: " << ptr
510 << " with integral " << ptr->Integral() << " and mean " << ptr->GetMean()
511 << std::endl;
512#endif
513
514
515 inFile->Close();
516
517 // Done
518 return ptr;
519
520}
#define e(i)
Definition: RSha256.hxx:103
unsigned int UInt_t
Definition: RtypesCore.h:42
This class encapsulates all information for the statistical interpretation of one experiment.
Definition: Channel.h:26
void Print(std::ostream &=std::cout)
Definition: Channel.cxx:73
TH1 * GetHistogram(std::string InputFile, std::string HistoPath, std::string HistoName)
Open a file and copy a histogram.
Definition: Channel.cxx:449
void PrintXML(std::string Directory, std::string Prefix="")
Definition: Channel.cxx:103
void SetData(const RooStats::HistFactory::Data &data)
set data object
Definition: Channel.h:50
void AddSample(RooStats::HistFactory::Sample sample)
Definition: Channel.cxx:65
void SetStatErrorConfig(double RelErrorThreshold, Constraint::Type ConstraintType)
Definition: Channel.cxx:201
void SetHisto(TH1 *Hist)
Definition: Data.h:50
std::string GetInputFile()
Definition: Data.h:37
std::string GetHistoName()
Definition: Data.h:40
std::string GetHistoPath()
Definition: Data.h:43
Configuration for an *un*constrained, coherent shape variation of affected samples.
Definition: Systematics.h:215
Configuration for a constrained, coherent shape variation of affected samples.
Definition: Systematics.h:205
const std::string & GetHistoNameHigh() const
Definition: Systematics.h:175
const std::string & GetHistoNameLow() const
Definition: Systematics.h:174
const std::string & GetHistoPathLow() const
Definition: Systematics.h:180
const std::string & GetInputFileHigh() const
Definition: Systematics.h:169
const std::string & GetInputFileLow() const
Definition: Systematics.h:168
const std::string & GetHistoPathHigh() const
Definition: Systematics.h:181
std::string GetHistoName() const
get histogram name
Definition: Sample.h:93
std::string GetName() const
get name of sample
Definition: Sample.h:83
const TH1 * GetHisto() const
Definition: Sample.cxx:99
void SetChannelName(const std::string &ChannelName)
set name of associated channel
Definition: Sample.h:105
void SetHisto(TH1 *histo)
Definition: Sample.h:46
std::string GetHistoPath() const
get histogram path
Definition: Sample.h:98
RooStats::HistFactory::StatError & GetStatError()
Definition: Sample.h:118
std::vector< RooStats::HistFactory::ShapeFactor > & GetShapeFactorList()
Definition: Sample.h:116
std::vector< RooStats::HistFactory::HistoFactor > & GetHistoFactorList()
Definition: Sample.h:113
std::string GetInputFile() const
get input ROOT file
Definition: Sample.h:88
std::vector< RooStats::HistFactory::HistoSys > & GetHistoSysList()
Definition: Sample.h:112
std::vector< RooStats::HistFactory::ShapeSys > & GetShapeSysList()
Definition: Sample.h:115
*Un*constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:267
const std::string & GetHistoPath() const
Definition: Systematics.h:306
const std::string & GetInputFile() const
Definition: Systematics.h:294
const std::string & GetHistoName() const
Definition: Systematics.h:300
Constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:225
std::string GetHistoPath() const
Definition: Systematics.h:243
const TH1 * GetErrorHist() const
Definition: Systematics.h:249
std::string GetHistoName() const
Definition: Systematics.h:240
std::string GetInputFile() const
Definition: Systematics.h:237
void SetErrorHist(TH1 *hError)
Definition: Systematics.h:252
const std::string & GetHistoPath() const
Definition: Systematics.h:346
const TH1 * GetErrorHist() const
Definition: Systematics.h:349
const std::string & GetInputFile() const
Definition: Systematics.h:340
const std::string & GetHistoName() const
Definition: Systematics.h:343
virtual TObject * Get(const char *namecycle)
Return pointer to object identified by namecycle.
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:48
virtual void Close(Option_t *option="")
Close a file.
Definition: TFile.cxx:914
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseGeneralPurpose, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3980
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
The TH1 histogram class.
Definition: TH1.h:56
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8351
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:7050
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2664
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7413
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:8635
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
The TTimeStamp encapsulates seconds and ns since EPOCH.
Definition: TTimeStamp.h:71
UInt_t GetDate(Bool_t inUTC=kTRUE, Int_t secOffset=0, UInt_t *year=0, UInt_t *month=0, UInt_t *day=0) const
Return date in form of 19971224 (i.e.
Definition: TTimeStamp.cxx:353
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:146
Type GetType(const std::string &Name)
Definition: Systematics.cxx:34
Namespace for the RooStats classes.
Definition: Asimov.h:20