From 9a024bf4c8991c802ca25dddd11f8a3b73a2bfdd Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Wed, 28 Aug 2013 16:36:46 -0400
Subject: [PATCH] Trace focal plane coordinates to target using COSY
 reconstruction coefficients.   Name of coefficient file in hcana.param   Some
 basic hms setting parameters added to test.database for run 50017   Traceback
 code in THcHalCSpectrometer

---
 examples/DBASE/test.database       |   5 +
 examples/PARAM/hcana.param         |   4 +
 examples/PARAM/hms_recon_coeff.dat | 407 +++++++++++++++++++++++++++++
 examples/PARAM/sos_recon_coeff.dat | 267 +++++++++++++++++++
 src/THcHallCSpectrometer.cxx       | 208 +++++++++++++--
 src/THcHallCSpectrometer.h         |  23 +-
 6 files changed, 895 insertions(+), 19 deletions(-)
 create mode 100644 examples/PARAM/hms_recon_coeff.dat
 create mode 100644 examples/PARAM/sos_recon_coeff.dat

diff --git a/examples/DBASE/test.database b/examples/DBASE/test.database
index 451a266..577f1a6 100644
--- a/examples/DBASE/test.database
+++ b/examples/DBASE/test.database
@@ -2,6 +2,11 @@
 50017
 g_ctp_parm_filename="PARAM/general.param"
 g_decode_map_filename="MAPS/jun04.map"
+# Values taken from Logbook.  Not necessarily official.
+gpbeam    = 4.015
+hpcentral = 2.074
+htheta_lab = 10.99
+hpartmass = 0.13957
 52949
 g_ctp_parm_filename="PARAM/52949/general.param"
 g_decode_map_filename="MAPS/july04.map"
diff --git a/examples/PARAM/hcana.param b/examples/PARAM/hcana.param
index 61d60c6..646ac74 100644
--- a/examples/PARAM/hcana.param
+++ b/examples/PARAM/hcana.param
@@ -15,3 +15,7 @@ hdc_plane_names = "1x1 1y1 1u1 1v1 1y2 1x2 2x1 2y1 2u1 2v1 2y2 2x2"
 hcal_layer_names = "1pr 2ta 3ta 4ta"
 
 hhodo_plane_names = "1x 1y 2x 2y"
