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 ]

Diff markup

Differences between /digits_hits/scorer/src/G4PSSphereSurfaceFlux.cc (Version 11.3.0) and /digits_hits/scorer/src/G4PSSphereSurfaceFlux.cc (Version 9.1)


  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 // $Id: G4PSSphereSurfaceFlux.cc,v 1.1 2007/07/11 01:31:03 asaim Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01 $
 27 //                                                 29 //
 28 // G4PSSphereSurfaceFlux                           30 // G4PSSphereSurfaceFlux
 29 #include "G4PSSphereSurfaceFlux.hh"                31 #include "G4PSSphereSurfaceFlux.hh"
 30                                                << 
 31 #include "G4SystemOfUnits.hh"                  << 
 32 #include "G4StepStatus.hh"                         32 #include "G4StepStatus.hh"
 33 #include "G4Track.hh"                              33 #include "G4Track.hh"
 34 #include "G4VSolid.hh"                         << 
 35 #include "G4VPhysicalVolume.hh"                << 
 36 #include "G4VPVParameterisation.hh"            << 
 37 #include "G4UnitsTable.hh"                         34 #include "G4UnitsTable.hh"
 38 #include "G4GeometryTolerance.hh"                  35 #include "G4GeometryTolerance.hh"
 39 //////////////////////////////////////////////     36 ////////////////////////////////////////////////////////////////////////////////
 40 // (Description)                                   37 // (Description)
 41 //   This is a primitive scorer class for scor     38 //   This is a primitive scorer class for scoring only Surface Flux.
 42 //  Flux version assumes only for G4Sphere sha <<  39 //  Flux version assumes only for G4Sphere shape. 
 43 //                                                 40 //
 44 // Surface is defined  at the inside of sphere     41 // Surface is defined  at the inside of sphere.
 45 // Direction                  -Rmin   +Rmax        42 // Direction                  -Rmin   +Rmax
 46 //   0  IN || OUT            ->|<-     |           43 //   0  IN || OUT            ->|<-     |
 47 //   1  IN                   ->|       |           44 //   1  IN                   ->|       |
 48 //   2  OUT                    |<-     |           45 //   2  OUT                    |<-     |
 49 //                                                 46 //
 50 // Created: 2005-11-14  Tsukasa ASO, Akinori K     47 // Created: 2005-11-14  Tsukasa ASO, Akinori Kimura.
 51 // 29-Mar-2007  T.Aso,  Bug fix for momentum d     48 // 29-Mar-2007  T.Aso,  Bug fix for momentum direction at outgoing flux.
 52 // 2010-07-22   Introduce Unit specification.  <<  49 // 
 53 // 2010-07-22   Add weighted and divideByAre o << 
 54 // 2011-02-21   Get correct momentum direction << 
 55 // 2011-09-09   Modify comment in PrintAll().  << 
 56 // 2014-03-03  T.Aso,  To use always positive  << 
 57 //////////////////////////////////////////////     50 ///////////////////////////////////////////////////////////////////////////////
 58                                                    51 
 59 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(c <<  52 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(G4String name, 
 60                                              G <<  53            G4int direction, G4int depth)
 61   : G4PSSphereSurfaceFlux(name, direction, "pe <<  54   :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction)
 62 {}                                             <<  55 {;}
 63                                                <<  56 
 64 G4PSSphereSurfaceFlux::G4PSSphereSurfaceFlux(c <<  57 G4PSSphereSurfaceFlux::~G4PSSphereSurfaceFlux()
 65                                              c <<  58 {;}
 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                                                    59 
 77 G4bool G4PSSphereSurfaceFlux::ProcessHits(G4St <<  60 G4bool G4PSSphereSurfaceFlux::ProcessHits(G4Step* aStep,G4TouchableHistory*)
 78 {                                                  61 {
 79   G4StepPoint* preStep = aStep->GetPreStepPoin     62   G4StepPoint* preStep = aStep->GetPreStepPoint();
 80                                                <<  63   G4VSolid * solid = 
 81   G4VPhysicalVolume* physVol       = preStep-> <<  64     preStep->GetPhysicalVolume()->GetLogicalVolume()->GetSolid();
 82   G4VPVParameterisation* physParam = physVol-> <<  65   if( solid->GetEntityType() != "G4Sphere" ){
 83   G4VSolid* solid                  = nullptr;  <<  66     G4Exception("G4PSSphereSurfaceFluxScorer. - Solid type is not supported.");
 84   if(physParam != nullptr)                     <<  67     return FALSE;
 85   {  // for parameterized volume               <<  68   }
 86     G4int idx =                                <<  69   G4Sphere* sphereSolid = (G4Sphere*)(solid);
 87       ((G4TouchableHistory*) (aStep->GetPreSte <<  70 
 88         ->GetReplicaNumber(indexDepth);        <<  71   G4int dirFlag =IsSelectedSurface(aStep,sphereSolid);
 89     solid = physParam->ComputeSolid(idx, physV <<  72   if ( dirFlag > 0 ) {
 90     solid->ComputeDimensions(physParam, idx, p <<  73     if ( fDirection == fFlux_InOut || fDirection == dirFlag ){
 91   }                                            <<  74 
 92   else                                         <<  75       G4StepPoint* thisStep=0;
 93   {  // for ordinary volume                    <<  76       if ( dirFlag == fFlux_In ){
 94     solid = physVol->GetLogicalVolume()->GetSo <<  77   thisStep = preStep;
 95   }                                            <<  78       }else if ( dirFlag == fFlux_Out ){
 96                                                <<  79   thisStep = aStep->GetPreStepPoint();
 97   auto  sphereSolid = (G4Sphere*) (solid);     <<  80       }else{
 98                                                <<  81   return FALSE;
 99   G4int dirFlag = IsSelectedSurface(aStep, sph << 
100   if(dirFlag > 0)                              << 
101   {                                            << 
102     if(fDirection == fFlux_InOut || fDirection << 
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       }                                            82       }
117                                                    83 
118       G4TouchableHandle theTouchable = thisSte     84       G4TouchableHandle theTouchable = thisStep->GetTouchableHandle();
119       G4ThreeVector pdirection       = thisSte <<  85       G4ThreeVector pdirection = thisStep->GetMomentumDirection();
120       G4ThreeVector localdir =                 <<  86       G4ThreeVector localdir  = 
121         theTouchable->GetHistory()->GetTopTran <<  87   theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection);
122       G4double localdirL2 = localdir.x() * loc <<  88       G4double localdirL2 = localdir.x()*localdir.x()
123                             localdir.y() * loc <<  89   +localdir.y()*localdir.y()
124                             localdir.z() * loc <<  90   +localdir.z()*localdir.z();
125       G4ThreeVector stppos1 = aStep->GetPreSte <<  91       G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
126       G4ThreeVector localpos1 =                <<  92       G4ThreeVector localpos1 = 
127         theTouchable->GetHistory()->GetTopTran <<  93   theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
128       G4double localR2 = localpos1.x() * local <<  94       G4double localR2 = localpos1.x()*localpos1.x()
129                          localpos1.y() * local <<  95   +localpos1.y()*localpos1.y()
130                          localpos1.z() * local <<  96   +localpos1.z()*localpos1.z();
131       G4double anglefactor =                   <<  97       G4double anglefactor = (localdir.x()*localpos1.x()
132         (localdir.x() * localpos1.x() + locald <<  98             +localdir.y()*localpos1.y()
133          localdir.z() * localpos1.z()) /       <<  99             +localdir.z()*localpos1.z())
134         std::sqrt(localdirL2) / std::sqrt(loca << 100   /std::sqrt(localdirL2)/std::sqrt(localR2);
135       if(anglefactor < 0.0)                    << 101 
136         anglefactor *= -1.0;                   << 102       G4double radi   = sphereSolid->GetInsideRadius();
137                                                << 103       G4double dph    = sphereSolid->GetDeltaPhiAngle()/radian;
138       G4double current = 1.0 / anglefactor;    << 104       G4double stth   = sphereSolid->GetStartThetaAngle()/radian;
139       if(weighted)                             << 105       G4double enth   = stth+sphereSolid->GetDeltaThetaAngle()/radian;
140         current *= thisStep->GetWeight();  //  << 106       G4double square = radi*radi*dph*( -std::cos(enth) + std::cos(stth) );
141       if(divideByArea)                     //  << 107 
142       {                                        << 108       G4double current = thisStep->GetWeight(); // Flux (Particle Weight)
143         G4double radi = sphereSolid->GetInnerR << 109       current = current/square;  // Flux with angle.
144         G4double dph  = sphereSolid->GetDeltaP << 110 
145         G4double stth = sphereSolid->GetStartT << 111       current /= anglefactor;
146         G4double enth = stth + sphereSolid->Ge << 
147         current /= radi * radi * dph * (-std:: << 
148       }                                        << 
149                                                   112 
150       G4int index = GetIndex(aStep);              113       G4int index = GetIndex(aStep);
151       EvtMap->add(index, current);             << 114       EvtMap->add(index,current);
152     }                                             115     }
153   }                                               116   }
154                                                   117 
155   return true;                                 << 118   return TRUE;
156 }                                                 119 }
157                                                   120 
158 G4int G4PSSphereSurfaceFlux::IsSelectedSurface << 121 G4int G4PSSphereSurfaceFlux::IsSelectedSurface(G4Step* aStep, G4Sphere* sphereSolid){
159                                                << 
160 {                                              << 
161   G4TouchableHandle theTouchable =             << 
162     aStep->GetPreStepPoint()->GetTouchableHand << 
163   G4double kCarTolerance =                     << 
164     G4GeometryTolerance::GetInstance()->GetSur << 
165                                                   122 
166   if(aStep->GetPreStepPoint()->GetStepStatus() << 123   G4TouchableHandle theTouchable = 
167   {                                            << 124     aStep->GetPreStepPoint()->GetTouchableHandle();
                                                   >> 125   G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
                                                   >> 126   
                                                   >> 127   if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
168     // Entering Geometry                          128     // Entering Geometry
169     G4ThreeVector stppos1 = aStep->GetPreStepP << 129     G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
170     G4ThreeVector localpos1 =                  << 130     G4ThreeVector localpos1 = 
171       theTouchable->GetHistory()->GetTopTransf    131       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
172     G4double localR2 = localpos1.x() * localpo << 132     G4double localR2 = localpos1.x()*localpos1.x()
173                        localpos1.y() * localpo << 133                       +localpos1.y()*localpos1.y()
174                        localpos1.z() * localpo << 134                       +localpos1.z()*localpos1.z();
175     // G4double InsideRadius2 =                << 135     G4double InsideRadius2 = 
176     //  sphereSolid->GetInsideRadius()*sphereS << 136       sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
177     // if(std::fabs( localR2 - InsideRadius2 ) << 137     if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
178     G4double InsideRadius = sphereSolid->GetIn << 
179     if(localR2 >                               << 
180          (InsideRadius - kCarTolerance) * (Ins << 
181        localR2 <                               << 
182          (InsideRadius + kCarTolerance) * (Ins << 
183     {                                          << 
184       return fFlux_In;                            138       return fFlux_In;
185     }                                             139     }
186   }                                               140   }
187                                                   141 
188   if(aStep->GetPostStepPoint()->GetStepStatus( << 142   if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
189   {                                            << 
190     // Exiting Geometry                           143     // Exiting Geometry
191     G4ThreeVector stppos2 = aStep->GetPostStep << 144     G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
192     G4ThreeVector localpos2 =                  << 145     G4ThreeVector localpos2 = 
193       theTouchable->GetHistory()->GetTopTransf    146       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
194     G4double localR2 = localpos2.x() * localpo << 147     G4double localR2 = localpos2.x()*localpos2.x()
195                        localpos2.y() * localpo << 148                       +localpos2.y()*localpos2.y()
196                        localpos2.z() * localpo << 149                       +localpos2.z()*localpos2.z();
197     // G4double InsideRadius2 =                << 150     G4double InsideRadius2 = 
198     //  sphereSolid->GetInsideRadius()*sphereS << 151       sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
199     // if(std::facb(localR2 - InsideRadius2) ) << 152     if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
200     G4double InsideRadius = sphereSolid->GetIn << 
201     if(localR2 >                               << 
202          (InsideRadius - kCarTolerance) * (Ins << 
203        localR2 <                               << 
204          (InsideRadius + kCarTolerance) * (Ins << 
205     {                                          << 
206       return fFlux_Out;                           153       return fFlux_Out;
207     }                                             154     }
208   }                                               155   }
209                                                   156 
210   return -1;                                      157   return -1;
211 }                                                 158 }
212                                                   159 
213 void G4PSSphereSurfaceFlux::Initialize(G4HCofT    160 void G4PSSphereSurfaceFlux::Initialize(G4HCofThisEvent* HCE)
214 {                                                 161 {
215   EvtMap = new G4THitsMap<G4double>(detector->    162   EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
216   if(HCID < 0)                                 << 163   if ( HCID < 0 ) HCID = GetCollectionID(0);
217     HCID = GetCollectionID(0);                 << 164   HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
218   HCE->AddHitsCollection(HCID, (G4VHitsCollect << 
219 }                                                 165 }
220                                                   166 
221 void G4PSSphereSurfaceFlux::clear() { EvtMap-> << 167 void G4PSSphereSurfaceFlux::EndOfEvent(G4HCofThisEvent*)
                                                   >> 168 {;}
                                                   >> 169 
                                                   >> 170 void G4PSSphereSurfaceFlux::clear(){
                                                   >> 171   EvtMap->clear();
                                                   >> 172 }
                                                   >> 173 
                                                   >> 174 void G4PSSphereSurfaceFlux::DrawAll()
                                                   >> 175 {;}
222                                                   176 
223 void G4PSSphereSurfaceFlux::PrintAll()            177 void G4PSSphereSurfaceFlux::PrintAll()
224 {                                                 178 {
225   G4cout << " MultiFunctionalDet  " << detecto    179   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
226   G4cout << " PrimitiveScorer " << GetName() < << 180   G4cout << " PrimitiveScorer " << GetName() <<G4endl; 
227   G4cout << " Number of entries " << EvtMap->e    181   G4cout << " Number of entries " << EvtMap->entries() << G4endl;
228   for(const auto& [copy, flux] : *(EvtMap->Get << 182   std::map<G4int,G4double*>::iterator itr = EvtMap->GetMap()->begin();
229   {                                            << 183   for(; itr != EvtMap->GetMap()->end(); itr++) {
230     G4cout << "  copy no.: " << copy           << 184     G4cout << "  copy no.: " << itr->first
231            << "  Flux  : " << *(flux) / GetUni << 185      << "  current  : " << *(itr->second)
232            << GetUnit() << "]" << G4endl;      << 186      << G4endl;
233   }                                            << 
234 }                                              << 
235                                                << 
236 void G4PSSphereSurfaceFlux::SetUnit(const G4St << 
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 + << 
252                      GetUnit() + "] ) for " +  << 
253       G4Exception("G4PSSphereSurfaceFlux::SetU << 
254                   msg);                        << 
255     }                                          << 
256   }                                               187   }
257 }                                                 188 }
258                                                   189 
259 void G4PSSphereSurfaceFlux::DefineUnitAndCateg << 
260 {                                              << 
261   // Per Unit Surface                          << 
262   new G4UnitDefinition("percentimeter2", "perc << 
263                        (1. / cm2));            << 
264   new G4UnitDefinition("permillimeter2", "perm << 
265                        (1. / mm2));            << 
266   new G4UnitDefinition("permeter2", "perm2", " << 
267 }                                              << 
268                                                   190