source: trunk/examples/advanced/Rich/src/PadHpdPhotoElectricEffect.cc@ 1194

Last change on this file since 1194 was 807, checked in by garnier, 17 years ago

update
File size: 9.7 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// Rich advanced example for Geant4
27// PadHpdPhotoElectricEffect.cc for Rich of LHCb
28// History:
29// Created: Sajan Easo (Sajan.Easo@cern.ch)
30// Revision: Patricia Mendez (Patricia.Mendez@cern.ch)
31/////////////////////////////////////////////////////////////////////////////
32#include "globals.hh"
33#include <cmath>
34
35#include "PadHpdPhotoElectricEffect.hh"
36#include "RichTbGeometryParameters.hh"
37#include "G4TransportationManager.hh"
38#include "G4TouchableHandle.hh"
39#include "G4GeometryTolerance.hh"
40#include "Randomize.hh"
41#include "RichTbAnalysisManager.hh"
42#include "RichTbRunConfig.hh"
43#include "RichTbMaterialParameters.hh"
44
45#include "RichTbAnalysisManager.hh"
46
47PadHpdPhotoElectricEffect::PadHpdPhotoElectricEffect(const G4String& processName ,
48 RichTbRunConfig* RConfig)
49 :G4VDiscreteProcess(processName),
50 DemagnificationFactor(std::vector<G4double>(NumHpdTot)),
51 DemagnificationQuadFactor(std::vector<G4double>(NumHpdTot)),
52 HpdQE(NumHpdTot, std::vector<G4double>( NumQEbins)),
53 HpdWabin(NumHpdTot, std::vector<G4double>( NumQEbins))
54{
55 rConfig=RConfig;
56 PrePhotoElectricVolName="PadHpdWindowQuartz";
57 PostPhotoElectricVolName="BiAlkaliPhCathode";
58 HpdPhElectronKE=(RConfig-> getHpdPhElectronEnergy())*keV;
59 PhCathodeToSilDetDist= HpdPhotoCathodeSiZdist;
60 PSFsigma=PadHpdPSFsigma;
61
62 for(G4int ihpdq=0; ihpdq<NumHpdTot; ihpdq++ ) {
63
64 if( HpdDemagLinearTerm[ihpdq] != 0.0 ) {
65 DemagnificationFactor[ihpdq]=1.0 / HpdDemagLinearTerm[ihpdq];
66 }else { DemagnificationFactor[ihpdq] =1.0 / 2.3;}
67
68 if( HpdDemagQuadraticTerm[ihpdq] != 0.0 ) {
69
70 DemagnificationQuadFactor[ihpdq]=1.0 / HpdDemagQuadraticTerm[ihpdq];
71
72 }else{DemagnificationQuadFactor[ihpdq]=0.0*(1.0/(1.0*mm)); }
73
74
75 // Now to apply the error on the HPD demag factor. SE 28-4-02
76 // for now a flat error is applied. the uniform number from -1.0 to
77 // 1.0 is obtained and then multiplied with the factor.
78
79 G4double DemagError= (HpdDemagErrorPercent/100.0)*(2.0*G4UniformRand()-1.0) ;
80
81 DemagnificationFactor[ihpdq] = DemagnificationFactor[ihpdq]*(1.0+DemagError);
82
83
84 std::vector<G4double>qeCurHpd = InitializeHpdQE(ihpdq);
85 std::vector<G4double>waCurHpd = InitializeHpdWaveL(ihpdq);
86 if(qeCurHpd.size() != waCurHpd.size() ) {
87 G4cout<<"Wrong size for Hpd QE "<<ihpdq<<" "<<qeCurHpd.size()
88 <<" "<< waCurHpd.size()<<G4endl;
89 }
90 for(size_t iqbin=0; iqbin < qeCurHpd.size(); iqbin++){
91 HpdQE[ihpdq][iqbin]=qeCurHpd[iqbin]/100;
92 HpdWabin[ihpdq][iqbin]=waCurHpd[iqbin];
93 }
94 }
95 G4cout<<GetProcessName() <<" is created "<<G4endl;
96
97
98 }
99PadHpdPhotoElectricEffect::~PadHpdPhotoElectricEffect() {; }
100
101G4bool PadHpdPhotoElectricEffect::IsApplicable(const G4ParticleDefinition&
102 aParticleType)
103{
104 return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
105}
106
107G4double PadHpdPhotoElectricEffect::GetMeanFreePath(const G4Track& ,
108 G4double ,
109 G4ForceCondition* condition)
110{
111 *condition = Forced;
112
113 return DBL_MAX;
114}
115
116G4VParticleChange* PadHpdPhotoElectricEffect::PostStepDoIt(const G4Track& aTrack,
117 const G4Step& aStep)
118{
119
120 aParticleChange.Initialize(aTrack);
121
122 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
123 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
124
125
126
127 if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
128
129 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
130
131 }
132
133
134 G4String PrePhName = pPreStepPoint -> GetPhysicalVolume() ->
135 GetLogicalVolume() -> GetMaterial()->GetName();
136 G4String PostPhName= pPostStepPoint -> GetPhysicalVolume() ->
137 GetLogicalVolume() -> GetMaterial() ->GetName();
138
139
140 if(( PrePhName == PrePhotoElectricVolName &&
141 PostPhName == PostPhotoElectricVolName) ||
142 ( PostPhName == PrePhotoElectricVolName &&
143 PrePhName == PostPhotoElectricVolName) ) {
144
145
146 }else {
147
148 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
149
150}
151
152 G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
153 if (aTrack.GetStepLength()<=kCarTolerance/2){
154 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
155 }
156
157 const G4DynamicParticle* aDynamicPhoton = aTrack.GetDynamicParticle();
158 G4double PhotonEnergy = aDynamicPhoton->GetKineticEnergy();
159
160 if(PhotonEnergy <= 0.0 ) {
161
162 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
163 }
164
165
166 //Now use the QE for the current HPD to determine if a
167 // photoelectron should be produced or not.
168
169 G4int currentHpdNumber= pPreStepPoint->GetTouchableHandle()
170 -> GetReplicaNumber(1);
171 if(currentHpdNumber >= NumHpdTot ){
172
173 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
174
175 }
176
177 G4double PhotWLength=PhotMomWaveConv/(PhotonEnergy/eV);
178
179 G4double PhCathodeQE = getHpdQEff(currentHpdNumber, PhotWLength);
180 G4double randomnum = G4UniformRand();
181
182
183 //the following three lines are copied from few lines later just
184 // for histogramming convenience.
185 G4ThreeVector GlobalElectronOrigin= pPostStepPoint->GetPosition();
186
187 G4Navigator* theNavigator =
188 G4TransportationManager::GetTransportationManager()->
189 GetNavigatorForTracking();
190 G4ThreeVector LocalElectronOrigin = theNavigator->
191 GetGlobalToLocalTransform().
192 TransformPoint(GlobalElectronOrigin);
193
194
195
196 // For the histogram of the radius of the cherenkov circle.
197 // This assumes that the beam is along 001 axis in the global
198 // coord system.
199 G4double GLx=GlobalElectronOrigin.x();
200 G4double GLy=GlobalElectronOrigin.y();
201 // G4double PhotCkvRad = std::pow((std::pow(GLx,2)+std::pow(GLy,2)),0.5);
202 G4double PhotCkvPhi = std::atan2(GLy,GLx)*180.0/pi;
203
204 if( PhotCkvPhi < - 180.0 )PhotCkvPhi+= 360.0;
205
206 if(randomnum < PhCathodeQE ) {
207
208
209 G4double CurDemagFactor=DemagnificationFactor[currentHpdNumber];
210 G4double CurDemagQuadFactor=DemagnificationQuadFactor[currentHpdNumber];
211
212 // now get the Point Spread function.
213
214 G4double PsfRandomAzimuth = twopi*G4UniformRand();
215 G4double PsfRandomRad= G4RandGauss::shoot(0.0,PSFsigma);
216 G4double PsfX= PsfRandomRad*std::cos( PsfRandomAzimuth);
217 G4double PsfY= PsfRandomRad*std::sin( PsfRandomAzimuth);
218
219
220 G4ThreeVector LocalElectronDirection(
221 (CurDemagFactor+CurDemagQuadFactor*LocalElectronOrigin.x()-1.0)*LocalElectronOrigin.x()+PsfX,
222 (CurDemagFactor+CurDemagQuadFactor*LocalElectronOrigin.y()-1.0)*LocalElectronOrigin.y()+PsfY,
223 -(PhCathodeToSilDetDist-
224 (HpdPhCathodeRInner-LocalElectronOrigin.z())));
225 //normalize this vector and then transform back to global coord system.
226 LocalElectronDirection = LocalElectronDirection.unit();
227
228 const G4ThreeVector GlobalElectronDirection = theNavigator->
229 GetLocalToGlobalTransform().
230 TransformAxis(LocalElectronDirection);
231
232 G4double ElecKineEnergy=getHpdPhElectronKE();
233
234 //create the electron
235 G4DynamicParticle* aElectron= new G4DynamicParticle (G4Electron::Electron(),
236 GlobalElectronDirection, ElecKineEnergy) ;
237
238 aParticleChange.SetNumberOfSecondaries(1) ;
239 // aParticleChange.AddSecondary( aElectron ) ;
240 aParticleChange.AddSecondary( aElectron,GlobalElectronOrigin,true ) ;
241
242
243 // Kill the incident photon when it has converted to photoelectron.
244
245 aParticleChange.ProposeLocalEnergyDeposit(PhotonEnergy);
246 aParticleChange.ProposeEnergy(0.);
247 aParticleChange.ProposeTrackStatus(fStopAndKill);
248 }
249 //photon is not killed if it is not converted to photoelectron
250 //SE 26-09-01.
251 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
252
253}
254
255G4double PadHpdPhotoElectricEffect::getHpdQEff(G4int HpdNum,
256 G4double PhotonWLength){
257
258 G4double hq1,hq2, wa1, wa2,aslope,aintc;
259 G4double qeff=0.0;
260 for (G4int ibinq=0 ; ibinq<NumQEbins-1 ; ibinq++ ){
261 wa1 = HpdWabin[HpdNum][ibinq];
262 wa2 = HpdWabin[HpdNum][ibinq+1];
263 if( PhotonWLength >= wa1 && PhotonWLength <= wa2 ) {
264 hq1 = HpdQE[HpdNum][ibinq];
265 hq2 = HpdQE[HpdNum][ibinq+1];
266 aslope = (hq2-hq1)/(wa2-wa1);
267 aintc = hq1 - (aslope * wa1 );
268 qeff= aintc + aslope * PhotonWLength ;
269 return qeff;
270 }
271
272 }
273 return qeff;
274}
275
276
277
278
Note: See TracBrowser for help on using the repository browser.