Geometry/src/TestHistoMgr.cc

0001 //#define StatTesting
0002 //#define PI121
0003
0004 #include "Validation/Geometry/interface/TestHistoMgr.h"
0005 #include <CLHEP/Units/SystemOfUnits.h>
0006 #ifdef StatTesting
0007 #include "Validation/SimG4GeometryValidation/interface/StatisticsComparator.h"
0008 #include "StatisticsTesting/Chi2ComparisonAlgorithm.h"
0009 #endif
0010
0011 #ifdef PI121
0012 #include "StatisticsTesting/ComparisonResult.h"
0013 #endif
0014 #include <iostream>
0015 #include <cstdlib>
0016
0017 //----------------------------------------------------------------------
0018 TestHistoMgr::TestHistoMgr() {
0019   //  pi_aida::Proxy_Selector::setHistogramType(pi_aida::Histogram_Native);
0020 }
0021
0022 //----------------------------------------------------------------------
0023 TestHistoMgr::~TestHistoMgr() {
0024   for (auto& it : theHistos1) {
0025     delete it.second;
0026   }
0027   for (auto& it : theHistos2) {
0028     delete it.second;
0029   }
0030   for (auto& it : theHistoProfs1) {
0031     delete it.second;
0032   }
0033   for (auto& it : theHistoProfs2) {
0034     delete it.second;
0035   }
0036 }
0037
0038 //----------------------------------------------------------------------
0039 void TestHistoMgr::save(const std::string& name) {
0040   edm::LogInfo("MaterialBudget") << "TestHistoMgr: Save user histos";
0041
0042   TFile fout(name.c_str(), "recreate");
0043
0044   // write out the histos
0045   mih1::const_iterator ite1;
0046   mih2::const_iterator ite2;
0047   mihp1::const_iterator itep1;
0048   mihp2::const_iterator itep2;
0049
0050   for (ite1 = theHistos1.begin(); ite1 != theHistos1.end(); ite1++) {
0051     ((*ite1).second)->Write();
0052   }
0053   for (ite2 = theHistos2.begin(); ite2 != theHistos2.end(); ite2++) {
0054     ((*ite2).second)->Write();
0055   }
0056   for (itep1 = theHistoProfs1.begin(); itep1 != theHistoProfs1.end(); itep1++) {
0057     ((*itep1).second)->Write();
0058   }
0059   for (itep2 = theHistoProfs2.begin(); itep2 != theHistoProfs2.end(); itep2++) {
0060     ((*itep2).second)->Write();
0061   }
0062 }
0063
0064 void TestHistoMgr::openSecondFile(const std::string& name) { theFileRef = std::make_unique<TFile>(name.c_str()); }
0065
0066 void TestHistoMgr::printComparisonResult(int ih) {
0067 #ifdef StatTesting
0068
0069   TH1F* histo2 = getHisto1FromSecondFile(histo1->GetName());
0070
0071   StatisticsTesting::StatisticsComparator<StatisticsTesting::Chi2ComparisonAlgorithm> comparator;
0072
0073 #ifdef PI121
0074
0075   qStatisticsTesting::ComparisonResult result = comparator.compare(*histo1, *histo2);
0076
0077 #else
0078
0079   double result = comparator.compare(*histo1, *histo2);
0080
0081 #endif
0082
0083   // ---------------------------------------------------------
0084   // Do something with (e.g. print out) the result of the test
0085   // ---------------------------------------------------------
0086   edm::LogInfo << "TestHistoMgr: Result of the Chi2 Statistical Test: " << histo1->GetName();
0087 #ifdef PI121
0088   << " distance = " << result.distance() << std::endl
0089   << " ndf = " << result.ndf() << std::endl
0090   << " p-value = " << result.quality() << std::endl
0091 #else
0092   << " p-value = " << result << std::endl
0093 #endif
0094   << "---------------- exampleReadXML  ENDS   ------------------ " << std::endl;
0095 #ifdef PI121
0096   std::cout << "[OVAL]: HISTO_QUALITY " << histo1->GetName() << " " << result.quality() << std::endl;
0097 #else
0098   std::cout << "[OVAL]: HISTO_QUALITY " << histo1->GetName() << " " << result << std::endl;
0099 #endif
0100   std::cout << std::endl
0101             << " mean= " << histo1->GetMean() << " " << histo1->GetName() << " " << histo1->GetTitle() << std::endl;
0102   std::cout << " rms= " << histo1->GetRMS() << " " << histo1->GetName() << " " << histo1->GetTitle() << std::endl;
0103
0104 #else
0105 #endif
0106 }
0107
0108 bool TestHistoMgr::addHisto1(TH1F* sih) {
0109   int ih = atoi(sih->GetName());
0110   theHistos1[ih] = sih;
0111   edm::LogInfo("MaterialBudget") << "TestHistoMgr: addHisto1: " << ih << " = " << sih->GetTitle();
0112   return true;
0113 }
0114
0115 bool TestHistoMgr::addHisto2(TH2F* sih) {
0116   int ih = atoi(sih->GetName());
0117   theHistos2[ih] = sih;
0118   return true;
0119 }
0120
0121 bool TestHistoMgr::addHistoProf1(TProfile* sih) {
0122   int ih = atoi(sih->GetName());
0123   theHistoProfs1[ih] = sih;
0124   edm::LogInfo("MaterialBudget") << "TestHistoMgr: addHistoProf1: " << ih << " = " << sih->GetTitle();
0125   return true;
0126 }
0127
0128 bool TestHistoMgr::addHistoProf2(TProfile2D* sih) {
0129   int ih = atoi(sih->GetName());
0130   theHistoProfs2[ih] = sih;
0131   edm::LogInfo("MaterialBudget") << "TestHistoMgr: addHistoProf2: " << ih << " = " << sih->GetTitle();
0132   return true;
0133 }
0134
0135 TH1F* TestHistoMgr::getHisto1(int ih) {
0136   TH1F* his = nullptr;
0137
0138   mih1::const_iterator ite = theHistos1.find(ih);
0139   if (ite != theHistos1.end()) {
0140     his = (*ite).second;
0141   } else {
0142     edm::LogError("MaterialBudget") << "TestHistoMgr: getHisto1 Histogram does not exist " << ih;
0143     std::exception();
0144   }
0145   return his;
0146 }
0147
0148 TH2F* TestHistoMgr::getHisto2(int ih) {
0149   TH2F* his = nullptr;
0150   mih2::const_iterator ite = theHistos2.find(ih);
0151   if (ite != theHistos2.end()) {
0152     his = (*ite).second;
0153   } else {
0154     edm::LogError("MaterialBudget") << "TestHistoMgr: getHisto2 Histogram does not exist " << ih;
0155     std::exception();
0156   }
0157   return his;
0158 }
0159
0160 TProfile* TestHistoMgr::getHistoProf1(int ih) {
0161   TProfile* his = nullptr;
0162   mihp1::const_iterator ite = theHistoProfs1.find(ih);
0163   if (ite != theHistoProfs1.end()) {
0164     his = (*ite).second;
0165   } else {
0166     edm::LogError("MaterialBudget") << "TestHistoMgr: Profile Histogram 1D does not exist " << ih;
0167     std::exception();
0168   }
0169   return his;
0170 }
0171
0172 TProfile2D* TestHistoMgr::getHistoProf2(int ih) {
0173   TProfile2D* his = nullptr;
0174   mihp2::const_iterator ite = theHistoProfs2.find(ih);
0175   if (ite != theHistoProfs2.end()) {
0176     his = (*ite).second;
0177   } else {
0178     edm::LogError("MaterialBudget") << "TestHistoMgr: Profile Histogram 2D does not exist " << ih;
0179     std::exception();
0180   }
0181   return his;
0182 }
0183
0184 TH1F* TestHistoMgr::getHisto1FromSecondFile(const char* hnam) {
0185   TH1F* his = new TH1F();
0186   if (!theFileRef) {
0187     edm::LogError("MaterialBudget") << "TestHistoMgr: Second file not yet opened ";
0188     std::exception();
0189   } else {
0190     his = (TH1F*)(*theFileRef).Get(hnam);
0191   }
0192
0193   if (!his) {
0194     edm::LogError("MaterialBudget") << "TestHistoMgr: FATAL ERROR Histogram does not exist in second file " << hnam;
0195     theFileRef->ls();
0196     std::exception();
0197   }
0198   return his;
0199 }