Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 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........oo 39 40 PulseAction::PulseAction() : G4UserTrackingAct 41 { 42 fpPulseInfo = std::make_unique<PulseInfo>(0) 43 fpMessenger = std::make_unique<PulseActionMe 44 Initialize(); 45 } 46 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 49 PulseInfo::PulseInfo(G4double delayedTime) : G 50 51 //....oooOO0OOooo........oooOO0OOooo........oo 52 53 G4double PulseInfo::GetDelayedTime() const 54 { 55 return fDelayedTime; 56 } 57 58 //....oooOO0OOooo........oooOO0OOooo........oo 59 60 PulseInfo::~PulseInfo() = default; 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 64 PulseAction::~PulseAction() = default; 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 68 void PulseAction::PreUserTrackingAction(const 69 { 70 if (fActivePulse) { 71 if (pTrack->GetParentID() == 0) { 72 fDelayedTime = RandomizeInPulse(); 73 fpPulseInfo = std::make_unique<PulseInfo 74 75 G4cout << "Particle comes at : " << G4Be 76 << G4endl; 77 if (fLonggestDelayedTime < fDelayedTime) 78 fLonggestDelayedTime = fDelayedTime; 79 } 80 } 81 auto pPulseInfo = new PulseInfo(*fpPulseIn 82 ((G4Track*)pTrack)->SetUserInformation(pPu 83 } 84 } 85 86 //....oooOO0OOooo........oooOO0OOooo........oo 87 88 G4double PulseAction::Interpolate(const std::a 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 98 G4double b = std::log10(xs2) - a * std::lo 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 - 107 } 108 return value; 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oo 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 foun 122 G4Exception("PulseAction::Initialize()", " 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........oo 140 141 G4double PulseAction::RandomizeInPulse() 142 { 143 const G4double minTime = 0.; 144 const G4double maxTime = fPulseLarger; // n 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 161 } while (G4UniformRand() * MaximumPulse > Pu 162 163 return selectedPulse; 164 } 165 166 //....oooOO0OOooo........oooOO0OOooo........oo 167 168 double PulseAction::PulseSpectrum(G4double tim 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.begi 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, valu 184 pulse = Interpolate(a); 185 } 186 return pulse; 187 } 188 189 //....oooOO0OOooo........oooOO0OOooo........oo 190 191 G4double PulseAction::GetLonggestDelayedTime() 192 { 193 return fLonggestDelayedTime; 194 } 195 196 //....oooOO0OOooo........oooOO0OOooo........oo 197