1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
/**
* Copyright (C) 2009
* The IceCube collaboration
* ID: $Id: G4TankIceSD.cxx 152814 2017-01-19 21:34:52Z jgonzalez $
*
* @file G4TankIceSD.cxx
* @version $Rev: 152814 $
* @date $Date: 2017-01-19 21:34:52 +0000 (Thu, 19 Jan 2017) $
* @author Tilo Waldenmaier, Thomas Melzig
*/
#include <g4-tankresponse/g4classes/G4TankIceSD.h>
#include <g4-tankresponse/g4classes/G4BeamTestTank.h>
#include <G4Step.hh>
#include <G4HCofThisEvent.hh>
#include <G4TouchableHistory.hh>
#include <G4ios.hh>
#include <iostream>
#include <ios>
#include <fstream>
#include <G4VProcess.hh>
#include <G4OpticalPhoton.hh>
#include "G4Poisson.hh"
#include <icetray/I3Units.h>
#include "G4Version.hh"
#if G4VERSION_NUMBER >= 950
// The G4MaterialPropertyVector is gone since 4.9.5.
// It has been typedef'd to a G4UnorderedPhysicsVector
// with a different interface. Try to support both
// versions with an ifdef.
#define MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
#endif
G4TankIceSD::G4TankIceSD(G4String name, const std::map<OMKey, G4ThreeVector>& domPositions)
: G4VSensitiveDetector(name), domPositions_(domPositions)
{}
G4TankIceSD::~G4TankIceSD() {}
void G4TankIceSD::Initialize(G4HCofThisEvent* HCE)
{
std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
{
sumEdep_[om_iter->first] = 0;
cogTime_[om_iter->first] = 0;
cherenkovCounter_[om_iter->first] = 0;
cherenkovCounterWeight_[om_iter->first] = 0;
}
}
G4bool G4TankIceSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
{
// Get energy deposit/time --> CHECK WHAT THESE QUANTITIES ACTUALLY MEAN !!!!!!!!!!!!!!!!!!!!!!!!
G4double edep = aStep->GetTotalEnergyDeposit();
G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
G4double cherenkovNumber = GetCerenkovNumber(aStep);
if(edep<=0 && cherenkovNumber<=0) return false;
// Get Position relative to ice center
G4ThreeVector preStepPoint = aStep->GetPreStepPoint()->GetPosition();
G4ThreeVector postStepPoint = aStep->GetPostStepPoint()->GetPosition();
G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle();
G4ThreeVector worldPosition = (preStepPoint + postStepPoint) / 2.0;
G4ThreeVector localPosition = theTouchable->GetHistory()->
GetTopTransform().TransformPoint(worldPosition);
std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
{
G4double radius = sqrt(pow(om_iter->second.x() - localPosition.x(), 2) +
pow(om_iter->second.y() - localPosition.y(), 2));
G4double height = (om_iter->second.z() - localPosition.z());
sumEdep_[om_iter->first] += edep;
cogTime_[om_iter->first] += edep*time;
cherenkovCounterWeight_[om_iter->first] += GetProbability(radius, height) * cherenkovNumber;
cherenkovCounter_[om_iter->first] += cherenkovNumber;
}
return true;
}
void G4TankIceSD::EndOfEvent(G4HCofThisEvent*)
{
std::map<OMKey, G4ThreeVector>::const_iterator om_iter;
for(om_iter=domPositions_.begin(); om_iter!=domPositions_.end(); ++om_iter)
{
if(sumEdep_[om_iter->first]>0)
{
cogTime_[om_iter->first] /= sumEdep_[om_iter->first];
}
}
}
G4double G4TankIceSD::GetCerenkovNumber(G4Step* aStep)
{
// same function as in "source/processes/electromagnetic/xrays/src/G4Cerenkov.cc" but without
// cerenkov angle and only for an energy independent refraction index
G4Track* aTrack = aStep->GetTrack();
const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
const G4Material* aMaterial = aTrack->GetMaterial();
G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
if (!aMaterialPropertiesTable) return 0.0;
G4MaterialPropertyVector* Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
if (!Rindex) return 0.0;
const G4double Rfact = 369.81 / (CLHEP::eV * CLHEP::cm);
const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
if (beta <= 0.0) return 0.0;
G4double BetaInverse = 1. / beta;
// Min and Max photon energies
#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
G4double Pmin = Rindex->GetMinLowEdgeEnergy();
G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
#else
G4double Pmin = Rindex->GetMinPhotonEnergy();
G4double Pmax = Rindex->GetMaxPhotonEnergy();
#endif
// Min and Max Refraction Indices
#ifdef MATERIAL_PROPERTY_VECTOR_IS_PHYSICS_VECTOR
G4double nMin = Rindex->GetMinValue();
G4double nMax = Rindex->GetMaxValue();
#else
G4double nMin = Rindex->GetMinProperty();
G4double nMax = Rindex->GetMaxProperty();
#endif
if (nMin!=nMax)
{
log_error("Support for energy dependent refraction index not yet implemented!");
return 0.0;
}
// If n(Pmax) < 1/Beta -- no photons generated
if (nMax < BetaInverse) return 0.0;
G4double MeanNumberOfPhotons = Rfact * charge/CLHEP::eplus * charge/CLHEP::eplus * (Pmax - Pmin) *
(1. - BetaInverse * BetaInverse / nMin / nMin);
G4double step_length = aStep->GetStepLength();
return MeanNumberOfPhotons * step_length;
}
G4double G4TankIceSD::GetProbability(G4double radius, G4double height)
{
height = 0.90 - height / CLHEP::m;
radius = radius / CLHEP::m;
G4double p0 = 0.0340 + 0.0120 * height;
G4double p1 = 0.0400 + 0.000622 * exp(-1.13 + 8.03 * height);
G4double p2 = 0.917 + 2.34 * exp(-1.82 + 3.12 * height);
G4double p3 = -0.00325 - 0.00199 * height - 0.0121 * height*height;
return (p0 + p3 * radius) + (p1 - p0) * exp(-p2*p2 * radius*radius);
}
|
