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

Last change on this file was 807, checked in by garnier, 16 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.