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 "G4Track.hh"
0033 #include "G4UnitsTable.hh"
0034 #include "Randomize.hh"
0035
0036 #include <memory>
0037
0038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0039
0040 PulseAction::PulseAction() : G4UserTrackingAction()
0041 {
0042   fpPulseInfo = std::make_unique<PulseInfo>(0);
0043   fpMessenger = std::make_unique<PulseActionMessenger>(this);
0044   Initialize();
0045 }
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 PulseInfo::PulseInfo(G4double delayedTime) : G4VUserPulseInfo(), fDelayedTime(delayedTime) {}
0050
0051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0052
0053 G4double PulseInfo::GetDelayedTime() const
0054 {
0055   return fDelayedTime;
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 PulseInfo::~PulseInfo() = default;
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 PulseAction::~PulseAction() = default;
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 void PulseAction::PreUserTrackingAction(const G4Track* pTrack)
0069 {
0070   if (fActivePulse) {
0071     if (pTrack->GetParentID() == 0) {
0072       fDelayedTime = RandomizeInPulse();
0073       fpPulseInfo = std::make_unique<PulseInfo>(fDelayedTime);
0074
0075       G4cout << "Particle comes at : " << G4BestUnit(fpPulseInfo->GetDelayedTime(), "Time")
0076              << G4endl;
0077       if (fLonggestDelayedTime < fDelayedTime) {
0078         fLonggestDelayedTime = fDelayedTime;
0079       }
0080     }
0081     auto pPulseInfo = new PulseInfo(*fpPulseInfo);
0082     ((G4Track*)pTrack)->SetUserInformation(pPulseInfo);
0083   }
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 G4double PulseAction::Interpolate(const std::array<G4double, 5>& data)
0089 {
0090   G4double e1 = data[0];
0091   G4double e2 = data[1];
0092   G4double e = data[2];
0093   G4double xs1 = data[3];
0094   G4double xs2 = data[4];
0095   G4double value = 0.;
0096   if ((std::log10(e2) - std::log10(e1)) != 0) {
0097     G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
0098     G4double b = std::log10(xs2) - a * std::log10(e2);
0099     G4double sigma = a * std::log10(e) + b;
0100     value = (std::pow(10., sigma));
0101   }
0102
0103   if ((e2 - e1) != 0) {
0104     G4double d1 = xs1;
0105     G4double d2 = xs2;
0106     value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
0107   }
0108   return value;
0109 }
0110
0111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0112
0113 void PulseAction::Initialize()
0114 {
0115   std::ostringstream FileName;
0116   FileName << "pulseShape.dat";
0117   std::ifstream input(FileName.str().c_str());
0118
0119   if (!input.is_open()) {
0120     G4ExceptionDescription exception;
0121     exception << "pulseShape.dat file not found. Please, provide";
0122     G4Exception("PulseAction::Initialize()", "PulseAction01", FatalException, exception);
0123   }
0124
0125   fPulseVector.clear();
0126   fPulseVector.push_back(0.);
0127   while (!input.eof()) {
0128     double aTDummy;
0129     double pTDummy;
0130     input >> aTDummy;
0131     if (aTDummy != fPulseVector.back()) {
0132       fPulseVector.push_back(aTDummy);
0133     }
0134     input >> pTDummy;
0135     fPulseData[aTDummy] = pTDummy;
0136   }
0137 }
0138
0139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0140
0141 G4double PulseAction::RandomizeInPulse()
0142 {
0143   const G4double minTime = 0.;
0144   const G4double maxTime = fPulseLarger;  // ns
0145
0146   G4double MaximumPulse = 0.;
0147   G4int nSteps = 50;
0148   G4double value(minTime);
0149
0150   for (G4int i = 0; i < nSteps; i++) {
0151     G4double PulseNumber = PulseSpectrum(value);
0152     if (PulseNumber >= MaximumPulse) {
0153       MaximumPulse = PulseNumber;
0154     }
0155     value += maxTime / nSteps;
0156   }
0157
0158   G4double selectedPulse = 0.;
0159   do {
0160     selectedPulse = G4UniformRand() * (maxTime - minTime);
0161   } while (G4UniformRand() * MaximumPulse > PulseSpectrum(selectedPulse));
0162
0163   return selectedPulse;
0164 }
0165
0166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0167
0168 double PulseAction::PulseSpectrum(G4double time)
0169 {
0170   G4double pulse = 0.;
0171   G4double valueT1 = 0;
0172   G4double valueT2 = 0;
0173   G4double xs1 = 0;
0174   G4double xs2 = 0;
0175   auto t2 = std::upper_bound(fPulseVector.begin(), fPulseVector.end(), time);
0176   auto t1 = t2 - 1;
0177   valueT1 = *t1;
0178   valueT2 = *t2;
0179   xs1 = fPulseData[valueT1];
0180   xs2 = fPulseData[valueT2];
0181   G4double xsProduct = xs1 * xs2;
0182   if (xsProduct != 0.) {
0183     std::array<G4double, 5> a = {valueT1, valueT2, time, xs1, xs2};
0184     pulse = Interpolate(a);
0185   }
0186   return pulse;
0187 }
0188
0189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0190
0191 G4double PulseAction::GetLonggestDelayedTime() const
0192 {
0193   return fLonggestDelayedTime;
0194 }
0195
0196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......