source: trunk/examples/advanced/Rich/src/RichTbMaterialParameters.cc@ 1277

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

update
File size: 7.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// RichTbMaterialParameters.cc for Rich of LHCb
28// History:
29// Created: Sajan Easo (Sajan.Easo@cern.ch)
30// Revision and changes: Patricia Mendez (Patricia.Mendez@cern.ch)
31/////////////////////////////////////////////////////////////////////////////
32#include <iostream>
33#include <fstream>
34#include "globals.hh"
35#include "RichTbGeometryParameters.hh"
36#include "RichTbMaterialParameters.hh"
37#include "FilterTrData.hh"
38#include "AerogelTypeSpec.hh"
39
40#include "RichTbAnalysisManager.hh"
41
42
43void InitializeRichTbMaterial(){
44
45}
46
47std::vector<G4double> InitializePhotonMomentumVector() {
48
49 G4double PhotonEnergyStep=(PhotonMaxEnergy-PhotonMinEnergy)/
50 NumPhotWaveLengthBins;
51 std::vector<G4double>PhotMomVect(NumPhotWaveLengthBins);
52 for (G4int ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
53 PhotMomVect[ibin]=PhotonMinEnergy+PhotonEnergyStep*ibin;
54 }
55 return PhotMomVect;
56}
57std::vector<G4double> InitN2RefIndex(G4double pressure, G4double temperature){
58
59 std::vector<G4double> PmV=InitN2RefPhotW();
60 std::vector<G4double> RefN2(NumPhotWaveLengthBins);
61 G4double GasRhoN2Cur=GasRhoN2atSTP*(GasTemperature_STP/temperature)*
62 (pressure/ GasPressure_STP);
63 G4double epho,pfe,cpfe;
64 for(G4int ibinwn =0; ibinwn<NumPhotWaveLengthBins ; ibinwn++ ){
65
66 epho = PmV[ibinwn]/eV;
67 pfe = SellN2F1/(SellN2E1*SellN2E1 - epho*epho ) +
68 SellN2F2/(SellN2E2*SellN2E2 - epho*epho );
69 cpfe=0.3738*(GasRhoN2Cur/GasMolWeightN2)*pfe;
70 RefN2[ibinwn]=std::pow((1.0+2*cpfe)/(1.0-cpfe),0.5);
71 }
72 return RefN2;
73}
74std::vector<G4double> InitN2RefPhotW() {
75 return InitializePhotonMomentumVector() ;
76}
77std::vector<G4double> InitAgelPhotW() {
78 return InitializePhotonMomentumVector() ;
79}
80std::vector<G4double> InitializeHpdQE(G4int ihpdqe) {
81 // Initialize the HPD QE
82 G4int iqb;
83 if(ihpdqe >= NumHpdTot ) {
84 G4cout<<"Wrong HPD Number for QE " <<ihpdqe<<" vs "
85 <<NumHpdTot <<G4endl;
86 }
87 std::vector<G4double>qeCurPerCent(NumQEbins);
88 if(ihpdqe == 0 ){
89 for(iqb=0; iqb<NumQEbins; iqb++){
90 qeCurPerCent[iqb] = Hpd0QEPerCent[iqb]* HpdQEReductionFactor;
91 }
92 }
93 if(ihpdqe == 1 ){
94 for(iqb=0; iqb<NumQEbins; iqb++){
95 qeCurPerCent[iqb] = Hpd1QEPerCent[iqb]* HpdQEReductionFactor;
96 }
97 }
98 if(ihpdqe == 2 ){
99 for(iqb=0; iqb<NumQEbins; iqb++){
100 qeCurPerCent[iqb] = Hpd2QEPerCent[iqb]* HpdQEReductionFactor;
101 }
102 }
103 if(ihpdqe == 3 ){
104 for(iqb=0; iqb<NumQEbins; iqb++){
105 qeCurPerCent[iqb] = Hpd3QEPerCent[iqb]* HpdQEReductionFactor;
106 }
107 }
108
109 return qeCurPerCent;
110}
111std::vector<G4double> InitializeHpdWaveL(G4int ihpdqe) {
112 G4int iqb;
113 if(ihpdqe >= NumHpdTot ) {
114 G4cout<<"Wrong HPD Number for QE wavelength " <<ihpdqe<<" vs "
115 <<NumHpdTot <<G4endl;
116 }
117 // for now all HPDs have the same wavelength bins.
118 std::vector<G4double>HpdQEW(NumQEbins);
119 for (iqb=0; iqb<NumQEbins; iqb++){
120 HpdQEW[iqb]= HpdQEWaveL[iqb];
121 }
122 return HpdQEW;
123}
124
125void HistoRichTbMaterialProperties(RichTbRunConfig* RConfig) {
126
127
128 G4int AerogelNum=0;
129 G4double waL=200;
130
131 G4double stepsize=7.0;
132
133 // G4double thickness=(GetCurAerogelLength(AerogelNum))/cm;
134 AerogelType CurAerogelType=RConfig-> GetCurAerogelType(AerogelNum);
135
136
137
138 G4double Aparam=0.;
139 G4double Cparam=0.;
140 if(CurAerogelType == AerogelTypeA ) {
141 Aparam = AerogelTypeATotTrans;
142 Cparam = AerogelTypeAClarity*cm/(micrometer*micrometer*micrometer*micrometer);
143 }
144
145
146
147 for(G4int Iabin=0; Iabin<100; Iabin ++ ) {
148
149 // G4double waLInmu = waL/1000.0;
150 // G4double Aetr = Aparam* std::exp(-Cparam * thickness / std::pow(waLInmu,4) );
151
152 waL += stepsize;
153 }
154
155
156
157
158
159
160 G4int ihpdqa;
161 ihpdqa=0;
162 std::vector<G4double>WaveL1 = InitializeHpdWaveL(ihpdqa);
163 std::vector<G4double>QEff1 = InitializeHpdQE(ihpdqa);
164
165
166
167}
168
169
170std::vector<G4int> getDeadPixelList(G4int ihpdNum, G4int){
171 std::vector<G4int>DeadPixelList;
172 // G4int isc,ipsc;
173
174
175 if(G4int(DeadPixelList.size()) > MaxNumDeadPixelPerHpdSect ){
176 G4cout<<" Too Many dead Pixels in Hpd "<<DeadPixelList.size()
177 <<" in Hpd "<<ihpdNum<<G4endl;
178 }
179
180 return DeadPixelList;
181}
182std::vector<G4double>GetAerogelRScatLength(AerogelType CurrentAerogelType) {
183
184 std::vector<G4double>AgelRayleighScatLength(NumPhotWaveLengthBins);
185 std::vector<G4double>AgelPhotW = InitAgelPhotW();
186 G4double aClarity=0.;
187 if(CurrentAerogelType == AerogelTypeA ) {
188 aClarity=AerogelTypeAClarity/(micrometer*micrometer*micrometer*micrometer);
189
190 }else if (CurrentAerogelType == AerogelTypeB ) {
191 aClarity=AerogelTypeBClarity/(micrometer*micrometer*micrometer*micrometer);
192 }else if (CurrentAerogelType == AerogelTypeC ) {
193 aClarity=AerogelTypeCClarity/(micrometer*micrometer*micrometer*micrometer);
194
195
196 }else if (CurrentAerogelType == AerogelTypeD ) {
197 aClarity=AerogelTypeDClarity/(micrometer*micrometer*micrometer*micrometer);
198
199 }else if (CurrentAerogelType == AerogelTypeE ) {
200 aClarity=AerogelTypeEClarity/(micrometer*micrometer*micrometer*micrometer);
201
202 }else {G4cout<<"Unknown Aerogel Type for Rayleigh Scat Length "<<G4endl; }
203
204 if(aClarity != 0.0 ) {
205 for(G4int ibinw=0; ibinw<NumPhotWaveLengthBins; ibinw++ ){
206 G4double ephoton=AgelPhotW[ibinw]/eV;
207 //In the following the 1000 is to convert form nm to micrometer
208 G4double wphoton=(PhotMomWaveConv/ephoton)/1000.0;
209 AgelRayleighScatLength[ibinw]=(std::pow(wphoton,4))/aClarity;
210
211 }
212 }
213
214 return AgelRayleighScatLength;
215}
216G4double GetCurrentBulkTrans(G4double currentMatRefIndex,
217 G4double currentNeighbourRefIndex,
218 G4double MaxTotMeasuredTransmission){
219 G4double ATrans=MaxTotMeasuredTransmission;
220 // G4double ePhot;
221 // in the following the energy of the photon is not used since
222 // it is only an approximate calulation.
223 G4double na= currentMatRefIndex;
224 G4double nb= currentNeighbourRefIndex;
225 G4double LossAtEntrance=std::pow(((na-nb)/(na+nb)),2.0);
226 G4double LossAtExit=std::pow(((nb-na)/(nb+na)),2.0);
227
228 G4double LightLossAtExternalSurface= LossAtEntrance+ LossAtExit;
229
230 ATrans += LightLossAtExternalSurface;
231 if(ATrans >= 1.0) ATrans=1.0;
232 return ATrans;
233}
234
235
236
237
238
Note: See TracBrowser for help on using the repository browser.