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 10.5)


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