Geant4 Cross Reference |
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 // ******************************************************************** 27 // 28 // CaTS (Calorimetry and Tracking Simulation) 29 // 30 // Authors : Hans Wenzel 31 // Soon Yung Jun 32 // (Fermi National Accelerator Laboratory) 33 // 34 // History 35 // October 18th, 2021 : first implementation 36 // 37 // ******************************************************************** 38 // 39 /// \file lArTPCSD.cc 40 /// \brief Implementation of the CaTS::lArTPCSD class 41 42 // Geant4 headers 43 #include "G4HCofThisEvent.hh" 44 #include "G4Step.hh" 45 #include "G4ThreeVector.hh" 46 #include "G4SDManager.hh" 47 #include "G4ios.hh" 48 #include "G4Track.hh" 49 #ifdef WITH_G4OPTICKS 50 # include "G4Opticks.hh" 51 # include "TrackInfo.hh" 52 # include "OpticksGenstep.h" 53 # include "OpticksFlags.hh" 54 # include "G4OpticksHit.hh" 55 # include "G4Cerenkov.hh" 56 # include "G4Event.hh" 57 # include "G4MaterialPropertiesTable.hh" 58 # include "G4PhysicalConstants.hh" 59 # include "G4RunManager.hh" 60 # include "G4SteppingManager.hh" 61 # include "G4SystemOfUnits.hh" 62 # include "G4UnitsTable.hh" 63 # include "G4VProcess.hh" 64 # include "G4VRestDiscreteProcess.hh" 65 # include "PhotonSD.hh" 66 # include "G4Cerenkov.hh" 67 # include "G4Scintillation.hh" 68 # include "G4Version.hh" 69 #endif 70 // project headers 71 #include "lArTPCSD.hh" 72 #include "ConfigurationManager.hh" 73 74 lArTPCSD::lArTPCSD(G4String name) 75 : G4VSensitiveDetector(name) 76 { 77 G4String HCname = name + "_HC"; 78 collectionName.insert(HCname); 79 verbose = ConfigurationManager::getInstance()->isEnable_verbose(); 80 if(verbose) 81 { 82 G4cout << collectionName.size() << " lArTPCSD name: " << name 83 << " collection Name: " << HCname << G4endl; 84 } 85 fHCID = -1; 86 first = true; 87 } 88 89 void lArTPCSD::Initialize(G4HCofThisEvent* hce) 90 { 91 flArTPCHitsCollection = 92 new lArTPCHitsCollection(SensitiveDetectorName, collectionName[0]); 93 if(fHCID < 0) 94 { 95 if(verbose) 96 { 97 G4cout << "lArTPCSD::Initialize: " << SensitiveDetectorName << " " 98 << collectionName[0] << G4endl; 99 } 100 fHCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); 101 } 102 hce->AddHitsCollection(fHCID, flArTPCHitsCollection); 103 } 104 105 G4bool lArTPCSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) 106 { 107 G4double edep = aStep->GetTotalEnergyDeposit(); 108 if(edep == 0.) 109 return false; 110 // only deal with charged particles 111 G4Track* aTrack = aStep->GetTrack(); 112 G4double charge = aTrack->GetDynamicParticle()->GetCharge(); 113 if(charge == 0) 114 return false; 115 G4double ds = aStep->GetStepLength(); 116 lArTPCHit* newHit = new lArTPCHit( 117 NumElectrons(edep, ds), aStep->GetPostStepPoint()->GetPosition().getX(), 118 aStep->GetPostStepPoint()->GetPosition().getY(), 119 aStep->GetPostStepPoint()->GetPosition().getZ()); 120 flArTPCHitsCollection->insert(newHit); 121 #ifdef WITH_G4OPTICKS 122 if(ConfigurationManager::getInstance()->isEnable_opticks()) 123 { 124 if(first) 125 { 126 aMaterial = aTrack->GetMaterial(); 127 materialIndex = aMaterial->GetIndex(); 128 if(verbose) 129 { 130 G4cout << "*******************************" << G4endl; 131 G4cout << "RadiatorSD::ProcessHits initializing Material: " 132 << aMaterial->GetName() << " " << G4endl; 133 G4cout << "RadiatorSD::ProcessHits: Name " 134 << aStep->GetPreStepPoint() 135 ->GetPhysicalVolume() 136 ->GetLogicalVolume() 137 ->GetName() 138 << G4endl; 139 } 140 aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable(); 141 if(verbose) 142 { 143 aMaterialPropertiesTable->DumpTable(); 144 } 145 // 146 // properties related to Scintillation 147 // 148 # if(G4VERSION_NUMBER > 1072) 149 YieldRatio = 150 aMaterialPropertiesTable->GetConstProperty(kSCINTILLATIONYIELD1) / 151 aMaterialPropertiesTable->GetConstProperty( 152 kSCINTILLATIONYIELD2); // slowerRatio, 153 FastTimeConstant = aMaterialPropertiesTable->GetConstProperty( 154 kSCINTILLATIONTIMECONSTANT1); // TimeConstant, 155 SlowTimeConstant = aMaterialPropertiesTable->GetConstProperty( 156 kSCINTILLATIONTIMECONSTANT2); // slowerTimeConstant, 157 # else 158 Fast_Intensity = aMaterialPropertiesTable->GetProperty(kFASTCOMPONENT); 159 Slow_Intensity = aMaterialPropertiesTable->GetProperty(kSLOWCOMPONENT); 160 YieldRatio = aMaterialPropertiesTable->GetConstProperty(kYIELDRATIO); 161 # endif 162 ScintillationType = Slow; 163 // 164 // properties related to Cerenkov 165 // 166 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX"); 167 # if(G4VERSION_NUMBER > 1072) 168 Pmin = Rindex->GetMinEnergy(); 169 Pmax = Rindex->GetMaxEnergy(); 170 # else 171 Pmin = Rindex->GetMinLowEdgeEnergy(); 172 Pmax = Rindex->GetMaxLowEdgeEnergy(); 173 # endif 174 dp = Pmax - Pmin; 175 if(verbose) 176 { 177 G4cout << "nMax: " << nMax << "Pmin: " << Pmin << "Pmax: " << Pmax 178 << "dp: " << dp << G4endl; 179 Rindex->DumpValues(); 180 } 181 // 182 first = false; 183 } 184 G4int Sphotons = 0; // number of scintillation photons this step 185 G4int Cphotons = 0; // number of Cerenkov photons this step 186 // 187 // info needed for generating Cerenkov photons on the GPU; 188 // 189 G4double maxCos = 0.0; 190 G4double maxSin2 = 0.0; 191 G4double beta = 0.0; 192 G4double beta1 = 0.0; 193 G4double beta2 = 0.0; 194 G4double BetaInverse = 0.0; 195 G4double MeanNumberOfPhotons1 = 0.0; 196 G4double MeanNumberOfPhotons2 = 0.0; 197 G4SteppingManager* fpSteppingManager = G4EventManager::GetEventManager() 198 ->GetTrackingManager() 199 ->GetSteppingManager(); 200 G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus(); 201 if(stepStatus != fAtRestDoItProc) 202 { 203 G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector(); 204 size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops(); 205 for(size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) 206 { 207 if((*procPost)[i3]->GetProcessName() == "Cerenkov") 208 { 209 G4Cerenkov* proc = (G4Cerenkov*) (*procPost)[i3]; 210 thePhysicsTable = proc->GetPhysicsTable(); 211 CerenkovAngleIntegrals = 212 (G4PhysicsOrderedFreeVector*) ((*thePhysicsTable)(materialIndex)); 213 Cphotons = proc->GetNumPhotons(); 214 if(Cphotons > 0) 215 { 216 beta1 = aStep->GetPreStepPoint()->GetBeta(); 217 beta2 = aStep->GetPostStepPoint()->GetBeta(); 218 beta = (beta1 + beta2) * 0.5; 219 BetaInverse = 1. / beta; 220 maxCos = BetaInverse / nMax; 221 maxSin2 = (1.0 - maxCos) * (1.0 + maxCos); 222 MeanNumberOfPhotons1 = 223 proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex); 224 MeanNumberOfPhotons2 = 225 proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex); 226 } 227 } 228 if((*procPost)[i3]->GetProcessName() == "Scintillation") 229 { 230 G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3]; 231 Sphotons = proc1->GetNumPhotons(); 232 } 233 } 234 } 235 tSphotons += Sphotons; 236 tCphotons += Cphotons; 237 G4ThreeVector deltaPosition = aStep->GetDeltaPosition(); 238 G4double ScintillationTime = 0. * ns; 239 G4int scntId = 1; 240 G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint(); 241 G4ThreeVector x0 = pPreStepPoint->GetPosition(); 242 G4ThreeVector p0 = aStep->GetDeltaPosition().unit(); 243 // 244 // harvest the Scintillation photon gensteps: 245 // 246 if(Sphotons > 0) 247 { 248 G4double ScintillationRiseTime = 0.0; 249 G4Opticks::Get()->collectGenstep_G4Scintillation_1042( 250 aTrack, aStep, Sphotons, scntId, ScintillationTime, 251 ScintillationRiseTime); 252 } 253 // 254 // harvest the Cerenkov photon gensteps: 255 // 256 if(Cphotons > 0) 257 { 258 G4Opticks::Get()->collectGenstep_G4Cerenkov_1042( 259 aTrack, aStep, Cphotons, BetaInverse, Pmin, Pmax, maxCos, maxSin2, 260 MeanNumberOfPhotons1, MeanNumberOfPhotons2); 261 } 262 G4Opticks* g4ok = G4Opticks::Get(); 263 G4RunManager* rm = G4RunManager::GetRunManager(); 264 const G4Event* event = rm->GetCurrentEvent(); 265 G4int eventid = event->GetEventID(); 266 G4OpticksHit hit; 267 unsigned num_photons = g4ok->getNumPhotons(); 268 if(num_photons > ConfigurationManager::getInstance()->getMaxPhotons()) 269 { 270 g4ok->propagateOpticalPhotons(eventid); 271 G4HCtable* hctable = G4SDManager::GetSDMpointer()->GetHCtable(); 272 for(G4int i = 0; i < hctable->entries(); ++i) 273 { 274 std::string sdn = hctable->GetSDname(i); 275 std::size_t found = sdn.find("Photondetector"); 276 if(found != std::string::npos) 277 { 278 PhotonSD* aSD = 279 (PhotonSD*) G4SDManager::GetSDMpointer()->FindSensitiveDetector( 280 sdn); 281 aSD->AddOpticksHits(); 282 } 283 } 284 g4ok->reset(); 285 } 286 } 287 #endif 288 return true; 289 } 290 291 void lArTPCSD::EndOfEvent(G4HCofThisEvent*) 292 { 293 tSphotons = 0; 294 tCphotons = 0; 295 G4int NbHits = flArTPCHitsCollection->entries(); 296 if(verbose) 297 { 298 G4cout << " Number of lArTPCHits: " << NbHits << G4endl; 299 } 300 } 301 302 G4double lArTPCSD::NumElectrons(G4double edep, G4double ds) 303 { 304 // Nucl.Instrum.Meth.A523:275-286,2004 305 G4double fGeVToElectrons = 4.237e+07; 306 G4double fModBoxA = 0.930; 307 G4double fModBoxB = 0.212; 308 G4double EFieldStep = 0.5; 309 G4double recomb = 0.0; 310 G4double dEdx = (ds <= 0.0) ? 0.0 : edep / ds; 311 if(dEdx < 1.) 312 dEdx = 1.; 313 if(ds > 0) 314 { 315 G4double Xi = fModBoxB * dEdx / EFieldStep; 316 recomb = std::log(fModBoxA + Xi) / Xi; 317 } 318 else 319 { 320 recomb = 0.0; 321 } 322 G4double fNumIonElectrons = fGeVToElectrons * 1.e-3 * edep * recomb; 323 return fNumIonElectrons; 324 } 325