drich-dev/scripts/extrapolate_material_tables.py

0001 #!/usr/bin/env python
0002 # Copyright 2023, Christopher Dilks
0003 # Subject to the terms in the LICENSE file found in the top-level directory.
0004
0005 # extrapolate some of the material property tables to a broader range
0006
0007 import ROOT as r
0008 import sys
0009 from numpy import linspace
0010
0011 ##### SETTINGS #########################
0012 FULL_WAVELENGTH_RANGE = [200, 1000]
0013 ########################################
0014
0015 root_file = r.TFile.Open('out/optical_materials_drich.root', 'READ')
0016
0017 # Material Propery Table class, to hold graphs, fit function, etc.
0018 # --------------------------------------------------------------------------
0019
0020 class MPT:
0021     def __init__(self, graph, func, fit_range, extrap_range, extrap_npoints, units='', sigfigs=[5,9,5]):
0022         self.graph          = graph
0023         self.func           = func
0024         self.fit_range      = fit_range
0025         self.extrap_range   = extrap_range
0026         self.extrap_npoints = extrap_npoints
0027         self.units          = units
0028         self.sigfigs        = sigfigs
0029         self.graph_range = [
0030             self.graph.GetPointX(0),
0031             self.graph.GetPointX(self.graph.GetN()-1)
0032         ]
0033
0034     def extrap(self):
0035         # perform the fit
0036         if self.func is not None:
0037             print(f'Fit graph "{self.graph.GetName()}" to function:')
0038             self.func.Print()
0039             self.graph.Fit(self.func, '', '', self.fit_range[0], self.fit_range[1])
0040         # extrapolate
0041         self.multi_gr = r.TMultiGraph()
0042         self.multi_gr.SetName(self.graph.GetName()+"_multi_gr")
0043         self.multi_gr.SetTitle(self.graph.GetTitle())
0044         self.multi_gr.Add(self.graph)
0045         self.graph_extrap = [ r.TGraphErrors(), r.TGraphErrors() ]
0046         self.graph_extrap[0].SetName(self.graph.GetName()+"_extrap_low")
0047         self.graph_extrap[1].SetName(self.graph.GetName()+"_extrap_high")
0048         self.graph_extrap[0].SetMarkerColor(r.kRed)
0049         self.graph_extrap[1].SetMarkerColor(r.kGreen+1)
0050         if self.func is not None:
0051             for i in range(2):
0052                 self.graph_extrap[i].SetTitle(self.graph.GetTitle())
0053                 self.graph_extrap[i].SetMarkerStyle(r.kStar)
0054                 if( ( i==0 and self.extrap_range[i]<self.graph_range[i] ) or ( i==1 and self.extrap_range[i]>self.graph_range[i] )):
0055                     self.multi_gr.Add(self.graph_extrap[i])
0056                     extrap_points = list(linspace(self.extrap_range[i], self.graph_range[i], self.extrap_npoints[i]+1))
0057                     del extrap_points[-1]
0058                     if( i==1 ):
0059                         extrap_points.reverse()
0060                     for x in extrap_points:
0061                         self.graph_extrap[i].AddPoint(x, self.func.Eval(x))
0062         # draw the results
0063         canv_name = f'{self.graph.GetName()}_canv'
0064         self.canv = r.TCanvas(canv_name, canv_name, 1000, 800)
0065         self.canv.SetGrid(1,1)
0066         self.multi_gr.Draw('APE')
0067         if self.func is not None:
0068             self.func.Draw('SAME')
0069         # print the table
0070         self.table = []
0071         for gr in [ self.graph_extrap[0], self.graph, self.graph_extrap[1] ]:
0072             for i in range(gr.GetN()):
0073                 energy = 1239.841875 / gr.GetPointX(i)
0074                 val = gr.GetPointY(i)
0075                 line = f'      {energy:<.{self.sigfigs[0]}f}*eV  {val:>{self.sigfigs[1]}.{self.sigfigs[2]}f}'
0076                 if(self.units!=''):
0077                     line += f'*{self.units}'
0078                 self.table.append(line)
0079         self.table.reverse()
0080         print('-'*80)
0081         print(f'TABLE: {self.graph.GetName()}')
0082         print('-'*80)
0083         for line in self.table:
0084             print(line)
0085         print('-'*80)
0086
0087     def set_fake_errors(self, err):
0088         for i in range(self.graph.GetN()):
0089             self.graph.SetPointError(i, 0.0, err)
0090
0091 tabs = {}
0092
0093 # fit function formula factories
0094 # --------------------------------------------------------------------------
0095
0096 def make_chebyshev(order):
0097     chebyshev = {
0098             0: "1",
0099             1: "x",
0100             2: "2*x^2-1",
0101             3: "4*x^3-3*x",
0102             4: "8*x^4-8*x^2+1",
0103             5: "16*x^5-20*x^3+5*x",
0104             6: "32*x^6-48*x^4+18*x^2-1",
0105             7: "64*x^7-112*x^5+56*x^3-7*x",
0106             }
0107     formula = []
0108     for i in range(order+1):
0109         formula.append(f'[{i}]*({chebyshev[i]})')
0110     return '+'.join(formula)
0111
0112 def make_sellmeier(order):
0113     formula = ["1"]
0114     p = 0
0115     for i in range(order):
0116         formula.append(f'[{p}]*x^2/(x^2-[{p+1}]^2)')
0117         p += 2
0118     return f'sqrt({"+".join(formula)})'
0119
0120 # replacements: replace any tables from `g4dRIChOptics`
0121 # --------------------------------------------------------------------------
0122
0123 ### aerogel ABSLENGTH: from transmittance analysis of Aerogel Factory samples
0124 ###                    used in the dRICH prototype
0125 aerogel_abslength_update_table = [  # [nm] [cm]
0126     [ 890, 58.661475  ],
0127     [ 880, 58.6551    ],
0128     [ 870, 58.64805   ],
0129     [ 860, 58.640225  ],
0130     [ 850, 58.631525  ],
0131     [ 840, 58.6219    ],
0132     [ 830, 58.611125  ],
0133     [ 820, 58.599175  ],
0134     [ 810, 58.585825  ],
0135     [ 800, 58.570925  ],
0136     [ 790, 58.554275  ],
0137     [ 780, 58.535575  ],
0138     [ 770, 58.51465   ],
0139     [ 760, 58.49115   ],
0140     [ 750, 58.464675  ],
0141     [ 740, 58.4349    ],
0142     [ 730, 58.4013    ],
0143     [ 720, 58.363375  ],
0144     [ 710, 58.320375  ],
0145     [ 700, 58.27175   ],
0146     [ 690, 58.216525  ],
0147     [ 680, 58.153725  ],
0148     [ 670, 58.0823    ],
0149     [ 660, 58.0008    ],
0150     [ 650, 57.907675  ],
0151     [ 640, 57.801175  ],
0152     [ 630, 57.67915   ],
0153     [ 620, 57.539075  ],
0154     [ 610, 57.378125  ],
0155     [ 600, 57.19285   ],
0156     [ 590, 56.979225  ],
0157     [ 580, 56.7327    ],
0158     [ 570, 56.447825  ],
0159     [ 560, 56.118275  ],
0160     [ 550, 55.7368    ],
0161     [ 540, 55.294975  ],
0162     [ 530, 54.78305   ],
0163     [ 520, 54.19      ],
0164     [ 510, 53.503475  ],
0165     [ 500, 52.709575  ],
0166     [ 490, 51.79315   ],
0167     [ 480, 50.73805   ],
0168     [ 470, 49.52745   ],
0169     [ 460, 48.14475   ],
0170     [ 450, 46.574425  ],
0171     [ 440, 44.8035    ],
0172     [ 430, 42.8232    ],
0173     [ 420, 40.630825  ],
0174     [ 410, 38.23185   ],
0175     [ 400, 35.641575  ],
0176     [ 390, 32.886125  ],
0177     [ 380, 30.00305   ],
0178     [ 370, 27.039875  ],
0179     [ 360, 24.05185   ],
0180     [ 350, 21.098575  ],
0181     [ 340, 18.239375  ],
0182     [ 330, 15.529475  ],
0183     [ 320, 13.0157    ],
0184     [ 310, 10.7339625 ],
0185     [ 300, 8.7075725  ],
0186     [ 290, 6.9466825  ],
0187     [ 280, 5.44934    ],
0188     [ 270, 4.2030425  ],
0189     [ 260, 3.1872325  ],
0190     [ 250, 2.3760525  ],
0191     [ 240, 1.7410525  ],
0192     [ 230, 1.2535325  ],
0193     [ 220, 0.88632675 ],
0194     [ 210, 0.61495675 ],
0195     [ 200, 0.4182277  ],
0196 ]
0197 aerogel_abslength_update_table.reverse()
0198 aerogel_abslength_update_graph = r.TGraphErrors()
0199 aerogel_abslength_update_graph.SetName("graph_Aerogel_ABSLENGTH__updated")
0200 aerogel_abslength_update_graph.SetTitle("Aerogel ABSLENGTH [cm]")
0201 aerogel_abslength_update_graph.SetMarkerColor(r.kBlue)
0202 aerogel_abslength_update_graph.SetMarkerStyle(r.kFullCircle)
0203 for wl, a in aerogel_abslength_update_table:
0204     aerogel_abslength_update_graph.AddPoint(wl, a)
0205
0206 # fits
0207 # --------------------------------------------------------------------------
0208
0209 ### aerogel - RINDEX
0210 ### - fit to 2nd order Sellmeier function
0211 tabs['aerogel'] = {}
0212 aerogel_rindex_fn = r.TF1("aerogel_rindex", make_sellmeier(2), *FULL_WAVELENGTH_RANGE)
0213 aerogel_rindex_fn.SetParLimits(1,0.01,400)
0214 tabs['aerogel']['rindex'] = MPT(
0215         root_file.Get('graph_Aerogel_RINDEX'),
0216         aerogel_rindex_fn,
0217         [200, 650],
0218         FULL_WAVELENGTH_RANGE,
0219         [10, 10],
0220         '',
0221         [5, 7, 5]
0222         )
0223 tabs['aerogel']['rindex'].set_fake_errors(3e-5)
0224
0225 ### aerogel - ABSLENGTH (OLD, CLAS12-based)
0226 ### - linear fit to 350 nm and above only
0227 # tabs['aerogel']['abslength_old'] = MPT(
0228 #         root_file.Get('graph_Aerogel_ABSLENGTH'),
0229 #         r.TF1("aerogel_abslength", "[0]+[1]*x", 350, FULL_WAVELENGTH_RANGE[-1]),
0230 #         [350, 600],
0231 #         FULL_WAVELENGTH_RANGE,
0232 #         [0, 10],
0233 #         'mm',
0234 #         [5, 7, 3]
0235 #         )
0236 # tabs['aerogel']['abslength_old'].set_fake_errors(0.5)
0237
0238 ### aerogel - ABSLENGTH (UPDATED, from above)
0239 ### - linear fit to the last few points
0240 tabs['aerogel']['abslength'] = MPT(
0241         aerogel_abslength_update_graph,
0242         r.TF1("aerogel_abslength", "[0]+[1]*x", 870, FULL_WAVELENGTH_RANGE[-1]),
0243         [870, 890],
0244         FULL_WAVELENGTH_RANGE,
0245         [0, 11],
0246         'cm',
0247         [5, 7, 3]
0248         )
0249 tabs['aerogel']['abslength'].set_fake_errors(0.5)
0250
0251 ### aerogel - RAYLEIGH
0252 ### - fit to lambda^4 dependence
0253 aerogel_rayleigh_fn = r.TF1("aerogel_rayleigh", "[0]+[1]*x^4", *FULL_WAVELENGTH_RANGE)
0254 tabs['aerogel']['rayleigh'] = MPT(
0255         root_file.Get('graph_Aerogel_RAYLEIGH'),
0256         aerogel_rayleigh_fn,
0257         [350, 600],
0258         FULL_WAVELENGTH_RANGE,
0259         [0, 10],
0260         'mm',
0261         [5, 7, 3]
0262         )
0263 tabs['aerogel']['rayleigh'].set_fake_errors(4)
0264
0265
0266 # extrapolate
0267 # --------------------------------------------------------------------------
0268 for obj_name, obj in tabs.items():
0269     for tab_name, tab in obj.items():
0270         tab.extrap()