ThreadsafeScorers/src/TSRunAction.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 TSRunAction.cc
0027 /// \brief Implementation of the TSRunAction class
0028
0029 #include "TSRunAction.hh"
0030
0031 #include "TSActionInitialization.hh"
0032 #include "TSDetectorConstruction.hh"
0033 #include "TSRun.hh"
0034
0035 #include "G4Run.hh"
0036 #include "G4RunManager.hh"
0037 #include "G4StatAnalysis.hh"
0038 #include "G4TaskRunManager.hh"
0039 #include "G4Timer.hh"
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042
0043 TSRunAction::TSRunAction()
0044   : fDetector(TSDetectorConstruction::Instance()), fName(fDetector->GetMFDName())
0045 {}
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 TSRunAction::~TSRunAction() {}
0050
0051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0052
0053 G4Run* TSRunAction::GenerateRun()
0054 {
0055   return new TSRun(fName);
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 void TSRunAction::BeginOfRunAction(const G4Run* aRun)
0061 {
0062   // G4int evts_to_process = aRun->GetNumberOfEventToBeProcessed();
0063   // G4RunManager::GetRunManager()->SetPrintProgress(
0064   //  (evts_to_process > 1000) ? evts_to_process / 1000 : 1);
0065   if (IsMaster() && aRun != nullptr) G4PrintEnv();
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 void TSRunAction::EndOfRunAction(const G4Run* aRun)
0071 {
0072   if (IsMaster()) {
0073     G4cout << " ###### EndOfTSRunAction ###### " << G4endl;
0074
0075     aRun->GetNumberOfEvent();
0076     std::ofstream fileout;
0077     G4String fname = "";
0078     std::stringstream separator;
0079
0080     separator << "============================================================";
0081
0082     typedef std::set<G4int> IDSet_t;
0083     IDSet_t IDs;
0084
0085     //- TSRun object.
0086     const TSRun* tsRun = static_cast<const TSRun*>(aRun);
0087     //--- Dump all scored quantities involved in TSRun.
0088
0089     //---------------------------------------------
0090     // Dump accumulated quantities for this RUN.
0091     //---------------------------------------------
0092     std::vector<G4String> primScorerNames{"EnergyDeposit", "NumberOfSteps"};
0093     std::vector<G4String> fnames{"mfd_tl", "mfd_tg"};
0094     std::vector<G4double> units{CLHEP::eV, CLHEP::keV, 1, 1};
0095     std::vector<G4String> unitstr{"keV", "steps"};
0096
0097     //----------------------------------------------------------------------//
0098     // lambda to print double value
0099     auto print = [](std::ostream& fout, G4int first, G4double second, G4double unit1,
0100                     G4double unit2, G4String unit2str) {
0101       if (fout) fout << first << "     " << second / unit1 << G4endl;
0102
0103       G4cout << "    " << std::setw(10) << first << "    " << std::setw(15) << std::setprecision(6)
0104              << std::fixed << second / unit2 << " " << unit2str << G4endl;
0105       G4cout.unsetf(std::ios::fixed);
0106     };
0107     //----------------------------------------------------------------------//
0108     // lambda to print statistics
0109     auto stat_print = [](std::ostream& fout, G4int first, G4StatAnalysis* stat,
0110                          G4ConvergenceTester* conv, G4double unit1, G4double unit2,
0111                          G4String unit2str) {
0112       if (!stat || !conv) return;
0113       auto fsecond = (*stat);
0114       auto psecond = (*stat);
0115       fsecond /= unit1;
0116       psecond /= unit2;
0117       if (fout) {
0118         fout << first << "     " << fsecond << G4endl;
0119         conv->ShowResult(fout);
0120       }
0121       std::stringstream ss;
0122       ss << "    " << std::setw(10) << first << "    " << std::setw(15) << std::setprecision(6)
0123          << std::fixed << psecond << " " << unit2str;
0124       // skip print of ConvergenceTester to stdout
0125       G4cout << ss.str() << G4endl;
0126     };
0127     //----------------------------------------------------------------------//
0128
0129     for (unsigned i = 0; i < primScorerNames.size(); ++i) {
0130       for (unsigned j = 0; j < fnames.size(); ++j) {
0131         fname = fnames.at(j) + "_" + primScorerNames.at(i) + ".out";
0132         fileout.open(fname);
0133         G4cout << separator.str() << G4endl;
0134         G4cout << " opened file " << fname << " for output" << G4endl;
0135         G4cout << separator.str() << G4endl;
0136
0137         G4bool valid = true;
0138         if (j == 0) {
0139           G4THitsMap<G4double>* hitmap = tsRun->GetHitsMap(fName + "/" + primScorerNames.at(i));
0140           G4StatContainer<G4StatAnalysis>* statmap =
0141             tsRun->GetStatMap(fName + "/" + primScorerNames.at(i));
0142           G4StatContainer<G4ConvergenceTester>* convmap =
0143             tsRun->GetConvMap(fName + "/" + primScorerNames.at(i));
0144
0145           if (hitmap && hitmap->size() != 0) {
0146             for (auto itr = hitmap->begin(); itr != hitmap->end(); itr++) {
0147               if (!hitmap->GetObject(itr)) continue;
0148               IDs.insert(itr->first);
0149               std::get<0>(fTypeCompare[primScorerNames.at(i)][itr->first]) =
0150                 *itr->second / units.at(i);
0151               print(fileout, itr->first, *itr->second, units.at(i), units.at(i + 1), unitstr.at(i));
0152             }
0153           }
0154           else {
0155             valid = false;
0156           }
0157
0158           if (statmap && statmap->size() != 0 && convmap && convmap->size() != 0) {
0159             auto stat_fname = "stat_" + fname;
0160             std::ofstream statout;
0161             statout.open(stat_fname);
0162             for (auto itr = statmap->begin(); itr != statmap->end(); itr++) {
0163               G4int _f = statmap->GetIndex(itr);
0164               G4StatAnalysis* _s = statmap->GetObject(itr);
0165               G4ConvergenceTester* _c = convmap->GetObject(_f);
0166               stat_print(statout, _f, _s, _c, units.at(i), units.at(i + 1), unitstr.at(i));
0167             }
0168             statout.close();
0169           }
0170           else {
0171             std::stringstream ss;
0172             ss << " StatMap/ConvMap is either not "
0173                << "created or the StatMap/ConvMap was empty";
0174             if (statmap) ss << " (StatMap size == " << statmap->size() << ")";
0175             if (convmap) ss << " (ConvMap size == " << convmap->size() << ")";
0176
0177             G4Exception("TSRunAction", "002", JustWarning,
0178                         G4String(primScorerNames.at(i) + ss.str()).c_str());
0179           }
0180
0181           if (!valid) {
0182             G4Exception("TSRunAction", "000", JustWarning,
0183                         G4String(primScorerNames.at(i)
0184                                  + " HitsMap is either not "
0185                                    "created or the HitsMap was empty")
0186                           .c_str());
0187           }
0188         }
0189         else {
0190           G4TAtomicHitsMap<G4double>* hitmap =
0191             tsRun->GetAtomicHitsMap(fName + "/" + primScorerNames.at(i));
0192           if (hitmap && hitmap->size() != 0) {
0193             for (auto itr = hitmap->begin(); itr != hitmap->end(); itr++) {
0194               IDs.insert(itr->first);
0195               std::get<1>(fTypeCompare[primScorerNames.at(i)][itr->first]) =
0196                 *itr->second / units.at(i);
0197               print(fileout, itr->first, *itr->second, units.at(i), units.at(i + 1), unitstr.at(i));
0198             }
0199           }
0200           else {
0201             valid = false;
0202           }
0203
0204           if (!valid) {
0205             G4Exception("TSRunAction", "001", JustWarning,
0206                         G4String(primScorerNames.at(i)
0207                                  + " HitsMap is either not "
0208                                    "created or the HitsMap was empty")
0209                           .c_str());
0210           }
0211         }
0212
0213         fileout.close();
0214         G4cout << separator.str() << G4endl;
0215         G4cout << " closed file " << fname << " for output" << G4endl;
0216       }
0217       // add the mutex data
0218       TSRun::MutexHitsMap_t* hitmap = tsRun->GetMutexHitsMap(fName + "/" + primScorerNames.at(i));
0219       if (hitmap && hitmap->size() != 0) {
0220         for (auto itr = hitmap->begin(); itr != hitmap->end(); itr++) {
0221           IDs.insert(itr->first);
0222           std::get<2>(fTypeCompare[primScorerNames.at(i)][itr->first]) = itr->second / units.at(i);
0223         }
0224       }
0225     }
0226
0227     //--------------------------------------------------------------------//
0228     // Check that the values are equivalent and there are no
0229     // IDs in one container that aren't in another
0230     //--------------------------------------------------------------------//
0231
0232     fname = "mfd_diff.out";
0233     fileout.open(fname);
0234
0235     G4cout << separator.str() << G4endl;
0236     G4cout << " opened file " << fname << " for difference output" << G4endl;
0237     G4cout << separator.str() << G4endl;
0238
0239     fileout << "    " << std::setw(10) << "ID"
0240             << "    " << std::setw(30) << std::setprecision(12) << std::fixed << "MFD value"
0241             << "    " << std::setw(30) << std::setprecision(12) << std::fixed
0242             << "Atomic Hits Map value"
0243             << "    " << std::setw(30) << std::setprecision(8) << std::scientific << "Difference"
0244             << "    " << std::setw(30) << std::setprecision(8) << std::scientific
0245             << "Diff (MFD - MUTEXED)"
0246             << "    " << std::setw(30) << std::setprecision(8) << std::scientific
0247             << "Diff (ATOM_HIT_MAP - MUTEXED)" << G4endl << G4endl;
0248
0249     //----------------------------------------------------------------------//
0250     //
0251     //    Example of using tasking in the user-application. Although this
0252     //    is sort of trivial case and might not result in any speed-up
0253     //    it is a good validation test because the order that strings
0254     //    are joined matters. Although the tasks are executed asynchronously
0255     //    and may complete at different times, the return values from
0256     //    the tasks are stored in futures and are "joined" in the order
0257     //    that they were submitted to the task-group
0258     //
0259     //----------------------------------------------------------------------//
0260     // do not directly call G4TaskManager::GetInstance() as this will generate
0261     // an instance
0262     auto tm = dynamic_cast<G4TaskRunManager*>(G4RunManager::GetRunManager());
0263     // Get the thread-pool if available
0264     auto tp = (tm) ? tm->GetThreadPool() : nullptr;
0265     // write a join algorithm which combines the strings from the tasks
0266     auto join_output = [](std::string& lhs, std::string&& rhs) {
0267       return (lhs += rhs);
0268     };
0269
0270     // this is the outer-loop of tasks
0271     auto report_type_comparison = [=](const G4String& id, const IDcompare_t& comp) {
0272       // the 'report_type_comparison' generates more tasks
0273       auto report_subtype_comparison = [](const G4int& idx, const Compare_t& value) {
0274         std::stringstream streamout;
0275         G4double d01 = std::fabs(std::get<0>(value) - std::get<1>(value));
0276         G4double d02 = std::fabs(std::get<0>(value) - std::get<2>(value));
0277         G4double d03 = std::fabs(std::get<1>(value) - std::get<2>(value));
0278
0279         auto _print_diff = [&](const G4double& _dval) {
0280           if (_dval > 0.0)
0281             streamout << std::setprecision(8) << std::scientific << std::setw(30) << _dval
0282                       << "    ";
0283           else
0284             streamout << std::setprecision(1) << std::fixed << std::setw(30) << _dval << "    ";
0285         };
0286
0287         streamout << "    " << std::setw(10) << idx << "    " << std::setw(30)
0288                   << std::setprecision(12) << std::fixed << std::get<0>(value) << "    "
0289                   << std::setw(30) << std::setprecision(12) << std::fixed << std::get<1>(value)
0290                   << "    ";
0291
0292         _print_diff(d01);
0293         _print_diff(d02);
0294         _print_diff(d03);
0295
0296         streamout << G4endl;
0297         return streamout.str();
0298       };
0299
0300       std::stringstream streamout;
0301       streamout << "\n\nType = " << id << "\n" << G4endl;
0302       if (tp) {
0303         // create a task group (nested inside the 'report_type_comparison' task)
0304         G4TaskGroup<std::string> tg(join_output, tp);
0305         // create the tasks in the task-group
0306         for (auto titr = comp.begin(); titr != comp.end(); ++titr)
0307           tg.exec(report_subtype_comparison, titr->first, titr->second);
0308         // wait on the tasks to finish and execute the join function
0309         // this will block the outer task from completing until all the inner
0310         // tasks have been completed
0311         streamout << tg.join();
0312       }
0313       else {
0314         // if there isn't a tasking thread-pool then we make traditional
0315         // function call on this thread
0316         for (auto titr = comp.begin(); titr != comp.end(); ++titr)
0317           streamout << report_subtype_comparison(titr->first, titr->second);
0318       }
0319       // this is the completion of the outer tasks
0320       return streamout.str();
0321     };
0322
0323     G4String tasking_result = "";
0324     if (tp) {
0325       G4cout << "\n\nGenerating diff output via tasking... ";
0326       // create a task group to
0327       G4TaskGroup<std::string> tg(join_output, tp);
0328       for (auto itr = fTypeCompare.begin(); itr != fTypeCompare.end(); ++itr)
0329         tg.exec(report_type_comparison, itr->first, itr->second);
0330       // wait on the tasks to finish and execute the join function
0331       tasking_result = tg.join();
0332     }
0333
0334     // if thread-pool was available, lets validate that tasking did what was
0335     // expected
0336     if (tp) {
0337       // generate the output serially
0338       G4String serial_result = "";
0339       for (auto itr = fTypeCompare.begin(); itr != fTypeCompare.end(); ++itr)
0340         serial_result += report_type_comparison(itr->first, itr->second);
0341
0342       // write the tasking result even if it was bad so that it can viewed
0343       fileout << tasking_result;
0344
0345       // compare the strings -- should be the same
0346       if (serial_result != tasking_result) {
0347         G4Exception("TSRunAction", "003", JustWarning,
0348                     "Output written via tasking did not match output written "
0349                     "serially. Appending serial result to output file");
0350         fileout << "\n\n#================CORRECT_SERIAL_OUTPUT================#\n\n";
0351         fileout << serial_result;
0352       }
0353     }
0354     else {
0355       // if thread-pool was not available, then just write serially
0356       for (auto itr = fTypeCompare.begin(); itr != fTypeCompare.end(); ++itr)
0357         fileout << report_type_comparison(itr->first, itr->second);
0358     }
0359
0360     fileout.close();
0361     G4cout << " closed file " << fname << " for difference output" << G4endl;
0362     G4cout << separator.str() << G4endl;
0363   }
0364 }
0365
0366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......