ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4dRICHOptics.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4dRICHOptics.hh
1 #ifndef G4E_CI_DRICH_MODEL_HH
2 #define G4E_CI_DRICH_MODEL_HH
3 
4 #include <fun4all/Fun4AllBase.h>
5 #include <fun4all/SubsysReco.h>
6 
7 #include <G4Color.hh>
8 #include <G4Element.hh>
9 #include <G4LogicalSkinSurface.hh>
10 #include <G4Material.hh>
11 #include <G4MaterialPropertiesTable.hh>
12 #include <G4OpticalSurface.hh>
13 #include <G4RotationMatrix.hh>
14 #include <G4SystemOfUnits.hh>
15 #include <G4VisAttributes.hh>
16 #include <G4tgbMaterialMgr.hh>
17 #include <G4tgbVolumeMgr.hh>
18 #include <G4PVDivision.hh>
19 
20 #include <iostream>
21 
22 /*
23  * Service Classes
24  */
25 
26 //
27 // Generic optical parameters class
28 //
29 
30 class EICG4dRICHOptics : public SubsysReco
31 {
32  public:
33  double *scaledE; // photon energies
34  double *scaledN; // real refractive index
35  double *scaledA; // absorption length
36  double *scaledS; // scattering length
37 
38  double *scaledSE; // surface efficiency
39  double *scaledSR; // surface reflectivity
40  double *scaledIN; // imaginary refractive index
41 
42  // Add optical properties either to material or volume surface
43  // matName: material name
44  // logVolName: logical volume name
45  // if name is "_NA_", properties are not applied to the corresponding
46  // component
47  EICG4dRICHOptics(const G4String matName, const G4String logVolName)
48  {
49  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
50  {
51  std::cout << "#=======================================================================" << std::endl;
52  std::cout << "# Set Optical Properties" << std::endl;
53  }
54 
55  materialName = matName;
56  logicalVName = logVolName;
57 
58  scaledE = NULL;
59  scaledN = NULL;
60  scaledA = NULL;
61  scaledS = NULL;
62 
63  scaledSE = NULL;
64  scaledSR = NULL;
65  scaledIN = NULL;
66 
67  if (matName != "_NA_")
68  {
69  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "# Material " << matName.data() << std::endl;
70 
71  mat = G4tgbMaterialMgr::GetInstance()->FindBuiltG4Material(matName);
72 
73  if (mat == NULL)
74  {
75  pTable = NULL;
76  std::cout << "# ERROR: Cannot retrieve " << matName.data() << " material in EICG4dRICH" << std::endl;
77  // handle error
78  }
79  else
80  {
81  pTable = mat->GetMaterialPropertiesTable();
82  }
83 
84  if (pTable == NULL)
85  {
86  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
87  {
88  std::cout << "# No properties table available for " << matName.data() << ", allocated a new one\n" << std::endl;
89  }
90  pTable = new G4MaterialPropertiesTable();
91  }
92  else
93  {
94  pTable->DumpTable();
95  }
96  }
97 
98  if (logVolName != "_NA_")
99  {
100  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE) std::cout << "# Logical Volume " << logVolName.data() << std::endl;
101  logVolume = G4tgbVolumeMgr::GetInstance()->FindG4LogVol(logVolName, 0);
102  if (logVolume == NULL)
103  {
104  std::cout << "# ERROR: Cannot retrieve " << logVolName.data() << " logical volume in EICG4dRICH" << std::endl;
105  // handle error
106  }
107  }
108  }
109 
110  ~EICG4dRICHOptics()
111  {
112  if (scaledE != NULL) delete[] scaledE;
113  if (scaledN != NULL) delete[] scaledN;
114  if (scaledA != NULL) delete[] scaledA;
115  if (scaledS != NULL) delete[] scaledS;
116 
117  if (scaledSE != NULL) delete[] scaledSE;
118  if (scaledSR != NULL) delete[] scaledSR;
119  if (scaledIN != NULL) delete[] scaledIN;
120  }
121 
122  // dvalue, ivalue, svalue may represent different quantities depending on
123  // implementation
124  virtual int setOpticalParams() { return -1; }
125  virtual int setOpticalParams(double dvalue) { return -1; }
126  virtual int setOpticalParams(int ivalue) { return -1; }
127  virtual int setOpticalParams(int ivalue, double dvalue) { return -1; }
128  virtual int setOpticalParams(G4String svalue) { return -1; }
129 
130  protected:
131  G4String materialName, logicalVName;
132  G4Material *mat;
133  G4LogicalVolume *logVolume; // used for skin surface
134 
135  G4MaterialPropertiesTable *pTable;
136 
137  double wl2e(double wl) { return 1239.84193 * eV / (wl / nm); } // wavelength to energy
138  double e2wl(double e) { return 1239.84193 * nm / (e / eV); } // energy to wavelength
139 
140  // add properties to the MaterialPropertiesTable and link to material
141  void setMatPropTable(int nEntries)
142  {
143  if (scaledN != NULL) pTable->AddProperty("RINDEX", scaledE, scaledN, nEntries)->SetSpline(true);
144  if (scaledA != NULL) pTable->AddProperty("ABSLENGTH", scaledE, scaledA, nEntries)->SetSpline(true);
145  if (scaledS != NULL) pTable->AddProperty("RAYLEIGH", scaledE, scaledS, nEntries)->SetSpline(true);
146  // pTable->AddConstProperty("SCINTILLATIONYIELD", 0. / MeV); // @@@ TBC
147  // @@@ pTable->AddConstProperty("RESOLUTIONSCALE", 1.0); // @@@ TBC @@@
148 
149  mat->SetMaterialPropertiesTable(pTable);
150  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
151  {
152  std::cout << "# Optical Table for material " << materialName.data() << " with " << nEntries << " points:" << std::endl;
153  pTable->DumpTable();
154  }
155  }
156 
157  // allocate and add properties to the MaterialPropertiesTable
158  G4MaterialPropertiesTable *addSkinPropTable(int nE)
159  {
160  G4MaterialPropertiesTable *pTab = new G4MaterialPropertiesTable();
161  if (scaledSE != NULL) pTab->AddProperty("EFFICIENCY", scaledE, scaledSE, nE);
162  if (scaledSR != NULL) pTab->AddProperty("REFLECTIVITY", scaledE, scaledSR, nE);
163  if (scaledN != NULL) pTab->AddProperty("REALRINDEX", scaledE, scaledN, nE);
164  if (scaledIN != NULL) pTab->AddProperty("IMAGINARYRINDEX", scaledE, scaledIN, nE);
165  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
166  {
167  std::cout << "# Optical Table for volume " << logicalVName.data() << " with " << nE << " points:" << std::endl;
168  pTab->DumpTable();
169  }
170  return pTab;
171  }
172 
173  // Linear Interpolation method
174  double linint(double val, int n, const double *x, const double *y)
175  {
176  if (val <= x[0]) return y[0];
177  if (val >= x[n - 1]) return y[n - 1];
178  for (int i = 0; i < (n - 1); i++)
179  {
180  if ((val >= x[i]) && (val < x[i + 1]))
181  {
182  return (y[i + 1] - y[i]) / (x[i + 1] - x[i]) * (val - x[i]) + y[i];
183  }
184  }
185  return 0.;
186  }
187 };
188 
189 //
190 // Aerogel
191 //
192 class EICG4dRICHAerogel : public EICG4dRICHOptics
193 {
194  public:
195  EICG4dRICHAerogel(const G4String matName) : EICG4dRICHOptics(matName, "_NA_"){}
196  //
197  // Compute the refractive index, absorption length, scattering length for
198  // different energies points
199  //
200  // different methods are available for the refractive index:
201  // 0 - Vorobiev
202  // 1 - Sellmeier - CLAS12
203  // 2 - Sellmeier - LHCB
204  // 3 - CLAS12 experimental points rescaled by Alessio/GEMC (the same used
205  // for the scattering and absorption length)
206  //
207  // data are scaled according to the input density of the aerogel
208  //
209  int setOpticalParams(int mode) override
210  {
211  const double aeroE[] = // energy : wavelenth 660 nm -> 200 nm
212  {1.87855 * eV, 1.96673 * eV, 2.05490 * eV, 2.14308 * eV, 2.23126 * eV,
213  2.31943 * eV, 2.40761 * eV, 2.49579 * eV, 2.58396 * eV, 2.67214 * eV,
214  2.76032 * eV, 2.84849 * eV, 2.93667 * eV, 3.02485 * eV, 3.11302 * eV,
215  3.20120 * eV, 3.28938 * eV, 3.37755 * eV, 3.46573 * eV, 3.55391 * eV,
216  3.64208 * eV, 3.73026 * eV, 3.81844 * eV, 3.90661 * eV, 3.99479 * eV,
217  4.08297 * eV, 4.17114 * eV, 4.25932 * eV, 4.34750 * eV, 4.43567 * eV,
218  4.52385 * eV, 4.61203 * eV, 4.70020 * eV, 4.78838 * eV, 4.87656 * eV,
219  4.96473 * eV, 5.05291 * eV, 5.14109 * eV, 5.22927 * eV, 5.31744 * eV,
220  5.40562 * eV, 5.49380 * eV, 5.58197 * eV, 5.67015 * eV, 5.75833 * eV,
221  5.84650 * eV, 5.93468 * eV, 6.02286 * eV, 6.11103 * eV, 6.19921 * eV};
222 
223  const double aeroN[] = // refractive index - scaled from CLAS12 RICH to
224  // n(300) = 1.02, rho ~ 0.088, n(400)~1.019
225  {1.01825, 1.01829, 1.01834, 1.01839, 1.01844, 1.01849, 1.01854, 1.01860,
226  1.01866, 1.01872, 1.01878, 1.01885, 1.01892, 1.01899, 1.01906, 1.01914,
227  1.01921, 1.01929, 1.01938, 1.01946, 1.01955, 1.01964, 1.01974, 1.01983,
228  1.01993, 1.02003, 1.02014, 1.02025, 1.02036, 1.02047, 1.02059, 1.02071,
229  1.02084, 1.02096, 1.02109, 1.02123, 1.02137, 1.02151, 1.02166, 1.02181,
230  1.02196, 1.02212, 1.02228, 1.02244, 1.02261, 1.02279, 1.02297, 1.02315,
231  1.02334, 1.02354};
232 
233  const double aeroA[] = // absorption length - scaled from CLAS12 RICH to
234  // n(300) = 1.02, rho ~ 0.088
235  {17.5000 * cm, 17.7466 * cm, 17.9720 * cm, 18.1789 * cm, 18.3694 * cm,
236  18.5455 * cm, 18.7086 * cm, 18.8602 * cm, 19.0015 * cm, 19.1334 * cm,
237  19.2569 * cm, 19.3728 * cm, 19.4817 * cm, 19.5843 * cm, 19.6810 * cm,
238  19.7725 * cm, 19.8590 * cm, 19.9410 * cm, 20.0188 * cm, 20.0928 * cm,
239  18.4895 * cm, 16.0174 * cm, 13.9223 * cm, 12.1401 * cm, 10.6185 * cm,
240  9.3147 * cm, 8.1940 * cm, 7.2274 * cm, 6.3913 * cm, 5.6659 * cm,
241  5.0347 * cm, 4.4841 * cm, 4.0024 * cm, 3.5801 * cm, 3.2088 * cm,
242  2.8817 * cm, 2.5928 * cm, 2.3372 * cm, 2.1105 * cm, 1.9090 * cm,
243  1.7296 * cm, 1.5696 * cm, 1.4266 * cm, 1.2986 * cm, 1.1837 * cm,
244  1.0806 * cm, 0.9877 * cm, 0.9041 * cm, 0.8286 * cm, 0.7603 * cm};
245 
246  const double aeroS[] = // rayleigh - scaled from CLAS12 RICH to n(300)
247  // = 1.02, rho ~ 0.088
248  {35.1384 * cm, 29.24805 * cm, 24.5418 * cm, 20.7453 * cm, 17.6553 * cm,
249  15.1197 * cm, 13.02345 * cm, 11.2782 * cm, 9.81585 * cm, 8.58285 * cm,
250  7.53765 * cm, 6.6468 * cm, 5.88375 * cm, 5.22705 * cm, 4.6596 * cm,
251  4.167 * cm, 3.73785 * cm, 3.36255 * cm, 3.03315 * cm, 2.7432 * cm,
252  2.48700 * cm, 2.26005 * cm, 2.05845 * cm, 1.87875 * cm, 1.71825 * cm,
253  1.57455 * cm, 1.44555 * cm, 1.3296 * cm, 1.2249 * cm, 1.1304 * cm,
254  1.04475 * cm, 0.9672 * cm, 0.89655 * cm, 0.83235 * cm, 0.77385 * cm,
255  0.7203 * cm, 0.67125 * cm, 0.6264 * cm, 0.58515 * cm, 0.54735 * cm,
256  0.51255 * cm, 0.48045 * cm, 0.45075 * cm, 0.4233 * cm, 0.39795 * cm,
257  0.37455 * cm, 0.3528 * cm, 0.33255 * cm, 0.3138 * cm, 0.29625 * cm};
258 
259  const int nEntries = sizeof(aeroE) / sizeof(double);
260 
261  double density = mat->GetDensity();
262  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
263  {
264  std::cout << "# Aerogel Density : " << density / (g / cm3) << " g/cm3" << std::endl;
265  }
266 
267  double refn = density2refIndex(density); // use a n vs rho formula with provide n at 400 nm
268  double refwl = 400 * nm;
269 
270  double refee = wl2e(refwl) / eV; // [eV] reference energy
271  double an0 = linint(refee, nEntries, aeroE, aeroN);
272  // double aa0 = linint(refee, nEntries, aeroE, aeroA);
273  // double as0 = linint(refee, nEntries, aeroE, aeroS);
274 
275  double aa;
276  double nn;
277  // double ri;
278  double a0;
279  // double wl0;
280  double rnscale;
281  double rho = 0.0;
282 
283  if (scaledE == NULL)
284  {
285  scaledE = new double[nEntries];
286  scaledN = new double[nEntries];
287  scaledS = new double[nEntries];
288  scaledA = new double[nEntries];
289  }
290 
291  for (int i = 0; i < nEntries; i++)
292  {
293  double ee = aeroE[i];
294  double wl = e2wl(ee) / nm; // from Energy to nm
295 
296  switch (mode)
297  {
298  case 0: // --- Vorobiev model
299  aa = airFraction(refn, refwl);
300  nn = aa * riAir(wl) + (1. - aa) * riQuartz(wl);
301  break;
302  case 1: // --- Sellmeier, 1 pole from (CLAS12/RICH EPJ A (2016) 52: 23)
303  // ri = 1.0494; // 400 nm
304  rho = 0.230; // g/cm3
305  a0 = 0.09683;
306  // wl0 = 84.13;
307  rnscale = sqrt(1. + (a0 * refwl * refwl) / (refwl * refwl - a0 * a0));
308  nn = sqrt(1. + (a0 * wl * wl) / (wl * wl - a0 * a0)) * refn / rnscale;
309  break;
310  case 2: // --- Sellmeier, 1 pole from T. Bellunato et al. Eur. Phys. J.
311  // C52, 759-764 (2007)
312  // ri = 1.03; // 400 nm
313  rho = 0.149; // g/cm3
314  a0 = 0.05639;
315  // wl0 = 84.22;
316  rnscale = sqrt(1. + (a0 * refwl * refwl) / (refwl * refwl - a0 * a0));
317  nn = sqrt(1. + (a0 * wl * wl) / (wl * wl - a0 * a0)) * refn / rnscale;
318  break;
319  case 3: // --- experimental points
320  rho = 0.088; // g/cm3
321  nn = aeroN[i] * refn / an0; // scale refractive index
322  break;
323  default:
324  nn = refn;
325  break;
326  }
327 
328  scaledE[i] = ee;
329  scaledN[i] = nn;
330  scaledA[i] = aeroA[i] * (rho * g / cm3) / density; // approx. larger the density, smaller the abs. length
331  scaledS[i] = aeroS[i] * (rho * g / cm3) / density; // approx. larger the density, smaller the abs. length
332  }
333 
334  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
335  {
336  std::cout << "# Aerogel Refractive Index, Absorption and Scattering Lengths rescaled to density ";
337  std::cout << density / g * cm3 << " g/cm3, method: " << mode << std::endl;
338  }
339  setMatPropTable(nEntries);
340 
341  return nEntries;
342  }
343 
344  private:
345  // Quartz (SiO2) Refractive index:
346  // https://refractiveindex.info/?shelf=main&book=SiO2&page=Malitson
347  double riQuartz(double wl)
348  { // wavelength
349  double x = wl / um;
350  double nn = sqrt(1 + 0.6961663 * x * x / (x * x - pow(0.0684043, 2)) +
351  0.4079426 * x * x / (x * x - pow(0.1162414, 2)) +
352  0.8974794 * x * x / (x * x - pow(9.896161, 2)));
353  if (nn < 1.)
354  {
355  std::cout << "# WARNING: estimated quartz refractive index is " << nn;
356  std::cout << "at wavelength " << x << " nm -> set to 1" << std::endl;
357  nn = 1.;
358  }
359  return nn;
360  }
361 
362  // Air Refractive index:
363  // https://refractiveindex.info/?shelf=other&book=air&page=Ciddor
364  double riAir(double wl)
365  { // wavelength
366  double x = wl / um;
367  double nn = 1.0 + (0.05792105 / (238.0185 - 1.0 / x / x) +
368  0.00167917 / (57.362 - 1.0 / x / x));
369  if (nn < 1.)
370  {
371  std::cout << "# WARNING: estimated air refractive index is " << nn;
372  std::cout << "at wavelength " << x << " nm -> set to 1" << std::endl;
373  nn = 1.;
374  }
375  return nn;
376  }
377 
378  /*
379  * n_aer = A * n_air + (1-A) * n_quartz : compute air weight-fraction given a
380  * reference n_air,lambda Vorobiev Model rn : reference refractive index
381  * of the aerogel rlambda: wavelength of the reference refractive index
382  */
383  double airFraction(double rn, double rlambda)
384  {
385  double rnq = riQuartz(rlambda);
386  double rna = riAir(rlambda);
387  double a = (rnq - rn) / (rnq - rna);
388  return a;
389  }
390 
391  // density of aerogel from air and quartz densities
392  // rho = (n-1)/0.21 g/cm3 (Bellunato) -> rho proportional to (n-1)
393  // rho_air = 0.0011939 g/cm3 T=25 deg
394  // rho_quartz = 2.196 g/cm3 (amorphous ?)
395  double Density(double rn, double rlambda)
396  {
397  double aa = airFraction(rn, rlambda);
398  return (aa * 1.1939 * mg / cm3 + (1. - aa) * 2.32 * g / cm3);
399  }
400 
401  // at 400 nm - (CLAS12/RICH EPJ A (2016) 52: 23)
402  double density2refIndex(double rho)
403  {
404  double nn2 = 1. + 0.438 * rho / g * cm3;
405  return sqrt(nn2);
406  }
407 };
408 
409 //
410 // Acrylic Filter
411 //
412 
413 class EICG4dRICHFilter : public EICG4dRICHOptics
414 {
415  public:
416  EICG4dRICHFilter(const G4String matName) : EICG4dRICHOptics(matName, "_NA_"){}
417 
418  // wlthr: threshold wavelength for low pass filter
419  // mode currently not used
420  int setOpticalParams(double wlthr) override
421  {
422  const double acryE[] = // energy : wavelenth 660 nm -> 200 nm
423  {1.87855 * eV, 1.96673 * eV, 2.05490 * eV, 2.14308 * eV, 2.23126 * eV,
424  2.31943 * eV, 2.40761 * eV, 2.49579 * eV, 2.58396 * eV, 2.67214 * eV,
425  2.76032 * eV, 2.84849 * eV, 2.93667 * eV, 3.02485 * eV, 3.11302 * eV,
426  3.20120 * eV, 3.28938 * eV, 3.37755 * eV, 3.46573 * eV, 3.55391 * eV,
427  3.64208 * eV, 3.73026 * eV, 3.81844 * eV, 3.90661 * eV, 3.99479 * eV,
428  4.08297 * eV, 4.17114 * eV, 4.25932 * eV, 4.34750 * eV, 4.43567 * eV,
429  4.52385 * eV, 4.61203 * eV, 4.70020 * eV, 4.78838 * eV, 4.87656 * eV,
430  4.96473 * eV, 5.05291 * eV, 5.14109 * eV, 5.22927 * eV, 5.31744 * eV,
431  5.40562 * eV, 5.49380 * eV, 5.58197 * eV, 5.67015 * eV, 5.75833 * eV,
432  5.84650 * eV, 5.93468 * eV, 6.02286 * eV, 6.11103 * eV, 6.19921 * eV};
433 
434  const double acryN[] = // refractive index
435  {1.4902, 1.4907, 1.4913, 1.4918, 1.4924, 1.4930, 1.4936, 1.4942, 1.4948,
436  1.4954, 1.4960, 1.4965, 1.4971, 1.4977, 1.4983, 1.4991, 1.5002, 1.5017,
437  1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017,
438  1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017,
439  1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017, 1.5017,
440  1.5017, 1.5017, 1.5017, 1.5017, 1.5017};
441 
442  const double acryA[] = // absorption length
443  {14.8495 * cm, 14.8495 * cm, 14.8495 * cm, 14.8495 * cm,
444  14.8495 * cm, 14.8495 * cm, 14.8495 * cm, 14.8495 * cm,
445  14.8495 * cm, 14.8495 * cm, 14.8495 * cm, 14.8495 * cm,
446  14.8495 * cm, 14.8495 * cm, 14.8495 * cm, 14.8495 * cm,
447  14.8495 * cm, 14.8494 * cm, 14.8486 * cm, 14.844 * cm,
448  14.8198 * cm, 14.7023 * cm, 14.1905 * cm, 12.3674 * cm,
449  8.20704 * cm, 3.69138 * cm, 1.33325 * cm, 0.503627 * cm,
450  0.23393 * cm, 0.136177 * cm, 0.0933192 * cm, 0.0708268 * cm,
451  0.0573082 * cm, 0.0483641 * cm, 0.0420282 * cm, 0.0373102 * cm,
452  0.033662 * cm, 0.0307572 * cm, 0.0283899 * cm, 0.0264235 * cm,
453  0.0247641 * cm, 0.0233451 * cm, 0.0221177 * cm, 0.0210456 * cm,
454  0.0201011 * cm, 0.0192627 * cm, 0.0185134 * cm, 0.0178398 * cm,
455  0.0172309 * cm, 0.0166779 * cm};
456 
457  const int nEntries = sizeof(acryE) / sizeof(double);
458 
459  if (scaledE == NULL)
460  {
461  scaledE = new double[nEntries];
462  scaledN = new double[nEntries];
463  scaledS = new double[nEntries];
464  scaledA = new double[nEntries];
465  }
466 
467  double e0 = acryE[24]; // wavelength corresponding to the above data
468  // threshold at about Half Maximum
469  double ethr = wl2e(wlthr);
470 
471  int ithr = -1;
472  double eshift = 0;
473 
474  for (int i = 0; i < (nEntries - 1); i++)
475  { // find closest
476  double d1 = ethr - acryE[i];
477  if (d1 >= 0)
478  { // good bin
479  ithr = i;
480  eshift = d1;
481  }
482  break;
483  }
484  if (ithr == -1)
485  {
486  std::cerr << "# ERROR filter: wavelength threshold " << wlthr / nm << " nm is out of range" << std::endl;
487  return 0;
488  }
489 
490  for (int i = 0; i < nEntries; i++)
491  {
492  scaledE[i] = acryE[i] + eshift; // to maky ethr corresponds to a sampled point
493  scaledN[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryN);
494  scaledA[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryA);
495  scaledS[i] = 100000. * cm; // @@@@
496  }
497 
498  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
499  {
500  std::cout << "# Acrylic Filter Refractive Index, Absorption and Scattering Lengths";
501  std::cout << " rescaled to wavelength threshold " << wlthr / nm << " nm" << std::endl;
502  }
503  setMatPropTable(nEntries);
504 
505  return nEntries;
506  }
507 
508  private:
509 };
510 
511 //
512 // gas
513 //
514 
515 class EICG4dRICHGas : public EICG4dRICHOptics
516 {
517  public:
518  EICG4dRICHGas(const G4String matName) : EICG4dRICHOptics(matName, "_NA_")
519  {
520 
521  int nel = mat->GetNumberOfElements();
522  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
523  {
524  std::cout << "# Gas material number of elements " << nel << std::endl;
525 
526  chemFormula = "";
527 
528  for (int i = 0; i < nel; i++)
529  { // extract chemical formula from gas material
530  auto ele = mat->GetElement(i);
531  std::cout << "# Element " << i;
532  std::cout << " : Z " << ele->GetZ();
533  std::cout << " Name " << ele->GetSymbol().data();
534  std::cout << " Atoms " << mat->GetAtomsVector()[i] << std::endl;
535 
536  chemFormula = chemFormula + ele->GetSymbol() + std::to_string(mat->GetAtomsVector()[i]);
537  }
538  std::cout << "# Chemical Formula : " << chemFormula.data() << std::endl;
539  }
540  }
541 
542  int setOpticalParams() override
543  {
544 
545  // very approximate values
546  /*
547  const double gasE[] =
548  { 2*eV, 2.5*eV, 3*eV, 3.5*eV, 4*eV, 4.5*eV, 5*eV, 5.5*eV, 6*eV, 6.5*eV,
549  7*eV } const double gasN[] = // (n-1)*10^6 { 823., 829., 835., 843., 852.,
550  863., 875., 889., 905., 923., 943. } const double gasA[] = { 100*cm,
551  100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
552  100*cm }
553  */
554  // different gas types parameters
555  G4String gasType[] = {"C2F6", "CF4"};
556 
557  // A.W. Burner and W. K. Goad - Measurement of the Specific Refractivities
558  // of CF4 and C2F6 for gases: n-1 = K*rho : K=specific refractivity or
559  // Gladstone-Dale constant C2F6: rho = 5.7 kg/m^3, K=0.131 cm^3/g +/- 0.0009
560  // cm^3/g at 300 K, lambda=633 nm CF4: rho = 7.2 kg/m^3, K=0.122 cm^3/g +/-
561  // 0.0009 cm^3/g at 300 K, lambda=633 nm
562  double Ksr[] = {0.131 * cm3 / g, 0.122 * cm3 / g};
563 
564  // One term Sellmeier formula: n-1 = A*10^-6 / (l0^-2 - l^-2)
565  // C2F6: A=0.18994, l0 = 65.47 [nm] (wavelength, E0=18.82 eV) : NIMA 354
566  // (1995) 417-418 CF4: A=0.124523, l0=61.88 nm (E0=20.04 eV) : NIMA 292
567  // (1990) 593-594
568  double Asel[] = {0.18994, 0.124523};
569  double L0sel[] = {65.47 * nm, 61.88 * nm};
570 
571  int igas = 0;
572 
573  for (int i = 0; i < 2; i++)
574  {
575  if (chemFormula == gasType[i]) igas = i;
576  }
577 
578  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
579  {
580  std::cout << "# Selected gas index " << igas << " for gas " << chemFormula.data() << std::endl;
581  }
582 
583  double density = mat->GetDensity();
584  double refn = Ksr[igas] * density + 1.;
585  double wlref = 633 * nm; // for density vs refractive index
586 
587  int nEntries = 10;
588  double wl0 = 200. * nm;
589  double wl1 = 700. * nm;
590  double dwl = (wl1 - wl0) / (nEntries - 1.);
591 
592  if (scaledE == NULL)
593  {
594  scaledE = new double[nEntries];
595  scaledN = new double[nEntries];
596  scaledS = new double[nEntries];
597  scaledA = new double[nEntries];
598  }
599 
600  double l02 = 1. / (L0sel[igas] / nm) / (L0sel[igas] / nm);
601 
602  double rnscale = Asel[igas] / 1e6 / (l02 - 1. / (wlref / nm) / (wlref / nm)) + 1.;
603 
604  for (int i = 0; i < nEntries; i++)
605  {
606  double wl = wl1 - i * dwl; // to get increasing energy
607  double ee = wl2e(wl);
608 
609  scaledE[i] = ee;
610  scaledN[i] = (Asel[igas] / 1e6 / (l02 - 1. / (wl / nm) / (wl / nm)) + 1.) * refn / rnscale;
611  scaledA[i] = 100. * cm; // @@@@
612  scaledS[i] = 100000. * cm; // @@@@
613  }
614 
615  if (Verbosity() >= Fun4AllBase::VERBOSITY_MORE)
616  {
617  std::cout << "# Gas Refractive Index, Absorption and Scattering Lengths rescaled ";
618  std::cout << "to density " << density / g * cm3 << " g/cm3, gas index: " << igas << std::endl;
619  }
620  setMatPropTable(nEntries);
621 
622  return nEntries;
623  }
624 
625  private:
626  G4String chemFormula; // chemical formula of the gas material
627 };
628 
629 //
630 // Mirror
631 //
632 
633 class EICG4dRICHMirror : public EICG4dRICHOptics
634 {
635  public:
636  EICG4dRICHMirror(const G4String logName) : EICG4dRICHOptics("_NA_", logName){}
637 
638  // pSurfName: prefix used to generate surface names inserted in optical table
639 
640  int setOpticalParams(G4String pSurfName) override
641  {
642 
643  G4String surfaceName = pSurfName + "mirrorSurf";
644  G4String skinSurfaceName = pSurfName + "mirrorSkinSurf";
645 
646  const double mirrorE[] = {
647  2.04358 * eV, 2.0664 * eV, 2.09046 * eV, 2.14023 * eV, 2.16601 * eV,
648  2.20587 * eV, 2.23327 * eV, 2.26137 * eV, 2.31972 * eV, 2.35005 * eV,
649  2.38116 * eV, 2.41313 * eV, 2.44598 * eV, 2.47968 * eV, 2.53081 * eV,
650  2.58354 * eV, 2.6194 * eV, 2.69589 * eV, 2.73515 * eV, 2.79685 * eV,
651  2.86139 * eV, 2.95271 * eV, 3.04884 * eV, 3.12665 * eV, 3.2393 * eV,
652  3.39218 * eV, 3.52508 * eV, 3.66893 * eV, 3.82396 * eV, 3.99949 * eV,
653  4.13281 * eV, 4.27679 * eV, 4.48244 * eV, 4.65057 * eV, 4.89476 * eV,
654  5.02774 * eV, 5.16816 * eV, 5.31437 * eV, 5.63821 * eV, 5.90401 * eV,
655  6.19921};
656 
657  const double mirrorR[] = // Reflectivity of AlMgF2 coated on thermally shaped acrylic
658  // sheets, measured by AJRP, 10/01/2012:
659  {0.8678125, 0.8651562, 0.8639063, 0.8637500, 0.8640625, 0.8645313,
660  0.8643750, 0.8656250, 0.8653125, 0.8650000, 0.8648437, 0.8638281,
661  0.8635156, 0.8631250, 0.8621875, 0.8617188, 0.8613281, 0.8610156,
662  0.8610938, 0.8616016, 0.8623047, 0.8637500, 0.8655859, 0.8673828,
663  0.8700586, 0.8741992, 0.8781055, 0.8825195, 0.8876172, 0.8937207,
664  0.8981836, 0.9027441, 0.9078369, 0.9102002, 0.9093164, 0.9061743,
665  0.9004223, 0.8915210, 0.8599536, 0.8208313, 0.7625024};
666 
667  const int nEntries = sizeof(mirrorE) / sizeof(double);
668 
669  if (scaledE == NULL)
670  {
671  scaledE = new double[nEntries];
672  scaledSR = new double[nEntries];
673  }
674 
675  for (int i = 0; i < nEntries; i++)
676  {
677  scaledE[i] = mirrorE[i];
678  scaledSR[i] = mirrorR[i];
679  }
680 
681  G4MaterialPropertiesTable *pT = addSkinPropTable(nEntries);
682 
683  G4OpticalSurface *pOps = new G4OpticalSurface(surfaceName, unified, polishedfrontpainted, dielectric_dielectric);
684  pOps->SetMaterialPropertiesTable(pT);
685 
686  new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
687 
688  return nEntries;
689 
690  /* from original Alessio GEMC code:
691  $mir{"name"} = "spherical_mirror";
692  $mir{"description"} = "reflective mirrors for eic rich";
693  $mir{"type"} = "dielectric_dielectric";
694  $mir{"finish"} = "polishedfrontpainted";
695  $mir{"model"} = "unified";
696  $mir{"border"} = "SkinSurface";
697  $mir{"photonEnergy"} = arrayToString(@PhotonEnergyBin1);
698  $mir{"reflectivity"} = arrayToString(@Reflectivity1);
699  print_mir(\%configuration, \%mir);
700  */
701  }
702 };
703 
704 //
705 // photo sensor
706 //
707 
708 class EICG4dRICHPhotosensor : public EICG4dRICHOptics
709 {
710  public:
711  EICG4dRICHPhotosensor(const G4String logName) : EICG4dRICHOptics("_NA_", logName){}
712 
713  // pSurfName: prefix used to generate surface names inserted in optical table
714 
715  int setOpticalParams(G4String pSurfName) override
716  {
717 
718  G4String surfaceName = pSurfName + "phseSurf";
719  G4String skinSurfaceName = pSurfName + "phseSkinSurf";
720 
721  double E[] = {1. * eV, 4. * eV, 7. * eV};
722  double SE[] = {1.0, 1.0, 1.0};
723  double N[] = {1.92, 1.92, 1.92};
724  double IN[] = {1.69, 1.69, 1.69};
725 
726  scaledE = new double[3];
727  scaledSE = new double[3];
728  scaledN = new double[3];
729  scaledIN = new double[3];
730 
731  for (int i = 0; i < 3; i++)
732  {
733  scaledE[i] = E[i];
734  scaledSE[i] = SE[i];
735  scaledN[i] = N[i];
736  scaledIN[i] = IN[i];
737  }
738 
739  G4MaterialPropertiesTable *pT = addSkinPropTable(3);
740 
741  G4OpticalSurface *pOps = new G4OpticalSurface(surfaceName, glisur, polished, dielectric_metal);
742  pOps->SetMaterialPropertiesTable(pT);
743 
744  new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
745 
746  return 2;
747  }
748 };
749 
750 #endif // G4E_CI_DRICH_MODEL_HH