Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/digits_hits/scorer/src/G4PSFlatSurfaceFlux.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/G4PSFlatSurfaceFlux.cc (Version 11.3.0) and /digits_hits/scorer/src/G4PSFlatSurfaceFlux.cc (Version 10.7.p1)


  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 // G4PSFlatSurfaceFlux                             28 // G4PSFlatSurfaceFlux
 29 #include "G4PSFlatSurfaceFlux.hh"                  29 #include "G4PSFlatSurfaceFlux.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 #include "G4VScoreHistFiller.hh"                   39 #include "G4VScoreHistFiller.hh"
 40                                                    40 
 41 //////////////////////////////////////////////     41 ////////////////////////////////////////////////////////////////////////////////
 42 // (Description)                                   42 // (Description)
 43 //   This is a primitive scorer class for scor     43 //   This is a primitive scorer class for scoring Surface Flux.
 44 //  Current version assumes only for G4Box sha     44 //  Current version assumes only for G4Box shape, and the surface
 45 //  is fixed on -Z plane.                          45 //  is fixed on -Z plane.
 46 //                                                 46 //
 47 // Surface is defined at the -Z surface.           47 // Surface is defined at the -Z surface.
 48 // Direction                  -Z   +Z              48 // Direction                  -Z   +Z
 49 //   0  IN || OUT            ->|<-  |              49 //   0  IN || OUT            ->|<-  |
 50 //   1  IN                   ->|    |              50 //   1  IN                   ->|    |
 51 //   2  OUT                    |<-  |              51 //   2  OUT                    |<-  |
 52 //                                                 52 //
 53 // Created: 2005-11-14  Tsukasa ASO, Akinori K     53 // Created: 2005-11-14  Tsukasa ASO, Akinori Kimura.
 54 //                                             <<  54 // 
 55 // 18-Nov-2005  T.Aso,  To use always positive     55 // 18-Nov-2005  T.Aso,  To use always positive value for anglefactor.
 56 // 29-Mar-2007  T.Aso,  Bug fix for momentum d     56 // 29-Mar-2007  T.Aso,  Bug fix for momentum direction at outgoing flux.
 57 // 2010-07-22   Introduce Unit specification.      57 // 2010-07-22   Introduce Unit specification.
 58 // 2010-07-22   Add weighted and divideByAre o     58 // 2010-07-22   Add weighted and divideByAre options
 59 // 2020-10-06   Use G4VPrimitivePlotter and fi     59 // 2020-10-06   Use G4VPrimitivePlotter and fill 1-D histo of kinetic energy (x)
 60 //              vs. cell flux * track weight       60 //              vs. cell flux * track weight             (Makoto Asai)
 61 //////////////////////////////////////////////     61 ///////////////////////////////////////////////////////////////////////////////
 62                                                    62 
 63 G4PSFlatSurfaceFlux::G4PSFlatSurfaceFlux(const <<  63 G4PSFlatSurfaceFlux::G4PSFlatSurfaceFlux(G4String name, 
 64                                          G4int <<  64            G4int direction, G4int depth)
 65   : G4PSFlatSurfaceFlux(name, direction, "perc <<  65   : G4VPrimitivePlotter(name,depth),HCID(-1),fDirection(direction),EvtMap(0),
 66 {}                                             <<  66     weighted(true),divideByArea(true)
 67                                                << 
 68 G4PSFlatSurfaceFlux::G4PSFlatSurfaceFlux(const << 
 69                                          const << 
 70   : G4VPrimitivePlotter(name, depth)           << 
 71   , HCID(-1)                                   << 
 72   , fDirection(direction)                      << 
 73   , EvtMap(nullptr)                            << 
 74   , weighted(true)                             << 
 75   , divideByArea(true)                         << 
 76 {                                                  67 {
 77   DefineUnitAndCategory();                     <<  68     DefineUnitAndCategory();
 78   SetUnit(unit);                               <<  69     SetUnit("percm2");
 79 }                                                  70 }
 80                                                    71 
 81 G4bool G4PSFlatSurfaceFlux::ProcessHits(G4Step <<  72 G4PSFlatSurfaceFlux::G4PSFlatSurfaceFlux(G4String name, 
                                                   >>  73            G4int direction, 
                                                   >>  74            const G4String& unit,
                                                   >>  75            G4int depth)
                                                   >>  76   : G4VPrimitivePlotter(name,depth),HCID(-1),fDirection(direction),EvtMap(0),
                                                   >>  77     weighted(true),divideByArea(true)
 82 {                                                  78 {
 83   G4StepPoint* preStep             = aStep->Ge <<  79     DefineUnitAndCategory();
 84   G4VPhysicalVolume* physVol       = preStep-> <<  80     SetUnit(unit);
                                                   >>  81 }
                                                   >>  82 
                                                   >>  83 G4PSFlatSurfaceFlux::~G4PSFlatSurfaceFlux()
                                                   >>  84 {;}
                                                   >>  85 
                                                   >>  86 G4bool G4PSFlatSurfaceFlux::ProcessHits(G4Step* aStep,G4TouchableHistory*)
                                                   >>  87 {
                                                   >>  88   G4StepPoint* preStep = aStep->GetPreStepPoint();
                                                   >>  89   G4VPhysicalVolume* physVol = preStep->GetPhysicalVolume();
 85   G4VPVParameterisation* physParam = physVol->     90   G4VPVParameterisation* physParam = physVol->GetParameterisation();
 86   G4VSolid* solid                  = nullptr;  <<  91   G4VSolid * solid = 0;
 87   if(physParam != nullptr)                     <<  92   if(physParam)
 88   {  // for parameterized volume               <<  93   { // for parameterized volume
 89     G4int idx =                                <<  94     G4int idx = ((G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable()))
 90       ((G4TouchableHistory*) (aStep->GetPreSte <<  95                 ->GetReplicaNumber(indexDepth);
 91         ->GetReplicaNumber(indexDepth);        << 
 92     solid = physParam->ComputeSolid(idx, physV     96     solid = physParam->ComputeSolid(idx, physVol);
 93     solid->ComputeDimensions(physParam, idx, p <<  97     solid->ComputeDimensions(physParam,idx,physVol);
 94   }                                                98   }
 95   else                                             99   else
 96   {  // for ordinary volume                    << 100   { // for ordinary volume
 97     solid = physVol->GetLogicalVolume()->GetSo    101     solid = physVol->GetLogicalVolume()->GetSolid();
 98   }                                               102   }
 99                                                   103 
100   auto  boxSolid = (G4Box*) (solid);           << 104   G4Box* boxSolid = (G4Box*)(solid);
101                                                   105 
102   G4int dirFlag = IsSelectedSurface(aStep, box << 106   G4int dirFlag =IsSelectedSurface(aStep,boxSolid);
103   if(dirFlag > 0)                              << 107   if ( dirFlag > 0 ) {
104   {                                            << 108     if ( fDirection == fFlux_InOut || fDirection == dirFlag ){
105     if(fDirection == fFlux_InOut || fDirection << 109 
106     {                                          << 110       G4StepPoint* thisStep=0;
107       G4StepPoint* thisStep = nullptr;         << 111       if ( dirFlag == fFlux_In ){
108       if(dirFlag == fFlux_In)                  << 112   thisStep = preStep;
109       {                                        << 113       }else if ( dirFlag == fFlux_Out ){
110         thisStep = preStep;                    << 114   thisStep = aStep->GetPostStepPoint();
111       }                                        << 115       }else{
112       else if(dirFlag == fFlux_Out)            << 116   return FALSE;
113       {                                        << 
114         thisStep = aStep->GetPostStepPoint();  << 
115       }                                        << 
116       else                                     << 
117       {                                        << 
118         return false;                          << 
119       }                                           117       }
120                                                   118 
121       G4TouchableHandle theTouchable = thisSte    119       G4TouchableHandle theTouchable = thisStep->GetTouchableHandle();
122       G4ThreeVector pdirection       = thisSte << 120       G4ThreeVector pdirection = thisStep->GetMomentumDirection();
123       G4ThreeVector localdir =                 << 121       G4ThreeVector localdir  = 
124         theTouchable->GetHistory()->GetTopTran << 122   theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection);
125       //                                          123       //
126       G4double angleFactor = localdir.z();        124       G4double angleFactor = localdir.z();
127       if(angleFactor < 0)                      << 125       if ( angleFactor < 0 ) angleFactor *= -1.;
128         angleFactor *= -1.;                    << 
129       G4double flux = 1.0;                        126       G4double flux = 1.0;
130       if(weighted)                             << 127       if ( weighted ) flux *=preStep->GetWeight(); // Current (Particle Weight)
131         flux *= preStep->GetWeight();  // Curr << 
132       //                                          128       //
133       G4double square =                        << 129       G4double square = 4.*boxSolid->GetXHalfLength()*boxSolid->GetYHalfLength();
134         4. * boxSolid->GetXHalfLength() * boxS << 
135       //                                          130       //
136       flux = flux / angleFactor;  // Flux with << 131       flux = flux/angleFactor;  // Flux with angle.
137       if(divideByArea)                         << 132       if ( divideByArea ) flux /= square;
138         flux /= square;                        << 
139       //                                          133       //
140       G4int index = GetIndex(aStep);              134       G4int index = GetIndex(aStep);
141       EvtMap->add(index, flux);                << 135       EvtMap->add(index,flux);
142                                                   136 
143       if(!hitIDMap.empty() && hitIDMap.find(in << 137       if(hitIDMap.size()>0 && hitIDMap.find(index)!=hitIDMap.end())
144       {                                           138       {
145         auto filler = G4VScoreHistFiller::Inst    139         auto filler = G4VScoreHistFiller::Instance();
146         if(filler == nullptr)                  << 140         if(!filler)
147         {                                         141         {
148           G4Exception("G4PSFlatSurfaceFlux::Pr << 142           G4Exception("G4PSFlatSurfaceFlux::ProcessHits","SCORER0123",JustWarning,
149                       JustWarning,             << 143              "G4TScoreHistFiller is not instantiated!! Histogram is not filled.");
150                       "G4TScoreHistFiller is n << 
151                       "not filled.");          << 
152         }                                         144         }
153         else                                      145         else
154         {                                         146         {
155           filler->FillH1(hitIDMap[index], preS << 147           filler->FillH1(hitIDMap[index],preStep->GetKineticEnergy(),flux);
156         }                                         148         }
157       }                                           149       }
                                                   >> 150 
158     }                                             151     }
159   }                                               152   }
160 #ifdef debug                                      153 #ifdef debug
161   G4cout << " PASSED vol " << index << " trk " << 154     G4cout << " PASSED vol " 
162          << fFlatSurfaceFlux << G4endl;        << 155      << index << " trk "<<trkid<<" len " << fFlatSurfaceFlux<<G4endl;
163 #endif                                            156 #endif
164                                                   157 
165   return true;                                 << 158   return TRUE;
166 }                                                 159 }
167                                                   160 
168 G4int G4PSFlatSurfaceFlux::IsSelectedSurface(G << 161 G4int G4PSFlatSurfaceFlux::IsSelectedSurface(G4Step* aStep, G4Box* boxSolid){
169 {                                              << 162 
170   G4TouchableHandle theTouchable =             << 163   G4TouchableHandle theTouchable = 
171     aStep->GetPreStepPoint()->GetTouchableHand    164     aStep->GetPreStepPoint()->GetTouchableHandle();
172   G4double kCarTolerance =                     << 165   G4double kCarTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
173     G4GeometryTolerance::GetInstance()->GetSur << 
174                                                   166 
175   if(aStep->GetPreStepPoint()->GetStepStatus() << 167   if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
176   {                                            << 
177     // Entering Geometry                          168     // Entering Geometry
178     G4ThreeVector stppos1 = aStep->GetPreStepP << 169     G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
179     G4ThreeVector localpos1 =                  << 170     G4ThreeVector localpos1 = 
180       theTouchable->GetHistory()->GetTopTransf    171       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
181     if(std::fabs(localpos1.z() + boxSolid->Get << 172     if(std::fabs( localpos1.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
182     {                                          << 
183       return fFlux_In;                            173       return fFlux_In;
184     }                                             174     }
185   }                                               175   }
186                                                   176 
187   if(aStep->GetPostStepPoint()->GetStepStatus( << 177   if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
188   {                                            << 
189     // Exiting Geometry                           178     // Exiting Geometry
190     G4ThreeVector stppos2 = aStep->GetPostStep << 179     G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
191     G4ThreeVector localpos2 =                  << 180     G4ThreeVector localpos2 = 
192       theTouchable->GetHistory()->GetTopTransf    181       theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
193     if(std::fabs(localpos2.z() + boxSolid->Get << 182     if(std::fabs( localpos2.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
194     {                                          << 
195       return fFlux_Out;                           183       return fFlux_Out;
196     }                                             184     }
197   }                                               185   }
198                                                   186 
199   return -1;                                      187   return -1;
200 }                                                 188 }
201                                                   189 
202 void G4PSFlatSurfaceFlux::Initialize(G4HCofThi    190 void G4PSFlatSurfaceFlux::Initialize(G4HCofThisEvent* HCE)
203 {                                                 191 {
204   EvtMap = new G4THitsMap<G4double>(GetMultiFu    192   EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(),
205                                     GetName()) << 193             GetName());
206   if(HCID < 0)                                 << 194   if ( HCID < 0 ) HCID = GetCollectionID(0);
207     HCID = GetCollectionID(0);                 << 195   HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
208   HCE->AddHitsCollection(HCID, (G4VHitsCollect << 196 }
                                                   >> 197 
                                                   >> 198 void G4PSFlatSurfaceFlux::EndOfEvent(G4HCofThisEvent*)
                                                   >> 199 {;}
                                                   >> 200 
                                                   >> 201 void G4PSFlatSurfaceFlux::clear(){
                                                   >> 202   EvtMap->clear();
209 }                                                 203 }
210                                                   204 
211 void G4PSFlatSurfaceFlux::clear() { EvtMap->cl << 205 void G4PSFlatSurfaceFlux::DrawAll()
                                                   >> 206 {;}
212                                                   207 
213 void G4PSFlatSurfaceFlux::PrintAll()              208 void G4PSFlatSurfaceFlux::PrintAll()
214 {                                                 209 {
215   G4cout << " MultiFunctionalDet  " << detecto    210   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
216   G4cout << " PrimitiveScorer" << GetName() << << 211   G4cout << " PrimitiveScorer" << GetName() <<G4endl; 
217   G4cout << " Number of entries " << EvtMap->e    212   G4cout << " Number of entries " << EvtMap->entries() << G4endl;
218   for(const auto& [copy, flux] : *(EvtMap->Get << 213   std::map<G4int,G4double*>::iterator itr = EvtMap->GetMap()->begin();
219   {                                            << 214   for(; itr != EvtMap->GetMap()->end(); itr++) {
220     G4cout << "  copy no.: " << copy           << 215     G4cout << "  copy no.: " << itr->first
221            << "  flux  : " << *(flux) / GetUni << 216      << "  flux  : " << *(itr->second)/GetUnitValue() 
222            << GetUnit() << "]" << G4endl;      << 217      << " [" << GetUnit() <<"]"
                                                   >> 218      << G4endl;
223   }                                               219   }
224 }                                                 220 }
225                                                   221 
226 void G4PSFlatSurfaceFlux::SetUnit(const G4Stri    222 void G4PSFlatSurfaceFlux::SetUnit(const G4String& unit)
227 {                                                 223 {
228   if(divideByArea)                             << 224     if ( divideByArea ) {
229   {                                            << 225   CheckAndSetUnit(unit,"Per Unit Surface");
230     CheckAndSetUnit(unit, "Per Unit Surface"); << 226     } else {
231   }                                            << 227   if (unit == "" ){
232   else                                         << 228       unitName = unit;
233   {                                            << 229       unitValue = 1.0;
234     if(unit.empty())                           << 230   }else{
235     {                                          << 231       G4String msg = "Invalid unit ["+unit+"] (Current  unit is [" +GetUnit()+"] ) for " + GetName();
236       unitName  = unit;                        << 232       G4Exception("G4PSFlatSurfaceFlux::SetUnit","DetPS0008",JustWarning,msg);
237       unitValue = 1.0;                         << 233   }
238     }                                             234     }
239     else                                       << 
240     {                                          << 
241       G4String msg = "Invalid unit [" + unit + << 
242                      GetUnit() + "] ) for " +  << 
243       G4Exception("G4PSFlatSurfaceFlux::SetUni << 
244                   msg);                        << 
245     }                                          << 
246   }                                            << 
247 }                                                 235 }
248                                                   236 
249 void G4PSFlatSurfaceFlux::DefineUnitAndCategor << 237 void G4PSFlatSurfaceFlux::DefineUnitAndCategory(){
250 {                                              << 238    // Per Unit Surface
251   // Per Unit Surface                          << 239    new G4UnitDefinition("percentimeter2","percm2","Per Unit Surface",(1./cm2));
252   new G4UnitDefinition("percentimeter2", "perc << 240    new G4UnitDefinition("permillimeter2","permm2","Per Unit Surface",(1./mm2));
253                        (1. / cm2));            << 241    new G4UnitDefinition("permeter2","perm2","Per Unit Surface",(1./m2));
254   new G4UnitDefinition("permillimeter2", "perm << 
255                        (1. / mm2));            << 
256   new G4UnitDefinition("permeter2", "perm2", " << 
257 }                                                 242 }
                                                   >> 243 
                                                   >> 244 
258                                                   245