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