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 // // G4PSCylinderSurfaceFlux 29 #include "G4PSCylinderSurfaceFlux.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 #include "G4VScoreHistFiller.hh" 40 41 // /////////////////////////////////////////// 42 // (Description) 43 // This is a primitive scorer class for scor 44 // Current version assumes only for G4Tubs sh 45 // is fixed on inner plane of the tube. 46 // 47 // Surface is defined at the innner surface of 48 // Direction R R+dR 49 // 0 IN || OUT ->|<- | 50 // 1 IN ->| | 51 // 2 OUT |<- | 52 // 53 // Created: 2007-03-29 Tsukasa ASO 54 // 2010-07-22 Introduce Unit specification. 55 // 2010-07-22 Add weighted and divideByArea 56 // 2011-02-21 Get correct momentum direction 57 // 2020-10-06 Use G4VPrimitivePlotter and fi 58 // vs. surface flux * track weigh 59 // 60 ////////////////////////////////////////////// 61 62 G4PSCylinderSurfaceFlux::G4PSCylinderSurfaceFl 63 64 : G4PSCylinderSurfaceFlux(name, direction, " 65 {} 66 67 G4PSCylinderSurfaceFlux::G4PSCylinderSurfaceFl 68 69 70 : G4VPrimitivePlotter(name, depth) 71 , HCID(-1) 72 , fDirection(direction) 73 , EvtMap(nullptr) 74 , weighted(true) 75 , divideByArea(true) 76 { 77 DefineUnitAndCategory(); 78 SetUnit(unit); 79 } 80 81 G4bool G4PSCylinderSurfaceFlux::ProcessHits(G4 82 { 83 G4StepPoint* preStep = aStep->GetPreStepPoin 84 G4VSolid* solid = ComputeCurrentSolid(a 85 assert(dynamic_cast<G4Tubs*>(solid)); 86 87 auto tubsSolid = static_cast<G4Tubs*>(solid 88 89 G4int dirFlag = IsSelectedSurface(aStep, tub 90 91 if(dirFlag > 0) 92 { 93 if(fDirection == fFlux_InOut || dirFlag == 94 { 95 G4StepPoint* thisStep = nullptr; 96 if(dirFlag == fFlux_In) 97 { 98 thisStep = preStep; 99 } 100 else if(dirFlag == fFlux_Out) 101 { 102 thisStep = aStep->GetPostStepPoint(); 103 } 104 else 105 { 106 return false; 107 } 108 109 G4TouchableHandle theTouchable = thisSte 110 G4ThreeVector pdirection = thisSte 111 G4ThreeVector localdir = 112 theTouchable->GetHistory()->GetTopTran 113 G4ThreeVector position = thisStep->GetPo 114 G4ThreeVector localpos = 115 theTouchable->GetHistory()->GetTopTran 116 G4double angleFactor = 117 (localdir.x() * localpos.x() + localdi 118 std::sqrt(localdir.x() * localdir.x() 119 localdir.z() * localdir.z()) 120 std::sqrt(localpos.x() * localpos.x() 121 122 if(angleFactor < 0) 123 angleFactor *= -1.; 124 G4double square = 2. * tubsSolid->GetZHa 125 tubsSolid->GetInnerRad 126 tubsSolid->GetDeltaPhi 127 128 G4double flux = 1.0; 129 if(weighted) 130 flux *= preStep->GetWeight(); 131 // Current (Particle Weight) 132 133 flux = flux / angleFactor; 134 if(divideByArea) 135 flux /= square; 136 // Flux with angle. 137 G4int index = GetIndex(aStep); 138 EvtMap->add(index, flux); 139 140 if(!hitIDMap.empty() && hitIDMap.find(in 141 { 142 auto filler = G4VScoreHistFiller::Inst 143 if(filler == nullptr) 144 { 145 G4Exception("G4PSCylinderSurfaceFlux 146 JustWarning, 147 "G4TScoreHistFiller is n 148 "not filled."); 149 } 150 else 151 { 152 filler->FillH1(hitIDMap[index], this 153 } 154 } 155 156 return true; 157 } 158 159 return false; 160 } 161 162 return false; 163 } 164 165 G4int G4PSCylinderSurfaceFlux::IsSelectedSurfa 166 167 { 168 G4TouchableHandle theTouchable = 169 aStep->GetPreStepPoint()->GetTouchableHand 170 G4double kCarTolerance = 171 G4GeometryTolerance::GetInstance()->GetSur 172 173 if(aStep->GetPreStepPoint()->GetStepStatus() 174 { 175 // Entering Geometry 176 G4ThreeVector stppos1 = aStep->GetPreStepP 177 G4ThreeVector localpos1 = 178 theTouchable->GetHistory()->GetTopTransf 179 if(std::fabs(localpos1.z()) > tubsSolid->G 180 return -1; 181 // if(std::fabs( localpos1.x()*localpos1.x 182 // - (tubsSolid->GetInnerRadius()*tubsSo 183 // <kCarTolerance ){ 184 G4double localR2 = 185 localpos1.x() * localpos1.x() + localpos 186 G4double InsideRadius = tubsSolid->GetInne 187 if(localR2 > 188 (InsideRadius - kCarTolerance) * (Ins 189 localR2 < 190 (InsideRadius + kCarTolerance) * (Ins 191 { 192 return fFlux_In; 193 } 194 } 195 196 if(aStep->GetPostStepPoint()->GetStepStatus( 197 { 198 // Exiting Geometry 199 G4ThreeVector stppos2 = aStep->GetPostStep 200 G4ThreeVector localpos2 = 201 theTouchable->GetHistory()->GetTopTransf 202 if(std::fabs(localpos2.z()) > tubsSolid->G 203 return -1; 204 // if(std::fabs( localpos2.x()*localpos2.x 205 // - (tubsSolid->GetInnerRadius()*tubsSo 206 // <kCarTolerance ){ 207 G4double localR2 = 208 localpos2.x() * localpos2.x() + localpos 209 G4double InsideRadius = tubsSolid->GetInne 210 if(localR2 > 211 (InsideRadius - kCarTolerance) * (Ins 212 localR2 < 213 (InsideRadius + kCarTolerance) * (Ins 214 { 215 return fFlux_Out; 216 } 217 } 218 219 return -1; 220 } 221 222 void G4PSCylinderSurfaceFlux::Initialize(G4HCo 223 { 224 EvtMap = new G4THitsMap<G4double>(GetMultiFu 225 GetName()) 226 if(HCID < 0) 227 HCID = GetCollectionID(0); 228 HCE->AddHitsCollection(HCID, (G4VHitsCollect 229 } 230 231 void G4PSCylinderSurfaceFlux::clear() { EvtMap 232 233 void G4PSCylinderSurfaceFlux::PrintAll() 234 { 235 G4cout << " MultiFunctionalDet " << detecto 236 G4cout << " PrimitiveScorer" << GetName() << 237 G4cout << " Number of entries " << EvtMap->e 238 for(const auto& [copy, flux] : *(EvtMap->Get 239 { 240 G4cout << " copy no.: " << copy 241 << " flux : " << *(flux) / GetUni 242 << GetUnit() << "]" << G4endl; 243 } 244 } 245 246 void G4PSCylinderSurfaceFlux::SetUnit(const G4 247 { 248 if(divideByArea) 249 { 250 CheckAndSetUnit(unit, "Per Unit Surface"); 251 } 252 else 253 { 254 if(unit.empty()) 255 { 256 unitName = unit; 257 unitValue = 1.0; 258 } 259 else 260 { 261 G4String msg = "Invalid unit [" + unit + 262 GetUnit() + "] ) for " + 263 G4Exception("G4PSCylinderSurfaceFlux::Se 264 msg); 265 } 266 } 267 } 268 269 void G4PSCylinderSurfaceFlux::DefineUnitAndCat 270 { 271 // Per Unit Surface 272 new G4UnitDefinition("percentimeter2", "perc 273 (1. / cm2)); 274 new G4UnitDefinition("permillimeter2", "perm 275 (1. / mm2)); 276 new G4UnitDefinition("permeter2", "perm2", " 277 } 278