Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/digits_hits/scorer/src/G4PSCylinderSurfaceFlux.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 // // G4PSCylinderSurfaceFlux
 29 #include "G4PSCylinderSurfaceFlux.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 #include "G4VScoreHistFiller.hh"
 40 
 41 // ////////////////////////////////////////////////////////////////////////////////
 42 // (Description)
 43 //   This is a primitive scorer class for scoring Surface Flux.
 44 //  Current version assumes only for G4Tubs shape, and the surface
 45 //  is fixed on inner plane of the tube.
 46 //
 47 // Surface is defined at the innner surface of the tube.
 48 // Direction                   R    R+dR
 49 //   0  IN || OUT            ->|<-  |
 50 //   1  IN                   ->|    |
 51 //   2  OUT                    |<-  |
 52 //
 53 // Created: 2007-03-29  Tsukasa ASO
 54 // 2010-07-22   Introduce Unit specification.
 55 // 2010-07-22   Add weighted and divideByArea options
 56 // 2011-02-21   Get correct momentum direction in Flux_Out.
 57 // 2020-10-06   Use G4VPrimitivePlotter and fill 1-D histo of kinetic energy (x)
 58 //              vs. surface flux * track weight (y)        (Makoto Asai)
 59 //
 60 ///////////////////////////////////////////////////////////////////////////////
 61 
 62 G4PSCylinderSurfaceFlux::G4PSCylinderSurfaceFlux(const G4String& name, G4int direction,
 63                                                  G4int depth)
 64   : G4PSCylinderSurfaceFlux(name, direction, "percm2", depth) 
 65 {}
 66 
 67 G4PSCylinderSurfaceFlux::G4PSCylinderSurfaceFlux(const G4String& name, G4int direction,
 68                                                  const G4String& unit,
 69                                                  G4int depth)
 70   : G4VPrimitivePlotter(name, depth)
 71   , HCID(-1)
 72   , fDirection(direction)
 73   , EvtMap(nullptr)
 74   , weighted(true)
 75   , divideByArea(true)
 76 {
 77   DefineUnitAndCategory();
 78   SetUnit(unit);
 79 }
 80 
 81 G4bool G4PSCylinderSurfaceFlux::ProcessHits(G4Step* aStep, G4TouchableHistory*)
 82 {
 83   G4StepPoint* preStep = aStep->GetPreStepPoint();
 84   G4VSolid* solid      = ComputeCurrentSolid(aStep);
 85   assert(dynamic_cast<G4Tubs*>(solid));
 86 
 87   auto  tubsSolid = static_cast<G4Tubs*>(solid);
 88 
 89   G4int dirFlag = IsSelectedSurface(aStep, tubsSolid);
 90 
 91   if(dirFlag > 0)
 92   {
 93     if(fDirection == fFlux_InOut || dirFlag == fDirection)
 94     {
 95       G4StepPoint* thisStep = nullptr;
 96       if(dirFlag == fFlux_In)
 97       {
 98         thisStep = preStep;
 99       }
100       else if(dirFlag == fFlux_Out)
101       {
102         thisStep = aStep->GetPostStepPoint();
103       }
104       else
105       {
106         return false;
107       }
108 
109       G4TouchableHandle theTouchable = thisStep->GetTouchableHandle();
110       G4ThreeVector pdirection       = thisStep->GetMomentumDirection();
111       G4ThreeVector localdir =
112         theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection);
113       G4ThreeVector position = thisStep->GetPosition();
114       G4ThreeVector localpos =
115         theTouchable->GetHistory()->GetTopTransform().TransformAxis(position);
116       G4double angleFactor =
117         (localdir.x() * localpos.x() + localdir.y() * localpos.y()) /
118         std::sqrt(localdir.x() * localdir.x() + localdir.y() * localdir.y() +
119                   localdir.z() * localdir.z()) /
120         std::sqrt(localpos.x() * localpos.x() + localpos.y() * localpos.y());
121 
122       if(angleFactor < 0)
123         angleFactor *= -1.;
124       G4double square = 2. * tubsSolid->GetZHalfLength() *
125                         tubsSolid->GetInnerRadius() *
126                         tubsSolid->GetDeltaPhiAngle() / radian;
127 
128       G4double flux = 1.0;
129       if(weighted)
130         flux *= preStep->GetWeight();
131       // Current (Particle Weight)
132 
133       flux = flux / angleFactor;
134       if(divideByArea)
135         flux /= square;
136       // Flux with angle.
137       G4int index = GetIndex(aStep);
138       EvtMap->add(index, flux);
139 
140       if(!hitIDMap.empty() && hitIDMap.find(index) != hitIDMap.cend())
141       {
142         auto filler = G4VScoreHistFiller::Instance();
143         if(filler == nullptr)
144         {
145           G4Exception("G4PSCylinderSurfaceFlux::ProcessHits", "SCORER0123",
146                       JustWarning,
147                       "G4TScoreHistFiller is not instantiated!! Histogram is "
148                       "not filled.");
149         }
150         else
151         {
152           filler->FillH1(hitIDMap[index], thisStep->GetKineticEnergy(), flux);
153         }
154       }
155 
156       return true;
157     }
158     
159     return false;
160   }
161   
162   return false;
163 }
164 
165 G4int G4PSCylinderSurfaceFlux::IsSelectedSurface(G4Step* aStep,
166                                                  G4Tubs* tubsSolid)
167 {
168   G4TouchableHandle theTouchable =
169     aStep->GetPreStepPoint()->GetTouchableHandle();
170   G4double kCarTolerance =
171     G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
172 
173   if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
174   {
175     // Entering Geometry
176     G4ThreeVector stppos1 = aStep->GetPreStepPoint()->GetPosition();
177     G4ThreeVector localpos1 =
178       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
179     if(std::fabs(localpos1.z()) > tubsSolid->GetZHalfLength())
180       return -1;
181     // if(std::fabs( localpos1.x()*localpos1.x()+localpos1.y()*localpos1.y()
182     //   - (tubsSolid->GetInnerRadius()*tubsSolid->GetInnerRadius()))
183     //   <kCarTolerance ){
184     G4double localR2 =
185       localpos1.x() * localpos1.x() + localpos1.y() * localpos1.y();
186     G4double InsideRadius = tubsSolid->GetInnerRadius();
187     if(localR2 >
188          (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
189        localR2 <
190          (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
191     {
192       return fFlux_In;
193     }
194   }
195 
196   if(aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
197   {
198     // Exiting Geometry
199     G4ThreeVector stppos2 = aStep->GetPostStepPoint()->GetPosition();
200     G4ThreeVector localpos2 =
201       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
202     if(std::fabs(localpos2.z()) > tubsSolid->GetZHalfLength())
203       return -1;
204     // if(std::fabs( localpos2.x()*localpos2.x()+localpos2.y()*localpos2.y()
205     //   - (tubsSolid->GetInnerRadius()*tubsSolid->GetInnerRadius()))
206     // <kCarTolerance ){
207     G4double localR2 =
208       localpos2.x() * localpos2.x() + localpos2.y() * localpos2.y();
209     G4double InsideRadius = tubsSolid->GetInnerRadius();
210     if(localR2 >
211          (InsideRadius - kCarTolerance) * (InsideRadius - kCarTolerance) &&
212        localR2 <
213          (InsideRadius + kCarTolerance) * (InsideRadius + kCarTolerance))
214     {
215       return fFlux_Out;
216     }
217   }
218 
219   return -1;
220 }
221 
222 void G4PSCylinderSurfaceFlux::Initialize(G4HCofThisEvent* HCE)
223 {
224   EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(),
225                                     GetName());
226   if(HCID < 0)
227     HCID = GetCollectionID(0);
228   HCE->AddHitsCollection(HCID, (G4VHitsCollection*) EvtMap);
229 }
230 
231 void G4PSCylinderSurfaceFlux::clear() { EvtMap->clear(); }
232 
233 void G4PSCylinderSurfaceFlux::PrintAll()
234 {
235   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
236   G4cout << " PrimitiveScorer" << GetName() << G4endl;
237   G4cout << " Number of entries " << EvtMap->entries() << G4endl;
238   for(const auto& [copy, flux] : *(EvtMap->GetMap()))
239   {
240     G4cout << "  copy no.: " << copy
241            << "  flux  : " << *(flux) / GetUnitValue() << " ["
242            << GetUnit() << "]" << G4endl;
243   }
244 }
245 
246 void G4PSCylinderSurfaceFlux::SetUnit(const G4String& unit)
247 {
248   if(divideByArea)
249   {
250     CheckAndSetUnit(unit, "Per Unit Surface");
251   }
252   else
253   {
254     if(unit.empty())
255     {
256       unitName  = unit;
257       unitValue = 1.0;
258     }
259     else
260     {
261       G4String msg = "Invalid unit [" + unit + "] (Current  unit is [" +
262                      GetUnit() + "] ) for " + GetName();
263       G4Exception("G4PSCylinderSurfaceFlux::SetUnit", "DetPS0003", JustWarning,
264                   msg);
265     }
266   }
267 }
268 
269 void G4PSCylinderSurfaceFlux::DefineUnitAndCategory()
270 {
271   // Per Unit Surface
272   new G4UnitDefinition("percentimeter2", "percm2", "Per Unit Surface",
273                        (1. / cm2));
274   new G4UnitDefinition("permillimeter2", "permm2", "Per Unit Surface",
275                        (1. / mm2));
276   new G4UnitDefinition("permeter2", "perm2", "Per Unit Surface", (1. / m2));
277 }
278