drich-dev/src/g4dRIChOptics.hh

0001 #ifndef G4E_CI_DRICH_MODEL_HH
0002 #define G4E_CI_DRICH_MODEL_HH
0003
0004 /*
0005  * g4dRIChOptics class hierarchy
0006  * -----------------------------
0007  * original authors: E. Cisbani, A. Del Dotto, C. Fanelli
0008  * source: git@github.com:cisbani/dRICh.git
0009  * -> adapted for usage in EPIC
0010  */
0011
0012 #include <iostream>
0013 #include <functional>
0014 #include <vector>
0015 #include <algorithm>
0016 #include "Geant4/G4Material.hh"
0017 #include "Geant4/G4SystemOfUnits.hh"
0018 #include "Geant4/G4MaterialPropertiesTable.hh"
0019 #include "Geant4/G4OpticalSurface.hh"
0020 #include "Geant4/G4LogicalSkinSurface.hh"
0021 #include "Geant4/G4tgbMaterialMgr.hh"
0022 #include "Geant4/G4tgbVolumeMgr.hh"
0023
0024 /*
0025  * Service Classes
0026  */
0027
0028 //
0029 // Generic optical parameters class
0030 //
0031
0032 class g4dRIChOptics {
0033
0034 public:
0035
0036   double *scaledE; // photon energies
0037   double *scaledN; // real refractive index
0038   double *scaledA; // absorption length
0039   double *scaledS; // scattering length
0040
0041   double *scaledSE; // surface efficiency
0042   double *scaledSR;  // surface reflectivity
0043   double *scaledIN; // imaginary refractive index
0044
0045   // Add optical properties either to material or volume surface
0046   // matName: material name
0047   // logVolName: logical volume name
0048   // if name is "_NA_", properties are not applied to the corresponding component
0049   g4dRIChOptics(const G4String matName, const G4String logVolName) {
0050
0051     printf("#=======================================================================\n");
0052     printf("# Set Optical Properties\n");
0053
0054     materialName = matName;
0055     logicalVName = logVolName;
0056
0057     scaledE=NULL;
0058     scaledN=NULL;
0059     scaledA=NULL;
0060     scaledS=NULL;
0061
0062     scaledSE = NULL;
0063     scaledSR = NULL;
0064     scaledIN = NULL;
0065
0066     if (matName != "_NA_") {
0067
0068       printf("# Material %s\n",matName.data());
0069
0070       mat = G4tgbMaterialMgr::GetInstance()->FindBuiltG4Material(matName);
0071
0072       if (mat == NULL) {
0073     pTable = NULL;
0074     printf("# ERROR: Cannot retrieve %s material in ci_DRICH\n",matName.data());
0075     // handle error
0076       } else {
0077     pTable = mat->GetMaterialPropertiesTable();
0078       }
0079
0080       if (pTable == NULL) {
0081     printf("# No properties table available for %s, allocated a new one\n",matName.data());
0082     pTable = new G4MaterialPropertiesTable();
0083       } else {
0084     pTable->DumpTable();
0085       }
0086     }
0087
0088     if (logVolName != "_NA_") {
0089       printf("# Logical Volume %s\n",logVolName.data());
0090       logVolume = G4tgbVolumeMgr::GetInstance()->FindG4LogVol(logVolName, 0);
0091       if (logVolume == NULL) {
0092     printf("# ERROR: Cannot retrieve %s logical volume in ci_DRICH\n",logVolName.data());
0093     // handle error
0094       }
0095     }
0096
0097   };
0098
0099   ~g4dRIChOptics() {
0100
0101     if (scaledE !=NULL) delete[] scaledE;
0102     if (scaledN !=NULL) delete[] scaledN;
0103     if (scaledA !=NULL) delete[] scaledA;
0104     if (scaledS !=NULL) delete[] scaledS;
0105
0106     if (scaledSE !=NULL) delete[] scaledSE;
0107     if (scaledSR !=NULL) delete[] scaledSR;
0108     if (scaledIN !=NULL) delete[] scaledIN;
0109
0110   };
0111
0112   // dvalue, ivalue, svalue may represent different quantities depending on implementation
0113   virtual int setOpticalParams() { return -1; };
0114   virtual int setOpticalParams(double dvalue) { return -1; };
0115   virtual int setOpticalParams(int ivalue) { return -1; };
0116   virtual int setOpticalParams(int ivalue, double dvalue) { return -1; };
0117   virtual int setOpticalParams(G4String svalue) { return -1; };
0118
0119   // accessors
0120   G4String getMaterialName() { return materialName; };
0121   G4String getLogicalVName() { return logicalVName; };
0122   G4OpticalSurface *getSurface() { return pOps; };
0123
0124   // get the appropriate material property table
0125   G4MaterialPropertiesTable *getMaterialPropertyTable() {
0126     auto t = pTable;
0127     if(t==nullptr && pOps!=nullptr) t = pOps->GetMaterialPropertiesTable();
0128     if(t==nullptr) fmt::print(stderr,"ERROR: cannot find material property table for ({},{})\n",materialName,logicalVName);
0129     return t;
0130   }
0131
0132   // get list of material property tables (that are in use)
0133   std::vector<G4String> getMaterialPropertyNames() {
0134     std::vector<G4String> validProps;
0135     auto tab = getMaterialPropertyTable();
0136     if(tab!=nullptr) {
0137       for(const auto& propName : tab->GetMaterialPropertyNames()) {
0138         if(tab->GetProperty(propName)!=nullptr) validProps.push_back(propName);
0139       }
0140     }
0141     return validProps;
0142   }
0143
0144   // get the size of the material property table
0145   int getMaterialPropertyTableSize() { return int(getMaterialPropertyNames().size()); }
0146
0147   // execute lambda function `block(energy,value)` for each entry (energy,value) in the material property table
0148   void loopMaterialPropertyTable(G4String propName, std::function<void(G4double,G4double)> block, bool reverseOrder=false) {
0149     auto tab = getMaterialPropertyTable();
0150     if(tab==nullptr) return;
0151     auto prop = tab->GetProperty(propName);
0152     if(prop==nullptr) return;
0153     if(reverseOrder) {
0154       for(int i=prop->GetVectorLength()-1; i>=0; i--) {
0155         auto energy = prop->Energy(i);
0156         auto value  = prop->operator[](i);
0157         block(energy,value);
0158       }
0159     } else {
0160       for(int i=0; i<prop->GetVectorLength(); i++) {
0161         auto energy = prop->Energy(i);
0162         auto value  = prop->operator[](i);
0163         block(energy,value);
0164       }
0165     }
0166   }
0167
0168   // converters
0169   static double wl2e(double wl) { // wavelength to energy
0170     return 1239.84193 * eV / (wl/nm);
0171   };
0172
0173   static double e2wl(double e) { // energy to wavelength
0174     return 1239.84193 *nm / (e / eV);
0175   };
0176
0177
0178 protected:
0179
0180   G4String materialName, logicalVName;
0181   G4Material *mat;
0182   G4OpticalSurface *pOps;
0183   G4LogicalVolume *logVolume; // used for skin surface
0184
0185   G4MaterialPropertiesTable *pTable;
0186
0187   // add properties to the MaterialPropertiesTable and link to material
0188   void setMatPropTable(int nEntries) {
0189
0190     if (scaledN!=NULL) pTable->AddProperty("RINDEX",    scaledE, scaledN, nEntries, false, true); // `true` replaced `SetSpline(true)`
0191     if (scaledA!=NULL) pTable->AddProperty("ABSLENGTH", scaledE, scaledA, nEntries, false, true);
0192     if (scaledS!=NULL) pTable->AddProperty("RAYLEIGH",  scaledE, scaledS, nEntries, false, true);
0193     //    pTable->AddConstProperty("SCINTILLATIONYIELD", 0. / MeV); // @@@ TBC @@@
0194     //    pTable->AddConstProperty("RESOLUTIONSCALE", 1.0); // @@@ TBC @@@
0195
0196     mat->SetMaterialPropertiesTable(pTable);
0197     printf("# Optical Table for material %s with %d points:\n",materialName.data(),nEntries);
0198     pTable->DumpTable();
0199
0200   };
0201
0202   // allocate and add properties to the MaterialPropertiesTable
0203   G4MaterialPropertiesTable *addSkinPropTable(int nE) {
0204
0205     G4MaterialPropertiesTable *pTab = new G4MaterialPropertiesTable();
0206     if (scaledSE !=NULL) pTab->AddProperty("EFFICIENCY", scaledE, scaledSE, nE);
0207     if (scaledSR !=NULL) pTab->AddProperty("REFLECTIVITY", scaledE, scaledSR, nE);
0208     if (scaledN !=NULL) pTab->AddProperty("REALRINDEX", scaledE, scaledN, nE);
0209     if (scaledIN !=NULL) pTab->AddProperty("IMAGINARYRINDEX", scaledE, scaledIN, nE);
0210     printf("# Optical Table for volume %s with %d points:\n",logicalVName.data(),nE);
0211     pTab->DumpTable();
0212
0213     return pTab;
0214
0215   };
0216
0217   // Linear Interpolation method
0218   double linint(double val, int n, const double* x, const double* y) {
0219     if (val<=x[0]) return y[0];
0220     if (val>=x[n-1]) return y[n-1];
0221     for (int i=0;i<(n-1);i++) {
0222       if ((val>=x[i]) && (val<x[i+1])) {
0223     return (y[i+1]-y[i])/(x[i+1]-x[i])*(val-x[i])+y[i];
0224       }
0225     }
0226     return 0.;
0227   };
0228
0229 };
0230
0231 //
0232 // Aerogel
0233 //
0234 class g4dRIChAerogel : public g4dRIChOptics {
0235
0236 public:
0237
0238   g4dRIChAerogel(const G4String matName) : g4dRIChOptics(matName, "_NA_") {};
0239   //
0240   // Compute the refractive index, absorption length, scattering length for different energies points
0241   //
0242   // different methods are available for the refractive index:
0243   //   0 - Vorobiev
0244   //   1 - Sellmeier - CLAS12
0245   //   2 - Sellmeier - LHCB
0246   //   3 - CLAS12 experimental points rescaled by Alessio/GEMC (the same used for the scattering and absorption length)
0247   //
0248   // data are scaled according to the input density of the aerogel
0249   //
0250   int setOpticalParams(int mode) {
0251
0252     const double aeroE[]= // energy : wavelenth 660 nm -> 200 nm
0253       { 1.87855*eV,1.96673*eV,2.05490*eV,2.14308*eV,2.23126*eV, 2.31943*eV,2.40761*eV,2.49579*eV,2.58396*eV,2.67214*eV,
0254     2.76032*eV,2.84849*eV,2.93667*eV,3.02485*eV,3.11302*eV, 3.20120*eV,3.28938*eV,3.37755*eV,3.46573*eV,3.55391*eV,
0255     3.64208*eV,3.73026*eV,3.81844*eV,3.90661*eV,3.99479*eV, 4.08297*eV,4.17114*eV,4.25932*eV,4.34750*eV,4.43567*eV,
0256     4.52385*eV,4.61203*eV,4.70020*eV,4.78838*eV,4.87656*eV, 4.96473*eV,5.05291*eV,5.14109*eV,5.22927*eV,5.31744*eV,
0257     5.40562*eV,5.49380*eV,5.58197*eV,5.67015*eV,5.75833*eV, 5.84650*eV,5.93468*eV,6.02286*eV,6.11103*eV,6.19921*eV };
0258
0259     const double aeroN[]= // refractive index - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088, n(400)~1.019
0260       { 1.01825,1.01829,1.01834,1.01839,1.01844, 1.01849,1.01854,1.01860,1.01866,1.01872,
0261     1.01878,1.01885,1.01892,1.01899,1.01906, 1.01914,1.01921,1.01929,1.01938,1.01946,
0262     1.01955,1.01964,1.01974,1.01983,1.01993, 1.02003,1.02014,1.02025,1.02036,1.02047,
0263     1.02059,1.02071,1.02084,1.02096,1.02109, 1.02123,1.02137,1.02151,1.02166,1.02181,
0264     1.02196,1.02212,1.02228,1.02244,1.02261, 1.02279,1.02297,1.02315,1.02334,1.02354 };
0265
0266     const double aeroA[]= // absorption length - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088
0267       { 17.5000*cm,17.7466*cm,17.9720*cm,18.1789*cm,18.3694*cm, 18.5455*cm,18.7086*cm,18.8602*cm,19.0015*cm,19.1334*cm,
0268     19.2569*cm,19.3728*cm,19.4817*cm,19.5843*cm,19.6810*cm, 19.7725*cm,19.8590*cm,19.9410*cm,20.0188*cm,20.0928*cm,
0269     18.4895*cm,16.0174*cm,13.9223*cm,12.1401*cm,10.6185*cm, 9.3147*cm,8.1940*cm,7.2274*cm,6.3913*cm,5.6659*cm,
0270     5.0347*cm,4.4841*cm,4.0024*cm,3.5801*cm,3.2088*cm, 2.8817*cm,2.5928*cm,2.3372*cm,2.1105*cm,1.9090*cm,
0271     1.7296*cm,1.5696*cm,1.4266*cm,1.2986*cm,1.1837*cm, 1.0806*cm,0.9877*cm,0.9041*cm,0.8286*cm,0.7603*cm };
0272
0273     const double aeroS[] = // rayleigh - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088
0274       { 35.1384*cm, 29.24805*cm, 24.5418*cm, 20.7453*cm, 17.6553*cm, 15.1197*cm, 13.02345*cm, 11.2782*cm, 9.81585*cm, 8.58285*cm,
0275     7.53765*cm, 6.6468*cm, 5.88375*cm, 5.22705*cm, 4.6596*cm, 4.167*cm, 3.73785*cm, 3.36255*cm, 3.03315*cm, 2.7432*cm,
0276     2.48700*cm, 2.26005*cm, 2.05845*cm, 1.87875*cm, 1.71825*cm, 1.57455*cm, 1.44555*cm, 1.3296*cm, 1.2249*cm, 1.1304*cm,
0277     1.04475*cm, 0.9672*cm, 0.89655*cm, 0.83235*cm, 0.77385*cm, 0.7203*cm, 0.67125*cm, 0.6264*cm, 0.58515*cm, 0.54735*cm,
0278     0.51255*cm, 0.48045*cm, 0.45075*cm, 0.4233*cm, 0.39795*cm, 0.37455*cm, 0.3528*cm, 0.33255*cm, 0.3138*cm, 0.29625*cm };
0279
0280     const int nEntries = sizeof(aeroE)/sizeof(double);
0281
0282     double density = mat->GetDensity();
0283     printf("# Aerogel Density : %f g/cm3\n",density/(g/cm3));
0284
0285     double refn = density2refIndex(density); // use a n vs rho formula with provide n at 400 nm
0286     double refwl = 400*nm;
0287
0288     double refee = wl2e(refwl)/eV; // [eV] reference energy
0289     double an0 = linint(refee, nEntries, aeroE, aeroN);
0290     //    double aa0 = linint(refee, nEntries, aeroE, aeroA);
0291     //    double as0 = linint(refee, nEntries, aeroE, aeroS);
0292
0293     double aa;
0294     double nn;
0295     double ri, a0, wl0, rnscale;
0296     double rho = 0.0;
0297
0298     if (scaledE==NULL) {
0299       scaledE = new double[nEntries];
0300       scaledN = new double[nEntries];
0301       scaledS = new double[nEntries];
0302       scaledA = new double[nEntries];
0303     }
0304
0305     for (int i=0;i<nEntries;i++) {
0306       double ee = aeroE[i];
0307       double wl = e2wl(ee)/nm; // from Energy to nm
0308
0309       switch (mode) {
0310       case 0:     // --- Vorobiev model
0311     aa = airFraction(refn, refwl);
0312     nn = aa * riAir(wl) + (1. - aa)*riQuartz(wl);
0313     break;
0314       case 1:     // --- Sellmeier, 1 pole from (CLAS12/RICH EPJ A (2016) 52: 23)
0315     ri = 1.0494; // 400 nm
0316     rho = 0.230; // g/cm3
0317     a0 = 0.09683;
0318     wl0 = 84.13;
0319     rnscale =  sqrt(1.+ (a0*refwl*refwl)/(refwl*refwl-a0*a0));
0320     nn = sqrt(1.+ (a0*wl*wl)/(wl*wl-a0*a0))* refn / rnscale;
0321     break;
0322       case 2:    // --- Sellmeier, 1 pole from T. Bellunato et al. Eur. Phys. J. C52, 759-764 (2007)
0323     ri = 1.03; // 400 nm
0324     rho = 0.149; // g/cm3
0325     a0 = 0.05639;
0326     wl0 = 84.22;
0327     rnscale =  sqrt(1.+ (a0*refwl*refwl)/(refwl*refwl-a0*a0));
0328     nn = sqrt(1.+ (a0*wl*wl)/(wl*wl-a0*a0)) * refn / rnscale;
0329     break;
0330       case 3:    // --- experimental points
0331     rho = 0.088; // g/cm3
0332     nn = aeroN[i] * refn / an0; // scale refractive index
0333     break;
0334       default:
0335     nn = refn;
0336     break;
0337       }
0338
0339       scaledE[i] = ee;
0340       scaledN[i] = nn;
0341       scaledA[i] = aeroA[i] * (rho*g/cm3)/density; // approx. larger the density, smaller the abs. length
0342       scaledS[i] = aeroS[i] * (rho*g/cm3)/density; // approx. larger the density, smaller the abs. length
0343
0344     }
0345
0346     printf("# Aerogel Refractive Index, Absorption and Scattering Lengths rescaled to density %.5f g/cm3, method: %d\n", density/g*cm3, mode);
0347
0348     setMatPropTable(nEntries);
0349
0350     return nEntries;
0351
0352   }
0353
0354 private:
0355
0356   // Quartz (SiO2) Refractive index: https://refractiveindex.info/?shelf=main&book=SiO2&page=Malitson
0357   double riQuartz(double wl) { // wavelength
0358     double x = wl / um;
0359     double nn = sqrt(1+0.6961663*x*x/(x*x-pow(0.0684043,2))+0.4079426*x*x/(x*x-pow(0.1162414,2))+0.8974794*x*x/(x*x-pow(9.896161,2)));
0360     if (nn<1.) {
0361       printf("# WARNING: estimated quartz refractive index is %f at wavelength %f nm -> set to 1\n",nn,x);
0362       nn = 1.;
0363     }
0364     return nn;
0365   }
0366
0367   // Air Refractive index: https://refractiveindex.info/?shelf=other&book=air&page=Ciddor
0368   double riAir(double wl) { // wavelength
0369     double x = wl / um;
0370     double nn = 1.0+(0.05792105/(238.0185-1.0/x/x)+0.00167917/(57.362-1.0/x/x));
0371     if (nn<1.) {
0372       printf("# WARNING: estimated air refractive index is %f at wavelength %f nm -> set to 1\n",nn,x);
0373       nn = 1.;
0374     }
0375     return nn;
0376   }
0377
0378   /*
0379    * n_aer = A * n_air + (1-A) * n_quartz  : compute air weight-fraction given a reference n_air,lambda
0380    *  Vorobiev Model
0381    *    rn     : reference refractive index of the aerogel
0382    *    rlambda: wavelength of the reference refractive index
0383    */
0384   double airFraction(double rn, double rlambda) {
0385     double rnq = riQuartz(rlambda);
0386     double rna = riAir(rlambda);
0387     double a = (rnq - rn)/(rnq - rna);
0388     return a;
0389   }
0390
0391   // density of aerogel from air and quartz densities
0392   // rho = (n-1)/0.21 g/cm3 (Bellunato) -> rho proportional to (n-1)
0393   // rho_air    = 0.0011939 g/cm3   T=25 deg
0394   // rho_quartz = 2.196 g/cm3  (amorphous ?)
0395   double Density(double rn, double rlambda) {
0396     double aa = airFraction(rn, rlambda);
0397     return (aa * 1.1939 * mg/cm3 + (1.- aa ) * 2.32 * g/cm3);
0398   }
0399
0400   // at 400 nm - (CLAS12/RICH EPJ A (2016) 52: 23)
0401   double density2refIndex(double rho) {
0402     double nn2 = 1.+0.438*rho/g*cm3;
0403     return sqrt(nn2);
0404   }
0405
0406 };
0407
0408 //
0409 // Acrylic Filter
0410 //
0411
0412 class g4dRIChFilter : public g4dRIChOptics {
0413
0414 public:
0415
0416   g4dRIChFilter(const G4String matName) : g4dRIChOptics(matName, "_NA_") {};
0417
0418   // wlthr: threshold wavelength for low pass filter
0419   // mode currently not used
0420   int setOpticalParams(double wlthr) {
0421
0422     const double acryE[]= // energy : wavelenth 660 nm -> 200 nm
0423       { 1.87855*eV,1.96673*eV,2.05490*eV,2.14308*eV,2.23126*eV, 2.31943*eV,2.40761*eV,2.49579*eV,2.58396*eV,2.67214*eV,
0424     2.76032*eV,2.84849*eV,2.93667*eV,3.02485*eV,3.11302*eV, 3.20120*eV,3.28938*eV,3.37755*eV,3.46573*eV,3.55391*eV,
0425     3.64208*eV,3.73026*eV,3.81844*eV,3.90661*eV,3.99479*eV, 4.08297*eV,4.17114*eV,4.25932*eV,4.34750*eV,4.43567*eV,
0426     4.52385*eV,4.61203*eV,4.70020*eV,4.78838*eV,4.87656*eV, 4.96473*eV,5.05291*eV,5.14109*eV,5.22927*eV,5.31744*eV,
0427     5.40562*eV,5.49380*eV,5.58197*eV,5.67015*eV,5.75833*eV, 5.84650*eV,5.93468*eV,6.02286*eV,6.11103*eV,6.19921*eV };
0428
0429     const double acryN[]= // refractive index
0430       { 1.4902, 1.4907, 1.4913, 1.4918, 1.4924, 1.4930,  1.4936,  1.4942,  1.4948,  1.4954,
0431     1.4960,  1.4965,  1.4971,  1.4977,  1.4983, 1.4991,  1.5002,  1.5017,  1.5017,  1.5017,
0432     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017,
0433     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017,
0434     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017};
0435
0436     const double acryA[]= // absorption length
0437       { 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm,
0438     14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8495*cm, 14.8494*cm, 14.8486*cm, 14.844*cm,
0439     14.8198*cm, 14.7023*cm, 14.1905*cm, 12.3674*cm, 8.20704*cm, 3.69138*cm, 1.33325*cm, 0.503627*cm, 0.23393*cm, 0.136177*cm,
0440     0.0933192*cm, 0.0708268*cm, 0.0573082*cm, 0.0483641*cm, 0.0420282*cm, 0.0373102*cm, 0.033662*cm, 0.0307572*cm, 0.0283899*cm, 0.0264235*cm,
0441     0.0247641*cm, 0.0233451*cm, 0.0221177*cm, 0.0210456*cm, 0.0201011*cm, 0.0192627*cm, 0.0185134*cm, 0.0178398*cm, 0.0172309*cm, 0.0166779*cm };
0442
0443     const int nEntries = sizeof(acryE)/sizeof(double);
0444
0445     if (scaledE==NULL) {
0446       scaledE = new double[nEntries];
0447       scaledN = new double[nEntries];
0448       scaledS = new double[nEntries];
0449       scaledA = new double[nEntries];
0450     }
0451
0452     double e0 = acryE[24]; // wavelength corresponding to the above data threshold at about Half Maximum
0453     double ethr= wl2e(wlthr);
0454
0455     int ithr=-1;
0456     double eshift=0;
0457
0458     for (int i=0;i<(nEntries-1);i++) { // find closest
0459       double d1 = ethr-acryE[i];
0460       if (d1>=0) { // good bin
0461     ithr = i;
0462     eshift = d1;
0463       }
0464       break;
0465     }
0466     if (ithr==-1) {
0467       fprintf(stderr,"# ERROR filter: wavelength threshold %f nm is out of range\n",wlthr/nm);
0468       return 0;
0469     }
0470
0471     for (int i=0;i<nEntries;i++) {
0472       scaledE[i] = acryE[i]+eshift; // to maky ethr corresponds to a sampled point
0473       scaledN[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryN);
0474       scaledA[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryA);
0475       scaledS[i] = 100000.*cm; // @@@@
0476     }
0477
0478     printf("# Acrylic Filter Refractive Index, Absorption and Scattering Lengths rescaled to wavelength threshold %.1f nm\n", wlthr/nm);
0479
0480     setMatPropTable(nEntries);
0481
0482     return nEntries;
0483
0484   };
0485
0486 private:
0487
0488 };
0489
0490 //
0491 // gas
0492 //
0493
0494 class g4dRIChGas : public g4dRIChOptics {
0495
0496 public:
0497
0498   g4dRIChGas(const G4String matName) : g4dRIChOptics(matName, "_NA_") {
0499
0500     int nel = mat->GetNumberOfElements();
0501     printf("# Gas material number of elements %d\n", nel);
0502
0503     chemFormula="";
0504
0505     for (int i=0;i<nel;i++) {  // extract chemical formula from gas material
0506       auto ele = mat->GetElement(i);
0507       printf("# Element %d : Z %f  Name %s Atoms %d\n",i,ele->GetZ(),ele->GetSymbol().data(),mat->GetAtomsVector()[i]);
0508       chemFormula = chemFormula + ele->GetSymbol() + std::to_string(mat->GetAtomsVector()[i]);
0509     }
0510     printf("# Chemical Formula : %s\n",chemFormula.data());
0511
0512   };
0513
0514   int setOpticalParams() {
0515
0516     // different gas types parameters
0517     G4String gasType[] = { "C2F6", "CF4", "C4F10" };
0518
0519     // absorption lengths
0520     // C2F6 and CF4: assumed to be 10m
0521     // C4F10: 6m, for wavelength 200-700 nm (ATLAS note Simulation of ATLAS Luminosity Monitoring with LUCID)
0522     const double absLength[] = { 30.*m, 10.*m, 6.*m };
0523
0524     // A.W. Burner and W. K. Goad - Measurement of the Specific Refractivities of CF4 and C2F6
0525     // for gases: n-1 = K*rho : K=specific refractivity or Gladstone-Dale constant
0526     // C2F6: rho = 5.7 kg/m^3, K=0.131 cm^3/g +/- 0.0009 cm^3/g at 300 K, lambda=633 nm
0527     // CF4:  rho = 7.2 kg/m^3, K=0.122 cm^3/g +/- 0.0009 cm^3/g at 300 K, lambda=633 nm
0528     //
0529     // C4F10: rho = 9.935 kg/m3 at 25°C and 1 atm // NICNAS file NA/317 (1996), perfluorobutane
0530     //                                            // PFG-5040 3M
0531     //        K = ??? (TODO; density-dependence disabled for now)
0532     double Ksr[]={ 0.131*cm3/g, 0.122*cm3/g, 0.000*cm3/g };
0533
0534     // One term Sellmeier formula: n-1 = A*10^-6 / (l0^-2 - l^-2)
0535     // C2F6: A=0.18994, l0 = 65.47 [nm] (wavelength, E0=18.82 eV)  :  NIMA 354 (1995) 417-418
0536     // CF4: A=0.124523, l0=61.88 nm (E0=20.04 eV) : NIMA 292 (1990) 593-594
0537     // C4F10: A=0.2375, l0=73.63 nm (E0=16.84 eV) : NIMA 510 (2003) 262–272
0538     double Asel[]={0.18994, 0.124523, 0.2375};
0539     double L0sel[]={65.47*nm, 61.88*nm, 73.63*nm};
0540
0541     int igas=0;
0542
0543     for (int i=0;i<3;i++) {
0544       if (chemFormula == gasType[i]) igas=i;
0545     }
0546
0547     printf("# Selected gas index %d for gas %s\n",igas, chemFormula.data());
0548
0549     double density = mat->GetDensity();
0550     double refn = Ksr[igas] * density + 1.;
0551     double wlref = 633*nm; // for density vs refractive index
0552
0553     int nEntries = 16;
0554     double wl0 = 200.*nm;
0555     double wl1 = 1000.*nm;
0556     double dwl = (wl1-wl0)/(nEntries-1.);
0557
0558     if (scaledE==NULL) {
0559       scaledE = new double[nEntries];
0560       scaledN = new double[nEntries];
0561       scaledS = new double[nEntries];
0562       scaledA = new double[nEntries];
0563     }
0564
0565     double l02 = 1./(L0sel[igas]/nm)/(L0sel[igas]/nm);
0566
0567     double rnscale = Asel[igas]/1e6/(l02 - 1./(wlref/nm)/(wlref/nm))+1.;
0568     if(chemFormula=="C4F10") rnscale = refn; // disable density-dependent refractive index
0569
0570     for (int i=0;i<nEntries;i++) {
0571
0572       double wl = wl1 - i*dwl; // to get increasing energy
0573       double ee = wl2e(wl);
0574
0575       scaledE[i]=ee;
0576       scaledN[i]=(Asel[igas]/1e6/(l02 - 1./(wl/nm)/(wl/nm))+1.) * refn/rnscale;
0577       scaledA[i]=absLength[igas];    // @@@@
0578       scaledS[i]=100000.*cm; // @@@@
0579     }
0580
0581     printf("# Gas Refractive Index, Absorption and Scattering Lengths rescaled to density %f g/cm3, gas index: %d\n", density/g*cm3, igas);
0582
0583     setMatPropTable(nEntries);
0584
0585     return nEntries;
0586
0587   };
0588
0589 private:
0590
0591   G4String chemFormula; // chemical formula of the gas material
0592
0593 };
0594
0595 //
0596 // Mirror
0597 //
0598
0599 class g4dRIChMirror : public g4dRIChOptics {
0600
0601 public:
0602   g4dRIChMirror(const G4String logName) : g4dRIChOptics("_NA_", logName) { };
0603
0604   // pSurfName: prefix used to generate surface names inserted in optical table
0605
0606   int setOpticalParams(G4String pSurfName) {
0607
0608     G4String surfaceName = pSurfName + "mirrorSurf";
0609     G4String skinSurfaceName = pSurfName + "mirrorSkinSurf";
0610
0611     const double mirrorE[] =
0612       { 2.04358*eV, 2.0664*eV, 2.09046*eV, 2.14023*eV, 2.16601*eV, 2.20587*eV, 2.23327*eV, 2.26137*eV,
0613     2.31972*eV, 2.35005*eV, 2.38116*eV, 2.41313*eV, 2.44598*eV, 2.47968*eV, 2.53081*eV, 2.58354*eV,
0614     2.6194*eV, 2.69589*eV, 2.73515*eV, 2.79685*eV, 2.86139*eV, 2.95271*eV, 3.04884*eV, 3.12665*eV,
0615     3.2393*eV, 3.39218*eV, 3.52508*eV, 3.66893*eV, 3.82396*eV, 3.99949*eV, 4.13281*eV, 4.27679*eV,
0616     4.48244*eV, 4.65057*eV, 4.89476*eV, 5.02774*eV, 5.16816*eV, 5.31437*eV, 5.63821*eV, 5.90401*eV,
0617     6.19921*eV };
0618
0619     const double mirrorR[]= // Reflectivity of AlMgF2 coated on thermally shaped acrylic sheets, measured by AJRP, 10/01/2012:
0620       { 0.8678125, 0.8651562, 0.8639063, 0.8637500, 0.8640625, 0.8645313, 0.8643750, 0.8656250,
0621     0.8653125, 0.8650000, 0.8648437, 0.8638281, 0.8635156, 0.8631250, 0.8621875, 0.8617188,
0622     0.8613281, 0.8610156, 0.8610938, 0.8616016, 0.8623047, 0.8637500, 0.8655859, 0.8673828,
0623     0.8700586, 0.8741992, 0.8781055, 0.8825195, 0.8876172, 0.8937207, 0.8981836, 0.9027441,
0624     0.9078369, 0.9102002, 0.9093164, 0.9061743, 0.9004223, 0.8915210, 0.8599536, 0.8208313,
0625     0.7625024
0626       };
0627
0628     const int nEntries = sizeof(mirrorE)/sizeof(double);
0629
0630     if (scaledE==NULL) {
0631       scaledE = new double[nEntries];
0632       scaledSR = new double[nEntries];
0633     }
0634
0635     for (int i=0;i<nEntries;i++) {
0636       scaledE[i] = mirrorE[i];
0637       scaledSR[i] = mirrorR[i];
0638     }
0639
0640     G4MaterialPropertiesTable * pT = addSkinPropTable(nEntries);
0641
0642     pOps = new G4OpticalSurface(surfaceName, unified, polishedfrontpainted, dielectric_dielectric); // to be parametrized
0643     pOps->SetMaterialPropertiesTable(pT);
0644     printf("# Surface properties:\n");
0645     pOps->DumpInfo();
0646
0647     new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
0648
0649     return nEntries;
0650
0651     /* from original Alessio GEMC code:
0652       $mir{"name"}         = "spherical_mirror";
0653         $mir{"description"}  = "reflective mirrors for eic rich";
0654         $mir{"type"}         = "dielectric_dielectric";
0655         $mir{"finish"}       = "polishedfrontpainted";
0656         $mir{"model"}        = "unified";
0657         $mir{"border"}       = "SkinSurface";
0658         $mir{"photonEnergy"} = arrayToString(@PhotonEnergyBin1);
0659         $mir{"reflectivity"} = arrayToString(@Reflectivity1);
0660         print_mir(\%configuration, \%mir);
0661     */
0662
0663
0664   };
0665
0666 };
0667
0668 //
0669 // photo sensor
0670 //
0671
0672 class g4dRIChPhotosensor : public g4dRIChOptics {
0673
0674 public:
0675
0676   g4dRIChPhotosensor(const G4String logName) : g4dRIChOptics("_NA_", logName) { };
0677
0678   // pSurfName: prefix used to generate surface names inserted in optical table
0679
0680   int setOpticalParams(G4String pSurfName) {
0681
0682     G4String surfaceName = pSurfName + "phseSurf";
0683     G4String skinSurfaceName = pSurfName + "phseSkinSurf";
0684
0685     // quantum effiency, from SiPM model S13361-3050NE-08
0686     std::vector<std::pair<double,double>> QE = { // wavelength [nm], quantum efficiency
0687       {315*nm, 0.00},
0688       {325*nm, 0.04},
0689       {340*nm, 0.10},
0690       {350*nm, 0.20},
0691       {370*nm, 0.30},
0692       {400*nm, 0.35},
0693       {450*nm, 0.40},
0694       {500*nm, 0.38},
0695       {550*nm, 0.35},
0696       {600*nm, 0.27},
0697       {650*nm, 0.20},
0698       {700*nm, 0.15},
0699       {750*nm, 0.12},
0700       {800*nm, 0.08},
0701       {850*nm, 0.06},
0702       {900*nm, 0.04},
0703       {1000*nm, 0.00}
0704     };
0705     std::reverse(QE.begin(), QE.end()); // order in increasing energy
0706     const int N_POINTS = QE.size();
0707     double E[N_POINTS], SE[N_POINTS], N[N_POINTS], IN[N_POINTS];
0708     int i_QE = 0;
0709     for(auto [w,q] : QE) {
0710       E[i_QE]  = wl2e(w);
0711       SE[i_QE] = q;
0712       N[i_QE]  = 1.92;
0713       IN[i_QE] = 1.69;
0714       i_QE++;
0715     }
0716
0717     scaledE  = new double[N_POINTS];
0718     scaledSE = new double[N_POINTS];
0719     scaledN  = new double[N_POINTS];
0720     scaledIN = new double[N_POINTS];
0721
0722     for (int i=0;i<N_POINTS;i++) {
0723       scaledE[i] = E[i];
0724       scaledSE[i] = SE[i];
0725       scaledN[i] = N[i];
0726       scaledIN[i] = IN[i];
0727     }
0728
0729     G4MaterialPropertiesTable * pT = addSkinPropTable(N_POINTS);
0730
0731     pOps = new G4OpticalSurface(surfaceName, glisur, polished, dielectric_dielectric);
0732     pOps->SetMaterialPropertiesTable(pT);
0733     printf("# Surface properties:\n");
0734     pOps->DumpInfo();
0735
0736     new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
0737
0738     return 2;
0739
0740   };
0741
0742 };
0743
0744 #endif //G4E_CI_DRICH_MODEL_HH
0745