Reports until 15:44, Thursday 31 October 2019
H1 CAL
jeffrey.kissel@LIGO.ORG - posted 15:44, Thursday 31 October 2019 - last comment - 16:27, Friday 01 November 2019(52837)
On the effect (and direction) of PCAL Systematic Error
J. Kissel

Recently, the PCAL team identified that we hadn't updated the model of the force-to-displacement transfer function on the PCAL systems for the new O3 test masses, and in parallel, further refined their estimate of the power-to-force coefficient due to more carefully addressing ADC gains and incorporating the results for many more measurements of the transfer coefficients between the gold standard NIST reference and regularly-in-use end-station references (e.g. the RX and TX PDs) -- see more details in LHO aLOG 52828. 

This has been updated and fixed today, so this will no longer be a problem for O3B. However, this means there's been a systematic error in the PCAL displacement estimate that we've been using as an absolute reference for the IFO calibration for O3A.  

In this aLOG, I interpret the level of difference between "now" and "earlier" as described in LHO aLOG 52828 in terms how how this impacts our model of the IFO response function, and thus the impact on our estimate of h(t) -- i.e. recasting the "now" vs. "earlier" as a systematic error that we report in our transfer function that we typically report as the "uncertainty budget," e.g. LHO aLOG 52727.

In that LHO aLOG 52828, (and the subsequent comment LHO aLOG 52834), we see the definition of the "change in calibrated displacement," xi, as reported is 

           [   NOW               ]
     xi =  [ --------    -    1  ]  * 100                                                      (1)
           [ EARLIER             ]

where NOW is after today when the coefficients were installed (updated to be closer to the correct truth), and EARLIER is prior to today, during all of O3A (which therefore had a level systematic error proportional to xi).

We've dealt with systematic error in the PCAL systems in the past: back in the winter of O2 when the alignment of the PCAL beam decayed to the level some level of time-dependent optical clipping on the integrating sphere. LHO aLOG 34153 gives us the formalism to address such systematic error, so we should recast xi into that formalism, which called the systematic error eps (where, for reasons that will become clear later, I re-write the definition of eps in terms of dL_{PCAL} instead of x_pcal):

     eps is a real number that converts the "apparent" displacement, dL_pcal, (i.e. that which includes a known systematic error) into the real displacement, dL_pcal' (i.e. the "truth" that has no known systematic error):

          dL_pcal' = eps * dL_pcal                                                             (2)

     so, if "NOW" is the real, truth displacement estimate, dL_pcal',  and EARLIER is the "apparent" displacement, dL_pcal, then

                            NOW            TRUTH            dL_pcal'
           (1 + xi/100) =  ---------   =  ---------   =    --------    = eps.                  (3)
                           EARLIER        APPARENT          dL_pcal

     And thus, skipping a few of the steps shown in LHO aLOG 34153 for brevity, the "true" response function, R_pcal', (determined by the NOW, true displacement estimate by PCAL system) is related to the apparent, previously report response function, R_pcal, (determined by the the EARLIER PCAL system which had systematic error) is
      
           R_pcal' = eps * R_pcal = (1 + xi/100) * R_pcal                                      (4)


So how does this relate to the various metrics we have for calibration? 
    - How will this systematic error be included in uncertainty budgets? 
    - What does it mean for error in the amplitude of previously reported h(t)? 
    - What does that mean for the BNS range?

To do this, we have to invoke the nomenclature in T1900169 (which is why I switched from x_pcal to dL_pcal above):

     Equation 23 of T19000169 states that the uncertainty budget (e.g. LHO aLOG 52727)

             R_sample     ~      R_pcal            dL_pcal
             ---------    =   -----------    =   -----------                                   (5)
               R_MAP            R_pyDARM          h_GDS * L 

where R_MAP == R_pyDARM is a single response function calculated from the reference model parameters, and R_sample is the posterior distribution of numerically evaluated response functions based on sampling the previously determined uncertainty (and systematic error) on the model parameters, which should bee equivalent to the response function as measured directly by the PCAL, i.e. R_pcal. Via the math shown in Eq. (1) of T1900169, it can be shown that any ratio of response functions is equivalent to the ratio of displacement, and thus the last equivalence in Eq. (5) here.

