| 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 // G4PSCylinderSurfaceCurrent
29 #include "G4PSCylinderSurfaceCurrent.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 scoring only Surface Current.
44 // Current version assumes only for G4Tubs shape.
45 //
46 // Surface is defined at the inner surface of the tube.
47 // Direction R R+dR
48 // 0 IN || OUT ->|<- |
49 // 1 IN ->| |
50 // 2 OUT |<- |
51 //
52 // Created: 2007-03-21 Tsukasa ASO
53 // 2010-07-22 Introduce Unit specification.
54 // 2020-10-06 Use G4VPrimitivePlotter and fill 1-D histo of kinetic energy (x)
55 // vs. surface current * track weight (y) (Makoto Asai)
56 //
57 ///////////////////////////////////////////////////////////////////////////////
58
59 G4PSCylinderSurfaceCurrent::G4PSCylinderSurfaceCurrent(const G4String& name,
60 G4int direction,
61 G4int depth)
62 : G4PSCylinderSurfaceCurrent(name, direction, "percm2", depth)
63 {}
64
65 G4PSCylinderSurfaceCurrent::G4PSCylinderSurfaceCurrent(const G4String& name,
66 G4int direction,
67 const G4String& unit,
68 G4int depth)
69 : G4VPrimitivePlotter(name, depth)
70 , HCID(-1)
71 , fDirection(direction)
72 , EvtMap(nullptr)
73 , weighted(true)
74 , divideByArea(true)
75 {
76 DefineUnitAndCategory();
77 SetUnit(unit);
78 }
79
80 G4bool G4PSCylinderSurfaceCurrent::ProcessHits(G4Step* aStep,
81 G4TouchableHistory*)
82 {
83 G4StepPoint* preStep = aStep->GetPreStepPoint();
84 G4VSolid* solid = ComputeCurrentSolid(aStep);
85 auto tubsSolid = static_cast<G4Tubs*>(solid);
86
87 G4int dirFlag = IsSelectedSurface(aStep, tubsSolid);
88 // G4cout << " pos " << preStep->GetPosition() <<" dirFlag " << G4endl;
89 if(dirFlag > 0)
90 {
91 if(fDirection == fCurrent_InOut || fDirection == dirFlag)
92 {
93 G4TouchableHandle theTouchable = preStep->GetTouchableHandle();
94 //
95 G4double current = 1.0;
96 if(weighted)
97 current = preStep->GetWeight(); // Current (Particle Weight)
98 //
99 if(divideByArea)
100 {
101 G4double square = 2. * tubsSolid->GetZHalfLength() *
102 tubsSolid->GetInnerRadius() *
103 tubsSolid->GetDeltaPhiAngle() / radian;
104 current = current / square; // Current normalized by Area
105 }
106
107 G4int index = GetIndex(aStep);
108 EvtMap->add(index, current);
109
110 if(!hitIDMap.empty() && hitIDMap.find(index) != hitIDMap.cend())
111 {
112 auto filler = G4VScoreHistFiller::Instance();
113 if(filler == nullptr)
114 {
115 G4Exception("G4PSCylinderSurfaceCurrent::ProcessHits", "SCORER0123",
116 JustWarning,
117 "G4TScoreHistFiller is not instantiated!! Histogram is "
118 "not filled.");
119 }
120 else
121 {
122 filler->FillH1(hitIDMap[index], preStep->GetKineticEnergy(), current);
123 }
124 }
125 }
126 }
127
128 return true;
129 }
130
131 G4int G4PSCylinderSurfaceCurrent::IsSelectedSurface(G4Step* aStep,
132 G4Tubs* tubsSolid)
133 {
134 G4TouchableHandle theTouchable =
135 aStep->GetPreStepPoint()->GetTouchableHandle();
136 G4double kCarTolerance =
137 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
138
139 if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
140 {
141 // Entering Geometry
142 G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
143 G4ThreeVector localpos1 =
144 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
145 if(std::fabs(localpos1.z()) > tubsSolid->GetZHalfLength())
146 return -1;
147 G4double localR2 =
148 localpos1.x() * localpos1.x() + localpos1.y() * localpos1.y();
149 G4double InsideRadius = tubsSolid->GetInnerRadius();
150 if(localR2 >
151 (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
152 localR2 <
153 (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
154 {
155 return fCurrent_In;
156 }
157 }
158
159 if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
160 {
161 // Exiting Geometry
162 G4ThreeVector stppos2 = aStep->GetPostStepPoint()->GetPosition();
163 G4ThreeVector localpos2 =
164 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
165 if(std::fabs(localpos2.z()) > tubsSolid->GetZHalfLength())
166 return -1;
167 G4double localR2 =
168 localpos2.x() * localpos2.x() + localpos2.y() * localpos2.y();
169 G4double InsideRadius = tubsSolid->GetInnerRadius();
170 if(localR2 >
171 (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
172 localR2 <
173 (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
174 {
175 return fCurrent_Out;
176 }
177 }
178
179 return -1;
180 }
181
182 void G4PSCylinderSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
183 {
184 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
185 if(HCID < 0)
186 HCID = GetCollectionID(0);
187 HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap);
188 }
189
190 void G4PSCylinderSurfaceCurrent::clear() { EvtMap->clear(); }
191
192 void G4PSCylinderSurfaceCurrent::PrintAll()
193 {
194 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
195 G4cout << " PrimitiveScorer " << GetName() << G4endl;
196 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
197 for(const auto& [copy, current] : *(EvtMap->GetMap()))
198 {
199 G4cout << " copy no.: " << copy << " current : ";
200 if(divideByArea)
201 {
202 G4cout << *(current) / GetUnitValue() << " [" << GetUnit() << "]";
203 }
204 else
205 {
206 G4cout << *(current) << " [tracks]";
207 }
208 G4cout << G4endl;
209 }
210 }
211
212 void G4PSCylinderSurfaceCurrent::SetUnit(const G4String& unit)
213 {
214 if(divideByArea)
215 {
216 CheckAndSetUnit(unit, "Per Unit Surface");
217 }
218 else
219 {
220 if(unit.empty())
221 {
222 unitName = unit;
223 unitValue = 1.0;
224 }
225 else
226 {
227 G4String msg = "Invalid unit [" + unit + "] (Current unit is [" +
228 GetUnit() + "] ) for " + GetName();
229 G4Exception("G4PSCylinderSurfaceCurrent::SetUnit", "DetPS0002",
230 JustWarning, msg);
231 }
232 }
233 }
234
235 void G4PSCylinderSurfaceCurrent::DefineUnitAndCategory()
236 {
237 // Per Unit Surface
238 new G4UnitDefinition("percentimeter2", "percm2", "Per Unit Surface",
239 (1. / cm2));
240 new G4UnitDefinition("permillimeter2", "permm2", "Per Unit Surface",
241 (1. / mm2));
242 new G4UnitDefinition("permeter2", "perm2", "Per Unit Surface", (1. / m2));
243 }
244