g4irt/include/g4dRIChOptics.hh

0001 #ifndef G4E_CI_DRICH_MODEL_HH
0002 #define G4E_CI_DRICH_MODEL_HH
0003
0004 #include <iostream>
0005 #include <G4PVDivision.hh>
0006 #include "G4RotationMatrix.hh"
0007 #include "G4Material.hh"
0008 #include "G4Element.hh"
0009 #include "G4Color.hh"
0010 #include "G4VisAttributes.hh"
0011 #include "G4SystemOfUnits.hh"
0012 #include "G4MaterialPropertiesTable.hh"
0013 #include "G4OpticalSurface.hh"
0014 #include "G4LogicalSkinSurface.hh"
0015 #include "G4tgbMaterialMgr.hh"
0016 #include "G4tgbVolumeMgr.hh"
0017
0018 //#include <tuning.h>
0019 #include <share.h>
0020
0021 /*
0022  * Service Classes
0023  */
0024
0025 //
0026 // Generic optical parameters class
0027 //
0028
0029 class g4dRIChOptics {
0030
0031 public:
0032
0033   double *scaledE; // photon energies
0034   double *scaledN; // real refractive index
0035   double *scaledA; // absorption length
0036   double *scaledS; // scattering length
0037
0038   double *scaledSE; // surface efficiency
0039   double *scaledSR;  // surface reflectivity
0040   double *scaledIN; // imaginary refractive index
0041
0042   // Add optical properties either to material or volume surface
0043   // matName: material name
0044   // logVolName: logical volume name
0045   // if name is "_NA_", properties are not applied to the corresponding component
0046   g4dRIChOptics(G4Material *mptr): mat(mptr) {
0047     materialName = mptr->GetName();
0048     logicalVName = "_NA_";
0049
0050     scaledE=NULL;
0051     scaledN=NULL;
0052     scaledA=NULL;
0053     scaledS=NULL;
0054
0055     scaledSE = NULL;
0056     scaledSR = NULL;
0057     scaledIN = NULL;
0058
0059     pTable = mat->GetMaterialPropertiesTable();
0060
0061     if (pTable == NULL) {
0062       printf("# No properties table available for %s, allocated a new one\n", mptr->GetName().data());
0063       pTable = new G4MaterialPropertiesTable();
0064     } else {
0065       //+pTable->DumpTable();
0066     }
0067   };
0068   g4dRIChOptics(const G4String matName, const G4String logVolName) {
0069
0070     printf("#=======================================================================\n");
0071     printf("# Set Optical Properties\n");
0072
0073     materialName = matName;
0074     logicalVName = logVolName;
0075
0076     scaledE=NULL;
0077     scaledN=NULL;
0078     scaledA=NULL;
0079     scaledS=NULL;
0080
0081     scaledSE = NULL;
0082     scaledSR = NULL;
0083     scaledIN = NULL;
0084
0085     if (matName != "_NA_") {
0086
0087       printf("# Material %s\n",matName.data());
0088
0089       mat = G4tgbMaterialMgr::GetInstance()->FindBuiltG4Material(matName);
0090
0091       if (mat == NULL) {
0092     pTable = NULL;
0093     printf("# ERROR: Cannot retrieve %s material in ci_DRICH\n",matName.data());
0094     // handle error
0095       } else {
0096     pTable = mat->GetMaterialPropertiesTable();
0097       }
0098
0099       if (pTable == NULL) {
0100     printf("# No properties table available for %s, allocated a new one\n",matName.data());
0101     pTable = new G4MaterialPropertiesTable();
0102       } else {
0103     //+pTable->DumpTable();
0104       }
0105     }
0106
0107     if (logVolName != "_NA_") {
0108       printf("# Logical Volume %s\n",logVolName.data());
0109       logVolume = G4tgbVolumeMgr::GetInstance()->FindG4LogVol(logVolName, 0);
0110       if (logVolume == NULL) {
0111     printf("# ERROR: Cannot retrieve %s logical volume in ci_DRICH\n",logVolName.data());
0112     // handle error
0113       }
0114     }
0115
0116   };
0117
0118   ~g4dRIChOptics() {
0119
0120     if (scaledE !=NULL) delete[] scaledE;
0121     if (scaledN !=NULL) delete[] scaledN;
0122     if (scaledA !=NULL) delete[] scaledA;
0123     if (scaledS !=NULL) delete[] scaledS;
0124
0125     if (scaledSE !=NULL) delete[] scaledSE;
0126     if (scaledSR !=NULL) delete[] scaledSR;
0127     if (scaledIN !=NULL) delete[] scaledIN;
0128
0129   };
0130
0131   // dvalue, ivalue, svalue may represent different quantities depending on implementation
0132   virtual int setOpticalParams() { return -1; };
0133   virtual int setOpticalParams(double dvalue) { return -1; };
0134   virtual int setOpticalParams(int ivalue) { return -1; };
0135   virtual int setOpticalParams(int ivalue, double dvalue) { return -1; };
0136   virtual int setOpticalParams(G4String svalue) { return -1; };
0137
0138   G4MaterialPropertiesTable *pTable;
0139 protected:
0140
0141   G4String materialName, logicalVName;
0142   G4Material *mat;
0143   G4LogicalVolume *logVolume; // used for skin surface
0144
0145   //G4MaterialPropertiesTable *pTable;
0146
0147   double wl2e(double wl) { // wavelength to energy
0148     return 1239.84193 * eV / (wl/nm);
0149   };
0150
0151   double e2wl(double e) { // energy to wavelength
0152     return 1239.84193 *nm / (e / eV);
0153   };
0154
0155   // add properties to the MaterialPropertiesTable and link to material
0156   void setMatPropTable(int nEntries) {
0157
0158     if (scaledN!=NULL) pTable->AddProperty("RINDEX",    scaledE, scaledN, nEntries);//->SetSpline(true);
0159 #ifndef _DISABLE_ABSORPTION_
0160     if (scaledA!=NULL) pTable->AddProperty("ABSLENGTH", scaledE, scaledA, nEntries);//->SetSpline(true);
0161 #endif
0162 #ifndef _DISABLE_RAYLEIGH_SCATTERING_
0163     if (scaledS!=NULL) pTable->AddProperty("RAYLEIGH",  scaledE, scaledS, nEntries);//->SetSpline(true);
0164 #endif
0165     //    pTable->AddConstProperty("SCINTILLATIONYIELD", 0. / MeV); // @@@ TBC @@@
0166     //    pTable->AddConstProperty("RESOLUTIONSCALE", 1.0); // @@@ TBC @@@
0167
0168     mat->SetMaterialPropertiesTable(pTable);
0169     printf("# Optical Table for material %s with %d points:\n",materialName.data(),nEntries);
0170     //+pTable->DumpTable();
0171   };
0172
0173   // allocate and add properties to the MaterialPropertiesTable
0174   G4MaterialPropertiesTable *addSkinPropTable(int nE) {
0175
0176     G4MaterialPropertiesTable *pTab = new G4MaterialPropertiesTable();
0177     if (scaledSE !=NULL) pTab->AddProperty("EFFICIENCY", scaledE, scaledSE, nE);
0178     if (scaledSR !=NULL) pTab->AddProperty("REFLECTIVITY", scaledE, scaledSR, nE);
0179     if (scaledN !=NULL) pTab->AddProperty("REALRINDEX", scaledE, scaledN, nE);
0180     if (scaledIN !=NULL) pTab->AddProperty("IMAGINARYRINDEX", scaledE, scaledIN, nE);
0181     printf("# Optical Table for volume %s with %d points:\n",logicalVName.data(),nE);
0182     //+pTab->DumpTable();
0183
0184     return pTab;
0185
0186   };
0187
0188   // Linear Interpolation method
0189   double linint(double val, int n, const double* x, const double* y) {
0190     if (val<=x[0]) return y[0];
0191     if (val>=x[n-1]) return y[n-1];
0192     for (int i=0;i<(n-1);i++) {
0193       if ((val>=x[i]) && (val<x[i+1])) {
0194     return (y[i+1]-y[i])/(x[i+1]-x[i])*(val-x[i])+y[i];
0195       }
0196     }
0197     return 0.;
0198   };
0199
0200 };
0201
0202 //
0203 // Aerogel
0204 //
0205 class g4dRIChAerogel : public g4dRIChOptics {
0206
0207 public:
0208
0209   //g4dRIChAerogel(const G4String matName) : g4dRIChOptics(matName, "_NA_") {};
0210   g4dRIChAerogel(G4Material *material) : g4dRIChOptics(material) {};//Name, "_NA_") {};
0211   //
0212   // Compute the refractive index, absorption length, scattering length for different energies points
0213   //
0214   // different methods are available for the refractive index:
0215   //   0 - Vorobiev
0216   //   1 - Sellmeier - CLAS12
0217   //   2 - Sellmeier - LHCB
0218   //   3 - CLAS12 experimental points rescaled by Alessio/GEMC (the same used for the scattering and absorption length)
0219   //
0220   // data are scaled according to the input density of the aerogel
0221   //
0222   int setOpticalParams(int mode) {
0223
0224     const double aeroE[]= // energy : wavelenth 660 nm -> 200 nm
0225       { 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,
0226     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,
0227     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,
0228     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,
0229     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 };
0230
0231     const double aeroN[]= // refractive index - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088, n(400)~1.019
0232       { 1.01825,1.01829,1.01834,1.01839,1.01844, 1.01849,1.01854,1.01860,1.01866,1.01872,
0233     1.01878,1.01885,1.01892,1.01899,1.01906, 1.01914,1.01921,1.01929,1.01938,1.01946,
0234     1.01955,1.01964,1.01974,1.01983,1.01993, 1.02003,1.02014,1.02025,1.02036,1.02047,
0235     1.02059,1.02071,1.02084,1.02096,1.02109, 1.02123,1.02137,1.02151,1.02166,1.02181,
0236     1.02196,1.02212,1.02228,1.02244,1.02261, 1.02279,1.02297,1.02315,1.02334,1.02354 };
0237
0238     const double aeroA[]= // absorption length - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088
0239       { 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,
0240     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,
0241     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,
0242     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,
0243     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 };
0244
0245     const double aeroS[] = // rayleigh - scaled from CLAS12 RICH to n(300) = 1.02, rho ~ 0.088
0246       { 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,
0247     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,
0248     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,
0249     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,
0250     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 };
0251
0252     const int nEntries = sizeof(aeroE)/sizeof(double);
0253
0254     double density = mat->GetDensity();
0255     printf("# Aerogel Density : %f g/cm3\n",density/(g/cm3));
0256
0257     double refn = density2refIndex(density); // use a n vs rho formula with provide n at 400 nm
0258     printf("  @@@ --> n: %f\n", refn);
0259     double refwl = 400*nm;
0260
0261     double refee = wl2e(refwl)/eV; // [eV] reference energy
0262     double an0 = linint(refee, nEntries, aeroE, aeroN);
0263     //    double aa0 = linint(refee, nEntries, aeroE, aeroA);
0264     //    double as0 = linint(refee, nEntries, aeroE, aeroS);
0265
0266     double aa;
0267     double nn;
0268     double /*ri,*/ a0, /*wl0,*/ rnscale;
0269     double rho = 0.0;
0270
0271     if (scaledE==NULL) {
0272       scaledE = new double[nEntries];
0273       scaledN = new double[nEntries];
0274       scaledS = new double[nEntries];
0275       scaledA = new double[nEntries];
0276     }
0277
0278     for (int i=0;i<nEntries;i++) {
0279       double ee = aeroE[i];
0280       double wl = e2wl(ee)/nm; // from Energy to nm
0281
0282       switch (mode) {
0283       case 0:     // --- Vorobiev model
0284     aa = airFraction(refn, refwl);
0285     nn = aa * riAir(wl) + (1. - aa)*riQuartz(wl);
0286     break;
0287       case 1:     // --- Sellmeier, 1 pole from (CLAS12/RICH EPJ A (2016) 52: 23)
0288     //+ri = 1.0494; // 400 nm
0289     rho = 0.230; // g/cm3
0290     a0 = 0.09683;
0291     //+wl0 = 84.13;
0292     rnscale =  sqrt(1.+ (a0*refwl*refwl)/(refwl*refwl-a0*a0));
0293     nn = sqrt(1.+ (a0*wl*wl)/(wl*wl-a0*a0))* refn / rnscale;
0294     break;
0295       case 2:    // --- Sellmeier, 1 pole from T. Bellunato et al. Eur. Phys. J. C52, 759-764 (2007)
0296     //+ri = 1.03; // 400 nm
0297     rho = 0.149; // g/cm3
0298     a0 = 0.05639;
0299     //+wl0 = 84.22;
0300     rnscale =  sqrt(1.+ (a0*refwl*refwl)/(refwl*refwl-a0*a0));
0301     nn = sqrt(1.+ (a0*wl*wl)/(wl*wl-a0*a0)) * refn / rnscale;
0302     break;
0303       case 3:    // --- experimental points
0304     rho = 0.088; // g/cm3
0305     nn = aeroN[i] * refn / an0; // scale refractive index
0306     break;
0307       default:
0308     nn = refn;
0309     break;
0310       }
0311
0312       scaledE[i] = ee;
0313       scaledN[i] = nn;
0314       scaledA[i] = aeroA[i] * (rho*g/cm3)/density; // approx. larger the density, smaller the abs. length
0315       scaledS[i] = aeroS[i] * (rho*g/cm3)/density; // approx. larger the density, smaller the abs. length
0316
0317       //printf("      %7.5f*eV  %7.5f\n", ee*1E6, nn);
0318       //printf("      %7.5f*eV  %7.3f*mm\n", ee*1E6, scaledA[i]);
0319       //printf("      %7.5f*eV  %7.3f*mm\n", ee*1E6, scaledS[i]);
0320     }
0321
0322     printf("# Aerogel Refractive Index, Absorption and Scattering Lengths rescaled to density %.5f g/cm3, method: %d\n", density/g*cm3, mode);
0323
0324     setMatPropTable(nEntries);
0325
0326     return nEntries;
0327
0328   }
0329
0330 private:
0331
0332   // Quartz (SiO2) Refractive index: https://refractiveindex.info/?shelf=main&book=SiO2&page=Malitson
0333   double riQuartz(double wl) { // wavelength
0334     double x = wl / um;
0335     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)));
0336     if (nn<1.) {
0337       printf("# WARNING: estimated quartz refractive index is %f at wavelength %f nm -> set to 1\n",nn,x);
0338       nn = 1.;
0339     }
0340     return nn;
0341   }
0342
0343   // Air Refractive index: https://refractiveindex.info/?shelf=other&book=air&page=Ciddor
0344   double riAir(double wl) { // wavelength
0345     double x = wl / um;
0346     double nn = 1.0+(0.05792105/(238.0185-1.0/x/x)+0.00167917/(57.362-1.0/x/x));
0347     if (nn<1.) {
0348       printf("# WARNING: estimated air refractive index is %f at wavelength %f nm -> set to 1\n",nn,x);
0349       nn = 1.;
0350     }
0351     return nn;
0352   }
0353
0354   /*
0355    * n_aer = A * n_air + (1-A) * n_quartz  : compute air weight-fraction given a reference n_air,lambda
0356    *  Vorobiev Model
0357    *    rn     : reference refractive index of the aerogel
0358    *    rlambda: wavelength of the reference refractive index
0359    */
0360   double airFraction(double rn, double rlambda) {
0361     double rnq = riQuartz(rlambda);
0362     double rna = riAir(rlambda);
0363     double a = (rnq - rn)/(rnq - rna);
0364     return a;
0365   }
0366
0367   // density of aerogel from air and quartz densities
0368   // rho = (n-1)/0.21 g/cm3 (Bellunato) -> rho proportional to (n-1)
0369   // rho_air    = 0.0011939 g/cm3   T=25 deg
0370   // rho_quartz = 2.196 g/cm3  (amorphous ?)
0371   double Density(double rn, double rlambda) {
0372     double aa = airFraction(rn, rlambda);
0373     return (aa * 1.1939 * mg/cm3 + (1.- aa ) * 2.32 * g/cm3);
0374   }
0375
0376   // at 400 nm - (CLAS12/RICH EPJ A (2016) 52: 23)
0377   double density2refIndex(double rho) {
0378     double nn2 = 1.+0.438*rho/g*cm3;
0379     return sqrt(nn2);
0380   }
0381
0382 };
0383
0384 //
0385 // Acrylic Filter
0386 //
0387
0388 class g4dRIChFilter : public g4dRIChOptics {
0389
0390 public:
0391
0392   //g4dRIChFilter(const G4String matName) : g4dRIChOptics(matName, "_NA_") {};
0393   g4dRIChFilter(G4Material *material) : g4dRIChOptics(material) {};
0394
0395   // wlthr: threshold wavelength for low pass filter
0396   // mode currently not used
0397   int setOpticalParams(double wlthr) {
0398
0399     const double acryE[]= // energy : wavelenth 660 nm -> 200 nm
0400       { 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,
0401     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,
0402     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,
0403     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,
0404     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 };
0405
0406     const double acryN[]= // refractive index
0407       { 1.4902, 1.4907, 1.4913, 1.4918, 1.4924, 1.4930,  1.4936,  1.4942,  1.4948,  1.4954,
0408     1.4960,  1.4965,  1.4971,  1.4977,  1.4983, 1.4991,  1.5002,  1.5017,  1.5017,  1.5017,
0409     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017,
0410     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017,
0411     1.5017,  1.5017,  1.5017,  1.5017,  1.5017, 1.5017,  1.5017,  1.5017,  1.5017,  1.5017};
0412
0413     const double acryA[]= // absorption length
0414       { 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,
0415     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,
0416     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,
0417     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,
0418     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 };
0419
0420     const int nEntries = sizeof(acryE)/sizeof(double);
0421
0422     if (scaledE==NULL) {
0423       scaledE = new double[nEntries];
0424       scaledN = new double[nEntries];
0425       scaledS = new double[nEntries];
0426       scaledA = new double[nEntries];
0427     }
0428
0429     double e0 = acryE[24]; // wavelength corresponding to the above data threshold at about Half Maximum
0430     double ethr= wl2e(wlthr);
0431
0432     int ithr=-1;
0433     double eshift=0;
0434
0435     for (int i=0;i<(nEntries-1);i++) { // find closest
0436       double d1 = ethr-acryE[i];
0437       if (d1>=0) { // good bin
0438     ithr = i;
0439     eshift = d1;
0440       }
0441       break;
0442     }
0443     if (ithr==-1) {
0444       fprintf(stderr,"# ERROR filter: wavelength threshold %f nm is out of range\n",wlthr/nm);
0445       return 0;
0446     }
0447
0448     for (int i=0;i<nEntries;i++) {
0449       scaledE[i] = acryE[i]+eshift; // to maky ethr corresponds to a sampled point
0450       scaledN[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryN);
0451       scaledA[i] = linint((scaledE[i] - ethr + e0), nEntries, acryE, acryA);
0452       scaledS[i] = 100000.*cm; // @@@@
0453     }
0454
0455     printf("# Acrylic Filter Refractive Index, Absorption and Scattering Lengths rescaled to wavelength threshold %.1f nm\n", wlthr/nm);
0456
0457     setMatPropTable(nEntries);
0458
0459     return nEntries;
0460
0461   };
0462
0463 private:
0464
0465 };
0466
0467 //
0468 // gas
0469 //
0470
0471 class g4dRIChGas : public g4dRIChOptics {
0472
0473 public:
0474
0475   //g4dRIChGas(const G4String matName) : g4dRIChOptics(matName, "_NA_") {
0476   g4dRIChGas(G4Material *material) : g4dRIChOptics(material) {
0477     int nel = mat->GetNumberOfElements();
0478     printf("# Gas material number of elements %d\n", nel);
0479
0480     chemFormula="";
0481
0482     for (int i=0;i<nel;i++) {  // extract chemical formula from gas material
0483       auto ele = mat->GetElement(i);
0484       printf("# Element %d : Z %f  Name %s Atoms %d\n",i,ele->GetZ(),ele->GetSymbol().data(),mat->GetAtomsVector()[i]);
0485       chemFormula = chemFormula + ele->GetSymbol() + std::to_string(mat->GetAtomsVector()[i]);
0486     }
0487     printf("# Chemical Formula : %s\n",chemFormula.data());
0488
0489   };
0490
0491   int setOpticalParams() {
0492
0493     // very approximate values
0494     //const double gasE[] =
0495     //{ 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, 7*eV };
0496     //const double gasN[] = // (n-1)*10^6
0497     //{ 823., 829., 835., 843., 852., 863., 875., 889., 905., 923., 943. };
0498     //const double gasA[] =
0499     //{ 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm };
0500
0501     // different gas types parameters
0502     G4String gasType[] = { "C2F6", "CF4" };
0503
0504     // A.W. Burner and W. K. Goad - Measurement of the Specific Refractivities of CF4 and C2F6
0505     // for gases: n-1 = K*rho : K=specific refractivity or Gladstone-Dale constant
0506     // C2F6: rho = 5.7 kg/m^3, K=0.131 cm^3/g +/- 0.0009 cm^3/g at 300 K, lambda=633 nm
0507     // CF4:  rho = 7.2 kg/m^3, K=0.122 cm^3/g +/- 0.0009 cm^3/g at 300 K, lambda=633 nm
0508     double Ksr[]={ 0.131*cm3/g, 0.122*cm3/g };
0509
0510     // One term Sellmeier formula: n-1 = A*10^-6 / (l0^-2 - l^-2)
0511     // C2F6: A=0.18994, l0 = 65.47 [nm] (wavelength, E0=18.82 eV)  :  NIMA 354 (1995) 417-418
0512     // CF4: A=0.124523, l0=61.88 nm (E0=20.04 eV) : NIMA 292 (1990) 593-594
0513     double Asel[]={0.18994, 0.124523};
0514     double L0sel[]={65.47*nm, 61.88*nm};
0515
0516     int igas=0;
0517
0518     for (int i=0;i<2;i++) {
0519       if (chemFormula == gasType[i]) igas=i;
0520     }
0521
0522     printf("# Selected gas index %d for gas %s\n",igas, chemFormula.data());
0523
0524     double density = mat->GetDensity();
0525     double refn = Ksr[igas] * density + 1.;
0526     double wlref = 633*nm; // for density vs refractive index
0527
0528     int nEntries = 10;
0529     double wl0 = 200.*nm;
0530     double wl1 = 700.*nm;
0531     double dwl = (wl1-wl0)/(nEntries-1.);
0532
0533     if (scaledE==NULL) {
0534       scaledE = new double[nEntries];
0535       scaledN = new double[nEntries];
0536       scaledS = new double[nEntries];
0537       scaledA = new double[nEntries];
0538     }
0539
0540     double l02 = 1./(L0sel[igas]/nm)/(L0sel[igas]/nm);
0541
0542     double rnscale = Asel[igas]/1e6/(l02 - 1./(wlref/nm)/(wlref/nm))+1.;
0543
0544     for (int i=0;i<nEntries;i++) {
0545
0546       double wl = wl1 - i*dwl; // to get increasing energy
0547       double ee = wl2e(wl);
0548
0549       scaledE[i]=ee;
0550       scaledN[i]=(Asel[igas]/1e6/(l02 - 1./(wl/nm)/(wl/nm))+1.) * refn/rnscale;
0551       scaledA[i]=10.*m;    // @@@@
0552       //scaledA[i]=100.0*m;    // @@@@
0553       scaledS[i]=100000.*cm; // @@@@
0554     }
0555
0556     printf("# Gas Refractive Index, Absorption and Scattering Lengths rescaled to density %f g/cm3, gas index: %d\n", density/g*cm3, igas);
0557
0558     setMatPropTable(nEntries);
0559
0560     return nEntries;
0561
0562   };
0563
0564 private:
0565
0566   G4String chemFormula; // chemical formula of the gas material
0567
0568 };
0569
0570 //
0571 // Mirror
0572 //
0573
0574 class g4dRIChMirror : public g4dRIChOptics {
0575
0576 public:
0577   g4dRIChMirror(const G4String logName) : g4dRIChOptics("_NA_", logName) { };
0578
0579   // pSurfName: prefix used to generate surface names inserted in optical table
0580
0581   int setOpticalParams(G4String pSurfName) {
0582
0583     G4String surfaceName = pSurfName + "mirrorSurf";
0584     G4String skinSurfaceName = pSurfName + "mirrorSkinSurf";
0585
0586     const double mirrorE[] =
0587       { 2.04358*eV, 2.0664*eV, 2.09046*eV, 2.14023*eV, 2.16601*eV, 2.20587*eV, 2.23327*eV, 2.26137*eV,
0588     2.31972*eV, 2.35005*eV, 2.38116*eV, 2.41313*eV, 2.44598*eV, 2.47968*eV, 2.53081*eV, 2.58354*eV,
0589     2.6194*eV, 2.69589*eV, 2.73515*eV, 2.79685*eV, 2.86139*eV, 2.95271*eV, 3.04884*eV, 3.12665*eV,
0590     3.2393*eV, 3.39218*eV, 3.52508*eV, 3.66893*eV, 3.82396*eV, 3.99949*eV, 4.13281*eV, 4.27679*eV,
0591     4.48244*eV, 4.65057*eV, 4.89476*eV, 5.02774*eV, 5.16816*eV, 5.31437*eV, 5.63821*eV, 5.90401*eV,
0592     6.19921*eV };
0593
0594     const double mirrorR[]= // Reflectivity of AlMgF2 coated on thermally shaped acrylic sheets, measured by AJRP, 10/01/2012:
0595       { 0.8678125, 0.8651562, 0.8639063, 0.8637500, 0.8640625, 0.8645313, 0.8643750, 0.8656250,
0596     0.8653125, 0.8650000, 0.8648437, 0.8638281, 0.8635156, 0.8631250, 0.8621875, 0.8617188,
0597     0.8613281, 0.8610156, 0.8610938, 0.8616016, 0.8623047, 0.8637500, 0.8655859, 0.8673828,
0598     0.8700586, 0.8741992, 0.8781055, 0.8825195, 0.8876172, 0.8937207, 0.8981836, 0.9027441,
0599     0.9078369, 0.9102002, 0.9093164, 0.9061743, 0.9004223, 0.8915210, 0.8599536, 0.8208313,
0600     0.7625024
0601       };
0602
0603     const int nEntries = sizeof(mirrorE)/sizeof(double);
0604
0605     if (scaledE==NULL) {
0606       scaledE = new double[nEntries];
0607       scaledSR = new double[nEntries];
0608     }
0609
0610     for (int i=0;i<nEntries;i++) {
0611       scaledE[i] = mirrorE[i];
0612       scaledSR[i] = mirrorR[i];
0613     }
0614
0615     G4MaterialPropertiesTable * pT = addSkinPropTable(nEntries);
0616
0617     G4OpticalSurface * pOps = new G4OpticalSurface(surfaceName, unified, polishedfrontpainted, dielectric_dielectric); // to be parametrized
0618     pOps->SetMaterialPropertiesTable(pT);
0619
0620     //++new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
0621     m_MirrorSurface = pOps;
0622
0623     return nEntries;
0624
0625     /* from original Alessio GEMC code:
0626       $mir{"name"}         = "spherical_mirror";
0627         $mir{"description"}  = "reflective mirrors for eic rich";
0628         $mir{"type"}         = "dielectric_dielectric";
0629         $mir{"finish"}       = "polishedfrontpainted";
0630         $mir{"model"}        = "unified";
0631         $mir{"border"}       = "SkinSurface";
0632         $mir{"photonEnergy"} = arrayToString(@PhotonEnergyBin1);
0633         $mir{"reflectivity"} = arrayToString(@Reflectivity1);
0634         print_mir(\%configuration, \%mir);
0635     */
0636
0637   };
0638
0639   G4OpticalSurface *m_MirrorSurface;
0640 };
0641
0642 //
0643 // photo sensor
0644 //
0645
0646 class g4dRIChPhotosensor : public g4dRIChOptics {
0647
0648 public:
0649
0650   g4dRIChPhotosensor(const G4String logName) : g4dRIChOptics("_NA_", logName) { };
0651
0652   // pSurfName: prefix used to generate surface names inserted in optical table
0653
0654   int setOpticalParams(G4String pSurfName) {
0655
0656     G4String surfaceName = pSurfName + "phseSurf";
0657     G4String skinSurfaceName = pSurfName + "phseSkinSurf";
0658
0659     double E[] = {1.*eV, 4.*eV, 7.*eV };
0660     double SE[] = {1.0, 1.0, 1.0 };
0661     double N[] = {1.92, 1.92, 1.92 };
0662     double IN[] = {1.69, 1.69, 1.69 };
0663
0664     scaledE  = new double[3];
0665     scaledSE = new double[3];
0666     scaledN  = new double[3];
0667     scaledIN = new double[3];
0668
0669     for (int i=0;i<3;i++) {
0670       scaledE[i] = E[i];
0671       scaledSE[i] = SE[i];
0672       scaledN[i] = N[i];
0673       scaledIN[i] = IN[i];
0674     }
0675
0676     G4MaterialPropertiesTable * pT = addSkinPropTable(3);
0677
0678     G4OpticalSurface * pOps = new G4OpticalSurface(surfaceName, glisur, polished, dielectric_metal);
0679     pOps->SetMaterialPropertiesTable(pT);
0680
0681     new G4LogicalSkinSurface(skinSurfaceName, logVolume, pOps);
0682
0683     return 2;
0684
0685   };
0686
0687 };
0688
0689 #endif //G4E_CI_DRICH_MODEL_HH
0690