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