source: trunk/source/digits_hits/utils/test/test2/src/Test2PhantomSD.cc@ 1355

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

update to last version 4.9.4
File size: 10.3 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#include "Test2PhantomSD.hh"
28#include "Test2PhantomHit.hh"
29#include "G4VPhysicalVolume.hh"
30#include "G4LogicalVolume.hh"
31#include "G4Track.hh"
32#include "G4Step.hh"
33#include "G4ParticleDefinition.hh"
34#include "G4VTouchable.hh"
35#include "G4TouchableHistory.hh"
36#include "G4ios.hh"
37#include "G4Box.hh"
38
39#include "G4PSDirectionFlag.hh"
40
41
42Test2PhantomSD::Test2PhantomSD(G4String name, G4int segment[3])
43 :G4VSensitiveDetector(name) {
44
45 nSegment[0]=segment[0];
46 nSegment[1]=segment[1];
47 nSegment[2]=segment[2];
48 G4String HCname;
49 collectionName.insert(HCname = "PhantomCollection");
50
51}
52
53Test2PhantomSD::~Test2PhantomSD() {
54 ;
55}
56
57void Test2PhantomSD::Initialize(G4HCofThisEvent *) {
58
59 fPhantomCollection = new Test2PhantomHitsCollection(SensitiveDetectorName,
60 collectionName[0]);
61 verboseLevel = 0;
62
63}
64
65G4bool Test2PhantomSD::ProcessHits(G4Step * aStep, G4TouchableHistory *) {
66
67 G4double edep = aStep->GetTotalEnergyDeposit();
68 G4double trklen = aStep->GetStepLength();
69 if(edep > 0. || trklen > 0.) {
70 edep /= MeV;
71 trklen /= mm;
72 // total energy deposit
73 G4double weight = aStep->GetPreStepPoint()->GetWeight();
74 G4double density = aStep->GetPreStepPoint()->GetMaterial()->GetDensity();
75 G4double volume = GetVolume(aStep);
76 G4double area = GetArea(aStep);
77 G4double charge = aStep->GetPreStepPoint()->GetCharge();
78
79 G4double dose = edep/(density*volume)/gray;
80
81 G4double flatSurfCurr = 0.0;
82 G4double flatSurfFlux= 0.0;
83 G4int dirFlag = IsSelectedSurface(aStep);
84 if ( dirFlag > 0 ){
85 flatSurfCurr = 1.0*weight/area/(1./cm2);
86 G4double anglefac = GetAngleFactor(aStep,dirFlag);
87 flatSurfFlux = weight/anglefac/area/(1./cm2);
88 }
89
90 G4double passageCellCurr = 0.0;
91 G4double passageCellFlux = 0.0;
92 if ( IsPassed(aStep) ){
93 passageCellCurr = weight;
94 passageCellFlux = fCellTrack/volume/(1./cm2); // fCellTrack is calculated in IsPassed().
95 }
96
97 G4double cellFlux = trklen*weight/volume/(1./cm2);
98
99 G4double nOfSecondary=0.0;
100 if ( IsSecondary(aStep) ){
101 nOfSecondary = weight;
102 }
103
104 G4double cellCharge = 0.0;
105 if ( IsEnterOrFirstStep(aStep) ){
106 cellCharge = charge*weight;
107 }else if ( IsExit(aStep) ){
108 cellCharge = -1.*charge*weight;
109 }
110
111 if(verboseLevel > 1) G4cout << "Next step edep [MeV] = " << edep/MeV << G4endl;
112
113 G4TouchableHistory * hist = (G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable());
114 G4int copyIDinX = hist->GetReplicaNumber(2);
115 G4int copyIDinY = hist->GetReplicaNumber(1);
116 G4int copyIDinZ = hist->GetReplicaNumber(0);
117
118 Test2PhantomHit* phantomHit
119 = new Test2PhantomHit(copyIDinX, copyIDinY, copyIDinZ,nSegment);
120 phantomHit->SetEdep(edep);
121 phantomHit->SetTrackLength(trklen);
122 phantomHit->SetParticleName(aStep->GetTrack()->GetParticleDefinition()->GetParticleName());
123 phantomHit->SetDose(dose);
124
125 phantomHit->SetFlatSurfaceCurrent(flatSurfCurr);
126 phantomHit->SetPassageCellCurrent(passageCellCurr);
127 phantomHit->SetFlatSurfaceFlux(flatSurfFlux);
128 phantomHit->SetCellFlux(cellFlux);
129 phantomHit->SetPassageCellFlux(passageCellFlux);
130 phantomHit->SetNofSecondary(nOfSecondary);
131 phantomHit->SetCellCharge(cellCharge);
132
133 fPhantomCollection->insert(phantomHit);
134 if(verboseLevel > 0) {
135 G4cout << " A hit of Test2PhantomHit is created in copy-id ("
136 << copyIDinX << ", " << copyIDinY << ", " << copyIDinZ
137 << ") at " << GetFullPathName() << G4endl;
138 }
139 }
140
141 return true;
142}
143
144void Test2PhantomSD::EndOfEvent(G4HCofThisEvent * HCE) {
145
146 static G4int HCID = -1;
147 if(HCID < 0) { HCID = GetCollectionID(0); }
148 HCE->AddHitsCollection( HCID, fPhantomCollection );
149}
150
151void Test2PhantomSD::clear() {
152 ;
153}
154
155void Test2PhantomSD::DrawAll() {
156 ;
157}
158
159void Test2PhantomSD::PrintAll() {
160 ;
161}
162
163
164G4VSolid* Test2PhantomSD::GetSolid(G4Step* aStep){
165 G4VPhysicalVolume* physVol = aStep->GetPreStepPoint()->GetPhysicalVolume();
166 G4VPVParameterisation* physParam = physVol->GetParameterisation();
167 G4VSolid* solid = 0;
168 G4int indexDepth = 0;
169 if(physParam)
170 { // for parameterized volume
171 G4int idx = ((G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable()))
172 ->GetReplicaNumber(indexDepth);
173 if(idx<0)
174 {
175 G4Exception("G4PSDoseDeposit","G4PSDoseDeposit::ProcessHits",JustWarning,
176 "Incorrect replica number");
177 G4cerr << " --- GetReplicaNumber : " << idx << G4endl;
178 }
179 solid = physParam->ComputeSolid(idx, physVol);
180 solid->ComputeDimensions(physParam,idx,physVol);
181 }
182 else
183 { // for ordinary volume
184 solid = physVol->GetLogicalVolume()->GetSolid();
185 }
186
187 return solid;
188}
189
190G4double Test2PhantomSD::GetVolume(G4Step* aStep){
191 G4VSolid* solid = GetSolid(aStep);
192 return (solid->GetCubicVolume());
193}
194
195G4double Test2PhantomSD::GetArea(G4Step* aStep){
196 G4Box* boxSolid = (G4Box*)(GetSolid(aStep));
197 return (4.*boxSolid->GetXHalfLength()*boxSolid->GetYHalfLength());
198}
199
200G4int Test2PhantomSD::IsSelectedSurface(G4Step* aStep){
201 G4Box* boxSolid = (G4Box*)GetSolid(aStep);
202 G4TouchableHandle theTouchable =
203 aStep->GetPreStepPoint()->GetTouchableHandle();
204 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
205
206 if (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ){
207 // Entering Geometry
208 G4ThreeVector stppos1= aStep->GetPreStepPoint()->GetPosition();
209 G4ThreeVector localpos1 =
210 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos1);
211 if(std::fabs( localpos1.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
212 return fCurrent_In;
213 }
214 }
215
216 if (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary ){
217 // Exiting Geometry
218 G4ThreeVector stppos2= aStep->GetPostStepPoint()->GetPosition();
219 G4ThreeVector localpos2 =
220 theTouchable->GetHistory()->GetTopTransform().TransformPoint(stppos2);
221 if(std::fabs( localpos2.z() + boxSolid->GetZHalfLength())<kCarTolerance ){
222 return fCurrent_Out;
223 }
224 }
225
226 return -1;
227}
228
229G4bool Test2PhantomSD::IsPassed(G4Step* aStep){
230 G4bool Passed = FALSE;
231
232 G4bool IsEnter = aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary;
233 G4bool IsExit = aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary;
234 G4double trklength = aStep->GetStepLength();
235 trklength *= aStep->GetPreStepPoint()->GetWeight();
236 G4int trkid = aStep->GetTrack()->GetTrackID();
237
238 if ( IsEnter &&IsExit ){ // Passed at one step
239 fCellTrack = trklength;
240 Passed = TRUE;
241 }else if ( IsEnter ){ // Enter a new geometry
242 fCurrentTrkID = trkid; // Resetting the current track.
243 fCellTrack = trklength;
244 }else if ( IsExit ){ // Exit a current geometry
245 if ( fCurrentTrkID == trkid ) {
246 fCellTrack += trklength;
247 Passed = TRUE; // if the track is same as entered.
248 }
249 }else{ // Inside geometry
250 if ( fCurrentTrkID == trkid ){ // Adding the track length to current one ,
251 fCellTrack = trklength;
252 }
253 }
254 return Passed;
255}
256
257
258G4double Test2PhantomSD::GetAngleFactor(G4Step* aStep, G4int dirFlag){
259 G4StepPoint* thisStep=0;
260 if ( dirFlag == fFlux_In ){
261 thisStep = aStep->GetPreStepPoint();
262 }else if ( dirFlag == fFlux_Out ){
263 thisStep = aStep->GetPostStepPoint();
264 }else{
265 return FALSE;
266 }
267 G4TouchableHandle theTouchable = thisStep->GetTouchableHandle();
268 G4ThreeVector pdirection = thisStep->GetMomentumDirection();
269 G4ThreeVector localdir =
270 theTouchable->GetHistory()->GetTopTransform().TransformAxis(pdirection);
271 //
272 G4double angleFactor = localdir.z();
273 if ( angleFactor < 0 ) angleFactor *= -1.;
274 return angleFactor;
275}
276
277G4bool Test2PhantomSD::IsSecondary(G4Step* aStep){
278 if ( aStep->GetTrack()->GetCurrentStepNumber() != 1) return FALSE;
279 //- check for this is not a primary particle. e.g. ParentID > 0 .
280 if ( aStep->GetTrack()->GetParentID() == 0 ) return FALSE;
281 //- check the particle if the partifle definition is given.
282
283 return TRUE;
284}
285G4bool Test2PhantomSD::IsEnterOrFirstStep(G4Step* aStep){
286 // Enter or First step of primary.
287 return (aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary
288 || ( aStep->GetTrack()->GetParentID() == 0 &&
289 aStep->GetTrack()->GetCurrentStepNumber() == 1 ) );
290}
291G4bool Test2PhantomSD::IsExit(G4Step* aStep){
292 return (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
293}
Note: See TracBrowser for help on using the repository browser.