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