source: trunk/examples/advanced/Rich/src/RichTbMaterial.cc@ 1341

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

update
File size: 41.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// Rich advanced example for Geant4
27// RichTbMaterial.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 <cmath>
34#include "globals.hh"
35#include "RichTbMaterial.hh"
36#include "G4Isotope.hh"
37#include "G4Element.hh"
38#include "G4ElementTable.hh"
39#include "G4Material.hh"
40#include "G4MaterialTable.hh"
41#include "G4UnitsTable.hh"
42
43#include "G4OpticalSurface.hh"
44#include "G4LogicalBorderSurface.hh"
45#include "G4LogicalSkinSurface.hh"
46#include "G4OpBoundaryProcess.hh"
47#include "RichTbMaterialParameters.hh"
48#include "RichTbGeometryParameters.hh"
49#include "G4MaterialPropertyVector.hh"
50
51
52RichTbMaterial::RichTbMaterial(RichTbRunConfig* RConfig):
53 RichTbAerogelMaterial(std::vector<G4Material*> (MaxNumberOfAerogelTypes)),
54 RichTbFilterMaterial(std::vector<G4Material*>(MaxNumberOfFilterTypes)){
55
56 rConfig=RConfig;
57
58 G4double a,z,density; //a=mass of a mole;
59 // z=mean number of protons;
60 G4String name,symbol;
61 //G4int isz,isn; //isz= number of protons in an isotope;
62 //isn= number of nucleons in an isotope;
63
64
65 G4int numel,natoms; //numel=Number of elements constituting a material.
66 G4double fractionmass;
67 G4double temperature, pressure;
68 // G4double FactorOne=1.0;
69 G4UnitDefinition::BuildUnitsTable();
70
71 //PhotonEnergy
72 G4int ibin=0;
73 G4double PhotonEnergyStep=(PhotonMaxEnergy-PhotonMinEnergy)/
74 NumPhotWaveLengthBins;
75 G4double* PhotonMomentum=new G4double[NumPhotWaveLengthBins];
76 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
77 PhotonMomentum[ibin]=PhotonMinEnergy+PhotonEnergyStep*ibin;
78 }
79
80 G4cout << "\nNow Define Elements ..\n" <<G4endl;
81
82// Nitrogen
83
84 a=14.01*g/mole;
85 G4Element* elN = new G4Element(name="Nitrogen",
86 symbol="N", z=7., a);
87
88//Oxygen
89
90 a=16.00*g/mole;
91 G4Element* elO = new G4Element(name="Oxygen",
92 symbol="O", z=8., a);
93
94//Hydrogen
95
96 a=1.01*g/mole;
97 G4Element* elH = new G4Element(name="Hydrogen",
98 symbol="H",z=1.,a);
99
100//Carbon
101
102 a=12.01*g/mole;
103 G4Element* elC = new G4Element(name="Carbon",
104 symbol="C",z=6.,a);
105
106//Silicon
107
108 a=28.09*g/mole;
109 G4Element* elSi = new G4Element(name="Silicon",
110 symbol="Si",z=14.,a);
111 //Fluorine
112 a=18.998*g/mole;
113 G4Element* elF = new G4Element(name="Fluorine",
114 symbol="F",z=9.,a);
115 //Aluminum
116 a=26.98*g/mole;
117 G4Element* elAL =new G4Element(name="Aluminium",
118 symbol="Al",z=13.,a);
119
120 //Sodium
121 a=22.99*g/mole;
122 G4Element* elNa = new G4Element(name="Sodium",
123 symbol="Na",z=11.,a);
124
125 //Potassium
126 a=39.10*g/mole;
127 G4Element* elK = new G4Element(name="Potassium",
128 symbol="K",z=19.,a);
129
130 //Cesium
131
132 // a=132.91*g/mole;
133 // G4Element* elCs = new G4Element(name="Cesium",
134 // symbol="Cs",z=55.,a);
135
136 //Antimony
137
138 a=121.76*g/mole;
139 G4Element* elSb = new G4Element(name="Antimony",
140 symbol="Sb",z=51.,a);
141
142
143//Define Materials
144 G4cout << "\nNow Define Materials ..\n" <<G4endl;
145 //
146
147//Air at 20 degree C and 1 atm for the ambiet air.
148 // Also Air as a radiator material for inside the tubes.
149//--
150 density = 1.205e-03*g/cm3;
151 pressure=1.*atmosphere;
152 temperature=293.*kelvin;
153 G4Material* Air = new G4Material(name="Air ", density, numel=2,
154 kStateGas,temperature,pressure);
155 Air->AddElement(elN, fractionmass=0.7);
156 Air->AddElement(elO, fractionmass=0.3);
157
158 G4double* AirAbsorpLength=new G4double[NumPhotWaveLengthBins];
159 G4double* AirRindex=new G4double[NumPhotWaveLengthBins];
160
161 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
162 AirAbsorpLength[ibin]=1.E32*mm;
163 AirRindex[ibin]=1.000273;
164 }
165 G4MaterialPropertiesTable* AirMPT =
166 new G4MaterialPropertiesTable();
167
168 AirMPT->AddProperty("ABSLENGTH",PhotonMomentum,
169 AirAbsorpLength,NumPhotWaveLengthBins);
170
171 Air->SetMaterialPropertiesTable(AirMPT);
172 RichTbAmbientAir = Air;
173
174
175 density = 1.205e-03*g/cm3;
176 pressure=1.*atmosphere;
177 temperature=293.*kelvin;
178 G4Material* TAir = new G4Material(name="TAir ", density, numel=2,
179 kStateGas,temperature,pressure);
180 TAir->AddElement(elN, fractionmass=0.7);
181 TAir->AddElement(elO, fractionmass=0.3);
182
183 G4double* TAirAbsorpLength=new G4double[NumPhotWaveLengthBins];
184 G4double* TAirRindex=new G4double[NumPhotWaveLengthBins];
185
186 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
187 TAirAbsorpLength[ibin]=1.E32*mm;
188 TAirRindex[ibin]=1.000273;
189 }
190 G4MaterialPropertiesTable* TAirMPT =
191 new G4MaterialPropertiesTable();
192
193 TAirMPT->AddProperty("ABSLENGTH",PhotonMomentum,
194 TAirAbsorpLength,NumPhotWaveLengthBins);
195
196 TAirMPT->AddProperty("RINDEX", PhotonMomentum,
197 AirRindex,NumPhotWaveLengthBins);
198
199 TAir->SetMaterialPropertiesTable(TAirMPT);
200 RichTbTubeAir = TAir;
201
202 //Nitrogen gas.
203
204 density = 0.8073e-03*g/cm3;
205 pressure = RConfig -> getPressureN2();
206 temperature = RConfig ->getTemperatureN2();
207
208 G4Material* NitrogenGas = new G4Material(name="NitrogenGas ",
209 density, numel=1,
210 kStateGas,temperature,pressure);
211 NitrogenGas->AddElement(elN, natoms=2);
212
213 G4double* NitrogenGasAbsorpLength=new G4double[NumPhotWaveLengthBins];
214 G4double* NitrogenGasRindex=new G4double[NumPhotWaveLengthBins];
215 G4double* NitrogenGasPhotW=new G4double[NumPhotWaveLengthBins];
216
217
218 std::vector<G4double>N2RefInd= InitN2RefIndex(pressure,temperature);
219 std::vector<G4double>N2RefPhotW=InitN2RefPhotW();
220
221 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
222 NitrogenGasAbsorpLength[ibin]=1.E32*mm;
223
224 NitrogenGasRindex[ibin]=N2RefInd[ibin];
225 NitrogenGasPhotW[ibin]=N2RefPhotW[ibin];
226
227 }
228 G4MaterialPropertiesTable* NitrogenGasMPT =
229 new G4MaterialPropertiesTable();
230
231 NitrogenGasMPT->AddProperty("ABSLENGTH",NitrogenGasPhotW,
232 NitrogenGasAbsorpLength,NumPhotWaveLengthBins);
233
234 NitrogenGasMPT->AddProperty("RINDEX", NitrogenGasPhotW,
235 NitrogenGasRindex,NumPhotWaveLengthBins);
236
237 NitrogenGas->SetMaterialPropertiesTable(NitrogenGasMPT);
238 RichTbNitrogenGas = NitrogenGas;
239
240//Water
241 density=1.000*g/cm3;
242 G4Material* H2O = new G4Material(name="Water",density,numel=2);
243 H2O->AddElement(elH,natoms=2);
244 H2O->AddElement(elO,natoms=1);
245
246 G4double* H2OAbsorpLength=new G4double[NumPhotWaveLengthBins];
247 G4double* H2ORindex=new G4double[NumPhotWaveLengthBins];
248
249 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
250 H2OAbsorpLength[ibin]=1.E32*mm;
251 H2ORindex[ibin]=1.33;
252 }
253
254
255 G4MaterialPropertiesTable* H2OMPT =
256 new G4MaterialPropertiesTable();
257
258 H2OMPT->AddProperty("ABSLENGTH",PhotonMomentum,
259 H2OAbsorpLength,NumPhotWaveLengthBins);
260
261 H2OMPT->AddProperty("RINDEX", PhotonMomentum,
262 H2ORindex,NumPhotWaveLengthBins);
263
264 H2O->SetMaterialPropertiesTable(H2OMPT);
265
266
267 RichTbH2O=H2O;
268//Sio2
269//There is a quartz for the mirror and
270 //another quartz which is used in aerogel and
271 // yet another quartz used for the quartz window.
272 //Mirrorquartz
273
274 density=2.200*g/cm3;
275 G4Material* SiO2MirrorQuartz = new G4Material(name="MirrorQuartz",
276 density,numel=2);
277 SiO2MirrorQuartz->AddElement(elSi,natoms=1);
278 SiO2MirrorQuartz->AddElement(elO,natoms=2);
279
280 G4double* MirrorQuartzRindex=new G4double[NumPhotWaveLengthBins];
281 G4double* MirrorQuartzAbsorpLength=new G4double[NumPhotWaveLengthBins];
282 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
283 MirrorQuartzAbsorpLength[ibin]=0.01*mm;
284
285 }
286 G4MaterialPropertiesTable* MirrorQuartzMPT =
287 new G4MaterialPropertiesTable();
288
289
290 MirrorQuartzMPT->AddProperty("ABSLENGTH",PhotonMomentum,
291 MirrorQuartzAbsorpLength,NumPhotWaveLengthBins);
292
293
294 SiO2MirrorQuartz->SetMaterialPropertiesTable(MirrorQuartzMPT);
295 RichTbMirrorQuartz=SiO2MirrorQuartz;
296
297 density=2.200*g/cm3;
298 G4Material* SiO2AerogelQuartz = new G4Material(name="AerogelQuartz",
299 density,numel=2);
300 SiO2AerogelQuartz->AddElement(elSi,natoms=1);
301 SiO2AerogelQuartz->AddElement(elO,natoms=2);
302
303 // QuartzWindow Quartz
304 density=2.200*g/cm3;
305 G4Material* WindowQuartz = new G4Material(name="WindowQuartz",
306 density,numel=2);
307 WindowQuartz->AddElement(elSi,natoms=1);
308 WindowQuartz->AddElement(elO,natoms=2);
309 G4double* WindowQuartzRindex=new G4double[NumPhotWaveLengthBins];
310 G4double* WindowQuartzAbsorpLength=new G4double[NumPhotWaveLengthBins];
311 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
312 WindowQuartzAbsorpLength[ibin]=1.E32*mm;
313 WindowQuartzRindex[ibin]=1.4;
314 }
315 G4MaterialPropertiesTable* WindowQuartzMPT =
316 new G4MaterialPropertiesTable();
317
318 WindowQuartzMPT->AddProperty("ABSLENGTH",PhotonMomentum,
319 WindowQuartzAbsorpLength,NumPhotWaveLengthBins);
320
321 WindowQuartzMPT->AddProperty("RINDEX", PhotonMomentum,
322 WindowQuartzRindex,NumPhotWaveLengthBins);
323 WindowQuartz->SetMaterialPropertiesTable(WindowQuartzMPT);
324
325 RichTbQuartzWindowMaterial=WindowQuartz;
326 //for now this is kept to be same as the hpdquartz window.
327 density=2.200*g/cm3;
328 G4Material* HpdWindowQuartz = new G4Material(name="HpdWindowQuartz",
329 density,numel=2);
330 HpdWindowQuartz->AddElement(elSi,natoms=1);
331 HpdWindowQuartz->AddElement(elO,natoms=2);
332 G4double* HpdWindowQuartzRindex=new G4double[NumPhotWaveLengthBins];
333 G4double* HpdWindowQuartzAbsorpLength=new G4double[NumPhotWaveLengthBins];
334 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
335 HpdWindowQuartzAbsorpLength[ibin]=1.E32*mm;
336 HpdWindowQuartzRindex[ibin]=1.40;
337 }
338 G4MaterialPropertiesTable* HpdWindowQuartzMPT =
339 new G4MaterialPropertiesTable();
340
341 HpdWindowQuartzMPT->AddProperty("ABSLENGTH",PhotonMomentum,
342 HpdWindowQuartzAbsorpLength,NumPhotWaveLengthBins);
343
344 HpdWindowQuartzMPT->AddProperty("RINDEX", PhotonMomentum,
345 HpdWindowQuartzRindex,NumPhotWaveLengthBins);
346 HpdWindowQuartz->SetMaterialPropertiesTable(HpdWindowQuartzMPT);
347
348 HpdQuartzWindowMaterial=HpdWindowQuartz;
349
350 // Borosilcate window of the Pad Hpd
351 // for now kept same as the other Hpd.
352 density=2.200*g/cm3;
353 G4Material* PadHpdWindowQuartz = new G4Material(name="PadHpdWindowQuartz",
354 density,numel=2);
355 PadHpdWindowQuartz->AddElement(elSi,natoms=1);
356 PadHpdWindowQuartz->AddElement(elO,natoms=2);
357 G4double* PadHpdWindowQuartzRindex=new G4double[NumPhotWaveLengthBins];
358 G4double* PadHpdWindowQuartzAbsorpLength=new G4double[NumPhotWaveLengthBins];
359 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
360 PadHpdWindowQuartzAbsorpLength[ibin]=1.E32*mm;
361 PadHpdWindowQuartzRindex[ibin]=1.40;
362 }
363 G4MaterialPropertiesTable* PadHpdWindowQuartzMPT =
364 new G4MaterialPropertiesTable();
365
366 PadHpdWindowQuartzMPT->AddProperty("ABSLENGTH",PhotonMomentum,
367 PadHpdWindowQuartzAbsorpLength,NumPhotWaveLengthBins);
368
369 PadHpdWindowQuartzMPT->AddProperty("RINDEX", PhotonMomentum,
370 PadHpdWindowQuartzRindex,NumPhotWaveLengthBins);
371 PadHpdWindowQuartz->SetMaterialPropertiesTable(PadHpdWindowQuartzMPT);
372
373 PadHpdQuartzWindowMaterial=PadHpdWindowQuartz;
374 //
375 //
376
377 G4int filterNumberThisRun=RConfig->GetFilterTNumber();
378 // now for the filter material glass d263
379 density=2.200*g/cm3;
380 G4Material* GlassD263 = new G4Material(name= FilterTypeString[0],
381 density,numel=2);
382 GlassD263->AddElement(elSi,natoms=1);
383 GlassD263->AddElement(elO,natoms=2);
384
385 if(filterNumberThisRun >= 0 ) {
386
387 //in the following the +2 is to match the materialproperty bins
388 // for the various materials, to avoid the tons of printout from G4.
389 // Please see the explanation below for getting the abosorption
390 // length of aerogel. The same comments apply here as well.
391 // Essentially the measured transmission input here is a combination of
392 // the bulk absorption and the fresnel surface loss. One needs to
393 // decouple them. Here a partial attempt is made to avoid
394 // modifying the G4OpBoundary process. SE. 15-11-2002.
395 G4double* GlassD263Rindex=new G4double[NumPhotBinGlassD263Trans+2];
396 G4double* GlassD263AbsorpLength=new G4double[NumPhotBinGlassD263Trans+2];
397 G4double* GlassD263MomValue = new G4double[NumPhotBinGlassD263Trans+2];
398 G4double* currBulkTransFilter = new G4double[NumPhotBinGlassD263Trans+2];
399 FilterTrData* CurFil = RConfig->GetFilterTrData();
400 std::vector<G4double>GlassD263TransWL = CurFil-> GetTransWL();
401 std::vector<G4double>GlassD263Transmis = CurFil->GetTransTotValue();
402 G4double FilterHalfZ= CurFil->GetCurFilterThickness();
403
404 for (ibin=0; ibin<NumPhotBinGlassD263Trans+2; ibin++){
405
406
407 GlassD263Rindex[ibin]=RefIndexGlassD263;
408 if(ibin > 0 && ibin < NumPhotBinGlassD263Trans+1 ){
409 //now using the formula trans=std::exp(-thickness/absorplength).
410 G4int ibina=ibin-1;
411
412 if(GlassD263TransWL[ibina] > 0.0 ) {
413 GlassD263MomValue[ibin]= PhotMomWaveConv*eV/GlassD263TransWL[ibina];
414 }
415 if(GlassD263Transmis[ibina] >0.0 ) {
416 // G4double currentfilterRefIndex= GlassD263Rindex[ibin];
417 G4double currentAdjacentMediumRefIndex=NitrogenNominalRefIndex;
418 // the following needs to be improved in the future
419 // to have a binary search and
420 // interpolation between the adjacent
421 // array elements etc. SE 15-11-2002.
422 for(size_t ibinr=0; ibinr<N2RefPhotW.size()-1 ; ibinr++){
423 G4double currMomA=GlassD263MomValue[ibin];
424 if(currMomA >= N2RefPhotW[ibinr] && currMomA <= N2RefPhotW[ibinr+1]){
425 currentAdjacentMediumRefIndex=N2RefInd[ibinr];
426 }
427
428 }
429
430 if( GlassD263Transmis[ibina] > 0.01 ) {
431 currBulkTransFilter[ibin]=
432 GetCurrentBulkTrans(GlassD263Rindex[ibin],
433 currentAdjacentMediumRefIndex,
434 GlassD263Transmis[ibina]);
435 } else {
436 currBulkTransFilter[ibin]=GlassD263Transmis[ibina];
437
438 }
439 if(currBulkTransFilter[ibin] > 0.0 &&
440 currBulkTransFilter[ibin] < 0.9995 ) {
441 GlassD263AbsorpLength[ibin]=
442 -(2.0*FilterHalfZ)/(std::log(currBulkTransFilter[ibin]));
443 }else if (currBulkTransFilter[ibin]== 0.0 ) {
444 GlassD263AbsorpLength[ibin]=FilterHalfZ/1.0E32;
445 }else {
446 GlassD263AbsorpLength[ibin]=DBL_MAX;
447 }
448 }else {
449
450 GlassD263AbsorpLength[ibin]=FilterHalfZ/1.0E32;
451 }
452 }
453
454 }
455 GlassD263MomValue[0]=PhotonMaxEnergy;
456 GlassD263AbsorpLength[0]=GlassD263AbsorpLength[1];
457 currBulkTransFilter[0]=currBulkTransFilter[1];
458
459 G4int mbin=NumPhotBinGlassD263Trans+1;
460 GlassD263MomValue[mbin]=PhotonMinEnergy;
461 GlassD263AbsorpLength[mbin]=GlassD263AbsorpLength[mbin-1];
462 currBulkTransFilter[mbin]=currBulkTransFilter[mbin-1];
463
464 G4MaterialPropertiesTable* GlassD263MPT =
465 new G4MaterialPropertiesTable();
466
467 GlassD263MPT->AddProperty("ABSLENGTH",GlassD263MomValue,
468 GlassD263AbsorpLength,NumPhotBinGlassD263Trans+2);
469
470 GlassD263MPT->AddProperty("RINDEX",GlassD263MomValue,
471 GlassD263Rindex,NumPhotBinGlassD263Trans+2);
472
473 GlassD263->SetMaterialPropertiesTable(GlassD263MPT);
474 }
475
476 GlassD263FilterMaterial=GlassD263;
477 RichTbFilterMaterial[0]=GlassD263;
478 //for the G4Example only 1 filter type is used.
479 G4cout << " Now Define Aerogel .." <<G4endl;
480
481
482//Aerogel upto five types considered so far.
483 // in the G4example the same type is repeated 5 times.
484//Now for TypeA
485
486 density=0.200*g/cm3;
487
488 G4Material* AerogTypeA =
489 new G4Material(name=AerogelTypeString[0], density, numel=2);
490 AerogTypeA->AddMaterial(SiO2AerogelQuartz, fractionmass=97.0*perCent);
491 AerogTypeA->AddMaterial(H2O, fractionmass=3.0*perCent);
492
493
494 G4double* AerogTypeARindex=new G4double[NumPhotWaveLengthBins];
495 G4double* AerogTypeAAbsorpLength=new G4double[NumPhotWaveLengthBins];
496 G4double* AerogTypeARScatLength = new G4double[NumPhotWaveLengthBins];
497 G4double* currentAgelTrans = new G4double[NumPhotWaveLengthBins];
498
499 std::vector<G4double>AerogelTypeASLength = GetAerogelRScatLength(AerogelTypeA);
500 G4int AerogNumber=0;
501 G4double AerogelLength=GetCurAerogelLength(AerogNumber);
502 G4double MaxTotTransmission=AerogelTypeATotTrans;
503 // Unfortunately the transmission measurement values only give the
504 // total transmission which includes the loss within aerogel
505 // and the Fresnel loss at the surface. In order to
506 // partially decouple this, the approximate loss at the
507 // the surface is calculated using the ref index of the
508 // aerogel and its surroundings. Then this is added to the
509 // measured transmission to get the transmission in the bulk of
510 // aerogel. This is then converted to an absorption length.
511 // In a more accurate implementation the loss at the surface
512 // should be calculated using a more precise formula. It is
513 // difficult since we do not know the direction of the photons
514 // at this point.
515 // One possibility is to modify the G4opBoundaryProcess
516 // for this, since we do know the direction of the photons by then.
517 // This is not done for this G4example, but only in the LHCb implementation.
518 // SE 15-11-2002.
519 // The aerogel is inside a volume made of Nitrogen
520
521 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
522 AerogTypeARindex[ibin]= ConvertAgelRIndex(PhotonMomentum[ibin],0);
523 AerogTypeARScatLength[ibin]=AerogelTypeASLength[ibin];
524 // G4double photwl = PhotMomWaveConv/ (PhotonMomentum[ibin]/eV);
525
526 G4double currentAgelRefIndex= AerogTypeARindex[ibin];
527 G4double currentNeighbourRefIndex= N2RefInd[ibin];
528 currentAgelTrans[ibin]=
529 GetCurrentBulkTrans( currentAgelRefIndex,
530 currentNeighbourRefIndex,MaxTotTransmission);
531 //now using the formula trans=std::exp(-thickness/absorplength)
532 // to get the absorplength.
533
534 if( currentAgelTrans[ibin] > 0.0 && currentAgelTrans[ibin] < 0.9995) {
535 AerogTypeAAbsorpLength[ibin]=
536 -(AerogelLength)/(std::log( currentAgelTrans[ibin]));
537 }else if (currentAgelTrans[ibin] == 0.0) {
538
539 AerogTypeAAbsorpLength[ibin]=AerogelLength/1.0E32;
540 }else {
541
542 AerogTypeAAbsorpLength[ibin]=DBL_MAX;
543 }
544
545 }
546
547 G4MaterialPropertiesTable* AerogTypeAMPT =
548 new G4MaterialPropertiesTable();
549
550 AerogTypeAMPT->AddProperty("ABSLENGTH",PhotonMomentum,
551 AerogTypeAAbsorpLength,NumPhotWaveLengthBins);
552
553
554 AerogTypeAMPT->AddProperty("RAYLEIGH",PhotonMomentum,
555 AerogTypeARScatLength,NumPhotWaveLengthBins);
556
557 AerogTypeAMPT->AddProperty("RINDEX", PhotonMomentum,
558 AerogTypeARindex,NumPhotWaveLengthBins);
559
560 AerogTypeA->SetMaterialPropertiesTable(AerogTypeAMPT);
561
562
563 RichTbAerogelTypeA = AerogTypeA;
564 RichTbAerogelMaterial[0] = AerogTypeA;
565 // In the G4example the same type is repeated 5 times.
566 // in the LHCb implementation 5 types of aerogel materials used.
567 //Now for Aerogel TypeB
568
569 RichTbAerogelTypeB = AerogTypeA;
570 RichTbAerogelMaterial[1] = AerogTypeA;
571
572 //Now for aerogel TypeC
573
574 RichTbAerogelTypeC = AerogTypeA;
575 RichTbAerogelMaterial[2] = AerogTypeA;
576
577 //Now for aerogel TypeD
578
579 RichTbAerogelTypeD = AerogTypeA;
580 RichTbAerogelMaterial[3] = AerogTypeA;
581
582 //Now for aerogel Type E
583
584
585 RichTbAerogelTypeE = AerogTypeA;
586 RichTbAerogelMaterial[4] = AerogTypeA;
587
588
589
590 //Bialkali Photocathode
591
592 //the following numbers on the property of the BiAlkali Photocathode
593 // may not be accurate.
594 //Some number is is jut put in for initial program test purposes.
595 density=0.100*g/cm3;
596 G4Material* BiAlkaliPhCathode = new G4Material(name="BiAlkaliPhCathode",
597 density, numel=3);
598 BiAlkaliPhCathode->AddElement(elNa, fractionmass=37.5*perCent);
599 BiAlkaliPhCathode->AddElement(elK, fractionmass=37.5*perCent);
600 BiAlkaliPhCathode->AddElement(elSb, fractionmass=25.0*perCent);
601
602 //for now properties for the ph cathode material.
603
604 G4double* BiAlkaliPhCathodeRindex=new G4double[NumPhotWaveLengthBins];
605 G4double* BiAlkaliPhCathodeAbsorpLength=new G4double[NumPhotWaveLengthBins];
606 G4double CathLen=PhotoCathodeThickness;
607 G4double CathTrans=PhCathodeNominalTransmission;
608 G4double CathAbsorpLen;
609 if(CathTrans > 0.0 && CathTrans < 0.9995 ) {
610 CathAbsorpLen = -(CathLen)/(std::log(CathTrans));
611 }else if (CathTrans > 0.0) {
612 CathAbsorpLen = CathLen/1.0E32;
613 }else {
614 CathAbsorpLen = DBL_MAX;
615 }
616
617 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
618 BiAlkaliPhCathodeAbsorpLength[ibin]=CathAbsorpLen;
619 BiAlkaliPhCathodeRindex[ibin]=1.40;
620 }
621 G4MaterialPropertiesTable* BiAlkaliPhCathodeMPT =
622 new G4MaterialPropertiesTable();
623 BiAlkaliPhCathodeMPT->AddProperty("ABSLENGTH",PhotonMomentum,
624 BiAlkaliPhCathodeAbsorpLength,NumPhotWaveLengthBins);
625
626 BiAlkaliPhCathodeMPT->AddProperty("RINDEX", PhotonMomentum,
627 BiAlkaliPhCathodeRindex,NumPhotWaveLengthBins);
628 BiAlkaliPhCathode->SetMaterialPropertiesTable(BiAlkaliPhCathodeMPT);
629 PadHpdPhCathodeMaterial=BiAlkaliPhCathode;
630
631//CF4
632//no data available at room temp and pressure;
633 density=0.003884*g/cm3;
634 temperature=273.*kelvin;
635 pressure=1.0*atmosphere;
636 a=88.01*g/mole;
637
638 G4Material* CF4 =new G4Material(name="CF4",density,numel=2,
639 kStateGas,temperature,pressure);
640 CF4->AddElement(elC,natoms=1);
641 CF4->AddElement(elF,natoms=4);
642 // Sellmeir coef to be added.
643 RichTbCF4=CF4;
644
645 G4cout << "\nNowDefineVacuum ..\n" <<G4endl;
646
647//Vacuum
648//
649 density=universe_mean_density;
650 a=1.01*g/mole;
651 pressure=1.e-19*pascal;
652 temperature=0.1*kelvin;
653
654 G4Material* vacuum = new G4Material(name="Galactic",density,numel=1,
655 kStateGas,temperature,pressure);
656 vacuum->AddElement(elH,natoms=1);
657
658 G4double* VacAbsorpLength=new G4double[NumPhotWaveLengthBins];
659 G4double* VacRindex=new G4double[NumPhotWaveLengthBins];
660
661 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
662 VacAbsorpLength[ibin]=1.E32*mm;
663 // the following ref index is just artifical, just to
664 // avoid the refraction between nitrogen gas and hpd master.
665 VacRindex[ibin]=1.000273;
666 }
667 G4MaterialPropertiesTable* VacMPT =
668 new G4MaterialPropertiesTable();
669
670 VacMPT->AddProperty("ABSLENGTH",PhotonMomentum,
671 VacAbsorpLength,NumPhotWaveLengthBins);
672 VacMPT->AddProperty("RINDEX", PhotonMomentum,
673 VacRindex,NumPhotWaveLengthBins);
674 vacuum->SetMaterialPropertiesTable(VacMPT);
675
676 RichTbVacuum=vacuum;
677
678//beamgas
679//
680 density=1.e-5*g/cm3;
681 pressure=2.e-2*bar;
682 temperature=STP_Temperature;
683 G4Material* beamgas = new G4Material(name="Beamgas",density,numel=1,
684 kStateGas,temperature,pressure);
685 beamgas->AddMaterial(Air,fractionmass=1.); // beware that air is at 20 deg;
686
687//
688//Aluminium
689 density=2.7*g/cm3;
690 G4Material* Aluminium =new G4Material(name="Aluminium",density,numel=1);
691 Aluminium->AddElement(elAL,natoms=1);
692
693 G4double* AluminiumAbsorpLength=new G4double[NumPhotWaveLengthBins];
694
695 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
696 AluminiumAbsorpLength[ibin]=0.0*mm;
697 }
698
699 G4MaterialPropertiesTable* AluminiumMPT =
700 new G4MaterialPropertiesTable();
701
702 AluminiumMPT->AddProperty("ABSLENGTH",PhotonMomentum,
703 AluminiumAbsorpLength,NumPhotWaveLengthBins);
704
705 Aluminium->SetMaterialPropertiesTable(AluminiumMPT);
706 RichTbAluminium=Aluminium;
707//PlasticAg , this is used as a wrap of aerogel and as upstream holder
708 // for aerogel frame. For now use same properties as that of Aluminium.
709 // this is just an opaque material.
710
711 density=2.7*g/cm3;
712 G4Material* PlasticAg =new G4Material(name="PlasticAg",density,numel=1);
713 PlasticAg->AddElement(elAL,natoms=1);
714
715 G4double* PlasticAgAbsorpLength=new G4double[NumPhotWaveLengthBins];
716
717 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
718 PlasticAgAbsorpLength[ibin]=0.0*mm;
719 }
720
721 G4MaterialPropertiesTable* PlasticAgMPT =
722 new G4MaterialPropertiesTable();
723
724 PlasticAgMPT->AddProperty("ABSLENGTH",PhotonMomentum,
725 PlasticAgAbsorpLength,NumPhotWaveLengthBins);
726
727 PlasticAg->SetMaterialPropertiesTable(PlasticAgMPT);
728 RichTbPlasticAg=PlasticAg;
729// Kovar
730 density=2.7*g/cm3;
731 G4Material* Kovar =new G4Material(name="Kovar",density,numel=1);
732 Kovar->AddElement(elAL,natoms=1);
733
734 G4double* KovarAbsorpLength=new G4double[NumPhotWaveLengthBins];
735
736 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
737 KovarAbsorpLength[ibin]=0.0*mm;
738 }
739
740 G4MaterialPropertiesTable* KovarMPT =
741 new G4MaterialPropertiesTable();
742
743 KovarMPT->AddProperty("ABSLENGTH",PhotonMomentum,
744 KovarAbsorpLength,NumPhotWaveLengthBins);
745
746 Kovar->SetMaterialPropertiesTable(KovarMPT);
747 HpdTubeMaterial=Kovar;
748
749// Silicon
750
751 density=2.33*g/cm3;
752 G4Material* Silicon =new G4Material(name="Silicon",density,numel=1);
753 Silicon->AddElement(elSi,natoms=1);
754
755 G4double* SiliconAbsorpLength=new G4double[NumPhotWaveLengthBins];
756
757 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
758 SiliconAbsorpLength[ibin]=0.0*mm;
759 }
760
761 G4MaterialPropertiesTable* SiliconMPT =
762 new G4MaterialPropertiesTable();
763
764 SiliconMPT->AddProperty("ABSLENGTH",PhotonMomentum,
765 SiliconAbsorpLength,NumPhotWaveLengthBins);
766
767 Silicon->SetMaterialPropertiesTable(SiliconMPT);
768 HpdSiDetMaterial=Silicon;
769
770// Silicon coating made of Si02.
771
772 density=2.33*g/cm3;
773 G4Material* SiliconCoating =new G4Material(name="SilCoat",density,numel=2);
774 SiliconCoating->AddElement(elSi,natoms=1);
775 SiliconCoating->AddElement(elO,natoms=2);
776
777 G4double* SiliconCoatingAbsorpLength=new G4double[NumPhotWaveLengthBins];
778 // G4double* SiliconCoatingRindex=new G4double[NumPhotWaveLengthBins];
779
780 for (ibin=0; ibin<NumPhotWaveLengthBins; ibin++){
781 SiliconCoatingAbsorpLength[ibin]=0.0001*mm;
782
783 }
784
785 G4MaterialPropertiesTable* SiliconCoatingMPT =
786 new G4MaterialPropertiesTable();
787
788 SiliconCoatingMPT->AddProperty("ABSLENGTH",PhotonMomentum,
789 SiliconCoatingAbsorpLength,NumPhotWaveLengthBins);
790
791 SiliconCoating->SetMaterialPropertiesTable(SiliconCoatingMPT);
792 HpdSiCoatingMaterial=SiliconCoating;
793
794 //
795 // Now for the material properties of Surfaces
796 //
797 //
798 //
799 //Front (reflecting surface of RichTb Mirror)
800
801 // First define wavelength in nm.
802 //For now assume that all segments have the same reflectivity.
803 // Hence the reflectivity is defined outside the loop of the
804 // the number of segments.
805 //Only the front surface is created.
806 // The abosorption length is set to a small value just to
807 // avoid photons exiting from the back of the mirror.
808 // the efficiency is for the absorption process.
809
810
811 G4double* PhotonMomentumRefl
812 =new G4double[NumPhotonRichMirrorReflWaveLengthBins];
813 G4double* PhotWaveRefl =
814 new G4double[NumPhotonRichMirrorReflWaveLengthBins];
815 G4double* PhotReflEff =new G4double[NumPhotonRichMirrorReflWaveLengthBins];
816 G4double* MirrorQuRefIndex
817 =new G4double[NumPhotonRichMirrorReflWaveLengthBins];
818
819 for (ibin=0; ibin<NumPhotonRichMirrorReflWaveLengthBins; ibin++){
820 PhotonMomentumRefl[ibin]=PhotMomWaveConv*eV/ PhotonWavelengthRefl[ibin];
821 PhotWaveRefl[ibin]= RichTbMirrorReflectivity[ibin];
822 PhotReflEff[ibin]= RichTbMirrorEfficiency[ibin];
823 //the following lines to avoid reflection at the mirror.
824
825 MirrorQuRefIndex[ibin] = 1.40;
826 }
827
828 G4OpticalSurface * OpRichTbMirrorSurface =
829 new G4OpticalSurface("RichTbMirrorSurface");
830
831 OpRichTbMirrorSurface->SetType(dielectric_metal);
832 OpRichTbMirrorSurface->SetFinish(polished);
833 OpRichTbMirrorSurface->SetModel(glisur);
834 G4MaterialPropertiesTable* OpRichTbMirrorSurfaceMPT =
835 new G4MaterialPropertiesTable();
836
837 OpRichTbMirrorSurfaceMPT->AddProperty("REFLECTIVITY",
838 PhotonMomentumRefl,
839 PhotWaveRefl,
840 NumPhotonRichMirrorReflWaveLengthBins);
841 OpRichTbMirrorSurfaceMPT->AddProperty("EFFICIENCY",
842 PhotonMomentumRefl,
843 PhotReflEff,
844 NumPhotonRichMirrorReflWaveLengthBins);
845 OpRichTbMirrorSurfaceMPT->AddProperty("RINDEX",
846 PhotonMomentumRefl,
847 MirrorQuRefIndex,
848 NumPhotonRichMirrorReflWaveLengthBins);
849
850 OpRichTbMirrorSurface->SetMaterialPropertiesTable(OpRichTbMirrorSurfaceMPT);
851 RichTbOpticalMirrorSurface=OpRichTbMirrorSurface;
852
853
854 // OpRichTbMirrorSurface->DumpInfo();
855
856 // Now for the Surface of the Vessel Enclosure.
857
858
859 G4OpticalSurface * OpRichTbEnclosureSurface =
860 new G4OpticalSurface("RichTbEnclosureSurface");
861 OpRichTbEnclosureSurface->SetType(dielectric_metal);
862 OpRichTbEnclosureSurface->SetFinish(polished);
863 OpRichTbEnclosureSurface->SetModel(glisur);
864
865 G4double NumPhotonRichEnclosureSurfaceWaveLengthBins=10;
866 G4double RichTbEnclosureSurfaceReflectivity[]=
867 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
868
869 G4double RichTbEnclosureSurfaceEfficiency[]=
870 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
871 G4double RichEnclosureSurfacePhotMom[]=
872 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
873 9.0*eV,10.0*eV};
874
875 G4MaterialPropertiesTable* OpRichTbEnclosureSurfaceMPT =
876 new G4MaterialPropertiesTable();
877
878 OpRichTbEnclosureSurfaceMPT->AddProperty("REFLECTIVITY",
879 RichEnclosureSurfacePhotMom,
880 RichTbEnclosureSurfaceReflectivity,
881 static_cast<int>(NumPhotonRichEnclosureSurfaceWaveLengthBins));
882 OpRichTbEnclosureSurfaceMPT->AddProperty("EFFICIENCY",
883 RichEnclosureSurfacePhotMom,
884 RichTbEnclosureSurfaceEfficiency,
885 static_cast<int>(NumPhotonRichEnclosureSurfaceWaveLengthBins));
886
887 OpRichTbEnclosureSurface->
888 SetMaterialPropertiesTable(OpRichTbEnclosureSurfaceMPT);
889
890 RichTbOpticalEnclosureSurface=OpRichTbEnclosureSurface;
891
892 //Now for the surface between the TAir and Quartz Window of the HPD
893
894 G4OpticalSurface * OpHpdQuartzWTSurface =
895 new G4OpticalSurface("HpdQuartzWTSurface");
896 OpHpdQuartzWTSurface->SetType(dielectric_dielectric);
897 OpHpdQuartzWTSurface->SetFinish(polished);
898 OpHpdQuartzWTSurface->SetModel(glisur);
899 //OpHpdQuartzWTSurface->SetModel(unified);
900
901 G4double NumPhotonHpdQuartzWTSurfaceWaveLengthBins=10;
902 G4double HpdQuartzWTSurfaceReflectivity[]=
903 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
904
905 G4double HpdQuartzWTSurfaceEfficiency[]=
906 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
907
908 G4double HpdQuartzWTSurfacePhotMom[]=
909 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
910 9.0*eV,10.0*eV};
911
912 G4MaterialPropertiesTable* OpHpdQuartzWTSurfaceMPT =
913 new G4MaterialPropertiesTable();
914
915
916 OpHpdQuartzWTSurfaceMPT->AddProperty("REFLECTIVITY",
917 HpdQuartzWTSurfacePhotMom,
918 HpdQuartzWTSurfaceReflectivity,
919 static_cast<int>(NumPhotonHpdQuartzWTSurfaceWaveLengthBins));
920 OpHpdQuartzWTSurfaceMPT->AddProperty("EFFICIENCY",
921 HpdQuartzWTSurfacePhotMom,
922 HpdQuartzWTSurfaceEfficiency,
923 static_cast<int>(NumPhotonHpdQuartzWTSurfaceWaveLengthBins));
924
925 OpHpdQuartzWTSurface->
926 SetMaterialPropertiesTable(OpHpdQuartzWTSurfaceMPT);
927
928 HpdTQuartzWSurface=OpHpdQuartzWTSurface;
929
930
931
932 //Now for the surface between the Quartz Window and Ph cathode of the HPD
933
934 G4OpticalSurface * OpHpdQuartzWPSurface =
935 new G4OpticalSurface("HpdQuartzWPSurface");
936 OpHpdQuartzWPSurface->SetType(dielectric_dielectric);
937 OpHpdQuartzWPSurface->SetFinish(polished);
938 OpHpdQuartzWPSurface->SetModel(glisur);
939
940 G4double NumPhotonHpdQuartzWPSurfaceWaveLengthBins=10;
941 G4double HpdQuartzWPSurfaceReflectivity[]=
942 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
943
944
945 G4double HpdQuartzWPSurfaceEfficiency[]=
946 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
947
948 G4double HpdQuartzWPSurfacePhotMom[]=
949 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
950 9.0*eV,10.0*eV};
951
952 G4MaterialPropertiesTable* OpHpdQuartzWPSurfaceMPT =
953 new G4MaterialPropertiesTable();
954
955
956 OpHpdQuartzWPSurfaceMPT->AddProperty("REFLECTIVITY",
957 HpdQuartzWPSurfacePhotMom,
958 HpdQuartzWPSurfaceReflectivity,
959 static_cast<int>(NumPhotonHpdQuartzWPSurfaceWaveLengthBins));
960 OpHpdQuartzWPSurfaceMPT->AddProperty("EFFICIENCY",
961 HpdQuartzWPSurfacePhotMom,
962 HpdQuartzWPSurfaceEfficiency,
963 static_cast<int>(NumPhotonHpdQuartzWPSurfaceWaveLengthBins));
964
965 OpHpdQuartzWPSurface->
966 SetMaterialPropertiesTable(OpHpdQuartzWPSurfaceMPT);
967
968 HpdQuartzWPhCathodeSurface=OpHpdQuartzWPSurface;
969
970
971
972 //Now for the skin surface of the PhCathode so that photons do
973 // not come out of the Photocathode.
974 // Changed to dielectric-dielectric so that photons DO come out
975 // of the photocathode. SE 26-9-01.
976
977 G4OpticalSurface * OpPhCathodeSurface =
978 new G4OpticalSurface("PhCathodeSurface");
979
980 OpPhCathodeSurface->SetType(dielectric_dielectric);
981 OpPhCathodeSurface->SetFinish(polished);
982 OpPhCathodeSurface->SetModel(glisur);
983
984
985 G4double NumPhotonPhCathodeSurfaceWaveLengthBins=10;
986 G4double PhCathodeSurfaceReflectivity[]=
987 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
988
989
990 G4double PhCathodeSurfaceEfficiency[]=
991 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
992
993 G4double PhCathodeSurfacePhotMom[]=
994 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
995 9.0*eV,10.0*eV};
996
997 G4MaterialPropertiesTable* OpPhCathodeSurfaceMPT =
998 new G4MaterialPropertiesTable();
999
1000
1001 OpPhCathodeSurfaceMPT->AddProperty("REFLECTIVITY",
1002 PhCathodeSurfacePhotMom,
1003 PhCathodeSurfaceReflectivity,
1004 static_cast<int>(NumPhotonPhCathodeSurfaceWaveLengthBins));
1005 OpPhCathodeSurfaceMPT->AddProperty("EFFICIENCY",
1006 PhCathodeSurfacePhotMom,
1007 PhCathodeSurfaceEfficiency,
1008 static_cast<int>(NumPhotonPhCathodeSurfaceWaveLengthBins));
1009
1010 OpPhCathodeSurface->
1011 SetMaterialPropertiesTable(OpPhCathodeSurfaceMPT);
1012
1013 PhCathodeSkinSurface=OpPhCathodeSurface;
1014 PhCathodeBorderSurface=OpPhCathodeSurface;
1015
1016
1017
1018
1019 //Now for the surface between Interior of HPD and Silicon Coating.
1020
1021 G4OpticalSurface * OpHpdSiCoatSurface =
1022 new G4OpticalSurface("HpdSiCoatSurface");
1023 OpHpdSiCoatSurface->SetType(dielectric_metal);
1024 OpHpdSiCoatSurface->SetFinish(polished);
1025 OpHpdSiCoatSurface->SetModel(glisur);
1026
1027
1028 G4double NumPhotonHpdSiCoatSurfaceWaveLengthBins=10;
1029
1030 G4double HpdSiCoatSurfaceReflectivity[]=
1031 {0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.9};
1032
1033 G4double HpdSiCoatSurfaceEfficiency[]=
1034 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1035
1036 G4double HpdSiCoatSurfacePhotMom[]=
1037 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
1038 9.0*eV,10.0*eV};
1039 G4double HpdSiCoatSurfaceRefInd[]=
1040 {1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4};
1041
1042
1043 G4MaterialPropertiesTable* OpHpdSiCoatSurfaceMPT =
1044 new G4MaterialPropertiesTable();
1045
1046
1047 OpHpdSiCoatSurfaceMPT->AddProperty("REFLECTIVITY",
1048 HpdSiCoatSurfacePhotMom,
1049 HpdSiCoatSurfaceReflectivity,
1050 static_cast<int>(NumPhotonHpdSiCoatSurfaceWaveLengthBins));
1051 OpHpdSiCoatSurfaceMPT->AddProperty("EFFICIENCY",
1052 HpdSiCoatSurfacePhotMom,
1053 HpdSiCoatSurfaceEfficiency,
1054 static_cast<int>(NumPhotonHpdSiCoatSurfaceWaveLengthBins));
1055 OpHpdSiCoatSurfaceMPT->AddProperty("RINDEX",
1056 HpdSiCoatSurfacePhotMom,
1057 HpdSiCoatSurfaceRefInd,
1058 static_cast<int>(NumPhotonHpdSiCoatSurfaceWaveLengthBins));
1059
1060 OpHpdSiCoatSurface->
1061 SetMaterialPropertiesTable(OpHpdSiCoatSurfaceMPT);
1062
1063 HpdSiCoatSurface=OpHpdSiCoatSurface;
1064
1065
1066
1067 // Now for the Surface of the MetalTube of HPD.
1068
1069
1070 G4OpticalSurface * OpRichTbHpdMetalSurface =
1071 new G4OpticalSurface("RichTbHpdMetalSurface");
1072 OpRichTbHpdMetalSurface->SetType(dielectric_metal);
1073 OpRichTbHpdMetalSurface->SetFinish(polished);
1074 OpRichTbHpdMetalSurface->SetModel(glisur);
1075
1076 G4double NumPhotonHpdMetalSurfaceWaveLengthBins=10;
1077 G4double RichHpdMetalSurfaceReflectivity[]=
1078 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1079
1080 G4double RichHpdMetalSurfaceEfficiency[]=
1081 {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1082 G4double RichHpdMetalSurfacePhotMom[]=
1083 {1.0*eV,2.0*eV, 3.0*eV,4.0*eV,5.0*eV,6.0*eV,7.0*eV,8.0*eV,
1084 9.0*eV,10.0*eV};
1085
1086 G4MaterialPropertiesTable* OpRichTbHpdMetalSurfaceMPT =
1087 new G4MaterialPropertiesTable();
1088
1089 OpRichTbHpdMetalSurfaceMPT->AddProperty("REFLECTIVITY",
1090 RichHpdMetalSurfacePhotMom,
1091 RichHpdMetalSurfaceReflectivity,
1092 static_cast<int>(NumPhotonHpdMetalSurfaceWaveLengthBins));
1093 OpRichTbHpdMetalSurfaceMPT->AddProperty("EFFICIENCY",
1094 RichHpdMetalSurfacePhotMom,
1095 RichHpdMetalSurfaceEfficiency,
1096 static_cast<int>(NumPhotonHpdMetalSurfaceWaveLengthBins));
1097
1098 OpRichTbHpdMetalSurface->
1099 SetMaterialPropertiesTable(OpRichTbHpdMetalSurfaceMPT);
1100
1101 RichTbOpticalHpdMetalSurface=OpRichTbHpdMetalSurface;
1102
1103
1104 //Now for the surface of the Filter
1105
1106 G4OpticalSurface * OpRichTbFilterSurface =
1107 new G4OpticalSurface("RichTbFilterSurface");
1108 OpRichTbFilterSurface->SetType(dielectric_dielectric);
1109 OpRichTbFilterSurface->SetFinish(polished);
1110 OpRichTbFilterSurface->SetModel(glisur);
1111
1112
1113
1114 if(filterNumberThisRun >= 0 ) {
1115
1116 G4int FilterNumbins=NumPhotonRichTbFilterSurfaceWaveLengthBins;
1117
1118
1119 G4double* FilterReflectivity = new G4double(FilterNumbins);
1120 G4double* FilterEff =new G4double(FilterNumbins);
1121 G4double* FilterPhotMom =new G4double(FilterNumbins);
1122
1123
1124 for(G4int ibinf =0 ; ibinf < FilterNumbins; ibinf++ ){
1125 FilterReflectivity[ibinf]= RichTbFilterSurfaceReflectivity[ibinf];
1126 FilterEff[ibinf]= RichTbFilterSurfaceEfficiency[ibinf];
1127 FilterPhotMom[ibinf]= RichTbFilterSurfacePhotMom[ibinf];
1128
1129// G4MaterialPropertiesTable* OpRichTbFilterSurfaceMPT =
1130// new G4MaterialPropertiesTable();
1131
1132 }
1133 RichTbOpticalFilterSurface=OpRichTbFilterSurface;
1134
1135 }
1136
1137 delete [] PhotonMomentum;
1138 delete [] AirAbsorpLength;
1139 delete [] AirRindex;
1140 delete [] MirrorQuartzRindex;
1141 delete [] MirrorQuartzAbsorpLength;
1142 delete [] WindowQuartzRindex;
1143 delete [] WindowQuartzAbsorpLength;
1144 delete [] AluminiumAbsorpLength;
1145 delete [] KovarAbsorpLength;
1146 delete [] PhotonMomentumRefl;
1147
1148
1149}
1150G4double RichTbMaterial::ConvertAgelRIndex(G4double phmom, G4int AgelTnum ) {
1151 AerogelRefData* AgData= rConfig -> GetAerogelRefdata();
1152 //Now to convert and interpolate to get the same binning
1153 // as the other property vectors.
1154 G4double Refind=0.;
1155 G4int Numphbin=AgData-> GetNumberOfRefIndBins();
1156 G4double phm1,phm2;
1157 if(phmom < AgData->GetAerogelRefphotE(0) ){
1158 Refind=AgData->GetCurAerogelRefIndValue(0,AgelTnum ); }
1159 if(phmom >= AgData->GetAerogelRefphotE(Numphbin-1 ) ) {
1160 Refind=AgData->GetCurAerogelRefIndValue(Numphbin-1,AgelTnum ); }
1161
1162 for( G4int iba=0; iba<Numphbin-1 ; iba ++ ) {
1163
1164 phm1=AgData->GetAerogelRefphotE(iba);
1165 phm2=AgData->GetAerogelRefphotE(iba+1);
1166
1167 if(phmom >= phm1 && phmom < phm2 ) {
1168
1169 G4double ref1=AgData->GetCurAerogelRefIndValue(iba,AgelTnum );
1170 G4double ref2=AgData->GetCurAerogelRefIndValue(iba+1,AgelTnum );
1171
1172 G4double grad = (ref2-ref1)/(phm2-phm1);
1173 G4double aint = ref1- grad*phm1;
1174 Refind = grad*phmom + aint ;
1175 break;
1176 }
1177 }
1178
1179 return Refind;
1180}
1181RichTbMaterial::RichTbMaterial() { ; }
1182RichTbMaterial::~RichTbMaterial(){ ; }
1183
1184
1185
Note: See TracBrowser for help on using the repository browser.