utils/src/tools_histo_flair.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///  \file tools_histo_flair.cc
0027 ///  \brief Tools to dump any G4H1 into Flair-compatible format.
0028 //
0029 //  Author: G.Hugo, 08 December 2022
0030 //
0031 // ***************************************************************************
0032 //
0033 //      tools_histo_flair
0034 //
0035 ///  Tools to dump any G4H1 into Flair-compatible format.
0036 ///
0037 ///  These tools are fully application-agnostic,
0038 ///  and could also be placed outside of the G4 examples,
0039 ///  as a core G4 Analysis Manager extension.
0040 //
0041 // ***************************************************************************
0042
0043 #include "tools_histo_flair.hh"
0044
0045 #include "G4SystemOfUnits.hh"
0046 #include "G4UnitsTable.hh"
0047 #include "G4ios.hh"
0048
0049 #include <cstring>
0050 #include <fstream>
0051 #include <iomanip>
0052 #include <memory>
0053 #include <stdarg.h>
0054 #include <stdlib.h>
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 namespace tools
0059 {
0060 namespace histo
0061 {
0062 namespace flair
0063 {
0064
0065 // ***************************************************************************
0066 // Returns abscissa name.
0067 // ***************************************************************************
0068 G4String getAbscissaName(const Abscissa abscissaKind)
0069 {
0070   G4String abscissaName = "";
0071   if (abscissaKind == Abscissa::KineticEnergy) {
0072     abscissaName = "kE";
0073   }
0074   else if (abscissaKind == Abscissa::Z) {
0075     abscissaName = "Z";
0076   }
0077   else if (abscissaKind == Abscissa::A) {
0078     abscissaName = "A";
0079   }
0080   return abscissaName;
0081 }
0082
0083 // ***************************************************************************
0084 // Returns abscissa unit.
0085 // ***************************************************************************
0086 G4String getAbscissaUnit(const Abscissa abscissaKind)
0087 {
0088   G4String abscissaUnit = "";
0089   if (abscissaKind == Abscissa::KineticEnergy) {
0090     abscissaUnit = "GeV";
0091   }
0092   return abscissaUnit;
0093 }
0094
0095 // ***************************************************************************
0096 // Returns abscissa unit value.
0097 // ***************************************************************************
0098 G4double getAbscissaUnitValue(const Abscissa abscissaKind)
0099 {
0100   G4double abscissaUnitValue = 1.;
0101   if (abscissaKind == Abscissa::KineticEnergy) {
0102     abscissaUnitValue = G4UnitDefinition::GetValueOf(getAbscissaUnit(Abscissa::KineticEnergy));
0103   }
0104   return abscissaUnitValue;
0105 }
0106
0107 // ***************************************************************************
0108 // String format like printf() for strings.
0109 // This sformat function is extracted as-is from geoviewer,
0110 // and was written by V. Vlachoudis in 2010.
0111 // ***************************************************************************
0112 G4String sformat(const char* fmt, ...)
0113 {
0114   std::unique_ptr<char[]> formatted;
0115   G4int n = (G4int)strlen(fmt) * 2;  // Reserve two times as much as the length of the fmt_str.
0116   while (true) {
0117     // Wrap the plain char array into the unique_ptr.
0118     formatted.reset(new char[n]);
0119     strcpy(&formatted[0], fmt);
0120     va_list ap;
0121     va_start(ap, fmt);
0122     G4int final_n = vsnprintf(&formatted[0], n, fmt, ap);
0123     va_end(ap);
0124     if (final_n < 0 || final_n >= n) {
0125       n += std::abs(final_n - n + 1);
0126     }
0127     else {
0128       break;
0129     }
0130   }
0131   return G4String(formatted.get());
0132 }
0133
0134 // ***************************************************************************
0135 // Returns MC error.
0136 // ***************************************************************************
0137 G4double computeError(const G4double mean, const G4double sumSquares, const G4int numEvents)
0138 {
0139   const G4double var = std::max(0., sumSquares / numEvents - mean * mean);
0140   const G4double err =
0141     ((mean != 0. && numEvents > 1) ? std::sqrt(var / (numEvents - 1)) / mean : 1.);
0142
0143   return err * 100.;
0144 }
0145
0146 // ***************************************************************************
0147 // Dump G4H1 (histo mode) to a Flair-compatible format.
0148 // ***************************************************************************
0149 void dumpG4H1HistoInFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
0150                                 const G4String& histoName, G4H1* const histo,
0151                                 const Abscissa abscissaKind, const G4String& binSchemeName,
0152                                 const G4int numEvents, const G4double sumSquaredEventTotals,
0153                                 const G4double sumSquaredEventInRangeTotals)
0154 {
0155   const G4bool isProfile = false;
0156   dumpG4H1InFlairFormat(output, indexInOutputFile, histoName, histo, abscissaKind, binSchemeName,
0157                         numEvents, isProfile, sumSquaredEventTotals, sumSquaredEventInRangeTotals);
0158 }
0159
0160 // ***************************************************************************
0161 // Dump G4H1 (profile mode, ie no stats) to a Flair-compatible format.
0162 // ***************************************************************************
0163 void dumpG4H1ProfileInFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
0164                                   const G4String& histoName, G4H1* const histo,
0165                                   const Abscissa abscissaKind, const G4String& binSchemeName)
0166 {
0167   const G4bool isProfile = true;
0168   const G4int numEvents = 1;
0169   dumpG4H1InFlairFormat(output, indexInOutputFile, histoName, histo, abscissaKind, binSchemeName,
0170                         numEvents, isProfile);
0171 }
0172
0173 // ***************************************************************************
0174 // Dump G4H1 (profile + histo modes) to a Flair-compatible format.
0175 // ***************************************************************************
0176 void dumpG4H1InFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
0177                            const G4String& histoName, G4H1* const histo,
0178                            const Abscissa abscissaKind, const G4String& binSchemeName,
0179                            const G4int numEvents, const G4bool isProfile,
0180                            const G4double sumSquaredEventTotals,
0181                            const G4double sumSquaredEventInRangeTotals)
0182 {
0183   const G4int numBins = histo->axis().bins();
0184   const G4double minAbscissa = histo->axis().lower_edge();
0185   const G4double maxAbscissa = histo->axis().upper_edge();
0186
0187   const G4String abscissaName = getAbscissaName(abscissaKind);
0188   const G4String abscissaUnit = getAbscissaUnit(abscissaKind);
0189   const G4double abscissaUnitValue = getAbscissaUnitValue(abscissaKind);
0190
0191   // START HEADER.
0192   if (indexInOutputFile != 1) {
0193     output << G4endl << G4endl;
0194   }
0195   output << std::left << std::setw(13) << "# Detector: " << indexInOutputFile << " " << histoName
0196          << G4endl;
0197   output << std::left << std::setw(13) << "# Dim: " << histo->get_dimension() << G4endl;
0198   output << std::left << std::setw(13) << "# Entries: " << numEvents << G4endl;
0199
0200   if (!isProfile) {
0201     // Total integral
0202     const G4double total = static_cast<G4double>(histo->sum_all_bin_heights()) / numEvents;
0203     const G4double totalError = computeError(total, sumSquaredEventTotals, numEvents);
0204     output << std::left << std::setw(13)
0205            << "# Total: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", total, totalError);
0206
0207     // In range integral
0208     const G4double inRange = static_cast<G4double>(histo->sum_bin_heights()) / numEvents;
0209     const G4double inRangeError = computeError(inRange, sumSquaredEventInRangeTotals, numEvents);
0210     output << std::left << std::setw(13)
0211            << "# InRange: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", inRange, inRangeError);
0212
0213     // Underflow
0214     const G4double underflow = static_cast<G4double>(histo->bins_sum_w().at(0)) / numEvents;
0215     const G4double underflowError = computeError(underflow, histo->bins_sum_w2().at(0), numEvents);
0216     output << std::left << std::setw(13)
0217            << "# Under: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", underflow, underflowError);
0218
0219     // Overflow
0220     const G4double overflow =
0221       static_cast<G4double>(histo->bins_sum_w().at(numBins + 1)) / numEvents;
0222     const G4double overflowError =
0223       computeError(overflow, histo->bins_sum_w2().at(numBins + 1), numEvents);
0224     output << std::left << std::setw(13)
0225            << "# Over: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", overflow, overflowError);
0226   }
0227
0228   // END HEADER.
0229   output << std::left << std::setw(13) << "# Log: " << binSchemeName << G4endl;
0230   output << std::left << std::setw(13) << "# Histogram: " << numBins << " "
0231          << minAbscissa / abscissaUnitValue << " " << maxAbscissa / abscissaUnitValue << "  ["
0232          << abscissaUnit << "]" << G4endl;
0233
0234   output << std::left << std::setw(13) << "# Column: "
0235          << "1  " << abscissaName << "_min  "
0236          << "[" << abscissaUnit << "]" << G4endl;
0237   output << std::left << std::setw(13) << "# Column: "
0238          << "2  " << abscissaName << "_max  "
0239          << "[" << abscissaUnit << "]" << G4endl;
0240   output << std::left << std::setw(13) << "# Column: "
0241          << "3  dN/d(" << abscissaName << ")  "
0242          << "[1" << (abscissaUnit.empty() ? "" : ("/" + abscissaUnit)) << "/pr]" << G4endl;
0243   if (!isProfile) {
0244     output << std::left << std::setw(13) << "# Column: "
0245            << "4  error  "
0246            << "[%]" << G4endl;
0247   }
0248
0249   // HISTO / PROFILE DATA.
0250   for (G4int binIndex = 0; binIndex < numBins; ++binIndex) {
0251     const G4double mean = histo->bin_height(binIndex) / numEvents;
0252     const G4double err = computeError(mean, histo->bins_sum_w2().at(binIndex + 1), numEvents);
0253
0254     // histo mode
0255     if (!isProfile) {
0256       output << sformat("%15g %15g %15g %6.2f\n",
0257                         histo->axis().bin_lower_edge(binIndex) / abscissaUnitValue,
0258                         histo->axis().bin_upper_edge(binIndex) / abscissaUnitValue,
0259                         mean / histo->axis().bin_width(binIndex) * abscissaUnitValue, err);
0260     }
0261     // profile mode
0262     else {
0263       output << sformat("%15g %15g %15g\n",
0264                         histo->axis().bin_lower_edge(binIndex) / abscissaUnitValue,
0265                         histo->axis().bin_upper_edge(binIndex) / abscissaUnitValue,
0266                         mean / histo->axis().bin_width(binIndex) * abscissaUnitValue);
0267     }
0268   }  // loop on all bins
0269 }  // dumpG4H1InFlairFormat
0270
0271 }  // namespace flair
0272 }  // namespace histo
0273 }  // namespace tools
0274
0275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......