source: trunk/source/digits_hits/scorer/src/G4PSCylinderSurfaceCurrent.cc@ 1355

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

update ti head
File size: 7.9 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: G4PSCylinderSurfaceCurrent.cc,v 1.5 2010/07/23 04:35:38 taso Exp $
28// GEANT4 tag $Name: $
29//
30// G4PSCylinderSurfaceCurrent
31#include "G4PSCylinderSurfaceCurrent.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 G4Tubs shape.
43//
44// Surface is defined at the inner surface of the tube.
45// Direction R R+dR
46// 0 IN || OUT ->|<- |
47// 1 IN ->| |
48// 2 OUT |<- |
49//
50// Created: 2007-03-21 Tsukasa ASO
51// 2010-07-22 Introduce Unit specification.
52//
53///////////////////////////////////////////////////////////////////////////////
54
55
56G4PSCylinderSurfaceCurrent::G4PSCylinderSurfaceCurrent(G4String name,
57 G4int direction, G4int depth)
58 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),
59 weighted(true),divideByArea(true)
60{
61 DefineUnitAndCategory();
62 SetUnit("percm2");
63}
64
65G4PSCylinderSurfaceCurrent::G4PSCylinderSurfaceCurrent(G4String name,
66 G4int direction,
67 const G4String& unit,
68 G4int depth)
69 :G4VPrimitiveScorer(name,depth),HCID(-1),fDirection(direction),
70 weighted(true),divideByArea(true)
71{
72 DefineUnitAndCategory();
73 SetUnit(unit);
74}
75
76G4PSCylinderSurfaceCurrent::~G4PSCylinderSurfaceCurrent()
77{;}
78
79G4bool G4PSCylinderSurfaceCurrent::ProcessHits(G4Step* aStep,G4TouchableHistory*)
80{
81 G4StepPoint* preStep = aStep->GetPreStepPoint();
82 G4VPhysicalVolume* physVol = preStep->GetPhysicalVolume();
83 G4VPVParameterisation* physParam = physVol->GetParameterisation();
84 G4VSolid * solid = 0;
85 if(physParam)
86 { // for parameterized volume
87 G4int idx = ((G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable()))
88 ->GetReplicaNumber(indexDepth);
89 solid = physParam->ComputeSolid(idx, physVol);
90 solid->ComputeDimensions(physParam,idx,physVol);
91 }
92 else
93 { // for ordinary volume
94 solid = physVol->GetLogicalVolume()->GetSolid();
95 }
96
97 G4Tubs* tubsSolid = (G4Tubs*)(solid);
98
99 G4int dirFlag =IsSelectedSurface(aStep,tubsSolid);
100 G4cout << " pos " << preStep->GetPosition() <<" dirFlag " << G4endl;
101 if ( dirFlag > 0 ) {
102 if ( fDirection == fCurrent_InOut || fDirection == dirFlag ){
103 G4TouchableHandle theTouchable = preStep->GetTouchableHandle();
104 //
105 G4double current = 1.0;
106 if ( weighted ) current = preStep->GetWeight(); // Current (Particle Weight)
107 //
108 if ( divideByArea ){
109 G4double square = 2.*tubsSolid->GetZHalfLength()
110 *tubsSolid->GetInnerRadius()* tubsSolid->GetDeltaPhiAngle()/radian;
111 current = current/square; // Current normalized by Area
112 }
113
114 G4int index = GetIndex(aStep);
115 EvtMap->add(index,current);
116 }
117
118 }
119
120 return TRUE;
121}
122
123G4int G4PSCylinderSurfaceCurrent::IsSelectedSurface(G4Step* aStep, G4Tubs* tubsSolid){
124
125 G4TouchableHandle theTouchable =
126 aStep->GetPreStepPoint()->GetTouchableHandle();
127 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
128
129 if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
130 // Entering Geometry
131 G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
132 G4ThreeVector localpos1 =
133 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
134 if ( std::fabs(localpos1.z()) > tubsSolid->GetZHalfLength() ) return -1;
135 G4double localR2 = localpos1.x()*localpos1.x()+localpos1.y()*localpos1.y();
136 G4double InsideRadius = tubsSolid->GetInnerRadius();
137 if (localR2 > (InsideRadius-kCarTolerance)*(InsideRadius-kCarTolerance)
138 &&localR2 < (InsideRadius+kCarTolerance)*(InsideRadius+kCarTolerance)){
139 return fCurrent_In;
140 }
141 }
142
143 if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
144 // Exiting Geometry
145 G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
146 G4ThreeVector localpos2 =
147 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
148 if ( std::fabs(localpos2.z()) > tubsSolid->GetZHalfLength() ) return -1;
149 G4double localR2 = localpos2.x()*localpos2.x()+localpos2.y()*localpos2.y();
150 G4double InsideRadius = tubsSolid->GetInnerRadius();
151 if (localR2 > (InsideRadius-kCarTolerance)*(InsideRadius-kCarTolerance)
152 &&localR2 < (InsideRadius+kCarTolerance)*(InsideRadius+kCarTolerance)){
153 return fCurrent_Out;
154 }
155 }
156
157 return -1;
158}
159
160void G4PSCylinderSurfaceCurrent::Initialize(G4HCofThisEvent* HCE)
161{
162 EvtMap = new G4THitsMap<G4double>(detector->GetName(), GetName());
163 if ( HCID < 0 ) HCID = GetCollectionID(0);
164 HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
165}
166
167void G4PSCylinderSurfaceCurrent::EndOfEvent(G4HCofThisEvent*)
168{;}
169
170void G4PSCylinderSurfaceCurrent::clear(){
171 EvtMap->clear();
172}
173
174void G4PSCylinderSurfaceCurrent::DrawAll()
175{;}
176
177void G4PSCylinderSurfaceCurrent::PrintAll()
178{
179 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
180 G4cout << " PrimitiveScorer " << GetName() <<G4endl;
181 G4cout << " Number of entries " << EvtMap->entries() << G4endl;
182 std::map<G4int,G4double*>::iterator itr = EvtMap->GetMap()->begin();
183 for(; itr != EvtMap->GetMap()->end(); itr++) {
184 G4cout << " copy no.: " << itr->first
185 << " current : " ;
186 if ( divideByArea ) {
187 G4cout << *(itr->second)/GetUnitValue()
188 << " ["<<GetUnit()<<"]";
189 } else {
190 G4cout << *(itr->second) << " [tracks]";
191 }
192 G4cout << G4endl;
193 }
194}
195
196void G4PSCylinderSurfaceCurrent::SetUnit(const G4String& unit)
197{
198 if ( divideByArea ) {
199 CheckAndSetUnit(unit,"Per Unit Surface");
200 } else {
201 if (unit == "" ){
202 unitName = unit;
203 unitValue = 1.0;
204 }else{
205 G4String msg = "Invalid unit ["+unit+"] (Current unit is [" +GetUnit()+"] )";
206 G4Exception(GetName(),"DetScorer0000",JustWarning,msg);
207 }
208 }
209}
210
211void G4PSCylinderSurfaceCurrent::DefineUnitAndCategory(){
212 // Per Unit Surface
213 new G4UnitDefinition("percentimeter2","percm2","Per Unit Surface",(1./cm2));
214 new G4UnitDefinition("permillimeter2","permm2","Per Unit Surface",(1./mm2));
215 new G4UnitDefinition("permeter2","perm2","Per Unit Surface",(1./m2));
216}
217
Note: See TracBrowser for help on using the repository browser.