But Eq. (5) is generic, and now we need to fold in the systematic error in PCAL; eps, or xi. For O3A, the denominators R_MAP, R_pyDARM, and h_GDS*L are now fixed. So, in order to reflect an error in dL_pcal, or R_pcal in to the ratio that we display for the uncertainty, R_sample / R_MAP, then we should *multiply* that ratio by eps = (1 + 100*xi):

             dL_pcal '      R_pcal '    eps * R_pcal          R_sample                    R_sample 
             ---------  = --------   =  ------------  = eps * --------   = (1 + xi/100) * ---------             (6)
            h_GDS * L      R_pyDARM       R_pyDARM              R_MAP                      R_MAP

which means that the reported h(t), h_GDS, is different from the truth, dL_pcal', by the following uncertainty and systematic error.

           dL_pcal' = eps * (R_sample/R_MAP) * h_GDS * L  = (1 + xi/100) * (R_sample/R_MAP) * h_GDS * L         (7)

or, if we had known these systematic errors and uncertainties ahead of time, we might say that the "real" "truth" strain amplitude, h_GDS', can be computed from the reported strain stored in the frames, h_GDS,

           h_GDS' = eps * (R_sample/R_MAP) * h_GDS  =  (1 + xi/100) * (R_sample/R_MAP) * h_GDS                  (8)

And, assuming no signal, the amplitude spectral density of h_GDS, the sensitivity of the detector, call it ASD_GDS, can be used to define that the "truth" BNS range, r', is related to the "apparent" or "reported" BNS range, r, by the following:
                    { f_max
          r' =  C * |       [ ASD_GDS' ]^{-1} * f^{-7/6}  df
                    } f_min

                    { f_max
              = C * |        [ eps * (R_sample / R_MAP) * ASD_GDS ]^{-1} * f^{-7/6} df
                    } f_min

                 1         { f_max
              = ---- * C * |        [ (R_sample / R_MAP) * ASD_GDS ]^-1 * f^(-7/6) df
                eps        } f_min

                 1                1 
           r' = ----  * r = ------------ * r                                                                    (9)
                eps         (1 + xi/100)

where, for brevity, the integrated (R_sample / R_MAP) has been dropped on the last line, because in the case where we're only considering a scalar systematic error, the integrand of (R_sample / R_MAP) * ASD_GDS is the same for both r' and r. 
This would not be true for a frequency dependent systematic error (either in PCAL or in any other component of R_sample/R_MAP).


So, depending on whether eps is greater or less than 1, or if xi is positive or negative, then 
- the uncertainty budget will have a corresponding multiplicative "offset" from unity magnitude, i.e. the "uncertainty budget" plot of R_sample / R_MAP should show a magnitude wiggling around (1 + eps),
- the reported strain amplitude has been wrong by a multiplicative factor of eps,
- the reported sensitivity, i.e. the ASD of the reported strain has been wrong by a multiplicative factor of eps, and
- the reported range has been wrong by by a multiplicative factor of (1 / eps),

So: 
- if eps > 1, or xi > 0, then 
    - if we don't correct the data, then the median line in the uncertainty budget goes up in magnitude, or, 
    - if we were to correct the live data stream (which we will not for O3A), the strain (and noise) amplitude gets larger, the ASD gets worse (less sensitive), and the range goes down.  
- if eps < 1, or xi < 0, then
    - if we don't correct the data, then the median line in the uncertainty budget goes down in magnitude, or, 
    - if we were to correct the live data stream (which we will not for O3A), the strain (and noise) amplitude gets smaller, the ASD gets better (more sensitive), and the range goes up.  
Comments related to this report
ling.sun@LIGO.ORG - 16:27, Friday 01 November 2019 (52897)

According to Rick, Dripta, and Jeff, the numbers to be applied to O3a are:

                               unc        sys (xi)         sys (eps)
