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 // 27 // 28 // G4PSSphereSurfaceFlux 29 #include "G4PSSphereSurfaceFlux.hh" 30 31 #include "G4SystemOfUnits.hh" 32 #include "G4StepStatus.hh" 33 #include "G4Track.hh" 34 #include "G4VSolid.hh" 35 #include "G4VPhysicalVolume.hh" 36 #include "G4VPVParameterisation.hh" 37 #include "G4UnitsTable.hh" 38 #include "G4GeometryTolerance.hh" 39 ////////////////////////////////////////////// 40 // (Description) 41 // This is a primitive scorer class for scor 42 // Flux version assumes only for G4Sphere sha 43 // 44 // Surface is defined at the inside of sphere 45 // Direction -Rmin +Rmax 46 // 0 IN || OUT ->|<- | 47 // 1 IN ->| | 48 // 2 OUT |<- | 49 // 50 // Created: 2005-11-14 Tsukasa ASO, Akinori K 51 // 29-Mar-2007 T.Aso, Bug fix for momentum d 52 // 2010-07-22 Introduce Unit specification. 53 // 2010-07-22 Add weighted and divideByAre o 54 // 2011-02-21 Get correct momentum direction 55 // 2011-09-09 Modify comment in PrintAll(). 56 // 2014-03-03 T.Aso, To use always positive 57 ////////////////////////////////////////////// 58 59 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(c 60 G 61 : G4PSSphereSurfaceFlux(name, direction, "pe 62 {} 63 64 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(c 65 c 66 : G4VPrimitiveScorer(name, depth) 67 , HCID(-1) 68 , fDirection(direction) 69 , EvtMap(nullptr) 70 , weighted(true) 71 , divideByArea(true) 72 { 73 DefineUnitAndCategory(); 74 SetUnit(unit); 75 } 76 77 G4bool G4PSSphereSurfaceFlux::ProcessHits(G4St 78 { 79 G4StepPoint* preStep = aStep->GetPreStepPoin 80 81 G4VPhysicalVolume* physVol = preStep-> 82 G4VPVParameterisation* physParam = physVol-> 83 G4VSolid* solid = nullptr; 84 if(physParam != nullptr) 85 { // for parameterized volume 86 G4int idx = 87 ((G4TouchableHistory*) (aStep->GetPreSte 88 ->GetReplicaNumber(indexDepth); 89 solid = physParam->ComputeSolid(idx, physV 90 solid->ComputeDimensions(physParam, idx, p 91 } 92 else 93 { // for ordinary volume 94 solid = physVol->GetLogicalVolume()->GetSo 95 } 96 97 auto sphereSolid = (G4Sphere*) (solid); 98 99 G4int dirFlag = IsSelectedSurface(aStep, sph 100 if(dirFlag > 0) 101 { 102 if(fDirection == fFlux_InOut || fDirection 103 { 104 G4StepPoint* thisStep = nullptr; 105 if(dirFlag == fFlux_In) 106 { 107 thisStep = preStep; 108 } 109 else if(dirFlag == fFlux_Out) 110 { 111 thisStep = aStep->GetPostStepPoint(); 112 } 113 else 114 { 115 return false; 116 } 117 118 G4TouchableHandle theTouchable = thisSte 119 G4ThreeVector pdirection = thisSte 120 G4ThreeVector localdir = 121 theTouchable->GetHistory()->GetTopTran 122 G4double localdirL2 = localdir.x() * loc 123 localdir.y() * loc 124 localdir.z() * loc 125 G4ThreeVector stppos1 = aStep->GetPreSte 126 G4ThreeVector localpos1 = 127 theTouchable->GetHistory()->GetTopTran 128 G4double localR2 = localpos1.x() * local 129 localpos1.y() * local 130 localpos1.z() * local 131 G4double anglefactor = 132 (localdir.x() * localpos1.x() + locald 133 localdir.z() * localpos1.z()) / 134 std::sqrt(localdirL2) / std::sqrt(loca 135 if(anglefactor < 0.0) 136 anglefactor *= -1.0; 137 138 G4double current = 1.0 / anglefactor; 139 if(weighted) 140 current *= thisStep->GetWeight(); // 141 if(divideByArea) // 142 { 143 G4double radi = sphereSolid->GetInnerR 144 G4double dph = sphereSolid->GetDeltaP 145 G4double stth = sphereSolid->GetStartT 146 G4double enth = stth + sphereSolid->Ge 147 current /= radi * radi * dph * (-std:: 148 } 149 150 G4int index = GetIndex(aStep); 151 EvtMap->add(index, current); 152 } 153 } 154 155 return true; 156 } 157 158 G4int G4PSSphereSurfaceFlux::IsSelectedSurface 159 160 { 161 G4TouchableHandle theTouchable = 162 aStep->GetPreStepPoint()->GetTouchableHand 163 G4double kCarTolerance = 164 G4GeometryTolerance::GetInstance()->GetSur 165 166 if(aStep->GetPreStepPoint()->GetStepStatus() 167 { 168 // Entering Geometry 169 G4ThreeVector stppos1 = aStep->GetPreStepP 170 G4ThreeVector localpos1 = 171 theTouchable->GetHistory()->GetTopTransf 172 G4double localR2 = localpos1.x() * localpo 173 localpos1.y() * localpo 174 localpos1.z() * localpo 175 // G4double InsideRadius2 = 176 // sphereSolid->GetInsideRadius()*sphereS 177 // if(std::fabs( localR2 - InsideRadius2 ) 178 G4double InsideRadius = sphereSolid->GetIn 179 if(localR2 > 180 (InsideRadius - kCarTolerance) * (Ins 181 localR2 < 182 (InsideRadius + kCarTolerance) * (Ins 183 { 184 return fFlux_In; 185 } 186 } 187 188 if(aStep->GetPostStepPoint()->GetStepStatus( 189 { 190 // Exiting Geometry 191 G4ThreeVector stppos2 = aStep->GetPostStep 192 G4ThreeVector localpos2 = 193 theTouchable->GetHistory()->GetTopTransf 194 G4double localR2 = localpos2.x() * localpo 195 localpos2.y() * localpo 196 localpos2.z() * localpo 197 // G4double InsideRadius2 = 198 // sphereSolid->GetInsideRadius()*sphereS 199 // if(std::facb(localR2 - InsideRadius2) ) 200 G4double InsideRadius = sphereSolid->GetIn 201 if(localR2 > 202 (InsideRadius - kCarTolerance) * (Ins 203 localR2 < 204 (InsideRadius + kCarTolerance) * (Ins 205 { 206 return fFlux_Out; 207 } 208 } 209 210 return -1; 211 } 212 213 void G4PSSphereSurfaceFlux::Initialize(G4HCofT 214 { 215 EvtMap = new G4THitsMap<G4double>(detector-> 216 if(HCID < 0) 217 HCID = GetCollectionID(0); 218 HCE->AddHitsCollection(HCID, (G4VHitsCollect 219 } 220 221 void G4PSSphereSurfaceFlux::clear() { EvtMap-> 222 223 void G4PSSphereSurfaceFlux::PrintAll() 224 { 225 G4cout << " MultiFunctionalDet " << detecto 226 G4cout << " PrimitiveScorer " << GetName() < 227 G4cout << " Number of entries " << EvtMap->e 228 for(const auto& [copy, flux] : *(EvtMap->Get 229 { 230 G4cout << " copy no.: " << copy 231 << " Flux : " << *(flux) / GetUni 232 << GetUnit() << "]" << G4endl; 233 } 234 } 235 236 void G4PSSphereSurfaceFlux::SetUnit(const G4St 237 { 238 if(divideByArea) 239 { 240 CheckAndSetUnit(unit, "Per Unit Surface"); 241 } 242 else 243 { 244 if(unit.empty()) 245 { 246 unitName = unit; 247 unitValue = 1.0; 248 } 249 else 250 { 251 G4String msg = "Invalid unit [" + unit + 252 GetUnit() + "] ) for " + 253 G4Exception("G4PSSphereSurfaceFlux::SetU 254 msg); 255 } 256 } 257 } 258 259 void G4PSSphereSurfaceFlux::DefineUnitAndCateg 260 { 261 // Per Unit Surface 262 new G4UnitDefinition("percentimeter2", "perc 263 (1. / cm2)); 264 new G4UnitDefinition("permillimeter2", "perm 265 (1. / mm2)); 266 new G4UnitDefinition("permeter2", "perm2", " 267 } 268