Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/digits_hits/scorer/src/G4PSSphereSurfaceCurrent.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 ]

Diff markup

Differences between /digits_hits/scorer/src/G4PSSphereSurfaceCurrent.cc (Version 11.3.0) and /digits_hits/scorer/src/G4PSSphereSurfaceCurrent.cc (Version 11.0.p4)


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