| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // G4PSFlatSurfaceCurrent 28 // G4PSFlatSurfaceCurrent
29 #include "G4PSFlatSurfaceCurrent.hh" 29 #include "G4PSFlatSurfaceCurrent.hh"
30 30
31 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
32 #include "G4StepStatus.hh" 32 #include "G4StepStatus.hh"
33 #include "G4Track.hh" 33 #include "G4Track.hh"
34 #include "G4VSolid.hh" 34 #include "G4VSolid.hh"
35 #include "G4VPhysicalVolume.hh" 35 #include "G4VPhysicalVolume.hh"
36 #include "G4VPVParameterisation.hh" 36 #include "G4VPVParameterisation.hh"
37 #include "G4UnitsTable.hh" 37 #include "G4UnitsTable.hh"
38 #include "G4GeometryTolerance.hh" 38 #include "G4GeometryTolerance.hh"
39 #include "G4VScoreHistFiller.hh" 39 #include "G4VScoreHistFiller.hh"
40 40
41 ////////////////////////////////////////////// 41 ////////////////////////////////////////////////////////////////////////////////
42 // (Description) 42 // (Description)
43 // This is a primitive scorer class for scor 43 // This is a primitive scorer class for scoring only Surface Flux.
44 // Current version assumes only for G4Box sha 44 // Current version assumes only for G4Box shape.
45 // 45 //
46 // Surface is defined at the -Z surface. 46 // Surface is defined at the -Z surface.
47 // Direction -Z +Z 47 // Direction -Z +Z
48 // 0 IN || OUT ->|<- | 48 // 0 IN || OUT ->|<- |
49 // 1 IN ->| | 49 // 1 IN ->| |
50 // 2 OUT |<- | 50 // 2 OUT |<- |
51 // 51 //
52 // Created: 2005-11-14 Tsukasa ASO, Akinori K 52 // Created: 2005-11-14 Tsukasa ASO, Akinori Kimura.
53 // 17-Nov-2005 T.Aso, Bug fix for area definit 53 // 17-Nov-2005 T.Aso, Bug fix for area definition.
54 // 31-Mar-2007 T.Aso, Add option for normalizi 54 // 31-Mar-2007 T.Aso, Add option for normalizing by the area.
55 // 2010-07-22 Introduce Unit specification. 55 // 2010-07-22 Introduce Unit specification.
56 // 2020-10-06 Use G4VPrimitivePlotter and fi 56 // 2020-10-06 Use G4VPrimitivePlotter and fill 1-D histo of kinetic energy (x)
57 // vs. Surface Current * track w 57 // vs. Surface Current * track weight (y) (Makoto Asai)
58 // 58 //
59 ////////////////////////////////////////////// 59 ///////////////////////////////////////////////////////////////////////////////
60 60
61 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent << 61 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent(G4String name, G4int direction,
62 62 G4int depth)
63 : G4PSFlatSurfaceCurrent(name, direction, "p << 63 : G4VPrimitivePlotter(name, depth)
64 {} << 64 , HCID(-1)
>> 65 , fDirection(direction)
>> 66 , EvtMap(0)
>> 67 , weighted(true)
>> 68 , divideByArea(true)
>> 69 {
>> 70 DefineUnitAndCategory();
>> 71 SetUnit("percm2");
>> 72 }
65 73
66 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent << 74 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent(G4String name, G4int direction,
67 75 const G4String& unit,
68 76 G4int depth)
69 : G4VPrimitivePlotter(name, depth) 77 : G4VPrimitivePlotter(name, depth)
70 , HCID(-1) 78 , HCID(-1)
71 , fDirection(direction) 79 , fDirection(direction)
72 , EvtMap(nullptr) << 80 , EvtMap(0)
73 , weighted(true) 81 , weighted(true)
74 , divideByArea(true) 82 , divideByArea(true)
75 { 83 {
76 DefineUnitAndCategory(); 84 DefineUnitAndCategory();
77 SetUnit(unit); 85 SetUnit(unit);
78 } 86 }
79 87
>> 88 G4PSFlatSurfaceCurrent::~G4PSFlatSurfaceCurrent() { ; }
>> 89
80 G4bool G4PSFlatSurfaceCurrent::ProcessHits(G4S 90 G4bool G4PSFlatSurfaceCurrent::ProcessHits(G4Step* aStep, G4TouchableHistory*)
81 { 91 {
82 G4StepPoint* preStep = aStep->Ge 92 G4StepPoint* preStep = aStep->GetPreStepPoint();
83 G4VPhysicalVolume* physVol = preStep-> 93 G4VPhysicalVolume* physVol = preStep->GetPhysicalVolume();
84 G4VPVParameterisation* physParam = physVol-> 94 G4VPVParameterisation* physParam = physVol->GetParameterisation();
85 G4VSolid* solid = nullptr; << 95 G4VSolid* solid = 0;
86 if(physParam != nullptr) << 96 if(physParam)
87 { // for parameterized volume 97 { // for parameterized volume
88 G4int idx = 98 G4int idx =
89 ((G4TouchableHistory*) (aStep->GetPreSte 99 ((G4TouchableHistory*) (aStep->GetPreStepPoint()->GetTouchable()))
90 ->GetReplicaNumber(indexDepth); 100 ->GetReplicaNumber(indexDepth);
91 solid = physParam->ComputeSolid(idx, physV 101 solid = physParam->ComputeSolid(idx, physVol);
92 solid->ComputeDimensions(physParam, idx, p 102 solid->ComputeDimensions(physParam, idx, physVol);
93 } 103 }
94 else 104 else
95 { // for ordinary volume 105 { // for ordinary volume
96 solid = physVol->GetLogicalVolume()->GetSo 106 solid = physVol->GetLogicalVolume()->GetSolid();
97 } 107 }
98 108
99 auto boxSolid = (G4Box*) (solid); << 109 G4Box* boxSolid = (G4Box*) (solid);
100 110
101 G4int dirFlag = IsSelectedSurface(aStep, box 111 G4int dirFlag = IsSelectedSurface(aStep, boxSolid);
102 if(dirFlag > 0) 112 if(dirFlag > 0)
103 { 113 {
104 if(fDirection == fCurrent_InOut || fDirect 114 if(fDirection == fCurrent_InOut || fDirection == dirFlag)
105 { 115 {
106 G4int index = GetInde 116 G4int index = GetIndex(aStep);
107 G4TouchableHandle theTouchable = preStep 117 G4TouchableHandle theTouchable = preStep->GetTouchableHandle();
108 G4double current = 1.0; 118 G4double current = 1.0;
109 if(weighted) 119 if(weighted)
110 current = preStep->GetWeight(); // Cu 120 current = preStep->GetWeight(); // Current (Particle Weight)
111 if(divideByArea) 121 if(divideByArea)
112 { 122 {
113 G4double square = 123 G4double square =
114 4. * boxSolid->GetXHalfLength() * bo 124 4. * boxSolid->GetXHalfLength() * boxSolid->GetYHalfLength();
115 current = current / square; // Normal 125 current = current / square; // Normalized by Area
116 } 126 }
117 EvtMap->add(index, current); 127 EvtMap->add(index, current);
118 128
119 if(!hitIDMap.empty() && hitIDMap.find(in << 129 if(hitIDMap.size() > 0 && hitIDMap.find(index) != hitIDMap.end())
120 { 130 {
121 auto filler = G4VScoreHistFiller::Inst 131 auto filler = G4VScoreHistFiller::Instance();
122 if(filler == nullptr) << 132 if(!filler)
123 { 133 {
124 G4Exception("G4PSFlatSurfaceCurrent: 134 G4Exception("G4PSFlatSurfaceCurrent::ProcessHits", "SCORER0123",
125 JustWarning, 135 JustWarning,
126 "G4TScoreHistFiller is n 136 "G4TScoreHistFiller is not instantiated!! Histogram is "
127 "not filled."); 137 "not filled.");
128 } 138 }
129 else 139 else
130 { 140 {
131 filler->FillH1(hitIDMap[index], preS 141 filler->FillH1(hitIDMap[index], preStep->GetKineticEnergy(), current);
132 } 142 }
133 } 143 }
134 } 144 }
135 } 145 }
136 146
137 return true; << 147 return TRUE;
138 } 148 }
139 149
140 G4int G4PSFlatSurfaceCurrent::IsSelectedSurfac 150 G4int G4PSFlatSurfaceCurrent::IsSelectedSurface(G4Step* aStep, G4Box* boxSolid)
141 { 151 {
142 G4TouchableHandle theTouchable = 152 G4TouchableHandle theTouchable =
143 aStep->GetPreStepPoint()->GetTouchableHand 153 aStep->GetPreStepPoint()->GetTouchableHandle();
144 G4double kCarTolerance = 154 G4double kCarTolerance =
145 G4GeometryTolerance::GetInstance()->GetSur 155 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
146 156
147 if(aStep->GetPreStepPoint()->GetStepStatus() 157 if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
148 { 158 {
149 // Entering Geometry 159 // Entering Geometry
150 G4ThreeVector stppos1 = aStep->GetPreStepP 160 G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
151 G4ThreeVector localpos1 = 161 G4ThreeVector localpos1 =
152 theTouchable->GetHistory()->GetTopTransf 162 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
153 if(std::fabs(localpos1.z() + boxSolid->Get 163 if(std::fabs(localpos1.z() + boxSolid->GetZHalfLength()) < kCarTolerance)
154 { 164 {
155 return fCurrent_In; 165 return fCurrent_In;
156 } 166 }
157 } 167 }
158 168
159 if(aStep->GetPostStepPoint()->GetStepStatus( 169 if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
160 { 170 {
161 // Exiting Geometry 171 // Exiting Geometry
162 G4ThreeVector stppos2 = aStep->GetPostStep 172 G4ThreeVector stppos2 = aStep->GetPostStepPoint()->GetPosition();
163 G4ThreeVector localpos2 = 173 G4ThreeVector localpos2 =
164 theTouchable->GetHistory()->GetTopTransf 174 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
165 if(std::fabs(localpos2.z() + boxSolid->Get 175 if(std::fabs(localpos2.z() + boxSolid->GetZHalfLength()) < kCarTolerance)
166 { 176 {
167 return fCurrent_Out; 177 return fCurrent_Out;
168 } 178 }
169 } 179 }
170 180
171 return -1; 181 return -1;
172 } 182 }
173 183
174 void G4PSFlatSurfaceCurrent::Initialize(G4HCof 184 void G4PSFlatSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
175 { 185 {
176 EvtMap = new G4THitsMap<G4double>(detector-> 186 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
177 if(HCID < 0) 187 if(HCID < 0)
178 HCID = GetCollectionID(0); 188 HCID = GetCollectionID(0);
179 HCE->AddHitsCollection(HCID, (G4VHitsCollect 189 HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap);
180 } 190 }
181 191
>> 192 void G4PSFlatSurfaceCurrent::EndOfEvent(G4HCofThisEvent*) { ; }
>> 193
182 void G4PSFlatSurfaceCurrent::clear() { EvtMap- 194 void G4PSFlatSurfaceCurrent::clear() { EvtMap->clear(); }
183 195
>> 196 void G4PSFlatSurfaceCurrent::DrawAll() { ; }
>> 197
184 void G4PSFlatSurfaceCurrent::PrintAll() 198 void G4PSFlatSurfaceCurrent::PrintAll()
185 { 199 {
186 G4cout << " MultiFunctionalDet " << detecto 200 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
187 G4cout << " PrimitiveScorer " << GetName() < 201 G4cout << " PrimitiveScorer " << GetName() << G4endl;
188 G4cout << " Number of entries " << EvtMap->e 202 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
189 for(const auto& [copy, current] : *(EvtMap-> << 203 std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin();
>> 204 for(; itr != EvtMap->GetMap()->end(); itr++)
190 { 205 {
191 G4cout << " copy no.: " << copy << " curr << 206 G4cout << " copy no.: " << itr->first << " current : ";
192 if(divideByArea) 207 if(divideByArea)
193 { 208 {
194 G4cout << *(current) / GetUnitValue() << << 209 G4cout << *(itr->second) / GetUnitValue() << " [" << GetUnit() << "]";
195 } 210 }
196 else 211 else
197 { 212 {
198 G4cout << *(current) / GetUnitValue() << << 213 G4cout << *(itr->second) / GetUnitValue() << " [tracks]";
199 } 214 }
200 G4cout << G4endl; 215 G4cout << G4endl;
201 } 216 }
202 } 217 }
203 218
204 void G4PSFlatSurfaceCurrent::SetUnit(const G4S 219 void G4PSFlatSurfaceCurrent::SetUnit(const G4String& unit)
205 { 220 {
206 if(divideByArea) 221 if(divideByArea)
207 { 222 {
208 CheckAndSetUnit(unit, "Per Unit Surface"); 223 CheckAndSetUnit(unit, "Per Unit Surface");
209 } 224 }
210 else 225 else
211 { 226 {
212 if(unit.empty()) << 227 if(unit == "")
213 { 228 {
214 unitName = unit; 229 unitName = unit;
215 unitValue = 1.0; 230 unitValue = 1.0;
216 } 231 }
217 else 232 else
218 { 233 {
219 G4String msg = "Invalid unit [" + unit + 234 G4String msg = "Invalid unit [" + unit + "] (Current unit is [" +
220 GetUnit() + "] ) for " + 235 GetUnit() + "] ) for " + GetName();
221 G4Exception("G4PSFlatSurfaceCurrent::Set 236 G4Exception("G4PSFlatSurfaceCurrent::SetUnit", "DetPS0007", JustWarning,
222 msg); 237 msg);
223 } 238 }
224 } 239 }
225 } 240 }
226 241
227 void G4PSFlatSurfaceCurrent::DefineUnitAndCate 242 void G4PSFlatSurfaceCurrent::DefineUnitAndCategory()
228 { 243 {
229 // Per Unit Surface 244 // Per Unit Surface
230 new G4UnitDefinition("percentimeter2", "perc 245 new G4UnitDefinition("percentimeter2", "percm2", "Per Unit Surface",
231 (1. / cm2)); 246 (1. / cm2));
232 new G4UnitDefinition("permillimeter2", "perm 247 new G4UnitDefinition("permillimeter2", "permm2", "Per Unit Surface",
233 (1. / mm2)); 248 (1. / mm2));
234 new G4UnitDefinition("permeter2", "perm2", " 249 new G4UnitDefinition("permeter2", "perm2", "Per Unit Surface", (1. / m2));
235 } 250 }
236 251