+
+# The following were defined in REPLAY.PARAM
+h_recon_coeff_filename =    'PARAM/hms_recon_coeff.dat'  ;hms optics matrix
+s_recon_coeff_filename =    'PARAM/sos_recon_coeff.dat'  ;sos optics matrix
diff --git a/examples/PARAM/hms_recon_coeff.dat b/examples/PARAM/hms_recon_coeff.dat
new file mode 100644
index 0000000..7b1ca30
--- /dev/null
+++ b/examples/PARAM/hms_recon_coeff.dat
@@ -0,0 +1,407 @@
+! recostruction matrix elements for HMS, tune HMS-100
+! <theta y phi delta | nmpq>;(x**n xp**m y**p yp**q)
+! as of 02/20/97, no z offset included: this statement is obsolete
+! created: 06/28/98
+! zeroth order matrix elements from fitting:
+! <theta|> -4.7276368293E-03 rad
+! <y|> -1.1606458055E-04 m
+! <phi|>  4.8280702357E-05 rad 
+! <p|> -3.8725980576E-03 %
+! <theta y phi delta | nmpq>;(x**n xp**m y**p yp**q)
+! **** mkj ****
+!  aug 26, 2003 added coefficients correcting theta
+!   for  raster ( x at target). Previously delta was corrected
+!   raster ( with up to 4 order) . y and phi were corrected
+!   for coefficient 00012 .
+!   Got coefficients from SIMC hms/recon_cosy.dat
+!     which is labeled:
+!        HMS Recon coefficients calculated 02-Feb-00  16:32:32 (HMS100 tune, P=-1.576
+!    The correction is up to 6th order.
+!  **** mkj *****
+ ---------------------------------------------------------------------
+   .342704157      .183629749E-01 -.742670280E-02  .256021141     1000
+  -3.02138749     -.152502758      .472422058E-01  .130192976     0100
+  -.673199668E-02 -.407567405      .258501913      .278565265E-02 0010
+   .110287597     -.483378951     -2.16704420      .407874298E-01 0001
+  -.495779826     -.548491564E-01  .147936777E-01 -.318166729     2000
+   4.49999138      .451213193     -.151686093      3.07552684     1100
+   .796294361      1.48286554     -.303847305     -.297352223E-01 1010
+   .261524802     -20.7938893      5.00402934     -.249121889     1001
+  -1.65128983     -.995337651      .536837816      2.14776190     0200
+  -6.47649939      .226936328     -.123499047     -.266124530     0110
+  -3.98529962     -6.32858639      1.39991277      4.27708696     0101
+   .113653600E-01  .690088486E-01 -.437114986E-01 -.438953228E-01 0020
+   .881941403      .227473684     -.628899349E-01  .475605891     0011
+  -3.71740449      3.93769576      .113368630     -3.53478948     0002
+   1.95723299     -.820253354E-01  .103519807      .936834319     3000
+  -31.4466915      2.57527314     -2.51407447     -13.5405899     2100
+  -5.89608158     -3.50686592      .471876566     -.312542194E-01 2010
+   40.5129114      58.0458860     -12.5399800     -2.74488175     2001
+   156.932838     -23.8014414      21.1502247      51.8797642     1200
+   57.3250219      40.7718723     -7.13532706      .284008412     1110
+  -298.395506     -639.340672      173.523125      10.6819836     1101
+   1.88702326      .266847247      .271167866      .954315611     1020
+  -24.4326614     -33.6584609     -.114525239      .434687181     1011
+   180.539265      53.7698885      3.76822635      31.9077575     1002
+  -279.667300      64.3177661     -60.6190378     -20.4989850     0300
+  -64.0800162     -86.2289723      11.5806188      3.95564331     0210
+  -247.924962      1524.67217     -483.832639     -46.7762650     0201
+  -21.7150970     -13.4646889      .401043380     -7.06883817     0120
+   211.296709      323.699221     -6.72775908      25.6924121     0111
+  -1440.32950      361.536621     -384.751537     -467.898102     0102
+  -.656643228      3.62511407     -1.27960041     -.358205667     0030
+  -1.58376006      67.9025208     -22.9733471     -3.15593520     0021
+   23.7224929     -987.049261      560.285176      34.1790370     0012
+  -180.294278      5983.27702     -3518.42552     -71.5085362     0003
+  -4.16665084      1.49723155     -.759630854     -2.30278775     4000
+   69.1944845     -46.6812620      21.8666725      58.0248309     3100
+   18.5878217      3.69400020      .602106657E-01 -1.46212838     3010
+  -227.551196     -96.3179597      .489353421      20.4063547     3001
+  -243.511203      535.564222     -239.722086     -491.562374     2200
+  -230.021793     -40.7163893      .940804015      47.2957535     2110
+   3150.68103      1424.96519      38.5622615     -521.870646     2101
+  -10.2926426     -3.56149925     -1.10754304     -1.00472512     2020
+   138.770921      263.397128     -15.7671082     -26.1090659     2011
+  -211.394983     -1895.91500      234.605250      406.369321     2002
+  -832.967425     -2704.20708      1207.70554      1718.29721     1300
+   429.578554     -234.967423     -20.3231936     -432.372956     1210
+  -11612.7737     -2668.84972     -1176.05843      4021.85608     1201
+   130.736123      70.7440580      5.42914965      16.7434624     1120
+  -1905.84788     -3551.03419      224.098878      343.660595     1111
+   363.670026      16804.9474     -2673.42755     -4586.96491     1102
+   6.91955519     -9.98074257      7.57861652      4.98155457     1030
+  -114.693315      63.5563817     -65.8809200     -36.5101930     1021
+   710.699660      4980.54374     -780.431756      60.2646297     1012
+  -7520.04385     -13187.0100      2486.26739      536.066295     1003
+   4347.54210      5161.66425     -2427.09769     -1950.10403     0400
+   1156.75126      1725.44332      188.803395      1323.39948     0310
+   12047.3853     -8815.37699      3407.72589     -11081.6594     0301
+  -26.9308548      60.3544886     -47.8411214     -67.2731587     0220
+   826.904198      4520.84740      916.341424     -974.052221     0211
+   19983.2818      7167.34390     -4225.20780      9358.46286     0202
+  -26.0105902     -188.661888      54.6866678     -37.4367181     0130
+   648.226508      2521.63446     -1052.88888      292.436840     0121
+  -8498.59600     -31139.3846      8203.03704     -697.205971     0112
+   63967.4235      128236.246     -24237.5465     -2336.48364     0103
+  -1.07911254      1.22231630     -.173675833      3.83526317     0040
+   9.94134973      258.488500     -10.9579494     -84.9391998     0031
+   158.508425     -2684.84975     -3.07844063      726.975013     0022
+  -3017.64724      14872.1349     -469.458295     -1667.46714     0013
+   466.801936     -55258.4859      4661.52259      3912.11649     0004
+   1.88180503     -2.54978699      1.02249535      4.96781903     5000
+   27.9678707      85.1165188     -25.6500831     -214.635508     4100
+  -19.6018738     -1.90234997     -2.26475906      3.07907961     4010
+   338.908065     -3.13670480      84.7364876     -48.5819865     4001
+  -1595.43911     -986.185572      116.663362      3438.51266     3200
+   256.718721     -2.53553430      86.9081338     -71.7570318     3110
+  -5619.82737      2554.66160     -3220.55423      1952.97191     3101
+   5.63505956      5.65523837     -.178636671     -9.82929708     3020
+   28.4611788     -535.882022      112.200518      63.7792928     3011
+  -3335.14896      6481.34708     -841.916671     -2407.64911     3002
+   17748.8967      4446.27557      1781.86526     -27888.4078     2300
+  -365.812295      713.968459     -1184.61359      175.230779     2210
+   26149.2365     -61279.3626      44610.2113     -20776.5089     2201
+   82.8634399     -82.7911036      40.7995535      223.223395     2120
+  -4711.16007      7972.91446     -3042.34279     -1765.82660     2111
+   111837.314     -122751.891      25531.4588      55520.9916     2102
+  -13.8546171     -44.6122707     -1.49148927      7.53241271     2030
+   379.629516      2219.19103     -283.396915     -8.69875712     2021
+  -2940.77125     -39662.7132      6848.07035      332.501630     2012
+  -4833.05018      291262.775     -40868.2995     -930.899000     2003
+  -75870.6114     -3035.26190     -19402.1263      118850.529     1400
+  -1090.33337     -5565.51679      8155.50571      3660.70792     1310
+  -76601.8398      456263.377     -281440.652      79440.7749     1301
+  -2721.78279     -512.627944     -452.280825     -1934.55934     1220
+   83864.4298      5816.23957      23306.1137      11961.6529     1211
+  -1306713.82      400125.281     -216284.553     -418149.697     1202
+   100.956187      1309.56623     -226.290462     -173.118732     1130
+  -4032.82257     -56824.1834      10616.2930      1607.73072     1121
+   79840.4758      782713.236     -150889.930     -14254.2451     1112
+  -330377.607     -4233718.14      646182.662      39788.7678     1103
+   6.06681234      7.12938618     -.220192683     -1.11240314     1040
+  -243.665570     -948.255482      153.244584      165.982247     1031
+   7181.11150      25876.5416     -2553.92502     -3998.95843     1022
+  -75005.8412     -104271.766     -23488.8683      29493.1698     1013
+   124433.465     -149264.761      219376.542     -13196.6766     1004
+   113262.801     -18667.6651      53113.3093     -209797.800     0500
+  -11632.6603      7242.92619     -21110.0334     -18679.7717     0410
+   301044.558     -1035121.76      655646.504     -71841.5676     0401
+   8371.11263      2221.54910      1543.28892      6295.96049     0320
+  -258028.488     -135832.589     -66162.4672     -36062.5709     0311
+   4218622.48      128274.356      663799.998      1118597.09     0302
+  -167.926507     -3321.19532      237.754153      1077.44822     0230
+   15453.6273      205593.650     -24454.8426     -11675.9844     0221
+  -438724.451     -3267312.69      463392.237      52657.7588     0212
+   2249544.29      17936516.3     -2298746.48     -118752.747     0203
+  -25.5889867     -61.6539042      41.8282483     -107.686882     0140
+   1866.01089     -1606.56535     -3207.84127      946.869731     0131
+  -71480.9824     -31234.8373      79600.1247      25102.1443     0122
+   814034.645      434180.428     -606101.312     -259877.198     0113
+  -2102030.39     -1584569.00      1583206.11      303094.278     0104
+  -3.59866781     -19.7084031     -18.4827350     -1.36754912     0050
+   50.9893561      384.737475      1103.07001      164.395576     0041
+   1986.13677      18001.7932     -34498.4820      90.2668576     0032
+  -30037.4064     -584192.598      674006.856      1100.41524     0023
+   20801.9238      6212983.14     -6436686.94     -17333.4694     0014
+   408508.663     -23722308.6      22707746.7      30963.1531     0005
+   1.12584242     0.000000000E+00 0.000000000E+00 0.768177700     00001
+  -3.30929606     0.000000000E+00 0.000000000E+00 -1.24824954     10001
+   9.55821560     0.000000000E+00 0.000000000E+00  7.18980365     01001
+   -12.8075508    0.000000000E+00 0.000000000E+00 -1.79480024     00002
+   31.3282870     0.000000000E+00 0.000000000E+00  4.30979026     20001
+   -366.358594    0.000000000E+00 0.000000000E+00 -47.5141035     11001
+   1281.27264     0.000000000E+00 0.000000000E+00  127.404997     02001
+   8.46492946     0.000000000E+00 0.000000000E+00  1.52132362     00201
+   -212.181332    0.000000000E+00 0.000000000E+00 -44.8420242     00111
+   86.2168045     0.000000000E+00 0.000000000E+00  13.9035869     10002
+   -202.375577    0.000000000E+00 0.000000000E+00 -75.1770025     01002
+   1292.40575     0.000000000E+00 0.000000000E+00  332.945301     00021
+   152.465743      412.398418      -118.551680    0.000000000E+00 00012
+   56.8142343     0.000000000E+00 0.000000000E+00  15.5951865     00003
+  -111.090641     0.000000000E+00 0.000000000E+00 -15.7177358     30001
+   1301.23090     0.000000000E+00 0.000000000E+00  249.924995     21001
+  -3377.64137     0.000000000E+00 0.000000000E+00 -1205.61003     12001
+  -1080.57253     0.000000000E+00 0.000000000E+00  1650.70638     03001
+  -206.758367     0.000000000E+00 0.000000000E+00 -8.63945412     10201
+   1286.12280     0.000000000E+00 0.000000000E+00  51.7345340     01201
+   457.388990     0.000000000E+00 0.000000000E+00  246.328770     10111
+   315.521007     0.000000000E+00 0.000000000E+00 -1468.88025     01111
+  -2.86376084     0.000000000E+00 0.000000000E+00 -77.6962058     20002
+  -284.532504     0.000000000E+00 0.000000000E+00  808.267218     11002
+  -54.9686786     0.000000000E+00 0.000000000E+00 -1887.27393     02002
+  -894.903229     0.000000000E+00 0.000000000E+00 -13.1721627     00202
+  -58.1868148     0.000000000E+00 0.000000000E+00 -1758.72381     10021
+   28.1001103     0.000000000E+00 0.000000000E+00  10443.6142     01021
+   24.9869024     0.000000000E+00 0.000000000E+00  377.176889     00112
+  -4.71048927     0.000000000E+00 0.000000000E+00 -174.332974     10003
+  -0.94935661     0.000000000E+00 0.000000000E+00  893.518077     01003
+   5.91096156     0.000000000E+00 0.000000000E+00 -2708.57078     00022
+   0.10590745     0.000000000E+00 0.000000000E+00 -149.156349     00004
+  0.744204334E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00101
+   2.84521581     0.000000000E+00 0.000000000E+00 0.000000000E+00 10101
+  -25.4393730     0.000000000E+00 0.000000000E+00 0.000000000E+00 01101
+ -0.932044457     0.000000000E+00 0.000000000E+00 0.000000000E+00 00011
+   13.7051555     0.000000000E+00 0.000000000E+00 0.000000000E+00 10011
+  -72.5425087     0.000000000E+00 0.000000000E+00 0.000000000E+00 01011
+   11.1885519     0.000000000E+00 0.000000000E+00 0.000000000E+00 00102
+  -64.8809742     0.000000000E+00 0.000000000E+00 0.000000000E+00 00030
+  -29.8576040     0.000000000E+00 0.000000000E+00 0.000000000E+00 20101
+   257.607594     0.000000000E+00 0.000000000E+00 0.000000000E+00 11101
+   222.800622     0.000000000E+00 0.000000000E+00 0.000000000E+00 02101
+  -5.74703998     0.000000000E+00 0.000000000E+00 0.000000000E+00 00301
+  -58.3774171     0.000000000E+00 0.000000000E+00 0.000000000E+00 20011
+  -75.7986407     0.000000000E+00 0.000000000E+00 0.000000000E+00 11011
+  -20.0306669     0.000000000E+00 0.000000000E+00 0.000000000E+00 02011
+   81.5699850     0.000000000E+00 0.000000000E+00 0.000000000E+00 00211
+   39.9072357     0.000000000E+00 0.000000000E+00 0.000000000E+00 10102
+   26.1337898     0.000000000E+00 0.000000000E+00 0.000000000E+00 01102
+   11.1028012     0.000000000E+00 0.000000000E+00 0.000000000E+00 00121
+   15.8064184     0.000000000E+00 0.000000000E+00 0.000000000E+00 10012
+   2.52332790     0.000000000E+00 0.000000000E+00 0.000000000E+00 01012
+   1.49611639     0.000000000E+00 0.000000000E+00 0.000000000E+00 00103
+   2.94427172     0.000000000E+00 0.000000000E+00 0.000000000E+00 00031
+   5.91096156     0.000000000E+00 0.000000000E+00 0.000000000E+00 00022
+  0.196444306     0.000000000E+00 0.000000000E+00 0.000000000E+00 00013
+  0.105907450     0.000000000E+00 0.000000000E+00 0.000000000E+00 00004
+   169.573680     0.000000000E+00 0.000000000E+00 0.000000000E+00 40001
+  -1188.40309     0.000000000E+00 0.000000000E+00 0.000000000E+00 31001
+  -2564.00859     0.000000000E+00 0.000000000E+00 0.000000000E+00 22001
+  -796.126266     0.000000000E+00 0.000000000E+00 0.000000000E+00 13001
+  -170.135856     0.000000000E+00 0.000000000E+00 0.000000000E+00 04001
+   25.2737324     0.000000000E+00 0.000000000E+00 0.000000000E+00 30101
+  -302.117798     0.000000000E+00 0.000000000E+00 0.000000000E+00 21101
+   10.5088748     0.000000000E+00 0.000000000E+00 0.000000000E+00 12101
+   17.9545609     0.000000000E+00 0.000000000E+00 0.000000000E+00 03101
+   393.042869     0.000000000E+00 0.000000000E+00 0.000000000E+00 20201
+  -2654.79358     0.000000000E+00 0.000000000E+00 0.000000000E+00 11201
+  -675.207959     0.000000000E+00 0.000000000E+00 0.000000000E+00 02201
+   31.5374558     0.000000000E+00 0.000000000E+00 0.000000000E+00 10301
+  -66.2872749     0.000000000E+00 0.000000000E+00 0.000000000E+00 01301
+   11.0483668     0.000000000E+00 0.000000000E+00 0.000000000E+00 00401
+   53.9952733     0.000000000E+00 0.000000000E+00 0.000000000E+00 30011
+   19.8514223     0.000000000E+00 0.000000000E+00 0.000000000E+00 21011
+   6.99512779     0.000000000E+00 0.000000000E+00 0.000000000E+00 12011
+   2.03058969     0.000000000E+00 0.000000000E+00 0.000000000E+00 03011
+  -414.755722     0.000000000E+00 0.000000000E+00 0.000000000E+00 20111
+  -242.924295     0.000000000E+00 0.000000000E+00 0.000000000E+00 11111
+  -53.2831769     0.000000000E+00 0.000000000E+00 0.000000000E+00 02111
+   14.0854426     0.000000000E+00 0.000000000E+00 0.000000000E+00 10211
+  -7.18970025     0.000000000E+00 0.000000000E+00 0.000000000E+00 01211
+   46.6175282     0.000000000E+00 0.000000000E+00 0.000000000E+00 00311
+   327.023510     0.000000000E+00 0.000000000E+00 0.000000000E+00 30002
+  -46.7603764     0.000000000E+00 0.000000000E+00 0.000000000E+00 21002
+  -13.1130776     0.000000000E+00 0.000000000E+00 0.000000000E+00 12002
+  -2.14608775     0.000000000E+00 0.000000000E+00 0.000000000E+00 03002
+  -405.162634     0.000000000E+00 0.000000000E+00 0.000000000E+00 20102
+  -37.8327394     0.000000000E+00 0.000000000E+00 0.000000000E+00 11102
+  -3.28616433     0.000000000E+00 0.000000000E+00 0.000000000E+00 02102
+  -797.756623     0.000000000E+00 0.000000000E+00 0.000000000E+00 10202
+  -100.253463     0.000000000E+00 0.000000000E+00 0.000000000E+00 01202
+  -58.1494000     0.000000000E+00 0.000000000E+00 0.000000000E+00 00302
+  -49.7218825     0.000000000E+00 0.000000000E+00 0.000000000E+00 20021
+  -23.4940000     0.000000000E+00 0.000000000E+00 0.000000000E+00 11021
+  -4.63065295     0.000000000E+00 0.000000000E+00 0.000000000E+00 02021
+  -7.01252658     0.000000000E+00 0.000000000E+00 0.000000000E+00 10121
+  -1.66128414     0.000000000E+00 0.000000000E+00 0.000000000E+00 01121
+   8.99761159     0.000000000E+00 0.000000000E+00 0.000000000E+00 00221
+  -26.9158564     0.000000000E+00 0.000000000E+00 0.000000000E+00 20012
+  -2.93740577     0.000000000E+00 0.000000000E+00 0.000000000E+00 11012
+ -0.281699956     0.000000000E+00 0.000000000E+00 0.000000000E+00 02012
+  -44.2613806     0.000000000E+00 0.000000000E+00 0.000000000E+00 10112
+  -5.18873926     0.000000000E+00 0.000000000E+00 0.000000000E+00 01112
+  -2.30913015     0.000000000E+00 0.000000000E+00 0.000000000E+00 00212
+   12.9742811     0.000000000E+00 0.000000000E+00 0.000000000E+00 20003
+   1.29863794     0.000000000E+00 0.000000000E+00 0.000000000E+00 11003
+  0.179319821     0.000000000E+00 0.000000000E+00 0.000000000E+00 02003
+  0.359636043     0.000000000E+00 0.000000000E+00 0.000000000E+00 10103
+ -0.686000125E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 01103
+  0.479897726     0.000000000E+00 0.000000000E+00 0.000000000E+00 00203
+   2.12875475     0.000000000E+00 0.000000000E+00 0.000000000E+00 10040
+   1.00542982     0.000000000E+00 0.000000000E+00 0.000000000E+00 01040
+  0.574499599E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00140
+  -1.51837641     0.000000000E+00 0.000000000E+00 0.000000000E+00 10031
+ -0.259505376     0.000000000E+00 0.000000000E+00 0.000000000E+00 01031
+   1.21378871     0.000000000E+00 0.000000000E+00 0.000000000E+00 00131
+  -3.83394487     0.000000000E+00 0.000000000E+00 0.000000000E+00 10022
+ -0.448803101     0.000000000E+00 0.000000000E+00 0.000000000E+00 01022
+  0.432742246E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00122
+  0.433795780E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 10013
+  0.381017134E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 01013
+  0.717327023E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00113
+  0.404935152E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 10004
+ -0.704014248E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 01004
+  0.417717692E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00104
+ -0.108814119     0.000000000E+00 0.000000000E+00 0.000000000E+00 00050
+  0.265212468     0.000000000E+00 0.000000000E+00 0.000000000E+00 00041
+  0.342938618E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00032
+  0.274219030E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00023
+  0.382346774E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00014
+  0.146595810E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00005
+  -221.978975     0.000000000E+00 0.000000000E+00 0.000000000E+00 50001
+   2081.48433     0.000000000E+00 0.000000000E+00 0.000000000E+00 41001
+  -430.560936     0.000000000E+00 0.000000000E+00 0.000000000E+00 32001
+  -218.050288     0.000000000E+00 0.000000000E+00 0.000000000E+00 23001
+  -55.5844501     0.000000000E+00 0.000000000E+00 0.000000000E+00 14001
+  -11.3642568     0.000000000E+00 0.000000000E+00 0.000000000E+00 05001
+   108.432463     0.000000000E+00 0.000000000E+00 0.000000000E+00 40101
+  -931.864934     0.000000000E+00 0.000000000E+00 0.000000000E+00 31101
+  -192.977158     0.000000000E+00 0.000000000E+00 0.000000000E+00 22101
+  -28.6473383     0.000000000E+00 0.000000000E+00 0.000000000E+00 13101
+  -3.41017534     0.000000000E+00 0.000000000E+00 0.000000000E+00 04101
+   61.3278498     0.000000000E+00 0.000000000E+00 0.000000000E+00 30201
+  -43.6289186     0.000000000E+00 0.000000000E+00 0.000000000E+00 21201
+  -47.1521819     0.000000000E+00 0.000000000E+00 0.000000000E+00 12201
+  -14.9198830     0.000000000E+00 0.000000000E+00 0.000000000E+00 03201
+   58.8310901     0.000000000E+00 0.000000000E+00 0.000000000E+00 20301
+  -36.3378650     0.000000000E+00 0.000000000E+00 0.000000000E+00 11301
+  -8.01512328     0.000000000E+00 0.000000000E+00 0.000000000E+00 02301
+  -104.944287     0.000000000E+00 0.000000000E+00 0.000000000E+00 10401
+   34.7713353     0.000000000E+00 0.000000000E+00 0.000000000E+00 01401
+  -122.288125     0.000000000E+00 0.000000000E+00 0.000000000E+00 00501
+  -77.5453423     0.000000000E+00 0.000000000E+00 0.000000000E+00 40011
+  -64.1314449     0.000000000E+00 0.000000000E+00 0.000000000E+00 31011
+  -13.9857505     0.000000000E+00 0.000000000E+00 0.000000000E+00 22011
+  -2.33654082     0.000000000E+00 0.000000000E+00 0.000000000E+00 13011
+ -0.340596491     0.000000000E+00 0.000000000E+00 0.000000000E+00 04011
+   29.3394194     0.000000000E+00 0.000000000E+00 0.000000000E+00 30111
+   48.1149331     0.000000000E+00 0.000000000E+00 0.000000000E+00 21111
+   9.50174924     0.000000000E+00 0.000000000E+00 0.000000000E+00 12111
+   1.38837541     0.000000000E+00 0.000000000E+00 0.000000000E+00 03111
+   52.4792515     0.000000000E+00 0.000000000E+00 0.000000000E+00 20211
+   3.80058742     0.000000000E+00 0.000000000E+00 0.000000000E+00 11211
+  0.234555136     0.000000000E+00 0.000000000E+00 0.000000000E+00 02211
+   39.5808085     0.000000000E+00 0.000000000E+00 0.000000000E+00 10311
+   15.2563542     0.000000000E+00 0.000000000E+00 0.000000000E+00 01311
+  -7.28080327     0.000000000E+00 0.000000000E+00 0.000000000E+00 00411
+  -597.497346     0.000000000E+00 0.000000000E+00 0.000000000E+00 40002
+  -102.721749     0.000000000E+00 0.000000000E+00 0.000000000E+00 31002
+  -15.3681107     0.000000000E+00 0.000000000E+00 0.000000000E+00 22002
+  -2.17208316     0.000000000E+00 0.000000000E+00 0.000000000E+00 13002
+ -0.274609264     0.000000000E+00 0.000000000E+00 0.000000000E+00 04002
+   106.837983     0.000000000E+00 0.000000000E+00 0.000000000E+00 30102
+   12.3766158     0.000000000E+00 0.000000000E+00 0.000000000E+00 21102
+   1.46530330     0.000000000E+00 0.000000000E+00 0.000000000E+00 12102
+  0.185393293     0.000000000E+00 0.000000000E+00 0.000000000E+00 03102
+   11.8887131     0.000000000E+00 0.000000000E+00 0.000000000E+00 20202
+  -2.55669988     0.000000000E+00 0.000000000E+00 0.000000000E+00 11202
+ -0.863809017     0.000000000E+00 0.000000000E+00 0.000000000E+00 02202
+   18.8047014     0.000000000E+00 0.000000000E+00 0.000000000E+00 10302
+   2.49667580     0.000000000E+00 0.000000000E+00 0.000000000E+00 01302
+  -3.93088497     0.000000000E+00 0.000000000E+00 0.000000000E+00 00402
+   12.7034651     0.000000000E+00 0.000000000E+00 0.000000000E+00 30021
+   7.96187516     0.000000000E+00 0.000000000E+00 0.000000000E+00 21021
+   1.60982507     0.000000000E+00 0.000000000E+00 0.000000000E+00 12021
+  0.257565675     0.000000000E+00 0.000000000E+00 0.000000000E+00 03021
+   5.93574507     0.000000000E+00 0.000000000E+00 0.000000000E+00 20121
+  0.679348453     0.000000000E+00 0.000000000E+00 0.000000000E+00 11121
+  0.757550151E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 02121
+   7.36294581     0.000000000E+00 0.000000000E+00 0.000000000E+00 10221
+   2.30080135     0.000000000E+00 0.000000000E+00 0.000000000E+00 01221
+ -0.118383007     0.000000000E+00 0.000000000E+00 0.000000000E+00 00321
+   14.4258145     0.000000000E+00 0.000000000E+00 0.000000000E+00 30012
+   1.52658163     0.000000000E+00 0.000000000E+00 0.000000000E+00 21012
+  0.163028817     0.000000000E+00 0.000000000E+00 0.000000000E+00 12012
+  0.171689509E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 03012
+   15.0743734     0.000000000E+00 0.000000000E+00 0.000000000E+00 20112
+   1.60706583     0.000000000E+00 0.000000000E+00 0.000000000E+00 11112
+  0.166465023     0.000000000E+00 0.000000000E+00 0.000000000E+00 02112
+   2.28803297     0.000000000E+00 0.000000000E+00 0.000000000E+00 10212
+  0.265413449     0.000000000E+00 0.000000000E+00 0.000000000E+00 01212
+   3.26683930     0.000000000E+00 0.000000000E+00 0.000000000E+00 00312
+  -1.44694805     0.000000000E+00 0.000000000E+00 0.000000000E+00 30003
+ -0.301551821     0.000000000E+00 0.000000000E+00 0.000000000E+00 21003
+ -0.979405958E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 12003
+ -0.218874633E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 03003
+  0.213731277     0.000000000E+00 0.000000000E+00 0.000000000E+00 20103
+ -0.424871757E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 11103
+ -0.263981236E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 02103
+ -0.390336182     0.000000000E+00 0.000000000E+00 0.000000000E+00 10203
+ -0.220563449E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 01203
+ -0.252061070E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00303
+  0.835395726     0.000000000E+00 0.000000000E+00 0.000000000E+00 20031
+  0.101082424     0.000000000E+00 0.000000000E+00 0.000000000E+00 11031
+  0.113826515E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 02031
+  0.967539242     0.000000000E+00 0.000000000E+00 0.000000000E+00 10131
+  0.283940573     0.000000000E+00 0.000000000E+00 0.000000000E+00 01131
+  0.704094684E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00231
+   1.96621152     0.000000000E+00 0.000000000E+00 0.000000000E+00 20022
+  0.211929662     0.000000000E+00 0.000000000E+00 0.000000000E+00 11022
+  0.222297020E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 02022
+  0.192176728     0.000000000E+00 0.000000000E+00 0.000000000E+00 10122
+  0.208429359E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 01122
+  0.409197722     0.000000000E+00 0.000000000E+00 0.000000000E+00 00222
+  0.116208798E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 20013
+ -0.861751979E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 11013
+ -0.368034008E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 02013
+ -0.305595883E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 10113
+  0.792820773E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 01113
+ -0.103581312E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00213
+ -0.317665296E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 20004
+ -0.868692894E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 11004
+ -0.130336633E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 02004
+ -0.485930059E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 10104
+ -0.181322677E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 01104
+ -0.225313015E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00204
+  0.102736083     0.000000000E+00 0.000000000E+00 0.000000000E+00 10041
+  0.316791817E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 01041
+  0.177406989E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00141
+  0.128130316E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 10032
+  0.134309873E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 01032
+  0.442654916E-01 0.000000000E+00 0.000000000E+00 0.000000000E+00 00132
+ -0.500745668E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 10023
+ -0.362264901E-04 0.000000000E+00 0.000000000E+00 0.000000000E+00 01023
+ -0.100846774E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 00123
+ -0.259911565E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 10014
+ -0.865988929E-05 0.000000000E+00 0.000000000E+00 0.000000000E+00 01014
+  0.126433262E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 00114
+ -0.376405317E-03 0.000000000E+00 0.000000000E+00 0.000000000E+00 10005
+ -0.495262130E-04 0.000000000E+00 0.000000000E+00 0.000000000E+00 01005
+  0.505435025E-04 0.000000000E+00 0.000000000E+00 0.000000000E+00 00105
+ -0.126889588E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00060
+  0.318778677E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00051
+  0.474276817E-02 0.000000000E+00 0.000000000E+00 0.000000000E+00 00042
+  0.108549858E-04 0.000000000E+00 0.000000000E+00 0.000000000E+00 00033
+  0.908354900E-04 0.000000000E+00 0.000000000E+00 0.000000000E+00 00024
+  0.311879063E-05 0.000000000E+00 0.000000000E+00 0.000000000E+00 00015
+  0.232698271E-05 0.000000000E+00 0.000000000E+00 0.000000000E+00 00006
+ ---------------------------------------------------------------------
diff --git a/examples/PARAM/sos_recon_coeff.dat b/examples/PARAM/sos_recon_coeff.dat
new file mode 100644
index 0000000..dde306f
--- /dev/null
+++ b/examples/PARAM/sos_recon_coeff.dat
@@ -0,0 +1,267 @@
+ ---------------------------------------------------------------------
+   .379139198E-01  .842839975E-03 -.165349417E-01  1.09008330     1000
+  -.340475497     -.190557451E-03 -.303884217E-01  .864312703E-01 0100
+   .412883275E-04 -.183808587      .979733182      .254302566E-02 0010
+   .815995425E-03  .575826053E-01 -5.49469348     -.117673035     0001
+  -.528286695E-01 -.422038183E-01  .689674787E-01  .795969827     2000
+  -.280522018E-01 -.615710778E-01 -.407260139E-01  3.54822776     1100
+   .580273258E-01  2.36544241      1.35260940      .230035665     1010
+  -.248614039     -14.5575974      7.01872929     -1.25071096     1001
+  -.252140610     -.101163199E-01 -.156760594E-01 -1.48667569     0200
+  -.811528766E-01  .901049007     -4.92046279      .288221600     0110
+  -.820406238E-01 -5.72170059      4.90060419     -.587392000     0101
+   .591913926E-01  .442996896E-01 -.230355364      .173773836     0020
+  -.219573091      .138562767      1.13340742     -5.56476527     0011
+   .183157973     -.826984028     -1.15673466      16.2159664     0002
+   .141613116      .563557025E-01 -.132614085      .893431379     3000
+  -.736574959     -.239878097      .448683179      6.22480992     2100
+  -.147618189      5.49472968     -3.14205915     -1.26710492     2010
+   .697380983      1.62324987     -12.8412738      6.95062674     2001
+   .799836904     -.338942744E-01 -.339159381      4.84876509     1200
+   1.04463497      5.12456575      11.5182373      .161979673     1110
+  -5.73077452     -65.4895006      21.7592854     -5.62857441     1101
+  -.940144651E-01  .729997893     -1.23899720      3.93503678     1020
+   .445867546     -5.51357171      1.00751517     -28.1235302     1011
+   5.30156098      6.84633968      19.6403132      16.8571460     1002
+  -3.12917608      .584819316      1.22766131     -2.41644263     0300
+   .446613746     -21.0828990      9.59314583     -5.28850025     0210
+  -3.88969612      80.7762158     -64.0493585      22.3948037     0201
+  -1.65486895      .429698026      7.45203435     -10.9048792     0120
+   15.9989982     -8.21241443     -75.9609989      81.9559180     0111
+  -60.4402763      21.2298911      212.509960     -248.259132     0102
+   .604652321      1.53919563      5.35871564     -6.62271245     0030
+  -9.25519025     -31.3433274     -16.5768903      38.0695796     0021
+   38.5001423      192.329974      34.8941730     -88.1044751     0012
+  -43.0872876     -450.298401     -695.196371      496.318701     0003
+  -.280910732      .110007520     -1.29442669     -6.87528629     4000
+   1.70456640     -.790200026E-01  3.24383795     -7.35688444     3100
+  -.980466960     -11.2835109      15.9675705     -14.5026852     3010
+   2.57242009      38.4628243     -12.1845687      66.4941464     3001
+  -3.08843217      .385636778     -2.10129804      32.8640525     2200
+   5.90646560      48.5058838     -79.1614198     -9.58971005     2110
+  -10.4317985     -17.6586954      56.1703738      103.838657     2101
+  -.495166447      .291029779      8.65190371     -61.9746102     2020
+   11.1187517     -13.7564113     -44.1495212      750.712571     2011
+  -29.1394954      29.2851208      7.16099318     -1767.78840     2002
+  -9.07930174      .465978293      .933872266     -16.0084695     1300
+  -1.83572469     -52.4119547      134.467671     -26.5805950     1210
+  -10.6276201     -115.283001     -131.157302      95.6252426     1201
+  -11.6763151      3.96624253      25.5702701     -120.170295     1120
+   81.2559728      24.2679709     -225.807372      1137.51843     1111
+  -112.359178     -116.026197      479.835157     -3361.40683     1102
+  -5.29485259     -39.0977840      114.868020     -29.3504217     1030
+   52.9540035      435.131830     -1915.00058      456.800801     1021
+  -255.245852     -1998.72127      10430.0610     -3507.86935     1012
+   710.282918      3346.50103     -18642.5174      9294.88698     1003
+   5.47692943     -3.42404936      .398772011      26.3688033     0400
+  -7.02852980      36.0408590     -125.997347      12.0775052     0310
+   26.0411563     -59.5277011      195.059592     -.870074117     0301
+  -1.70642864     -8.05137099      66.2405899     -108.403205     0220
+   116.037347      25.1955957     -419.542715      685.197915     0211
+  -498.764979     -12.3929700      868.952277     -1583.11852     0202
+  -11.4654209      3.91337164      105.681673      .948774663     0130
+   182.564797      367.199616     -3314.03838      275.734707     0121
+  -785.464215     -1180.95553      21302.3634     -3567.09815     0112
+   1074.60409     -936.552421     -39399.4921      14185.0001     0103
+   27.3214997     -32.7467761      27.7875506     -31.0089464     0040
+  -72.0083764      499.365480      171.394594     -263.310809     0031
+  -1228.82936     -3038.62107     -5146.59971      3660.01403     0022
+   6990.44924      7629.29170      26216.2932     -19292.3817     0013
+  -12389.5540     -5952.19133     -40019.6818      27353.7886     0004
+   .233451694     -.725229200      4.49313708     -43.3437356     5000
+  -2.63479797      2.33506679     -12.0336095     -40.6088590     4100
+   2.80978052      56.6704274     -103.441626     -72.3865672     4010
+  -7.68051561     -132.315167      517.040649      601.038133     4001
+   1.55805319     -3.82640802      15.5718115      47.3538727     3200
+  -17.8565522     -346.988718      507.124617     -323.449119     3110
+   59.8473991      827.189754     -1801.08879      1569.31536     3101
+  -4.92016617     -7.94687178      12.2896473     -135.280605     3020
+   56.5503782      109.903033     -172.184206      2326.00503     3011
+  -184.441942     -120.544100      410.841450     -5124.02325     3002
+   34.2299151     -5.62614655     -19.6635668      28.5705216     2300
+   50.2036531      730.763870     -1384.99022      323.129619     2210
+  -92.7242514     -1185.38953      3546.92453     -1234.40563     2201
+  -21.1694556     -4.26971067     -187.353025     -263.232713     2120
+   50.8019739     -9.39109616      1961.55237      2703.80378     2111
+   407.045396      336.291695     -4962.65404     -7780.19411     2102
+  -32.8883563     -267.260855     -481.397420     -103.624057     2030
+   553.125925      3101.61079      6834.48349      6266.75033     2021
+  -2233.40649     -10995.3016     -40114.5660     -41964.2086     2012
+   1986.62069      12638.4725      95744.3925      70261.8242     2003
+  -157.435585      19.8341717      13.9972473      315.679803     1400
+  -52.8533936     -330.790519      1380.72305     -25.2066060     1310
+   43.1893277     -589.102339     -2937.18781      26.6981676     1301
+  -.131234530     -112.413423     -23.3238739      1664.16589     1220
+   284.871958      401.268662      200.871609     -12034.7631     1211
+  -1481.98521     -206.864693     -1872.19232      21158.0840     1202
+  -26.4458611     -371.408877     -1592.68708      147.661187     1130
+  -376.587514      3004.18852      14892.4013      5188.12933     1121
+   2666.07022     -20519.1906     -55737.5598     -46520.0024     1112
+  -2562.26184      61442.9024      72144.0869      108561.498     1103
+   15.2846925     -178.339099      98.7028679      778.606191     1040
+   259.514383      2089.90320      1620.09957     -20530.4149     1031
+  -972.532761     -10654.4336     -22167.5340      179631.580     1022
+  -6644.40422      31062.5766      95321.3704     -687066.885     1013
+   18138.9284     -39607.3611     -127382.262      803139.181     1004
+   260.033050     -9.37175174     -1.22854324     -149.131593     0500
+   12.9346160     -209.050875     -1828.13619      367.568826     0410
+   1.03860566      1642.89764      6259.32653     -1101.69817     0401
+   228.168555      47.9387217     -1007.66858     -1366.33588     0320
+  -2473.77471      401.516926      8723.22926      10463.1494     0311
+   7453.65344     -2459.90132     -20213.8613     -11893.8889     0302
+   43.0657501      2527.25046     -3579.69004     -265.543867     0230
+   85.9932457     -26732.0327      29217.0535      9454.46277     0221
+  -1310.73720      95335.7974     -74758.8642     -71048.3225     0212
+   1714.54526     -88924.4487      80562.0556      168080.435     0203
+   392.721684      363.565031     -671.427466     -301.828056     0140
+  -4902.70897     -5019.37580      3974.69466      12449.5401     0131
+   15252.5174      23287.3413      5749.87907     -48658.7614     0122
+   24169.2549     -36349.3754     -39136.7356     -186715.777     0113
+  -104121.773      895.479217     -3184.62532      792792.749     0104
+   41.5090986     -1148.66253     -1959.97263      651.990633     0050
+  -902.499745      24376.2669      13971.6905     -6281.70119     0041
+   8862.26009     -201261.260     -102332.755      6816.59750     0032
+  -33284.0983      822732.760      1107381.57      144256.913     0023
+   37274.2374     -1756990.04     -5432659.57     -988369.924     0014
+   3194.06668      1808484.40      9426156.81      2089893.33     0005
+  -.734005948      .000000000E+00  .000000000E+00  101.080958     6000
+   4.03546917      .000000000E+00  .000000000E+00 -235.474771     5100
+   .000000000E+00 -48.0176565      274.841961      .000000000E+00 5010
+   .000000000E+00  97.8419414     -1415.79591      .000000000E+00 5001
+   2.05839214      .000000000E+00  .000000000E+00 -174.906709     4200
+   .000000000E+00  811.417591     -1360.99969      .000000000E+00 4110
+   .000000000E+00 -2087.95015      5179.72150      .000000000E+00 4101
+  -24.6478776      .000000000E+00  .000000000E+00  1387.09849     4020
+   209.361971      .000000000E+00  .000000000E+00 -13347.2648     4011
+  -506.249544      .000000000E+00  .000000000E+00  25591.7786     4002
+  -35.9951898      .000000000E+00  .000000000E+00  1370.20393     3300
+   .000000000E+00 -3600.45042      3538.05250      .000000000E+00 3210
+   .000000000E+00  7524.67426     -10368.9029      .000000000E+00 3201
+   16.7701194      .000000000E+00  .000000000E+00  1785.11257     3120
+   448.015228      .000000000E+00  .000000000E+00 -12501.4175     3111
+  -2582.00272      .000000000E+00  .000000000E+00  70665.3189     3102
+   .000000000E+00  26.3558036     -3598.27809      .000000000E+00 3030
+   .000000000E+00 -1707.32364      32248.0060      .000000000E+00 3021
+   .000000000E+00  2366.24326     -95444.5145      .000000000E+00 3012
+   .000000000E+00  3260.93076      58043.6500      .000000000E+00 3003
+   85.1713624      .000000000E+00  .000000000E+00 -2956.22410     2400
+   .000000000E+00  7780.01535     -2907.83488      .000000000E+00 2310
+   .000000000E+00 -14636.9488      8134.53415      .000000000E+00 2301
+  -301.445311      .000000000E+00  .000000000E+00  4634.51958     2220
+   509.781116      .000000000E+00  .000000000E+00 -23255.2766     2211
+   3016.84963      .000000000E+00  .000000000E+00 -2992.69769     2202
+   .000000000E+00 -2765.24904      2365.34945      .000000000E+00 2130
+   .000000000E+00  29060.8651      30172.3857      .000000000E+00 2121
+   .000000000E+00 -105263.627     -343723.291      .000000000E+00 2112
+   .000000000E+00  164542.543      859596.652      .000000000E+00 2103
+   617.734883      .000000000E+00  .000000000E+00  13164.1544     2040
+  -9762.94511      .000000000E+00  .000000000E+00 -329312.117     2031
+   61777.0098      .000000000E+00  .000000000E+00  2755829.40     2022
+  -188522.811      .000000000E+00  .000000000E+00 -9589667.98     2013
+   293596.776      .000000000E+00  .000000000E+00  11777955.3     2004
+   155.742569      .000000000E+00  .000000000E+00  8004.70162     1500
+   .000000000E+00 -10710.7950     -3347.29371      .000000000E+00 1410
+   .000000000E+00  26613.4941      7824.54712      .000000000E+00 1401
+   1595.56179      .000000000E+00  .000000000E+00 -13614.5237     1320
+  -2346.18950      .000000000E+00  .000000000E+00  69708.0101     1311
+  -14502.8338      .000000000E+00  .000000000E+00 -54026.6988     1302
+   .000000000E+00  10576.5378      15937.5378      .000000000E+00 1230
+   .000000000E+00 -62086.0056     -76687.9044      .000000000E+00 1221
+   .000000000E+00 -55297.1950     -298675.928      .000000000E+00 1212
+   .000000000E+00  588445.113      1388085.99      .000000000E+00 1203
+  -597.751474      .000000000E+00  .000000000E+00  2913.83304     1140
+   6306.24490      .000000000E+00  .000000000E+00 -120786.591     1131
+   9249.80576      .000000000E+00  .000000000E+00  1516439.68     1122
+  -392504.112      .000000000E+00  .000000000E+00 -7016627.96     1113
+   1095950.97      .000000000E+00  .000000000E+00  9306653.35     1104
+   .000000000E+00 -8733.05877     -8552.74975      .000000000E+00 1050
+   .000000000E+00  137325.477      352699.311      .000000000E+00 1041
+   .000000000E+00 -910386.624     -3196889.05      .000000000E+00 1032
+   .000000000E+00  3234363.59      10796249.9      .000000000E+00 1023
+   .000000000E+00 -7252714.69     -10916643.6      .000000000E+00 1014
+   .000000000E+00  9543523.46     -6478540.83      .000000000E+00 1005
+   187.600021      .000000000E+00  .000000000E+00 -7467.84730     0600
+   .000000000E+00  8248.31060      9123.23763      .000000000E+00 0510
+   .000000000E+00 -26718.3257     -20637.9818      .000000000E+00 0501
+  -972.185544      .000000000E+00  .000000000E+00  28252.5400     0420
+  -5564.66121      .000000000E+00  .000000000E+00 -252365.927     0411
+   37465.5972      .000000000E+00  .000000000E+00  611723.395     0402
+   .000000000E+00 -22890.0920      22487.2028      .000000000E+00 0330
+   .000000000E+00  229636.659     -133458.844      .000000000E+00 0321
+   .000000000E+00 -795394.324      102279.096      .000000000E+00 0312
+   .000000000E+00  1010473.66      227935.761      .000000000E+00 0303
+  -824.302042      .000000000E+00  .000000000E+00  6197.39315     0240
+  -2787.07768      .000000000E+00  .000000000E+00 -40717.8144     0231
+   192262.614      .000000000E+00  .000000000E+00 -104938.059     0222
+  -1015988.95      .000000000E+00  .000000000E+00  1807893.56     0213
+   1581572.44      .000000000E+00  .000000000E+00 -3864158.57     0204
+   .000000000E+00  4929.32145      8058.93465      .000000000E+00 0150
+   .000000000E+00 -157973.959      292812.788      .000000000E+00 0141
+   .000000000E+00  1124967.46     -4749970.24      .000000000E+00 0132
+   .000000000E+00 -1039698.52      25070359.8      .000000000E+00 0123
+   .000000000E+00 -14297693.9     -62207004.0      .000000000E+00 0114
+   .000000000E+00  38481213.4      72831847.8      .000000000E+00 0105
+  -5527.79362      .000000000E+00  .000000000E+00 -5997.21048     0060
+   80926.6935      .000000000E+00  .000000000E+00  130411.524     0051
+  -538292.630      .000000000E+00  .000000000E+00 -303996.532     0042
+   1573640.14      .000000000E+00  .000000000E+00 -2290518.44     0033
+   257457.657      .000000000E+00  .000000000E+00  6087225.09     0024
+  -9176115.43      .000000000E+00  .000000000E+00  5616794.55     0015
+   12403312.1      .000000000E+00  .000000000E+00  2571808.44     0006
+ -0.432787425     0.000000000E+00 0.000000000E+00 0.391005873     00001
+  0.111255322     0.000000000E+00 0.000000000E+00 0.801847205     10001
+ -0.222099911     0.000000000E+00 0.000000000E+00  1.39223024     01001
+  0.000000000E+00  1.05404300    -0.472399367     0.000000000E+00 00101
+  0.000000000E+00 -6.27039977      3.78630828     0.000000000E+00 00011
+  0.187705863E-01 0.000000000E+00 0.000000000E+00 0.288248330     00002
+  0.466755273E-01 0.000000000E+00 0.000000000E+00 0.544095505     20001
+  0.255381022E-01 0.000000000E+00 0.000000000E+00  7.40084080     11001
+ -0.695597072     0.000000000E+00 0.000000000E+00  2.94559377     02001
+  0.000000000E+00  2.46445435     -2.54411388     0.000000000E+00 10101
+  0.000000000E+00-0.875373764      2.46186900     0.000000000E+00 01101
+  0.206789175     0.000000000E+00 0.000000000E+00-0.296304303     00201
+  0.000000000E+00 0.328318087      2.21842346     0.000000000E+00 10011
+  0.000000000E+00 -17.7113496      10.4387423     0.000000000E+00 01011
+  -2.83396800     0.000000000E+00 0.000000000E+00-0.810485111     00111
+  0.333321181E-01 0.000000000E+00 0.000000000E+00 0.727649980     10002
+ -0.106615205     0.000000000E+00 0.000000000E+00  2.24994616     01002
+  0.000000000E+00-0.440727542E-01 0.829650831E-03 0.000000000E+00 00102
+   9.05531754     0.000000000E+00 0.000000000E+00  4.35294783     00021
+  0.000000000E+00  1.17053477    -0.377328467     0.000000000E+00 00012
+  0.462025566E-02 0.000000000E+00 0.000000000E+00 0.224334624     00003
+ -0.213218941E-01 0.000000000E+00 0.000000000E+00-0.231138001     30001
+  0.671951131     0.000000000E+00 0.000000000E+00  7.83883863     21001
+  -2.48151318     0.000000000E+00 0.000000000E+00  51.8024010     12001
+  -1.53316700     0.000000000E+00 0.000000000E+00 -12.1339968     03001
+  0.000000000E+00 -1.93828614     -1.12464412     0.000000000E+00 20101
+  0.000000000E+00  14.7092776     -10.0892893     0.000000000E+00 11101
+  0.000000000E+00 -7.70961193      11.6264715     0.000000000E+00 02101
+   1.21564770     0.000000000E+00 0.000000000E+00-0.997171793     10201
+  -2.13666486     0.000000000E+00 0.000000000E+00 -8.63467712     01201
+  0.000000000E+00 -5.64716834      6.40534160     0.000000000E+00 00301
+  0.000000000E+00  10.3926000     -4.84593297     0.000000000E+00 20011
+  0.000000000E+00 -16.1788532      21.9577619     0.000000000E+00 11011
+  0.000000000E+00 -39.5563571      24.3791266     0.000000000E+00 02011
+  -8.29395208     0.000000000E+00 0.000000000E+00 -12.3262090     10111
+   14.3644302     0.000000000E+00 0.000000000E+00  66.1129198     01111
+  0.000000000E+00 -11.0477274     -7.23587036     0.000000000E+00 00211
+  0.254197599E-01 0.000000000E+00 0.000000000E+00 0.395279317     20002
+  0.219900371E-01 0.000000000E+00 0.000000000E+00  11.1361643     11002
+ -0.763757716     0.000000000E+00 0.000000000E+00  11.5920111     02002
+  0.000000000E+00 -1.23069700     0.143156389     0.000000000E+00 10102
+  0.000000000E+00  5.40559063     -2.75062054     0.000000000E+00 01102
+  0.489849474E-01 0.000000000E+00 0.000000000E+00-0.413153323     00202
+   3.71345858     0.000000000E+00 0.000000000E+00  35.6736779     10021
+  -3.93972567     0.000000000E+00 0.000000000E+00 -139.812832     01021
+  0.000000000E+00  271.617600     -130.343209     0.000000000E+00 00121
+  0.000000000E+00 0.685511203     0.670438936     0.000000000E+00 10012
+  0.000000000E+00 -12.9229381      10.2806595     0.000000000E+00 01012
+ -0.469630205E-01 0.000000000E+00 0.000000000E+00 0.929541757E-01 00112
+  0.165068696E-01 0.000000000E+00 0.000000000E+00 0.666969645     10003
+ -0.495401630E-01 0.000000000E+00 0.000000000E+00  3.03133290     01003
+  0.000000000E+00 0.265342622    -0.214067719     0.000000000E+00 00103
+  0.000000000E+00 -745.760209      400.637642     0.000000000E+00 00031
+  -1.44823908     0.000000000E+00 0.000000000E+00  1.73786094     00022
+  0.000000000E+00 -1.62597887      1.25528356     0.000000000E+00 00013
+  0.169392626E-02 0.000000000E+00 0.000000000E+00 0.195681174     00004
+ ---------------------------------------------------------------------
diff --git a/src/THcHallCSpectrometer.cxx b/src/THcHallCSpectrometer.cxx
index a23020b..fad3810 100644
--- a/src/THcHallCSpectrometer.cxx
+++ b/src/THcHallCSpectrometer.cxx
@@ -44,18 +44,16 @@
 
 #include "THcHallCSpectrometer.h"
 #include "THaTrackingDetector.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
 #include "THaTrack.h"
 #include "THaTrackProj.h"
 #include "THaTriggerTime.h"
 #include "TMath.h"
 #include "TList.h"
 
-#include "TList.h"
-#include "TMath.h"
-
-#ifdef WITH_DEBUG
 #include <iostream>
-#endif
+#include <fstream>
 
 using namespace std;
 
@@ -95,23 +93,194 @@ Bool_t THcHallCSpectrometer::GetTrSorting() const
 }
  
 //_____________________________________________________________________________
-Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
+void THcHallCSpectrometer::InitializeReconstruction()
+{
+  fNReconTerms = 0;
+  for(Int_t i=0;i<fNReconTerms;i++) {
+    for(Int_t j=0;j<4;j++) {
+      fReconCoeff[i][j] = 0.0;
+    }
+    for(Int_t j=0;j<5;j++) {
+      fReconExponents[i][j] = 0;
+    }
+  }
+  fAngSlope_x = 0.0;
+  fAngSlope_y = 0.0;
+  fAngOffset_x = 0.0;
+  fAngOffset_y = 0.0;
+  fDetOffset_x = 0.0;
+  fDetOffset_y = 0.0;
+  fZTrueFocus = 0.0;
+}
+//_____________________________________________________________________________
+Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
 {
-  // Reconstruct target coordinates for all tracks found in the focal plane.
 
-  TIter nextTrack( fTrackingDetectors );
+  static const char* const here = "THcHallCSpectrometer::ReadDatabase";
 
-  nextTrack.Reset();
-  while( THaTrackingDetector* theTrackDetector =
-	 static_cast<THaTrackingDetector*>( nextTrack() )) {
 #ifdef WITH_DEBUG
-    if( fDebug>1 ) cout << "Call FineTrack() for " 
-			<< theTrackDetector->GetName() << "... ";
+  cout << "In THcHallCSpectrometer::ReadDatabase()" << endl;
 #endif
-    theTrackDetector->FindVertices( tracks );
-#ifdef WITH_DEBUG
-    if( fDebug>1 ) cout << "done.\n";
+
+  // Get the matrix element filename from the variable store
+  // Read in the matrix
+
+  InitializeReconstruction();
+
+  char prefix[2];
+
+  cout << " GetName() " << GetName() << endl;
+
+  prefix[0]=tolower(GetName()[0]);
+  prefix[1]='\0';
+
+  string reconCoeffFilename;
+  DBRequest list[]={
+    {"_recon_coeff_filename",&reconCoeffFilename,kString},
+    {"theta_offset",&fThetaOffset,kDouble},
+    {"phi_offset",&fPhiOffset,kDouble},
+    {"delta_offset",&fDeltaOffset,kDouble},
+    {"pcentral",&fPCentral,kDouble},
+    {0}
+  };
+  gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+  ifstream ifile;
+  ifile.open(reconCoeffFilename.c_str());
+  if(!ifile.is_open()) {
+    Error(here, "error opening reconstruction coefficient file %s",reconCoeffFilename.c_str());
+    return kInitError; // Is this the right return code?
+  }
+  
+  string line="!";
+  int good=1;
+  while(good && line[0]=='!') {
+    good = getline(ifile,line).good();
+    cout << line << endl;
+  }
+  // Read in focal plane rotation coefficients
+  // Probably not used, so for now, just paste in fortran code as a comment
+  while(good && line.compare(0,4," ---")!=0) {
+    //  if(line(1:13).eq.'h_ang_slope_x')read(line,1201,err=94)h_ang_slope_x
+    //  if(line(1:13).eq.'h_ang_slope_y')read(line,1201,err=94)h_ang_slope_y
+    //  if(line(1:14).eq.'h_ang_offset_x')read(line,1201,err=94)h_ang_offset_x
+    //  if(line(1:14).eq.'h_ang_offset_y')read(line,1201,err=94)h_ang_offset_y
+    //  if(line(1:14).eq.'h_det_offset_x')read(line,1201,err=94)h_det_offset_x
+    //  if(line(1:14).eq.'h_det_offset_y')read(line,1201,err=94)h_det_offset_y
+    //  if(line(1:14).eq.'h_z_true_focus')read(line,1201,err=94)h_z_true_focus
+    good = getline(ifile,line).good();
+  }
+  // Read in reconstruction coefficients and exponents
+  line=" ";
+  good = getline(ifile,line).good();
+  cout << line << endl;
+  fNReconTerms = 0;
+  cout << "Reading matrix elements" << endl;
+  while(good && line.compare(0,4," ---")!=0) {
+    if(fNReconTerms >= fMaxReconElements) {
+      Error(here, "too much data in reconstruction coefficient file %s",reconCoeffFilename.c_str());
+      return kInitError; // Is this the right return code?
+    }
+#if 0
+    Double_t x,y,z,t;
+    Int_t a,b,c,d,e;
+    sscanf(line.c_str()," %le %le %le %le %1d%1d%1d%1d%1d",&x,&y,&z,&t,
+	   &a,&b,&c,&d,&e);
+    cout << x << " " << y << " " << z << " " << t << " "
+	 << a << b << c << d << endl;
 #endif
+    sscanf(line.c_str()," %le %le %le %le %1d%1d%1d%1d%1d"
+	   ,&fReconCoeff[fNReconTerms][0],&fReconCoeff[fNReconTerms][1]
+	   ,&fReconCoeff[fNReconTerms][2],&fReconCoeff[fNReconTerms][3]
+	   ,&fReconExponents[fNReconTerms][0]
+	   ,&fReconExponents[fNReconTerms][1]
+	   ,&fReconExponents[fNReconTerms][2]
+	   ,&fReconExponents[fNReconTerms][3]
+	   ,&fReconExponents[fNReconTerms][4]);
+    // Parse line into fReconCoeff[fNReconTerms][0 - 3] and
+    //                 fReconExponents[fNReconTerms][0 - 5]
+    fNReconTerms++;
+    good = getline(ifile,line).good();    
+  }
+#if 0
+  cout << fNReconTerms << " Reconstruction terms" << endl;
+  for (Int_t i=0;i<fNReconTerms;i++) {
+    cout << fReconCoeff[i][0] << " " << fReconCoeff[i][1] << " "
+	 << fReconCoeff[i][2] << " " << fReconCoeff[i][3] << " "
+	 << fReconExponents[i][0] << fReconExponents[i][1]
+	 << fReconExponents[i][2] << fReconExponents[i][3]
+	 << fReconExponents[i][4] << endl;
+  }
+#endif
+  if(!good) {
+    Error(here, "error processing reconstruction coefficient file %s",reconCoeffFilename.c_str());
+    return kInitError; // Is this the right return code?
+  }
+
+  return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
+{
+  // Reconstruct target coordinates for all tracks found in the focal plane.
+
+  // In Hall A, this is passed off to the tracking detectors.
+  // In Hall C, we do the target traceback here since the traceback should
+  // not depend on which tracking detectors are used.
+
+  cout << "In THcHallCSpectrometer::FindVertices" << endl;
+  for (Int_t it=0;it<tracks.GetLast()+1;it++) {
+    THaTrack* track = static_cast<THaTrack*>( tracks[it] );
+
+    Double_t hut[5];
+    Double_t hut_rot[5];
+    
+    hut[0] = track->GetX()/100.0 + fZTrueFocus*track->GetTheta() + fDetOffset_x;//m
+    hut[1] = track->GetTheta() + fAngOffset_x;//radians
+    hut[2] = track->GetY()/100.0 + fZTrueFocus*track->GetPhi() + fDetOffset_y;//m
+    hut[3] = track->GetPhi() + fAngOffset_y;//radians
+
+    Double_t gbeam_y = 0.0;// This will be the y position from the fast raster
+
+    hut[4] = -gbeam_y/100.0;
+
+    // Retrieve the focal plane coordnates
+    // Do the transpormation
+    // Stuff results into track
+    hut_rot[0] = hut[0];
+    hut_rot[1] = hut[1] + hut[0]*fAngSlope_x;
+    hut_rot[2] = hut[2];
+    hut_rot[3] = hut[3] + hut[2]*fAngSlope_y;
+    hut_rot[4] = hut[4];
+
+    // Compute COSY sums
+    Double_t sum[4];
+    for(Int_t k=0;k<4;k++) {
+      sum[k] = 0.0;
+    }
+    for(Int_t iterm=0;iterm<fNReconTerms;iterm++) {
+      Double_t term=1.0;
+      for(Int_t j=0;j<5;j++) {
+	if(fReconExponents[iterm][j]!=0) {
+	  term *= pow(hut_rot[j],fReconExponents[iterm][j]);
+	}
+      }
+      for(Int_t k=0;k<4;k++) {
+	sum[k] += term*fReconCoeff[iterm][k];
+      }
+    }
+    // Transfer results to track
+    // No beam raster yet
+    //; In transport coordinates phi = hyptar = dy/dz and theta = hxptar = dx/dz 
+    //;    but for unknown reasons the yp offset is named  htheta_offset
+    //;    and  the xp offset is named  hphi_offset
+
+    track->SetTarget(0.0, sum[1]*100.0, sum[0]+fPhiOffset, sum[2]+fThetaOffset);
+    track->SetDp(sum[3]*100.0+fDeltaOffset);	// Percent.  (Don't think podd cares if it is % or fraction)
+    // There is an hpcentral_offset that needs to be applied somewhere.
+    // (happly_offs)
+    track->SetMomentum(fPCentral*(1+track->GetDp()/100.0));
   }
 
   // If enabled, sort the tracks by chi2/ndof
@@ -157,16 +326,19 @@ Int_t THcHallCSpectrometer::TrackTimes( TClonesArray* Tracks ) {
   // To be useful, a meaningful timing resolution should be assigned
   // to each Scintillator object (part of the database).
   
+#if 0
   if ( !Tracks ) return -1;
-  
+
   THaTrack *track=0;
   Int_t ntrack = GetNTracks();
 
+  THaTrack *track=0;
+
   // linear regression to:  t = t0 + pathl/(beta*c)
   //   where t0 is the time of the track at the reference plane (sc_ref).
   //   t0 and beta are solved for.
   //
-#if 0
+
   for ( Int_t i=0; i < ntrack; i++ ) {
     track = static_cast<THaTrack*>(Tracks->At(i));
     THaTrackProj* tr_ref = static_cast<THaTrackProj*>
diff --git a/src/THcHallCSpectrometer.h b/src/THcHallCSpectrometer.h
index f560bbf..356e6cf 100644
--- a/src/THcHallCSpectrometer.h
+++ b/src/THcHallCSpectrometer.h
@@ -17,14 +17,35 @@ public:
   THcHallCSpectrometer( const char* name, const char* description );
   virtual ~THcHallCSpectrometer();
 
+  virtual Int_t   ReadDatabase( const TDatime& date );
   virtual Int_t   FindVertices( TClonesArray& tracks );
   virtual Int_t   TrackCalc();
   virtual Int_t   TrackTimes( TClonesArray* tracks );
 
   Bool_t SetTrSorting( Bool_t set = kFALSE );
   Bool_t GetTrSorting() const;
-  
+
 protected:
+  void InitializeReconstruction();
+
+  // Should look at the ThaMatrixElement class in THaVDC.h for better way
+  // to store matrix element data
+#define fMaxReconElements 1000
+  Int_t fNReconTerms;
+  Double_t fReconCoeff[fMaxReconElements][4];
+  Int_t fReconExponents[fMaxReconElements][5];
+  Double_t fAngSlope_x;
+  Double_t fAngSlope_y;
+  Double_t fAngOffset_x;
+  Double_t fAngOffset_y;
+  Double_t fDetOffset_x;
+  Double_t fDetOffset_y;
+  Double_t fZTrueFocus;
+  Double_t fThetaOffset;
+  Double_t fPhiOffset;
+  Double_t fDeltaOffset;
+  Double_t fPCentral;
+
   //  THaScintillator *sc_ref;  // calculate time track hits this plane
 
   // Flag for fProperties indicating that tracks are to be sorted by chi2
-- 
GitLab