source: trunk/source/digits_hits/scorer/src/G4PSSphereSurfaceCurrent.cc@ 1354

Last change on this file since 1354 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 8.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4PSSphereSurfaceCurrent.cc,v 1.6 2010/07/23 04:35:38 taso Exp $
28// GEANT4 tag $Name: $
29//
30// G4PSSphereSurfaceCurrent
31#include "G4PSSphereSurfaceCurrent.hh"
32#include "G4StepStatus.hh"
33#include "G4Track.hh"
34#include "G4VSolid.hh"
35#include "G4VPhysicalVolume.hh"
36#include "G4VPVParameterisation.hh"
37#include "G4UnitsTable.hh"
38#include "G4GeometryTolerance.hh"
39////////////////////////////////////////////////////////////////////////////////
40// (Description)
41// This is a primitive scorer class for scoring only Surface Current.
42// Current version assumes only for G4Sphere shape.
43//
44// Surface is defined at the inside of sphere.
45// Direction -Rmin +Rmax
46// 0 IN || OUT ->|<- |
47// 1 IN ->| |
48// 2 OUT |<- |
49//
50// Created: 2005-11-14 Tsukasa ASO, Akinori Kimura.
51// 2010-07-22 Introduce Unit specification.
52//
53///////////////////////////////////////////////////////////////////////////////
54
55G4PSSphereSurfaceCurrent::G4PSSphereSurfaceCurrent(G4String name,
56 G4int direction, G4int depth)
57 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),
58 weighted(true),divideByArea(true)
59{
60 DefineUnitAndCategory();
61 SetUnit("percm2");
62}
63
64G4PSSphereSurfaceCurrent::G4PSSphereSurfaceCurrent(G4String name,
65 G4int direction,
66 const G4String& unit,
67 G4int depth)
68 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),
69 weighted(true),divideByArea(true)
70{
71 DefineUnitAndCategory();
72 SetUnit(unit);
73}
74
75G4PSSphereSurfaceCurrent::~G4PSSphereSurfaceCurrent()
76{;}
77
78G4bool G4PSSphereSurfaceCurrent::ProcessHits(G4Step* aStep,G4TouchableHistory*)
79{
80 G4StepPoint* preStep = aStep->GetPreStepPoint();
81 G4VPhysicalVolume* physVol = preStep->GetPhysicalVolume();
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 }
95
96 G4Sphere* sphereSolid = (G4Sphere*)(solid);
97
98 G4int dirFlag =IsSelectedSurface(aStep,sphereSolid);
99 if ( dirFlag > 0 ) {
100 if ( fDirection == fCurrent_InOut || fDirection == dirFlag ){
101 G4double radi = sphereSolid->GetInsideRadius();
102 G4double dph = sphereSolid->GetDeltaPhiAngle()/radian;
103 G4double stth = sphereSolid->GetStartThetaAngle()/radian;
104 G4double enth = stth+sphereSolid->GetDeltaThetaAngle()/radian;
105 G4double current = 1.0;
106 if ( weighted) current = preStep->GetWeight(); // Current (Particle Weight)
107 if ( divideByArea ){
108 G4double square = radi*radi*dph*( -std::cos(enth) + std::cos(stth) );
109 current = current/square; // Current with angle.
110 }
111
112 G4int index = GetIndex(aStep);
113 EvtMap->add(index,current);
114 }
115 }
116
117 return TRUE;
118}
119
120G4int G4PSSphereSurfaceCurrent::IsSelectedSurface(G4Step* aStep, G4Sphere* sphereSolid){
121
122 G4TouchableHandle theTouchable =
123 aStep->GetPreStepPoint()->GetTouchableHandle();
124 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
125
126 if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
127 // Entering Geometry
128 G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
129 G4ThreeVector localpos1 =
130 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
131 G4double localR2 = localpos1.x()*localpos1.x()
132 +localpos1.y()*localpos1.y()
133 +localpos1.z()*localpos1.z();
134 //G4double InsideRadius2 =
135 // sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
136 //if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
137 G4double InsideRadius = sphereSolid->GetInsideRadius();
138 if ( localR2 > (InsideRadius-kCarTolerance)*(InsideRadius-kCarTolerance)
139 &&localR2 < (InsideRadius+kCarTolerance)*(InsideRadius+kCarTolerance)){
140 return fCurrent_In;
141 }
142 }
143
144 if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
145 // Exiting Geometry
146 G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
147 G4ThreeVector localpos2 =
148 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
149 G4double localR2 = localpos2.x()*localpos2.x()
150 +localpos2.y()*localpos2.y()
151 +localpos2.z()*localpos2.z();
152 //G4double InsideRadius2 =
153 // sphereSolid->GetInsideRadius()*sphereSolid->GetInsideRadius();
154 //if(std::fabs( localR2 - InsideRadius2 ) < kCarTolerance ){
155 G4double InsideRadius = sphereSolid->GetInsideRadius();
156 if ( localR2 > (InsideRadius-kCarTolerance)*(InsideRadius-kCarTolerance)
157 &&localR2 < (InsideRadius+kCarTolerance)*(InsideRadius+kCarTolerance)){
158 return fCurrent_Out;
159 }
160 }
161
162 return -1;
163}
164
165void G4PSSphereSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
166{
167 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
168 if ( HCID < 0 ) HCID = GetCollectionID(0);
169 HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
170}
171
172void G4PSSphereSurfaceCurrent::EndOfEvent(G4HCofThisEvent*)
173{;}
174
175void G4PSSphereSurfaceCurrent::clear(){
176 EvtMap->clear();
177}
178
179void G4PSSphereSurfaceCurrent::DrawAll()
180{;}
181
182void G4PSSphereSurfaceCurrent::PrintAll()
183{
184 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
185 G4cout << " PrimitiveScorer " << GetName() <<G4endl;
186 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
187 std::map<G4int,G4double*>::iterator itr = EvtMap->GetMap()->begin();
188 for(; itr != EvtMap->GetMap()->end(); itr++) {
189 G4cout << " copy no.: " << itr->first << " current : " ;
190 if ( divideByArea ) {
191 G4cout << *(itr->second)/GetUnitValue()
192 << " [" <<GetUnit()<<"]";
193 }else {
194 G4cout << *(itr->second) << " [tracks]" ;
195 }
196 G4cout << G4endl;
197 }
198}
199
200
201void G4PSSphereSurfaceCurrent::SetUnit(const G4String& unit)
202{
203 if ( divideByArea ) {
204 CheckAndSetUnit(unit,"Per Unit Surface");
205 } else {
206 if (unit == "" ){
207 unitName = unit;
208 unitValue = 1.0;
209 }else{
210 G4String msg = "Invalid unit ["+unit+"] (Current unit is [" +GetUnit()+"] )";
211 G4Exception(GetName(),"DetScorer0000",JustWarning,msg);
212 }
213 }
214}
215
216void G4PSSphereSurfaceCurrent::DefineUnitAndCategory(){
217 // Per Unit Surface
218 new G4UnitDefinition("percentimeter2","percm2","Per Unit Surface",(1./cm2));
219 new G4UnitDefinition("permillimeter2","permm2","Per Unit Surface",(1./mm2));
220 new G4UnitDefinition("permeter2","perm2","Per Unit Surface",(1./m2));
221}
Note: See TracBrowser for help on using the repository browser.