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 // Class G4PSVolumeFlux 26 // Class G4PSVolumeFlux 27 // 27 // 28 // Class description: 28 // Class description: 29 // Scorer that scores number of tracks coming 29 // Scorer that scores number of tracks coming into the associated volume. 30 // Optionally number of tracks can be divided 30 // Optionally number of tracks can be divided by the surface area to 31 // score the volume current and divided by co 31 // score the volume current and divided by cos(theta) for volume flux 32 // where theta is the incident angle 32 // where theta is the incident angle 33 // 33 // 34 // - Created M. Asai, Sept. 2020 34 // - Created M. Asai, Sept. 2020 35 // 35 // 36 // 36 // 37 #include "G4PSVolumeFlux.hh" 37 #include "G4PSVolumeFlux.hh" 38 38 39 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh" 40 #include "G4StepStatus.hh" 40 #include "G4StepStatus.hh" 41 #include "G4Track.hh" 41 #include "G4Track.hh" 42 #include "G4VSolid.hh" 42 #include "G4VSolid.hh" 43 #include "G4VPhysicalVolume.hh" 43 #include "G4VPhysicalVolume.hh" 44 #include "G4VPVParameterisation.hh" 44 #include "G4VPVParameterisation.hh" 45 #include "G4UnitsTable.hh" 45 #include "G4UnitsTable.hh" 46 #include "G4GeometryTolerance.hh" 46 #include "G4GeometryTolerance.hh" 47 #include "G4VScoreHistFiller.hh" 47 #include "G4VScoreHistFiller.hh" 48 48 49 G4PSVolumeFlux::G4PSVolumeFlux(const G4String& << 49 G4PSVolumeFlux::G4PSVolumeFlux(G4String name, G4int direction, G4int depth) 50 : G4VPrimitivePlotter(name, depth) 50 : G4VPrimitivePlotter(name, depth) 51 , fDirection(direction) 51 , fDirection(direction) 52 {} 52 {} 53 53 54 G4bool G4PSVolumeFlux::ProcessHits(G4Step* aSt 54 G4bool G4PSVolumeFlux::ProcessHits(G4Step* aStep, G4TouchableHistory*) 55 { 55 { 56 G4StepPoint* preStepPoint = aStep->GetPreSt 56 G4StepPoint* preStepPoint = aStep->GetPreStepPoint(); 57 G4StepPoint* postStepPoint = aStep->GetPostS 57 G4StepPoint* postStepPoint = aStep->GetPostStepPoint(); 58 G4StepPoint* thisStepPoint = nullptr; 58 G4StepPoint* thisStepPoint = nullptr; 59 if(fDirection == 1) // Score only the inwar 59 if(fDirection == 1) // Score only the inward particle 60 { 60 { 61 if(preStepPoint->GetStepStatus() == fGeomB 61 if(preStepPoint->GetStepStatus() == fGeomBoundary) 62 { 62 { 63 thisStepPoint = preStepPoint; 63 thisStepPoint = preStepPoint; 64 } 64 } 65 else 65 else 66 { 66 { 67 return false; 67 return false; 68 } 68 } 69 } 69 } 70 else if(fDirection == 2) // Score only the 70 else if(fDirection == 2) // Score only the outward particle 71 { 71 { 72 if(postStepPoint->GetStepStatus() == fGeom 72 if(postStepPoint->GetStepStatus() == fGeomBoundary) 73 { 73 { 74 thisStepPoint = postStepPoint; 74 thisStepPoint = postStepPoint; 75 } 75 } 76 else 76 else 77 { 77 { 78 return false; 78 return false; 79 } 79 } 80 } 80 } 81 81 82 G4double flux = preStepPoint->GetWeight(); 82 G4double flux = preStepPoint->GetWeight(); 83 83 84 if(divare || divcos) 84 if(divare || divcos) 85 { 85 { 86 G4VPhysicalVolume* physVol = preStep 86 G4VPhysicalVolume* physVol = preStepPoint->GetPhysicalVolume(); 87 G4VPVParameterisation* physParam = physVol 87 G4VPVParameterisation* physParam = physVol->GetParameterisation(); 88 G4VSolid* solid = nullptr 88 G4VSolid* solid = nullptr; 89 if(physParam != nullptr) 89 if(physParam != nullptr) 90 { // for parameterized volume 90 { // for parameterized volume 91 auto idx = ((G4TouchableHistory*) (preSt 91 auto idx = ((G4TouchableHistory*) (preStepPoint->GetTouchable())) 92 ->GetReplicaNumber(indexDep 92 ->GetReplicaNumber(indexDepth); 93 solid = physParam->ComputeSolid(idx, phy 93 solid = physParam->ComputeSolid(idx, physVol); 94 solid->ComputeDimensions(physParam, idx, 94 solid->ComputeDimensions(physParam, idx, physVol); 95 } 95 } 96 else 96 else 97 { // for ordinary volume 97 { // for ordinary volume 98 solid = physVol->GetLogicalVolume()->Get 98 solid = physVol->GetLogicalVolume()->GetSolid(); 99 } 99 } 100 100 101 if(divare) 101 if(divare) 102 { 102 { 103 flux /= solid->GetSurfaceArea(); 103 flux /= solid->GetSurfaceArea(); 104 } 104 } 105 105 106 if(divcos) 106 if(divcos) 107 { 107 { 108 G4TouchableHandle theTouchable = thisSte 108 G4TouchableHandle theTouchable = thisStepPoint->GetTouchableHandle(); 109 G4ThreeVector pdirection = thisSte 109 G4ThreeVector pdirection = thisStepPoint->GetMomentumDirection(); 110 G4ThreeVector localdir = 110 G4ThreeVector localdir = 111 theTouchable->GetHistory()->GetTopTran 111 theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection); 112 G4ThreeVector globalPos = thisStepPoint- 112 G4ThreeVector globalPos = thisStepPoint->GetPosition(); 113 G4ThreeVector localPos = 113 G4ThreeVector localPos = 114 theTouchable->GetHistory()->GetTopTran 114 theTouchable->GetHistory()->GetTopTransform().TransformPoint(globalPos); 115 G4ThreeVector surfNormal = solid->Surfac 115 G4ThreeVector surfNormal = solid->SurfaceNormal(localPos); 116 G4double cosT = surfNormal.co 116 G4double cosT = surfNormal.cosTheta(localdir); 117 if(cosT != 0.) 117 if(cosT != 0.) 118 flux /= std::abs(cosT); 118 flux /= std::abs(cosT); 119 } 119 } 120 } 120 } 121 121 122 G4int index = GetIndex(aStep); 122 G4int index = GetIndex(aStep); 123 EvtMap->add(index, flux); 123 EvtMap->add(index, flux); 124 124 125 if(!hitIDMap.empty() && hitIDMap.find(index) 125 if(!hitIDMap.empty() && hitIDMap.find(index) != hitIDMap.cend()) 126 { 126 { 127 auto filler = G4VScoreHistFiller::Instance 127 auto filler = G4VScoreHistFiller::Instance(); 128 if(filler == nullptr) 128 if(filler == nullptr) 129 { 129 { 130 G4Exception( 130 G4Exception( 131 "G4PSVolumeFlux::ProcessHits", "SCORER 131 "G4PSVolumeFlux::ProcessHits", "SCORER0123", JustWarning, 132 "G4TScoreHistFiller is not instantiate 132 "G4TScoreHistFiller is not instantiated!! Histogram is not filled."); 133 } 133 } 134 else 134 else 135 { 135 { 136 filler->FillH1(hitIDMap[index], thisStep 136 filler->FillH1(hitIDMap[index], thisStepPoint->GetKineticEnergy(), flux); 137 } 137 } 138 } 138 } 139 139 140 return true; 140 return true; 141 } 141 } 142 142 143 void G4PSVolumeFlux::Initialize(G4HCofThisEven 143 void G4PSVolumeFlux::Initialize(G4HCofThisEvent* HCE) 144 { 144 { 145 if(HCID < 0) 145 if(HCID < 0) 146 HCID = GetCollectionID(0); 146 HCID = GetCollectionID(0); 147 EvtMap = new G4THitsMap<G4double>(GetMultiFu 147 EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), 148 GetName()) 148 GetName()); 149 HCE->AddHitsCollection(HCID, (G4VHitsCollect 149 HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap); 150 } 150 } 151 151 152 void G4PSVolumeFlux::clear() { EvtMap->clear() 152 void G4PSVolumeFlux::clear() { EvtMap->clear(); } 153 153 154 void G4PSVolumeFlux::PrintAll() 154 void G4PSVolumeFlux::PrintAll() 155 { 155 { 156 G4cout << " MultiFunctionalDet " << detecto 156 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl; 157 G4cout << " PrimitiveScorer" << GetName() << 157 G4cout << " PrimitiveScorer" << GetName() << G4endl; 158 G4cout << " Number of entries " << EvtMap->e 158 G4cout << " Number of entries " << EvtMap->entries() << G4endl; 159 for(const auto& [copy, flux] : *(EvtMap->Get 159 for(const auto& [copy, flux] : *(EvtMap->GetMap())) 160 { 160 { 161 G4cout << " copy no.: " << copy << " flu 161 G4cout << " copy no.: " << copy << " flux : " << *(flux) 162 << G4endl; 162 << G4endl; 163 } 163 } 164 } 164 } 165 165