UHDR/src/PulseAction.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 PulseAction.cc
0027 /// \brief Implementation of the PulseAction class
0028
0029 // author: hoang tran
0030 // Important : This example is provided as a prototype. Some mistakes may be present in the code.
0031 // We would be happy to hear from you — please don’t hesitate to share your feedback and
0032 // let us know about any difficulties you may encounter.
0033
0034 #include "PulseAction.hh"
0035
0036 #include "PulseActionMessenger.hh"
0037
0038 #include "G4RunManager.hh"
0039 #include "G4Scheduler.hh"
0040 #include "G4String.hh"
0041 #include "G4SystemOfUnits.hh"
0042 #include "G4Track.hh"
0043 #include "G4UIcommand.hh"
0044 #include "G4UnitsTable.hh"
0045 #include "Randomize.hh"
0046
0047 #include <memory>
0048
0049 G4Mutex PulseAction::gUHDRMutex = G4MUTEX_INITIALIZER;  // Le Tuan Anh: for autolock in MT mode
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051
0052 PulseAction::PulseAction(const G4String& pulse, G4bool useHisto)
0053   : G4UserTrackingAction(), fFileName(pulse), fUseHistoInput(useHisto)
0054 {
0055   fpPulseInfo = std::make_unique<PulseInfo>(0);
0056   fpMessenger = std::make_unique<PulseActionMessenger>(this);
0057   if (fUseHistoInput) {
0058     InitializeForHistoInput();
0059   }
0060   else {
0061     Initialize();
0062   }
0063 }
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 PulseInfo::PulseInfo(G4double delayedTime) : G4VUserPulseInfo(), fDelayedTime(delayedTime) {}
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070
0071 G4double PulseInfo::GetDelayedTime() const
0072 {
0073   return fDelayedTime;
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 PulseInfo::~PulseInfo() = default;
0079
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081
0082 PulseAction::~PulseAction() = default;
0083
0084 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0085
0086 void PulseAction::PreUserTrackingAction(const G4Track* pTrack)
0087 {
0088   if (!fActivePulse || fFileName == "") {
0089     return;
0090   }
0091   if (fActivePulse && fFileName == "") {
0092     return;
0093   }
0094   if (fFileName != "" && !fActivePulse) {
0095     return;
0096   }
0097
0098   if (pTrack->GetParentID() == 0) {
0099     fDelayedTime = RandomizeInPulse();
0100     fpPulseInfo = std::make_unique<PulseInfo>(fDelayedTime);
0101     if (fVerbose > 1) {
0102       G4cout << "Particle comes at : " << G4BestUnit(fpPulseInfo->GetDelayedTime(), "Time")
0103              << G4endl;
0104     }
0105     if (fLonggestDelayedTime < fDelayedTime) {
0106       fLonggestDelayedTime = fDelayedTime;
0107     }
0108   }
0109   auto pPulseInfo = new PulseInfo(*fpPulseInfo);
0110   ((G4Track*)pTrack)->SetUserInformation(pPulseInfo);
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 G4double PulseAction::Interpolate(const std::array<G4double, 5>& data)
0116 {
0117   G4double e1 = data[0];
0118   G4double e2 = data[1];
0119   G4double e = data[2];
0120   G4double xs1 = data[3];
0121   G4double xs2 = data[4];
0122   G4double value = 0.;
0123
0124   if ((e2 - e1) != 0) {
0125     G4double d1 = xs1;
0126     G4double d2 = xs2;
0127     value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
0128   }
0129   return value;
0130 }
0131
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133
0134 void PulseAction::Initialize()
0135 {
0136   if (fFileName.empty()) {
0137     return;
0138   }
0139
0140   fPulseVector = {0.};  // Commence avec zéro comme première valeur
0141   fPulseData.clear();
0142   fPulseLarger = 0;
0143   G4MUTEXLOCK(&gUHDRMutex);  // Le Tuan Anh: for autolock in MT mode
0144
0145   std::ifstream inputFile(fFileName);
0146   if (!inputFile) {
0147     G4ExceptionDescription exception;
0148     exception << "pulse Shape file not found. Please, provide : " << fFileName;
0149     G4Exception("PulseAction::Initialize()", "PulseAction01", FatalException, exception);
0150   }
0151
0152   G4double time, pulseAmplitude;
0153   while (inputFile >> time) {
0154     time *= CLHEP::us;  // Conversion en unités
0155
0156     if (!(inputFile >> pulseAmplitude)) {
0157       break;
0158     }
0159     if (time != fPulseVector.back()) {
0160       fPulseVector.push_back(time);
0161     }
0162     fPulseData[time] = pulseAmplitude;
0163     fPulseLarger = std::max(fPulseLarger, time);
0164   }
0165   G4MUTEXUNLOCK(&gUHDRMutex);  // Le Tuan Anh: for autolock in MT mode
0166 }
0167
0168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0169 void PulseAction::InitializeForHistoInput()
0170 {
0171   // L.T. Anh: read histogram for Randomizing pulse using CLHEP::RandGeneral
0172   std::ostringstream FileName;
0173   FileName << fFileName;
0174   if (fFileName == "") {
0175     return;
0176   }
0177   fPulseVector.clear();
0178   G4int nbins = 0;
0179   G4double Tmin = 0, Tmax = 0, pdata;
0180   G4MUTEXLOCK(&gUHDRMutex);
0181   std::ifstream input(FileName.str().c_str());
0182   if (!input.is_open()) {
0183     G4ExceptionDescription exception;
0184     exception << "pulse Shape file " << fFileName
0185               << " not found. Please, provide correct file name!!!";
0186     G4Exception("PulseAction::InitializeForHistoInput()", "PulseAction01", FatalException,
0187                 exception);
0188   }
0189
0190   G4String aline;
0191   while (std::getline(input, aline)) {
0192     if (aline.empty()) continue;
0193     std::istringstream issLine(aline);
0194     G4String firstWord;
0195     issLine >> firstWord;
0196     if (firstWord == "#") continue;
0197     if (firstWord == "nbins") {
0198       issLine >> nbins;
0199     }
0200     else if (firstWord == "Tmin") {
0201       G4String units = "";
0202       issLine >> Tmin >> units;
0203       fTmin = Tmin * G4UIcommand::ValueOf(units);
0204     }
0205     else if (firstWord == "Tmax") {
0206       G4String units = "";
0207       issLine >> Tmax >> units;
0208       fTmax = Tmax * G4UIcommand::ValueOf(units);
0209     }
0210     else {
0211       pdata = std::stod(firstWord);
0212       fPulseVector.push_back(pdata);
0213     }
0214   }
0215   input.close();
0216   G4MUTEXUNLOCK(&gUHDRMutex);
0217   if (nbins != (G4int)fPulseVector.size() || nbins == 0) {
0218     G4ExceptionDescription exception;
0219     exception << "Nbins =  " << nbins;
0220     if (nbins != 0) exception << " not equal to data-size = " << fPulseVector.size();
0221     exception << "!!! \nPlease check the content/format of input file.";
0222     G4Exception("PulseAction::InitializeForHistoInput()", "PulseAction01", FatalException,
0223                 exception);
0224   }
0225   fRandGeneral = std::make_unique<CLHEP::RandGeneral>(&fPulseVector.at(0), nbins);
0226 }
0227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0228
0229 G4double PulseAction::RandomizeInPulse()
0230 {
0231   if (fUseHistoInput) {
0232     // L.T. Anh: Randomizing pulse using CLHEP::RandGeneral
0233     G4double r = fRandGeneral->shoot();
0234     G4double tp = (fTmax - fTmin) * r + fTmin;
0235     return tp;
0236   }
0237   else {
0238     const G4double minTime = 0.;
0239     const G4double maxTime = fPulseLarger;  // us
0240
0241     G4double MaximumPulse = 0.;
0242     G4int nSteps = 50;
0243     G4double value(minTime);
0244
0245     for (G4int i = 0; i < nSteps; i++) {
0246       G4double PulseNumber = PulseSpectrum(value);
0247       if (PulseNumber >= MaximumPulse) {
0248         MaximumPulse = PulseNumber;
0249       }
0250       value += maxTime / nSteps;
0251     }
0252
0253     G4double selectedPulse;
0254     do {
0255       selectedPulse = G4UniformRand() * (maxTime - minTime);
0256     } while (G4UniformRand() * MaximumPulse > PulseSpectrum(selectedPulse));
0257
0258     return selectedPulse;
0259   }
0260 }
0261
0262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0263
0264 double PulseAction::PulseSpectrum(G4double time)
0265 {
0266   G4double pulse = 0.;
0267   G4double valueT1, valueT2, xs1, xs2;
0268   auto t2 = std::upper_bound(fPulseVector.begin(), fPulseVector.end(), time);
0269   auto t1 = t2 - 1;
0270   valueT1 = *t1;
0271   valueT2 = *t2;
0272   xs1 = fPulseData[valueT1];
0273   xs2 = fPulseData[valueT2];
0274   auto xsProduct = xs1 * xs2;
0275   if (xsProduct != 0.) {
0276     std::array<G4double, 5> a = {valueT1, valueT2, time, xs1, xs2};
0277     pulse = Interpolate(a);
0278   }
0279   return pulse;
0280 }
0281
0282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0283
0284 G4double PulseAction::GetLonggestDelayedTime() const
0285 {
0286   return fLonggestDelayedTime;
0287 }
0288
0289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0290 G4double PulseAction::GetPulseLarger() const
0291 {
0292   return fPulseLarger;
0293 }
0294
0295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......