Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/digits_hits/scorer/src/G4PSSphereSurfaceFlux.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // G4PSSphereSurfaceFlux
 29 #include "G4PSSphereSurfaceFlux.hh"
 30 
 31 #include "G4SystemOfUnits.hh"
 32 #include "G4StepStatus.hh"
 33 #include "G4Track.hh"
 34 #include "G4VSolid.hh"
 35 #include "G4VPhysicalVolume.hh"
 36 #include "G4VPVParameterisation.hh"
 37 #include "G4UnitsTable.hh"
 38 #include "G4GeometryTolerance.hh"
 39 ////////////////////////////////////////////////////////////////////////////////
 40 // (Description)
 41 //   This is a primitive scorer class for scoring only Surface Flux.
 42 //  Flux version assumes only for G4Sphere shape.
 43 //
 44 // Surface is defined  at the inside of sphere.
 45 // Direction                  -Rmin   +Rmax
 46 //   0  IN || OUT            ->|<-     |
 47 //   1  IN                   ->|       |
 48 //   2  OUT                    |<-     |
 49 //
 50 // Created: 2005-11-14  Tsukasa ASO, Akinori Kimura.
 51 // 29-Mar-2007  T.Aso,  Bug fix for momentum direction at outgoing flux.
 52 // 2010-07-22   Introduce Unit specification.
 53 // 2010-07-22   Add weighted and divideByAre options
 54 // 2011-02-21   Get correct momentum direction in Flux_Out.
 55 // 2011-09-09   Modify comment in PrintAll().
 56 // 2014-03-03  T.Aso,  To use always positive value for anglefactor.
 57 ///////////////////////////////////////////////////////////////////////////////
 58 
 59 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(const G4String& name, G4int direction,
 60                                              G4int depth)
 61   : G4PSSphereSurfaceFlux(name, direction, "percm2", depth)
 62 {}
 63 
 64 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(const G4String& name, G4int direction,
 65                                              const G4String& unit, G4int depth)
 66   : G4VPrimitiveScorer(name, depth)
 67   , HCID(-1)
 68   , fDirection(direction)
 69   , EvtMap(nullptr)
 70   , weighted(true)
 71   , divideByArea(true)
 72 {
 73   DefineUnitAndCategory();
 74   SetUnit(unit);
 75 }
 76 
 77 G4bool G4PSSphereSurfaceFlux::ProcessHits(G4Step* aStep, G4TouchableHistory*)
 78 {
 79   G4StepPoint* preStep = aStep->GetPreStepPoint();
 80 
 81   G4VPhysicalVolume* physVol       = preStep->GetPhysicalVolume();
 82   G4VPVParameterisation* physParam = physVol->GetParameterisation();
 83   G4VSolid* solid                  = nullptr;
 84   if(physParam != nullptr)
 85   {  // for parameterized volume
 86     G4int idx =
 87       ((G4TouchableHistory*) (aStep->GetPreStepPoint()->GetTouchable()))
 88         ->GetReplicaNumber(indexDepth);
 89     solid = physParam->ComputeSolid(idx, physVol);
 90     solid->ComputeDimensions(physParam, idx, physVol);
 91   }
 92   else
 93   {  // for ordinary volume
 94     solid = physVol->GetLogicalVolume()->GetSolid();
 95   }
 96 
 97   auto  sphereSolid = (G4Sphere*) (solid);
 98 
 99   G4int dirFlag = IsSelectedSurface(aStep, sphereSolid);
100   if(dirFlag > 0)
101   {
102     if(fDirection == fFlux_InOut || fDirection == dirFlag)
103     {
104       G4StepPoint* thisStep = nullptr;
105       if(dirFlag == fFlux_In)
106       {
107         thisStep = preStep;
108       }
109       else if(dirFlag == fFlux_Out)
110       {
111         thisStep = aStep->GetPostStepPoint();
112       }
113       else
114       {
115         return false;
116       }
117 
118       G4TouchableHandle theTouchable = thisStep->GetTouchableHandle();
119       G4ThreeVector pdirection       = thisStep->GetMomentumDirection();
120       G4ThreeVector localdir =
121         theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection);
122       G4double localdirL2 = localdir.x() * localdir.x() +
123                             localdir.y() * localdir.y() +
124                             localdir.z() * localdir.z();
125       G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
126       G4ThreeVector localpos1 =
127         theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
128       G4double localR2 = localpos1.x() * localpos1.x() +
129                          localpos1.y() * localpos1.y() +
130                          localpos1.z() * localpos1.z();
131       G4double anglefactor =
132         (localdir.x() * localpos1.x() + localdir.y() * localpos1.y() +
133          localdir.z() * localpos1.z()) /
134         std::sqrt(localdirL2) / std::sqrt(localR2);
135       if(anglefactor < 0.0)
136         anglefactor *= -1.0;
137 
138       G4double current = 1.0 / anglefactor;
139       if(weighted)
140         current *= thisStep->GetWeight();  // Flux (Particle Weight)
141       if(divideByArea)                     // Flux with angle.
142       {
143         G4double radi = sphereSolid->GetInnerRadius();
144         G4double dph  = sphereSolid->GetDeltaPhiAngle() / radian;
145         G4double stth = sphereSolid->GetStartThetaAngle() / radian;
146         G4double enth = stth + sphereSolid->GetDeltaThetaAngle() / radian;
147         current /= radi * radi * dph * (-std::cos(enth) + std::cos(stth));
148       }
149 
150       G4int index = GetIndex(aStep);
151       EvtMap->add(index, current);
152     }
153   }
154 
155   return true;
156 }
157 
158 G4int G4PSSphereSurfaceFlux::IsSelectedSurface(G4Step* aStep,
159                                                G4Sphere* sphereSolid)
160 {
161   G4TouchableHandle theTouchable =
162     aStep->GetPreStepPoint()->GetTouchableHandle();
163   G4double kCarTolerance =
164     G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
165 
166   if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
167   {
168     // Entering Geometry
169     G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
170     G4ThreeVector localpos1 =
171       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
172     G4double localR2 = localpos1.x() * localpos1.x() +
173                        localpos1.y() * localpos1.y() +
174                        localpos1.z() * localpos1.z();
175     // G4double InsideRadius2 =
176     //  sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
177     // if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
178     G4double InsideRadius = sphereSolid->GetInnerRadius();
179     if(localR2 >
180          (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
181        localR2 <
182          (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
183     {
184       return fFlux_In;
185     }
186   }
187 
188   if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
189   {
190     // Exiting Geometry
191     G4ThreeVector stppos2 = aStep->GetPostStepPoint()->GetPosition();
192     G4ThreeVector localpos2 =
193       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
194     G4double localR2 = localpos2.x() * localpos2.x() +
195                        localpos2.y() * localpos2.y() +
196                        localpos2.z() * localpos2.z();
197     // G4double InsideRadius2 =
198     //  sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
199     // if(std::facb(localR2 - InsideRadius2) ) < kCarTolerance ){
200     G4double InsideRadius = sphereSolid->GetInnerRadius();
201     if(localR2 >
202          (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
203        localR2 <
204          (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
205     {
206       return fFlux_Out;
207     }
208   }
209 
210   return -1;
211 }
212 
213 void G4PSSphereSurfaceFlux::Initialize(G4HCofThisEvent* HCE)
214 {
215   EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
216   if(HCID < 0)
217     HCID = GetCollectionID(0);
218   HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap);
219 }
220 
221 void G4PSSphereSurfaceFlux::clear() { EvtMap->clear(); }
222 
223 void G4PSSphereSurfaceFlux::PrintAll()
224 {
225   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
226   G4cout << " PrimitiveScorer " << GetName() << G4endl;
227   G4cout << " Number of entries " << EvtMap->entries() << G4endl;
228   for(const auto& [copy, flux] : *(EvtMap->GetMap()))
229   {
230     G4cout << "  copy no.: " << copy
231            << "  Flux  : " << *(flux) / GetUnitValue() << " ["
232            << GetUnit() << "]" << G4endl;
233   }
234 }
235 
236 void G4PSSphereSurfaceFlux::SetUnit(const G4String& unit)
237 {
238   if(divideByArea)
239   {
240     CheckAndSetUnit(unit, "Per Unit Surface");
241   }
242   else
243   {
244     if(unit.empty())
245     {
246       unitName  = unit;
247       unitValue = 1.0;
248     }
249     else
250     {
251       G4String msg = "Invalid unit [" + unit + "] (Current  unit is [" +
252                      GetUnit() + "] ) for " + GetName();
253       G4Exception("G4PSSphereSurfaceFlux::SetUnit", "DetPS0016", JustWarning,
254                   msg);
255     }
256   }
257 }
258 
259 void G4PSSphereSurfaceFlux::DefineUnitAndCategory()
260 {
261   // Per Unit Surface
262   new G4UnitDefinition("percentimeter2", "percm2", "Per Unit Surface",
263                        (1. / cm2));
264   new G4UnitDefinition("permillimeter2", "permm2", "Per Unit Surface",
265                        (1. / mm2));
266   new G4UnitDefinition("permeter2", "perm2", "Per Unit Surface", (1. / m2));
267 }
268