gorad/src/GRScoreWriter.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 //  Gorad (Geant4 Open-source Radiation Analysis and Design)
0027 //
0028 //  Author : Makoto Asai (SLAC National Accelerator Laboratory)
0029 //
0030 //  Development of Gorad is funded by NASA Johnson Space Center (JSC)
0031 //  under the contract NNJ15HK11B.
0032 //
0033 // ********************************************************************
0034 //
0035 // GRScoreWriter.hh
0036 //   Defines the printout format of primitive scorer
0037 //
0038 // History
0039 //   September 8th, 2020 : first implementation
0040 //
0041 // ********************************************************************
0042
0043 #include "GRScoreWriter.hh"
0044
0045 #include <map>
0046 #include <fstream>
0047
0048 #include "G4MultiFunctionalDetector.hh"
0049 #include "G4SDParticleFilter.hh"
0050 #include "G4VPrimitiveScorer.hh"
0051 #include "G4VScoringMesh.hh"
0052
0053 GRScoreWriter::GRScoreWriter()
0054 {;}
0055
0056 GRScoreWriter::~GRScoreWriter()
0057 {;}
0058
0059 void GRScoreWriter::DumpQuantityToFile(const G4String& psName,
0060                                         const G4String& fileName,
0061                                         const G4String& option) {
0062
0063   // change the option string into lowercase to the case-insensitive.
0064   G4String opt = option;
0065   std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
0066
0067   // confirm the option
0068   if(opt.size() == 0) opt = "csv";
0069   if(opt.find("csv") == std::string::npos &&
0070      opt.find("sequence") == std::string::npos) {
0071     G4cerr << "ERROR : DumpToFile : Unknown option -> "
0072            << option << G4endl;
0073     return;
0074   }
0075
0076   // open the file
0077   std::ofstream ofile(fileName);
0078   if(!ofile) {
0079     G4cerr << "ERROR : DumpToFile : File open error -> "
0080            << fileName << G4endl;
0081     return;
0082   }
0083   ofile << "# Mesh or volume name: " << fScoringMesh->GetWorldName() << G4endl;
0084
0085   using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
0086   // retrieve the map
0087   MeshScoreMap fSMap = fScoringMesh->GetScoreMap();
0088
0089
0090   MeshScoreMap::const_iterator msMapItr = fSMap.find(psName);
0091   if(msMapItr == fSMap.end()) {
0092     G4cerr << "ERROR : DumpToFile : Unknown quantity, \""
0093            << psName << "\"." << G4endl;
0094     return;
0095   }
0096
0097   std::map<G4int, G4StatDouble*> * score = msMapItr->second->GetMap();
0098   ofile << "# Primitive scorer name: " << msMapItr->first << G4endl;
0099   if(fact!=1.0)
0100   { ofile << "# Multiplication factor : " << fact << G4endl; }
0101
0102   G4double unitValue = fScoringMesh->GetPSUnitValue(psName);
0103   G4String unit = fScoringMesh->GetPSUnit(psName);
0104   G4String divisionAxisNames[3];
0105   fScoringMesh->GetDivisionAxisNames(divisionAxisNames);
0106   ofile << "# First three integer entries: index of a cell of a mesh, just one cell for a volume tally." << G4endl;
0107   ofile << "# Forth entry: sum of scores." << G4endl;
0108   ofile << "# Fifth entry: sum of squared scores." << G4endl;
0109   ofile << "# Sixth entry: number of events with non-zero effect." << G4endl;
0110   ofile << "# Seventh entry: relative statistical error in %." << G4endl << G4endl;
0111   // index of the cell
0112   ofile << "# i" << divisionAxisNames[0]
0113         << ", i" << divisionAxisNames[1]
0114         << ", i" << divisionAxisNames[2];
0115   // unit of scored value
0116   ofile << ", total(value) ";
0117   if(unit.size() > 0) ofile << "[" << unit << "]";
0118   ofile << ", total(value^2), number of entries, relative error (%)" << G4endl;
0119
0120   // "sequence" option: write header info
0121   if(opt.find("sequence") != std::string::npos) {
0122     ofile << fNMeshSegments[0] << " " << fNMeshSegments[1] << " " << fNMeshSegments[2]
0123           << G4endl;
0124   }
0125
0126   // write quantity
0127   long count = 0;
0128   ofile << std::setprecision(16); // for double value with 8 bytes
0129   for(int x = 0; x < fNMeshSegments[0]; x++) {
0130     for(int y = 0; y < fNMeshSegments[1]; y++) {
0131       for(int z = 0; z < fNMeshSegments[2]; z++) {
0132         G4int idx = GetIndex(x, y, z);
0133
0134         if(opt.find("csv") != std::string::npos)
0135           ofile << x << "," << y << "," << z << ",";
0136
0137         std::map<G4int, G4StatDouble*>::iterator value = score->find(idx);
0138         if(value == score->end()) {
0139           ofile << 0. << "," << 0. << "," << 0;
0140         } else {
0141           G4double x1 = value->second->sum_wx()/unitValue*fact;
0142           G4double x2 = value->second->sum_wx2()/unitValue/unitValue*fact*fact;
0143           G4int n = value->second->n();
0144           // rms  is sigma = sqrt((x2-x1^2/n)/(n-1))
0145           // var = mean +/- rms/sqrt(n): means that the relative error of the sum = rms*sqrt(n)/x1
0146           G4double rms = value->second->rms()/unitValue*fact;
0147           // Relative error in %
0148           G4double relError = rms*std::sqrt(n)*100./x1;
0149           ofile << x1 << ", " << x2 << ", " << n << ", " << relError;
0150           ofile << G4endl;
0151           G4double factor = relError*relError/100.;
0152           if (factor > 1.)
0153           {
0154              G4cout << "# Mesh or volume name: " << fScoringMesh->GetWorldName()
0155                << "  -- # Primitive scorer name: " << msMapItr->first << G4endl
0156                << "  bin " << x << "," << y << "," << z << " : statistical error " << relError << "(%)" << G4endl
0157                << "   to reduce the statistical error below 10%, increase number of events approximately "
0158                << factor << " times." << G4endl;
0159           }
0160         }
0161
0162         if(opt.find("csv") != std::string::npos) {
0163           ofile << G4endl;
0164         } else if(opt.find("sequence") != std::string::npos) {
0165           ofile << " ";
0166           if(count++%5 == 4) ofile << G4endl;
0167         }
0168
0169       } // z
0170     } // y
0171   } // x
0172   ofile << std::setprecision(6);
0173
0174   // close the file
0175   ofile.close();
0176 }
0177
0178 void GRScoreWriter::DumpAllQuantitiesToFile(const G4String& fileName,
0179                                              const G4String& option) {
0180
0181   // change the option string into lowercase to the case-insensitive.
0182   G4String opt = option;
0183   std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
0184
0185   // confirm the option
0186   if(opt.size() == 0) opt = "csv";
0187   if(opt.find("csv") == std::string::npos &&
0188      opt.find("sequence") == std::string::npos) {
0189     G4cerr << "ERROR : DumpToFile : Unknown option -> "
0190            << option << G4endl;
0191     return;
0192   }
0193
0194   // open the file
0195   std::ofstream ofile(fileName);
0196   if(!ofile) {
0197     G4cerr << "ERROR : DumpToFile : File open error -> "
0198            << fileName << G4endl;
0199     return;
0200   }
0201   ofile << "# Mesh or volume name: " << fScoringMesh->GetWorldName() << G4endl;
0202   if(fact!=1.0)
0203   { ofile << "# Multiplication factor : " << fact << G4endl; }
0204   ofile << "# First three integer entries: index of a cell of a mesh, just one cell for a volume tally." << G4endl;
0205   ofile << "# Forth entry: sum of scores." << G4endl;
0206   ofile << "# Fifth entry: sum of squared scores." << G4endl;
0207   ofile << "# Sixth entry: number of events with non-zero effect." << G4endl;
0208   ofile << "# Seventh entry: relative statistical error in %." << G4endl << G4endl;
0209
0210   // retrieve the map
0211   using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
0212   MeshScoreMap fSMap = fScoringMesh->GetScoreMap();
0213   MeshScoreMap::const_iterator msMapItr = fSMap.begin();
0214   std::map<G4int, G4StatDouble*> * score;
0215   for(; msMapItr != fSMap.end(); msMapItr++) {
0216
0217     G4String psname = msMapItr->first;
0218
0219     score = msMapItr->second->GetMap();
0220     ofile << "# Primitive scorer name: " << msMapItr->first << G4endl;
0221
0222     G4double unitValue = fScoringMesh->GetPSUnitValue(psname);
0223     G4String unit = fScoringMesh->GetPSUnit(psname);
0224     G4String divisionAxisNames[3];
0225     fScoringMesh->GetDivisionAxisNames(divisionAxisNames);
0226     // index order
0227     ofile << "# i" << divisionAxisNames[0]
0228           << ", i" << divisionAxisNames[1]
0229           << ", i" << divisionAxisNames[2];
0230     // unit of scored value
0231     ofile << ", total(value) ";
0232     if(unit.size() > 0) ofile << "[" << unit << "]";
0233     ofile << ", total(value^2), number of entries, relative error (%)" << G4endl;
0234
0235
0236     // "sequence" option: write header info
0237     if(opt.find("sequence") != std::string::npos) {
0238       ofile << fNMeshSegments[0] << " " << fNMeshSegments[1] << " " << fNMeshSegments[2]
0239             << G4endl;
0240     }
0241
0242     // write quantity
0243     long count = 0;
0244     ofile << std::setprecision(16); // for double value with 8 bytes
0245     for(int x = 0; x < fNMeshSegments[0]; x++) {
0246       for(int y = 0; y < fNMeshSegments[1]; y++) {
0247         for(int z = 0; z < fNMeshSegments[2]; z++) {
0248           G4int idx = GetIndex(x, y, z);
0249
0250           if(opt.find("csv") != std::string::npos)
0251             ofile << x << "," << y << "," << z << ",";
0252
0253           std::map<G4int, G4StatDouble*>::iterator value = score->find(idx);
0254           if(value == score->end()) {
0255             ofile << 0. << "," << 0. << "," << 0;
0256           } else {
0257             G4double x1 = value->second->sum_wx()/unitValue*fact;
0258             G4double x2 = value->second->sum_wx2()/unitValue/unitValue*fact*fact;
0259             G4int n = value->second->n();
0260             // rms  is sigma = sqrt((x2-x1^2/n)/(n-1))
0261             // var = mean +/- rms/sqrt(n): means that the relative error of the sum = rms*sqrt(n)/x1
0262             G4double rms = value->second->rms()/unitValue*fact;
0263             // Relative error in %
0264             G4double relError = rms*std::sqrt(n)*100./x1;
0265             ofile << x1 << ", " << x2 << ", " << n << ", " << relError;
0266             ofile << G4endl;
0267             G4double factor = relError*relError/100.;
0268             if (factor > 1.)
0269             {
0270                G4cout << "# Mesh or volume name: " << fScoringMesh->GetWorldName()
0271                  << "  -- # Primitive scorer name: " << msMapItr->first << G4endl
0272                  << "  bin " << x << "," << y << "," << z << " : statistical error " << relError << "(%)" << G4endl
0273                  << "   to reduce the statistical error below 10%, increase number of events approximately "
0274                  << factor << " times." << G4endl;
0275             }
0276           }
0277
0278           if(opt.find("csv") != std::string::npos) {
0279             ofile << G4endl;
0280           } else if(opt.find("sequence") != std::string::npos) {
0281             ofile << " ";
0282             if(count++%5 == 4) ofile << G4endl;
0283           }
0284
0285         } // z
0286       } // y
0287     } // x
0288     ofile << std::setprecision(6);
0289
0290   } // for(; msMapItr ....)
0291
0292   // close the file
0293   ofile.close();
0294
0295 }
0296
0297