Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/UHDR/src/PulseAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // author: hoang tran
 27 
 28 #include "PulseAction.hh"
 29 
 30 #include "PulseActionMessenger.hh"
 31 
 32 #include "G4Track.hh"
 33 #include "G4UnitsTable.hh"
 34 #include "Randomize.hh"
 35 
 36 #include <memory>
 37 
 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 39 
 40 PulseAction::PulseAction() : G4UserTrackingAction()
 41 {
 42   fpPulseInfo = std::make_unique<PulseInfo>(0);
 43   fpMessenger = std::make_unique<PulseActionMessenger>(this);
 44   Initialize();
 45 }
 46 
 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 48 
 49 PulseInfo::PulseInfo(G4double delayedTime) : G4VUserPulseInfo(), fDelayedTime(delayedTime) {}
 50 
 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 52 
 53 G4double PulseInfo::GetDelayedTime() const
 54 {
 55   return fDelayedTime;
 56 }
 57 
 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 59 
 60 PulseInfo::~PulseInfo() = default;
 61 
 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 63 
 64 PulseAction::~PulseAction() = default;
 65 
 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67 
 68 void PulseAction::PreUserTrackingAction(const G4Track* pTrack)
 69 {
 70   if (fActivePulse) {
 71     if (pTrack->GetParentID() == 0) {
 72       fDelayedTime = RandomizeInPulse();
 73       fpPulseInfo = std::make_unique<PulseInfo>(fDelayedTime);
 74 
 75       G4cout << "Particle comes at : " << G4BestUnit(fpPulseInfo->GetDelayedTime(), "Time")
 76              << G4endl;
 77       if (fLonggestDelayedTime < fDelayedTime) {
 78         fLonggestDelayedTime = fDelayedTime;
 79       }
 80     }
 81     auto pPulseInfo = new PulseInfo(*fpPulseInfo);
 82     ((G4Track*)pTrack)->SetUserInformation(pPulseInfo);
 83   }
 84 }
 85 
 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 87 
 88 G4double PulseAction::Interpolate(const std::array<G4double, 5>& data)
 89 {
 90   G4double e1 = data[0];
 91   G4double e2 = data[1];
 92   G4double e = data[2];
 93   G4double xs1 = data[3];
 94   G4double xs2 = data[4];
 95   G4double value = 0.;
 96   if ((std::log10(e2) - std::log10(e1)) != 0) {
 97     G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
 98     G4double b = std::log10(xs2) - a * std::log10(e2);
 99     G4double sigma = a * std::log10(e) + b;
100     value = (std::pow(10., sigma));
101   }
102 
103   if ((e2 - e1) != 0) {
104     G4double d1 = xs1;
105     G4double d2 = xs2;
106     value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
107   }
108   return value;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 void PulseAction::Initialize()
114 {
115   std::ostringstream FileName;
116   FileName << "pulseShape.dat";
117   std::ifstream input(FileName.str().c_str());
118 
119   if (!input.is_open()) {
120     G4ExceptionDescription exception;
121     exception << "pulseShape.dat file not found. Please, provide";
122     G4Exception("PulseAction::Initialize()", "PulseAction01", FatalException, exception);
123   }
124 
125   fPulseVector.clear();
126   fPulseVector.push_back(0.);
127   while (!input.eof()) {
128     double aTDummy;
129     double pTDummy;
130     input >> aTDummy;
131     if (aTDummy != fPulseVector.back()) {
132       fPulseVector.push_back(aTDummy);
133     }
134     input >> pTDummy;
135     fPulseData[aTDummy] = pTDummy;
136   }
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
141 G4double PulseAction::RandomizeInPulse()
142 {
143   const G4double minTime = 0.;
144   const G4double maxTime = fPulseLarger;  // ns
145 
146   G4double MaximumPulse = 0.;
147   G4int nSteps = 50;
148   G4double value(minTime);
149 
150   for (G4int i = 0; i < nSteps; i++) {
151     G4double PulseNumber = PulseSpectrum(value);
152     if (PulseNumber >= MaximumPulse) {
153       MaximumPulse = PulseNumber;
154     }
155     value += maxTime / nSteps;
156   }
157 
158   G4double selectedPulse = 0.;
159   do {
160     selectedPulse = G4UniformRand() * (maxTime - minTime);
161   } while (G4UniformRand() * MaximumPulse > PulseSpectrum(selectedPulse));
162 
163   return selectedPulse;
164 }
165 
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167 
168 double PulseAction::PulseSpectrum(G4double time)
169 {
170   G4double pulse = 0.;
171   G4double valueT1 = 0;
172   G4double valueT2 = 0;
173   G4double xs1 = 0;
174   G4double xs2 = 0;
175   auto t2 = std::upper_bound(fPulseVector.begin(), fPulseVector.end(), time);
176   auto t1 = t2 - 1;
177   valueT1 = *t1;
178   valueT2 = *t2;
179   xs1 = fPulseData[valueT1];
180   xs2 = fPulseData[valueT2];
181   G4double xsProduct = xs1 * xs2;
182   if (xsProduct != 0.) {
183     std::array<G4double, 5> a = {valueT1, valueT2, time, xs1, xs2};
184     pulse = Interpolate(a);
185   }
186   return pulse;
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
191 G4double PulseAction::GetLonggestDelayedTime() const
192 {
193   return fLonggestDelayedTime;
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197