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