| 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 // $Id: G4PSFlatSurfaceCurrent.cc 81087 2014-05-20 15:44:27Z gcosmo $
27 // 28 //
28 // G4PSFlatSurfaceCurrent 29 // G4PSFlatSurfaceCurrent
29 #include "G4PSFlatSurfaceCurrent.hh" 30 #include "G4PSFlatSurfaceCurrent.hh"
30 31
31 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
32 #include "G4StepStatus.hh" 33 #include "G4StepStatus.hh"
33 #include "G4Track.hh" 34 #include "G4Track.hh"
34 #include "G4VSolid.hh" 35 #include "G4VSolid.hh"
35 #include "G4VPhysicalVolume.hh" 36 #include "G4VPhysicalVolume.hh"
36 #include "G4VPVParameterisation.hh" 37 #include "G4VPVParameterisation.hh"
37 #include "G4UnitsTable.hh" 38 #include "G4UnitsTable.hh"
38 #include "G4GeometryTolerance.hh" 39 #include "G4GeometryTolerance.hh"
39 #include "G4VScoreHistFiller.hh" <<
40 <<
41 ////////////////////////////////////////////// 40 ////////////////////////////////////////////////////////////////////////////////
42 // (Description) 41 // (Description)
43 // This is a primitive scorer class for scor 42 // This is a primitive scorer class for scoring only Surface Flux.
44 // Current version assumes only for G4Box sha << 43 // Current version assumes only for G4Box shape.
45 // 44 //
46 // Surface is defined at the -Z surface. 45 // Surface is defined at the -Z surface.
47 // Direction -Z +Z 46 // Direction -Z +Z
48 // 0 IN || OUT ->|<- | 47 // 0 IN || OUT ->|<- |
49 // 1 IN ->| | 48 // 1 IN ->| |
50 // 2 OUT |<- | 49 // 2 OUT |<- |
51 // 50 //
52 // Created: 2005-11-14 Tsukasa ASO, Akinori K 51 // Created: 2005-11-14 Tsukasa ASO, Akinori Kimura.
53 // 17-Nov-2005 T.Aso, Bug fix for area definit 52 // 17-Nov-2005 T.Aso, Bug fix for area definition.
54 // 31-Mar-2007 T.Aso, Add option for normalizi 53 // 31-Mar-2007 T.Aso, Add option for normalizing by the area.
55 // 2010-07-22 Introduce Unit specification. 54 // 2010-07-22 Introduce Unit specification.
56 // 2020-10-06 Use G4VPrimitivePlotter and fi << 55 //
57 // vs. Surface Current * track w <<
58 // <<
59 ////////////////////////////////////////////// 56 ///////////////////////////////////////////////////////////////////////////////
60 57
61 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent << 58
62 << 59 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent(G4String name,
63 : G4PSFlatSurfaceCurrent(name, direction, "p << 60 G4int direction, G4int depth)
64 {} << 61 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),EvtMap(0),
65 << 62 weighted(true),divideByArea(true)
66 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent <<
67 <<
68 <<
69 : G4VPrimitivePlotter(name, depth) <<
70 , HCID(-1) <<
71 , fDirection(direction) <<
72 , EvtMap(nullptr) <<
73 , weighted(true) <<
74 , divideByArea(true) <<
75 { 63 {
76 DefineUnitAndCategory(); << 64 DefineUnitAndCategory();
77 SetUnit(unit); << 65 SetUnit("percm2");
78 } 66 }
79 67
80 G4bool G4PSFlatSurfaceCurrent::ProcessHits(G4S << 68 G4PSFlatSurfaceCurrent::G4PSFlatSurfaceCurrent(G4String name,
>> 69 G4int direction,
>> 70 const G4String& unit,
>> 71 G4int depth)
>> 72 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),EvtMap(0),
>> 73 weighted(true),divideByArea(true)
>> 74 {
>> 75 DefineUnitAndCategory();
>> 76 SetUnit(unit);
>> 77 }
>> 78
>> 79 G4PSFlatSurfaceCurrent::~G4PSFlatSurfaceCurrent()
>> 80 {;}
>> 81
>> 82 G4bool G4PSFlatSurfaceCurrent::ProcessHits(G4Step* aStep,G4TouchableHistory*)
81 { 83 {
82 G4StepPoint* preStep = aStep->Ge << 84 G4StepPoint* preStep = aStep->GetPreStepPoint();
83 G4VPhysicalVolume* physVol = preStep-> << 85 G4VPhysicalVolume* physVol = preStep->GetPhysicalVolume();
84 G4VPVParameterisation* physParam = physVol-> 86 G4VPVParameterisation* physParam = physVol->GetParameterisation();
85 G4VSolid* solid = nullptr; << 87 G4VSolid * solid = 0;
86 if(physParam != nullptr) << 88 if(physParam)
87 { // for parameterized volume << 89 { // for parameterized volume
88 G4int idx = << 90 G4int idx = ((G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable()))
89 ((G4TouchableHistory*) (aStep->GetPreSte << 91 ->GetReplicaNumber(indexDepth);
90 ->GetReplicaNumber(indexDepth); <<
91 solid = physParam->ComputeSolid(idx, physV 92 solid = physParam->ComputeSolid(idx, physVol);
92 solid->ComputeDimensions(physParam, idx, p << 93 solid->ComputeDimensions(physParam,idx,physVol);
93 } 94 }
94 else 95 else
95 { // for ordinary volume << 96 { // for ordinary volume
96 solid = physVol->GetLogicalVolume()->GetSo 97 solid = physVol->GetLogicalVolume()->GetSolid();
97 } 98 }
98 99
99 auto boxSolid = (G4Box*) (solid); << 100 G4Box* boxSolid = (G4Box*)(solid);
100 101
101 G4int dirFlag = IsSelectedSurface(aStep, box << 102 G4int dirFlag =IsSelectedSurface(aStep,boxSolid);
102 if(dirFlag > 0) << 103 if ( dirFlag > 0 ) {
103 { << 104 if ( fDirection == fCurrent_InOut || fDirection == dirFlag ){
104 if(fDirection == fCurrent_InOut || fDirect << 105 G4int index = GetIndex(aStep);
105 { << 106 G4TouchableHandle theTouchable = preStep->GetTouchableHandle();
106 G4int index = GetInde << 107 G4double current = 1.0;
107 G4TouchableHandle theTouchable = preStep << 108 if ( weighted ) current=preStep->GetWeight(); // Current (Particle Weight)
108 G4double current = 1.0; << 109 if ( divideByArea ){
109 if(weighted) << 110 G4double square = 4.*boxSolid->GetXHalfLength()*boxSolid->GetYHalfLength();
110 current = preStep->GetWeight(); // Cu << 111 current = current/square; // Normalized by Area
111 if(divideByArea) << 112 }
112 { << 113 EvtMap->add(index,current);
113 G4double square = << 114 }
114 4. * boxSolid->GetXHalfLength() * bo <<
115 current = current / square; // Normal <<
116 } <<
117 EvtMap->add(index, current); <<
118 <<
119 if(!hitIDMap.empty() && hitIDMap.find(in <<
120 { <<
121 auto filler = G4VScoreHistFiller::Inst <<
122 if(filler == nullptr) <<
123 { <<
124 G4Exception("G4PSFlatSurfaceCurrent: <<
125 JustWarning, <<
126 "G4TScoreHistFiller is n <<
127 "not filled."); <<
128 } <<
129 else <<
130 { <<
131 filler->FillH1(hitIDMap[index], preS <<
132 } <<
133 } <<
134 } <<
135 } 115 }
136 116
137 return true; << 117 return TRUE;
138 } 118 }
139 119
140 G4int G4PSFlatSurfaceCurrent::IsSelectedSurfac << 120 G4int G4PSFlatSurfaceCurrent::IsSelectedSurface(G4Step* aStep, G4Box* boxSolid){
141 { << 121
142 G4TouchableHandle theTouchable = << 122 G4TouchableHandle theTouchable =
143 aStep->GetPreStepPoint()->GetTouchableHand 123 aStep->GetPreStepPoint()->GetTouchableHandle();
144 G4double kCarTolerance = << 124 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
145 G4GeometryTolerance::GetInstance()->GetSur <<
146 125
147 if(aStep->GetPreStepPoint()->GetStepStatus() << 126 if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
148 { <<
149 // Entering Geometry 127 // Entering Geometry
150 G4ThreeVector stppos1 = aStep->GetPreStepP << 128 G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
151 G4ThreeVector localpos1 = << 129 G4ThreeVector localpos1 =
152 theTouchable->GetHistory()->GetTopTransf 130 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
153 if(std::fabs(localpos1.z() + boxSolid->Get << 131 if(std::fabs( localpos1.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
154 { <<
155 return fCurrent_In; 132 return fCurrent_In;
156 } 133 }
157 } 134 }
158 135
159 if(aStep->GetPostStepPoint()->GetStepStatus( << 136 if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
160 { <<
161 // Exiting Geometry 137 // Exiting Geometry
162 G4ThreeVector stppos2 = aStep->GetPostStep << 138 G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
163 G4ThreeVector localpos2 = << 139 G4ThreeVector localpos2 =
164 theTouchable->GetHistory()->GetTopTransf 140 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
165 if(std::fabs(localpos2.z() + boxSolid->Get << 141 if(std::fabs( localpos2.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
166 { <<
167 return fCurrent_Out; 142 return fCurrent_Out;
168 } 143 }
169 } 144 }
170 145
171 return -1; 146 return -1;
172 } 147 }
173 148
174 void G4PSFlatSurfaceCurrent::Initialize(G4HCof 149 void G4PSFlatSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
175 { 150 {
176 EvtMap = new G4THitsMap<G4double>(detector-> 151 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
177 if(HCID < 0) << 152 if ( HCID < 0 ) HCID = GetCollectionID(0);
178 HCID = GetCollectionID(0); << 153 HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
179 HCE->AddHitsCollection(HCID, (G4VHitsCollect <<
180 } 154 }
181 155
182 void G4PSFlatSurfaceCurrent::clear() { EvtMap- << 156 void G4PSFlatSurfaceCurrent::EndOfEvent(G4HCofThisEvent*)
>> 157 {;}
>> 158
>> 159 void G4PSFlatSurfaceCurrent::clear(){
>> 160 EvtMap->clear();
>> 161 }
>> 162
>> 163 void G4PSFlatSurfaceCurrent::DrawAll()
>> 164 {;}
183 165
184 void G4PSFlatSurfaceCurrent::PrintAll() 166 void G4PSFlatSurfaceCurrent::PrintAll()
185 { 167 {
186 G4cout << " MultiFunctionalDet " << detecto 168 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
187 G4cout << " PrimitiveScorer " << GetName() < << 169 G4cout << " PrimitiveScorer " << GetName() <<G4endl;
188 G4cout << " Number of entries " << EvtMap->e 170 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
189 for(const auto& [copy, current] : *(EvtMap-> << 171 std::map<G4int,G4double*>::iterator itr = EvtMap->GetMap()->begin();
190 { << 172 for(; itr != EvtMap->GetMap()->end(); itr++) {
191 G4cout << " copy no.: " << copy << " curr << 173 G4cout << " copy no.: " << itr->first << " current : " ;
192 if(divideByArea) << 174 if ( divideByArea ) {
193 { << 175 G4cout << *(itr->second)/GetUnitValue()
194 G4cout << *(current) / GetUnitValue() << << 176 << " ["<<GetUnit()<<"]";
195 } << 177 }else {
196 else << 178 G4cout << *(itr->second)/GetUnitValue() << " [tracks]";
197 { << 179 }
198 G4cout << *(current) / GetUnitValue() << << 180 G4cout << G4endl;
199 } <<
200 G4cout << G4endl; <<
201 } 181 }
202 } 182 }
203 183
204 void G4PSFlatSurfaceCurrent::SetUnit(const G4S 184 void G4PSFlatSurfaceCurrent::SetUnit(const G4String& unit)
205 { 185 {
206 if(divideByArea) << 186 if ( divideByArea ) {
207 { << 187 CheckAndSetUnit(unit,"Per Unit Surface");
208 CheckAndSetUnit(unit, "Per Unit Surface"); << 188 } else {
209 } << 189 if (unit == "" ){
210 else << 190 unitName = unit;
211 { << 191 unitValue = 1.0;
212 if(unit.empty()) << 192 }else{
213 { << 193 G4String msg = "Invalid unit ["+unit+"] (Current unit is [" +GetUnit()+"] ) for " + GetName();
214 unitName = unit; << 194 G4Exception("G4PSFlatSurfaceCurrent::SetUnit","DetPS0007",JustWarning,msg);
215 unitValue = 1.0; << 195 }
216 } <<
217 else <<
218 { <<
219 G4String msg = "Invalid unit [" + unit + <<
220 GetUnit() + "] ) for " + <<
221 G4Exception("G4PSFlatSurfaceCurrent::Set <<
222 msg); <<
223 } 196 }
224 } <<
225 } 197 }
226 198
227 void G4PSFlatSurfaceCurrent::DefineUnitAndCate << 199 void G4PSFlatSurfaceCurrent::DefineUnitAndCategory(){
228 { << 200 // Per Unit Surface
229 // Per Unit Surface << 201 new G4UnitDefinition("percentimeter2","percm2","Per Unit Surface",(1./cm2));
230 new G4UnitDefinition("percentimeter2", "perc << 202 new G4UnitDefinition("permillimeter2","permm2","Per Unit Surface",(1./mm2));
231 (1. / cm2)); << 203 new G4UnitDefinition("permeter2","perm2","Per Unit Surface",(1./m2));
232 new G4UnitDefinition("permillimeter2", "perm <<
233 (1. / mm2)); <<
234 new G4UnitDefinition("permeter2", "perm2", " <<
235 } 204 }
>> 205
236 206