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