| 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 // G4PSSphereSurfaceCurrent
29 #include "G4PSSphereSurfaceCurrent.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 scoring only Surface Current.
42 // Current version assumes only for G4Sphere shape.
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 Kimura.
51 // 2010-07-22 Introduce Unit specification.
52 //
53 ///////////////////////////////////////////////////////////////////////////////
54
55 G4PSSphereSurfaceCurrent::G4PSSphereSurfaceCurrent(const G4String& name,
56 G4int direction, G4int depth)
57 : G4PSSphereSurfaceCurrent(name, direction, "percm2", depth)
58 {}
59
60 G4PSSphereSurfaceCurrent::G4PSSphereSurfaceCurrent(const G4String& name,
61 G4int direction,
62 const G4String& unit,
63 G4int depth)
64 : G4VPrimitiveScorer(name, depth)
65 , HCID(-1)
66 , fDirection(direction)
67 , EvtMap(nullptr)
68 , weighted(true)
69 , divideByArea(true)
70 {
71 DefineUnitAndCategory();
72 SetUnit(unit);
73 }
74
75 G4bool G4PSSphereSurfaceCurrent::ProcessHits(G4Step* aStep, G4TouchableHistory*)
76 {
77 G4StepPoint* preStep = aStep->GetPreStepPoint();
78 G4VSolid* solid = ComputeCurrentSolid(aStep);
79 assert(dynamic_cast<G4Sphere*>(solid) != nullptr);
80
81 auto sphereSolid = static_cast<G4Sphere*>(solid);
82
83 G4int dirFlag = IsSelectedSurface(aStep, sphereSolid);
84 if(dirFlag > 0)
85 {
86 if(fDirection == fCurrent_InOut || fDirection == dirFlag)
87 {
88 G4double radi = sphereSolid->GetInnerRadius();
89 G4double dph = sphereSolid->GetDeltaPhiAngle() / radian;
90 G4double stth = sphereSolid->GetStartThetaAngle() / radian;
91 G4double enth = stth + sphereSolid->GetDeltaThetaAngle() / radian;
92 G4double current = 1.0;
93 if(weighted)
94 current = preStep->GetWeight(); // Current (Particle Weight)
95 if(divideByArea)
96 {
97 G4double square =
98 radi * radi * dph * (-std::cos(enth) + std::cos(stth));
99 current = current / square; // Current with angle.
100 }
101
102 G4int index = GetIndex(aStep);
103 EvtMap->add(index, current);
104 }
105 }
106
107 return true;
108 }
109
110 G4int G4PSSphereSurfaceCurrent::IsSelectedSurface(G4Step* aStep,
111 G4Sphere* sphereSolid)
112 {
113 G4TouchableHandle theTouchable =
114 aStep->GetPreStepPoint()->GetTouchableHandle();
115 G4double kCarTolerance =
116 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
117
118 if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
119 {
120 // Entering Geometry
121 G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
122 G4ThreeVector localpos1 =
123 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
124 G4double localR2 = localpos1.x() * localpos1.x() +
125 localpos1.y() * localpos1.y() +
126 localpos1.z() * localpos1.z();
127 // G4double InsideRadius2 =
128 // sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
129 // if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
130 G4double InsideRadius = sphereSolid->GetInnerRadius();
131 if(localR2 >
132 (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
133 localR2 <
134 (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
135 {
136 return fCurrent_In;
137 }
138 }
139
140 if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
141 {
142 // Exiting Geometry
143 G4ThreeVector stppos2 = aStep->GetPostStepPoint()->GetPosition();
144 G4ThreeVector localpos2 =
145 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
146 G4double localR2 = localpos2.x() * localpos2.x() +
147 localpos2.y() * localpos2.y() +
148 localpos2.z() * localpos2.z();
149 // G4double InsideRadius2 =
150 // sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
151 // if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
152 G4double InsideRadius = sphereSolid->GetInnerRadius();
153 if(localR2 >
154 (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
155 localR2 <
156 (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
157 {
158 return fCurrent_Out;
159 }
160 }
161
162 return -1;
163 }
164
165 void G4PSSphereSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
166 {
167 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
168 if(HCID < 0)
169 HCID = GetCollectionID(0);
170 HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap);
171 }
172
173 void G4PSSphereSurfaceCurrent::clear() { EvtMap->clear(); }
174
175 void G4PSSphereSurfaceCurrent::PrintAll()
176 {
177 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
178 G4cout << " PrimitiveScorer " << GetName() << G4endl;
179 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
180 for(const auto& [copy, current] : *(EvtMap->GetMap()))
181 {
182 G4cout << " copy no.: " << copy << " current : ";
183 if(divideByArea)
184 {
185 G4cout << *(current) / GetUnitValue() << " [" << GetUnit() << "]";
186 }
187 else
188 {
189 G4cout << *(current) << " [tracks]";
190 }
191 G4cout << G4endl;
192 }
193 }
194
195 void G4PSSphereSurfaceCurrent::SetUnit(const G4String& unit)
196 {
197 if(divideByArea)
198 {
199 CheckAndSetUnit(unit, "Per Unit Surface");
200 }
201 else
202 {
203 if(unit.empty())
204 {
205 unitName = unit;
206 unitValue = 1.0;
207 }
208 else
209 {
210 G4String msg = "Invalid unit [" + unit + "] (Current unit is [" +
211 GetUnit() + "] ) for " + GetName();
212 G4Exception("G4PSSphereSurfaceCurrent::SetUnit", "DetPS0015", JustWarning,
213 msg);
214 }
215 }
216 }
217
218 void G4PSSphereSurfaceCurrent::DefineUnitAndCategory()
219 {
220 // Per Unit Surface
221 new G4UnitDefinition("percentimeter2", "percm2", "Per Unit Surface",
222 (1. / cm2));
223 new G4UnitDefinition("permillimeter2", "permm2", "Per Unit Surface",
224 (1. / mm2));
225 new G4UnitDefinition("permeter2", "perm2", "Per Unit Surface", (1. / m2));
226 }
227