| Geant4 Cross Reference |
1 // ******************************************* 1 // ********************************************************************
2 // * License and Disclaimer 2 // * License and Disclaimer *
3 // * 3 // * *
4 // * The Geant4 software is copyright of th 4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided 5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License 6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/ 7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. 8 // * include a list of copyright holders. *
9 // * 9 // * *
10 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file 14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitatio 15 // * for the full disclaimer and the limitation of liability. *
16 // * 16 // * *
17 // * This code implementation is the result 17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboratio 18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distri 19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you ag 20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publicati 21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Sof 22 // * acceptance of all terms of the Geant4 Software license. *
23 // ******************************************* 23 // ********************************************************************
24 // 24 //
25 // ******************************************* 25 // ********************************************************************
26 // 26 //
27 // CaTS (Calorimetry and Tracking Simulation) 27 // CaTS (Calorimetry and Tracking Simulation)
28 // 28 //
29 // Authors : Hans Wenzel 29 // Authors : Hans Wenzel
30 // Soon Yung Jun 30 // Soon Yung Jun
31 // (Fermi National Accelerator Labo 31 // (Fermi National Accelerator Laboratory)
32 // 32 //
33 // History 33 // History
34 // October 18th, 2021 : first implementation 34 // October 18th, 2021 : first implementation
35 // 35 //
36 // ******************************************* 36 // ********************************************************************
37 // 37 //
38 /// \file RadiatorSD.cc 38 /// \file RadiatorSD.cc
39 /// \brief Implementation of the CaTS::Radiato 39 /// \brief Implementation of the CaTS::RadiatorSD class
40 40
41 // Geant4 headers 41 // Geant4 headers
42 #include "G4Step.hh" 42 #include "G4Step.hh"
43 #include "G4Track.hh" 43 #include "G4Track.hh"
44 #ifdef WITH_G4OPTICKS 44 #ifdef WITH_G4OPTICKS
45 # include "G4ThreeVector.hh" 45 # include "G4ThreeVector.hh"
46 # include "G4ios.hh" 46 # include "G4ios.hh"
47 # include "G4UnitsTable.hh" 47 # include "G4UnitsTable.hh"
48 # include "G4SystemOfUnits.hh" 48 # include "G4SystemOfUnits.hh"
49 # include "G4VProcess.hh" 49 # include "G4VProcess.hh"
50 # include "G4VRestDiscreteProcess.hh" 50 # include "G4VRestDiscreteProcess.hh"
51 # include "G4SDManager.hh" 51 # include "G4SDManager.hh"
52 # include "G4HCofThisEvent.hh" 52 # include "G4HCofThisEvent.hh"
53 # include "G4Opticks.hh" 53 # include "G4Opticks.hh"
54 # include "TrackInfo.hh" 54 # include "TrackInfo.hh"
55 # include "OpticksGenstep.h" 55 # include "OpticksGenstep.h"
56 # include "OpticksFlags.hh" 56 # include "OpticksFlags.hh"
57 # include "G4OpticksHit.hh" 57 # include "G4OpticksHit.hh"
58 # include "G4EventManager.hh" 58 # include "G4EventManager.hh"
59 # include "G4Event.hh" 59 # include "G4Event.hh"
60 # include "G4RunManager.hh" 60 # include "G4RunManager.hh"
61 # include "G4Version.hh" 61 # include "G4Version.hh"
62 # include "PhotonSD.hh" 62 # include "PhotonSD.hh"
63 # include "G4Cerenkov.hh" 63 # include "G4Cerenkov.hh"
64 # include "G4Scintillation.hh" 64 # include "G4Scintillation.hh"
65 # include <string> 65 # include <string>
66 #endif 66 #endif
67 // project headers 67 // project headers
68 #include "RadiatorSD.hh" 68 #include "RadiatorSD.hh"
69 #include "ConfigurationManager.hh" 69 #include "ConfigurationManager.hh"
70 70
71 RadiatorSD::RadiatorSD(G4String name) 71 RadiatorSD::RadiatorSD(G4String name)
72 : G4VSensitiveDetector(name) 72 : G4VSensitiveDetector(name)
73 { 73 {
74 verbose = ConfigurationManager::getInstance( 74 verbose = ConfigurationManager::getInstance()->isEnable_verbose();
75 } 75 }
76 76
77 void RadiatorSD::Initialize(G4HCofThisEvent*) 77 void RadiatorSD::Initialize(G4HCofThisEvent*) {}
78 78
79 G4bool RadiatorSD::ProcessHits(G4Step* aStep, 79 G4bool RadiatorSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
80 { 80 {
81 G4double edep = aStep->GetTotalEnergyDeposit 81 G4double edep = aStep->GetTotalEnergyDeposit();
82 if(edep == 0.) 82 if(edep == 0.)
83 return false; 83 return false;
84 // only deal with charged particles 84 // only deal with charged particles
85 G4Track* aTrack = aStep->GetTrack(); 85 G4Track* aTrack = aStep->GetTrack();
86 G4double charge = aTrack->GetDynamicParticle 86 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
87 if(charge == 0) 87 if(charge == 0)
88 return false; 88 return false;
89 #ifdef WITH_G4OPTICKS 89 #ifdef WITH_G4OPTICKS
90 if(ConfigurationManager::getInstance()->isEn 90 if(ConfigurationManager::getInstance()->isEnable_opticks())
91 { 91 {
92 G4int materialIndex = 0; 92 G4int materialIndex = 0;
93 if(first) 93 if(first)
94 { 94 {
95 aMaterial = aTrack->GetMaterial(); 95 aMaterial = aTrack->GetMaterial();
96 materialIndex = aMaterial->GetIndex(); 96 materialIndex = aMaterial->GetIndex();
97 if(verbose) 97 if(verbose)
98 { 98 {
99 G4cout << "*************************** 99 G4cout << "*******************************" << G4endl;
100 G4cout << "RadiatorSD::ProcessHits ini 100 G4cout << "RadiatorSD::ProcessHits initializing Material: "
101 << aMaterial->GetName() << " " 101 << aMaterial->GetName() << " " << G4endl;
102 G4cout << "RadiatorSD::ProcessHits: Na 102 G4cout << "RadiatorSD::ProcessHits: Name "
103 << aStep->GetPreStepPoint() 103 << aStep->GetPreStepPoint()
104 ->GetPhysicalVolume() 104 ->GetPhysicalVolume()
105 ->GetLogicalVolume() 105 ->GetLogicalVolume()
106 ->GetName() 106 ->GetName()
107 << G4endl; 107 << G4endl;
108 } 108 }
109 aMaterialPropertiesTable = aMaterial->Ge 109 aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
110 if(verbose) 110 if(verbose)
111 { 111 {
112 aMaterialPropertiesTable->DumpTable(); 112 aMaterialPropertiesTable->DumpTable();
113 } 113 }
114 // 114 //
115 // properties related to Scintillation 115 // properties related to Scintillation
116 // 116 //
117 # if(G4VERSION_NUMBER > 1072) 117 # if(G4VERSION_NUMBER > 1072)
118 YieldRatio = 118 YieldRatio =
119 aMaterialPropertiesTable->GetConstProp 119 aMaterialPropertiesTable->GetConstProperty(kSCINTILLATIONYIELD1) /
120 aMaterialPropertiesTable->GetConstProp 120 aMaterialPropertiesTable->GetConstProperty(
121 kSCINTILLATIONYIELD2); // slowerRat 121 kSCINTILLATIONYIELD2); // slowerRatio,
122 FastTimeConstant = aMaterialPropertiesTa 122 FastTimeConstant = aMaterialPropertiesTable->GetConstProperty(
123 kSCINTILLATIONTIMECONSTANT1); // Time 123 kSCINTILLATIONTIMECONSTANT1); // TimeConstant,
124 SlowTimeConstant = aMaterialPropertiesTa 124 SlowTimeConstant = aMaterialPropertiesTable->GetConstProperty(
125 kSCINTILLATIONTIMECONSTANT2); // slow 125 kSCINTILLATIONTIMECONSTANT2); // slowerTimeConstant,
126 # else 126 # else
127 Fast_Intensity = aMaterialPropertiesTabl 127 Fast_Intensity = aMaterialPropertiesTable->GetProperty(kFASTCOMPONENT);
128 Slow_Intensity = aMaterialPropertiesTabl 128 Slow_Intensity = aMaterialPropertiesTable->GetProperty(kSLOWCOMPONENT);
129 YieldRatio = aMaterialPropertiesTabl 129 YieldRatio = aMaterialPropertiesTable->GetConstProperty(kYIELDRATIO);
130 # endif 130 # endif
131 ScintillationType = Slow; 131 ScintillationType = Slow;
132 // 132 //
133 // properties related to Cerenkov 133 // properties related to Cerenkov
134 // 134 //
135 Rindex = aMaterialPropertiesTable->GetPr 135 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
136 # if(G4VERSION_NUMBER > 1072) 136 # if(G4VERSION_NUMBER > 1072)
137 Pmin = Rindex->GetMinEnergy(); 137 Pmin = Rindex->GetMinEnergy();
138 Pmax = Rindex->GetMaxEnergy(); 138 Pmax = Rindex->GetMaxEnergy();
139 # else 139 # else
140 Pmin = Rindex->GetMinLowEdgeEn 140 Pmin = Rindex->GetMinLowEdgeEnergy();
141 Pmax = Rindex->GetMaxLowEdgeEn 141 Pmax = Rindex->GetMaxLowEdgeEnergy();
142 # endif 142 # endif
143 dp = Pmax - Pmin; 143 dp = Pmax - Pmin;
144 nMax = Rindex->GetMaxValue(); 144 nMax = Rindex->GetMaxValue();
145 if(verbose) 145 if(verbose)
146 { 146 {
147 G4cout << "nMax: " << nMax << " Pmin: 147 G4cout << "nMax: " << nMax << " Pmin: " << Pmin << " Pmax: " << Pmax
148 << " dp: " << dp << G4endl; 148 << " dp: " << dp << G4endl;
149 Rindex->DumpValues(); 149 Rindex->DumpValues();
150 } 150 }
151 // 151 //
152 first = false; 152 first = false;
153 } // end if first 153 } // end if first
154 G4int Sphotons = 0; // number of scintill 154 G4int Sphotons = 0; // number of scintillation photons this step
155 G4int Cphotons = 0; // number of Cerenkov 155 G4int Cphotons = 0; // number of Cerenkov photons this step
156 // 156 //
157 // info needed for generating Cerenkov pho 157 // info needed for generating Cerenkov photons on the GPU;
158 // 158 //
159 G4double maxCos = 0.0 159 G4double maxCos = 0.0;
160 G4double maxSin2 = 0.0 160 G4double maxSin2 = 0.0;
161 G4double beta = 0.0 161 G4double beta = 0.0;
162 G4double beta1 = 0.0 162 G4double beta1 = 0.0;
163 G4double beta2 = 0.0 163 G4double beta2 = 0.0;
164 G4double BetaInverse = 0.0 164 G4double BetaInverse = 0.0;
165 G4double MeanNumberOfPhotons1 = 0.0 165 G4double MeanNumberOfPhotons1 = 0.0;
166 G4double MeanNumberOfPhotons2 = 0.0 166 G4double MeanNumberOfPhotons2 = 0.0;
167 G4SteppingManager* fpSteppingManager = G4E 167 G4SteppingManager* fpSteppingManager = G4EventManager::GetEventManager()
168 - 168 ->GetTrackingManager()
169 - 169 ->GetSteppingManager();
170 G4StepStatus stepStatus = fpSteppingManage 170 G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
171 if(stepStatus != fAtRestDoItProc) 171 if(stepStatus != fAtRestDoItProc)
172 { 172 {
173 G4ProcessVector* procPost = fpSteppingMa 173 G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector();
174 size_t MAXofPostStepLoops = fpSteppingMa 174 size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops();
175 for(size_t i3 = 0; i3 < MAXofPostStepLoo 175 for(size_t i3 = 0; i3 < MAXofPostStepLoops; i3++)
176 { 176 {
177 if((*procPost)[i3]->GetProcessName() = 177 if((*procPost)[i3]->GetProcessName() == "Cerenkov")
178 { 178 {
179 G4Cerenkov* proc = (G4Cerenkov*) (*p 179 G4Cerenkov* proc = (G4Cerenkov*) (*procPost)[i3];
180 thePhysicsTable = proc->GetPhysicsT 180 thePhysicsTable = proc->GetPhysicsTable();
181 CerenkovAngleIntegrals = 181 CerenkovAngleIntegrals =
182 (G4PhysicsOrderedFreeVector*) ((*t 182 (G4PhysicsOrderedFreeVector*) ((*thePhysicsTable)(materialIndex));
183 Cphotons = proc->GetNumPhotons(); 183 Cphotons = proc->GetNumPhotons();
184 if(Cphotons > 0) 184 if(Cphotons > 0)
185 { 185 {
186 beta1 = aStep->GetPreStepPoi 186 beta1 = aStep->GetPreStepPoint()->GetBeta();
187 beta2 = aStep->GetPostStepPo 187 beta2 = aStep->GetPostStepPoint()->GetBeta();
188 beta = (beta1 + beta2) * 0. 188 beta = (beta1 + beta2) * 0.5;
189 BetaInverse = 1. / beta; 189 BetaInverse = 1. / beta;
190 maxCos = BetaInverse / nMax; 190 maxCos = BetaInverse / nMax;
191 maxSin2 = (1.0 - maxCos) * (1. 191 maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
192 MeanNumberOfPhotons1 = 192 MeanNumberOfPhotons1 =
193 proc->GetAverageNumberOfPhotons( 193 proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex);
194 MeanNumberOfPhotons2 = 194 MeanNumberOfPhotons2 =
195 proc->GetAverageNumberOfPhotons( 195 proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex);
196 } 196 }
197 } 197 }
198 if((*procPost)[i3]->GetProcessName() = 198 if((*procPost)[i3]->GetProcessName() == "Scintillation")
199 { 199 {
200 G4Scintillation* proc1 = (G4Scintill 200 G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
201 Sphotons = proc1->GetN 201 Sphotons = proc1->GetNumPhotons();
202 } 202 }
203 } 203 }
204 } 204 }
205 tSphotons += Sphotons; 205 tSphotons += Sphotons;
206 tCphotons += Cphotons; 206 tCphotons += Cphotons;
207 G4ThreeVector deltaPosition = aStep->GetDe 207 G4ThreeVector deltaPosition = aStep->GetDeltaPosition();
208 G4double ScintillationTime = 0. * ns; 208 G4double ScintillationTime = 0. * ns;
209 G4int scntId = 1; 209 G4int scntId = 1;
210 G4StepPoint* pPreStepPoint = aStep->GetPr 210 G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
211 G4ThreeVector x0 = pPreStepPoin 211 G4ThreeVector x0 = pPreStepPoint->GetPosition();
212 G4ThreeVector p0 = aStep->GetDe 212 G4ThreeVector p0 = aStep->GetDeltaPosition().unit();
213 // 213 //
214 // harvest the Scintillation photon genste 214 // harvest the Scintillation photon gensteps:
215 // 215 //
216 if(Sphotons > 0) 216 if(Sphotons > 0)
217 { 217 {
218 G4double ScintillationRiseTime = 0.0; 218 G4double ScintillationRiseTime = 0.0;
219 G4Opticks::Get()->collectGenstep_G4Scint 219 G4Opticks::Get()->collectGenstep_G4Scintillation_1042(
220 aTrack, aStep, Sphotons, scntId, Scint 220 aTrack, aStep, Sphotons, scntId, ScintillationTime,
221 ScintillationRiseTime); 221 ScintillationRiseTime);
222 } 222 }
223 // 223 //
224 // harvest the Cerenkov photon gensteps: 224 // harvest the Cerenkov photon gensteps:
225 // 225 //
226 if(Cphotons > 0) 226 if(Cphotons > 0)
227 { 227 {
228 G4Opticks::Get()->collectGenstep_G4Ceren 228 G4Opticks::Get()->collectGenstep_G4Cerenkov_1042(
229 aTrack, aStep, Cphotons, BetaInverse, 229 aTrack, aStep, Cphotons, BetaInverse, Pmin, Pmax, maxCos, maxSin2,
230 MeanNumberOfPhotons1, MeanNumberOfPhot 230 MeanNumberOfPhotons1, MeanNumberOfPhotons2);
231 } 231 }
232 G4Opticks* g4ok = G4Opticks::Get(); 232 G4Opticks* g4ok = G4Opticks::Get();
233 G4RunManager* rm = G4RunManager::GetRu 233 G4RunManager* rm = G4RunManager::GetRunManager();
234 const G4Event* event = rm->GetCurrentEvent 234 const G4Event* event = rm->GetCurrentEvent();
235 G4int eventid = event->GetEventID() 235 G4int eventid = event->GetEventID();
236 G4OpticksHit hit; 236 G4OpticksHit hit;
237 unsigned num_photons = g4ok->getNumPhotons 237 unsigned num_photons = g4ok->getNumPhotons();
238 if(num_photons > ConfigurationManager::get 238 if(num_photons > ConfigurationManager::getInstance()->getMaxPhotons())
239 { 239 {
240 g4ok->propagateOpticalPhotons(eventid); 240 g4ok->propagateOpticalPhotons(eventid);
241 G4HCtable* hctable = G4SDManager::GetSDM 241 G4HCtable* hctable = G4SDManager::GetSDMpointer()->GetHCtable();
242 for(G4int i = 0; i < hctable->entries(); 242 for(G4int i = 0; i < hctable->entries(); ++i)
243 { 243 {
244 std::string sdn = hctable->GetSDname 244 std::string sdn = hctable->GetSDname(i);
245 std::size_t found = sdn.find("Photonde 245 std::size_t found = sdn.find("Photondetector");
246 if(found != std::string::npos) 246 if(found != std::string::npos)
247 { 247 {
248 PhotonSD* aSD = 248 PhotonSD* aSD =
249 (PhotonSD*) G4SDManager::GetSDMpoi 249 (PhotonSD*) G4SDManager::GetSDMpointer()->FindSensitiveDetector(
250 sdn); 250 sdn);
251 aSD->AddOpticksHits(); 251 aSD->AddOpticksHits();
252 } 252 }
253 } 253 }
254 g4ok->reset(); 254 g4ok->reset();
255 } 255 }
256 } 256 }
257 #endif 257 #endif
258 return true; 258 return true;
259 } 259 }
260 260
261 void RadiatorSD::EndOfEvent(G4HCofThisEvent*) 261 void RadiatorSD::EndOfEvent(G4HCofThisEvent*)
262 { 262 {
263 tSphotons = 0; 263 tSphotons = 0;
264 tCphotons = 0; 264 tCphotons = 0;
265 } 265 }
266 266