LHO EY RXPD      0.54%      +0.42%          1.0042
LLO EY RXPD      0.54%      +0.31%          1.0031

The RRNom changes are as below. (Note: valid GPS time for this correction is set to be chunk 2+3)

Index: RRNom.py
===================================================================
--- RRNom.py    (revision 8656)
+++ RRNom.py    (working copy)
@@ -172,11 +172,23 @@
 # FIXME: Check clipping in O1/O2 (e.g., alog LHO 36884)
 # clipCorrection = hfile['Clipping/ClipFactor'].value
 # We believe in O3 there's no clipping, so set the correction to 1.
-clipCorrection = 1.0
+# O3a: Jun 11 2019 16:57:18 UTC -- Oct 01 2019 16:23:10 UTC, add Pcal sys error
+# It should apply to all O3a, but we do not correct for the released chunk1 before Jun 11
+if run == 'O3' and args.gpsTime >= 1244307456 and args.gpsTime <= 1253982208:
+    if args.IFO == 'LHO':
+        clipCorrection = 1.0042
+    else:
+        clipCorrection = 1.0031
+else:
+    clipCorrection = 1.0

 # TODO: update the value when more measurements are taken
 # In O3, the optimistic lower limit is 0.5% (G1900395 Page 14)
-PCALSystUnc = 0.0079
+# It should apply to all O3, but we do not correct for the released chunk1 before Jun 11
+if run == 'O3' and args.gpsTime >= 1244307456:
+    PCALSystUnc = 0.0054
+else:
+    PCALSystUnc = 0.0079

 # -------------------------------------------------------
 # Create the nominal model using model file
@@ -670,7 +682,7 @@
     thisRR = responseFunc(thisSenseTheta, thisActUIMTheta, thisActPUMTheta, thisActTSTTheta, \
                           thisCCSyst, thisAUSyst, thisAPSyst, thisATSyst, ff, True, \
                           thisKU, thisKP, thisKT, thisKC, thisFC, thisSenseTheta[2], thisSenseTheta[3], \
-                          'all', thisPcalSystErr, clipCorrection)
+                          'all', thisPcalSystErr, 1.0)

 

The two attached plots shows the comparison before (old.png) and after (new.png) the PCAL corrections. The commands for generating these plots are:

OLD:

python /home/ling.sun/aligocalibration/trunk/Common/pyDARM/RRNom.py --IFO=LHO --HDF5_A_MCMCresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_A_MCMC_20190404.hdf5 --HDF5_A_GPRresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_A_GPR_20190416_meas0327-0925.hdf5 --HDF5_C_MCMCresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_C_MCMC_20190404.hdf5 --HDF5_C_GPRresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_C_GPR_20190416_meas0328-0828_withHFafter0611_nofsQcorr_fmin30Hz_lengthscaleZp5.hdf5 --modelPath=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/params/ --modelFilename=modelparams_H1_20190416 --IFOmodel=modelPars --sampleNumber=1000 --seed=1111 --version=ref --outDir=/home/ling.sun/tmpUnc/  --plot1SigmaUncs --saveSummaries --gpsTime=1244307455

NEW:

python /home/ling.sun/aligocalibration/trunk/Common/pyDARM/RRNom.py --IFO=LHO --HDF5_A_MCMCresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_A_MCMC_20190404.hdf5 --HDF5_A_GPRresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_A_GPR_20190416_meas0327-0925.hdf5 --HDF5_C_MCMCresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_C_MCMC_20190404.hdf5 --HDF5_C_GPRresults=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/Results/Uncertainty/O3_H1_C_GPR_20190416_meas0328-0828_withHFafter0611_nofsQcorr_fmin30Hz_lengthscaleZp5.hdf5 --modelPath=/home/ling.sun/aligocalibration/trunk/Runs/O3/H1/params/ --modelFilename=modelparams_H1_20190416 --IFOmodel=modelPars --sampleNumber=1000 --seed=1111 --version=ref --outDir=/home/ling.sun/tmpUnc/  --plot1SigmaUncs --saveSummaries --gpsTime=1244307456

Images attached to this comment