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 ]

Diff markup

Differences between /examples/extended/medical/dna/UHDR/src/PulseAction.cc (Version 11.3.0) and /examples/extended/medical/dna/UHDR/src/PulseAction.cc (Version 11.2.2)


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