• [X] arraymancer NN DSL PR

• [X] write mail to CAST PC about talk at Patras.
• [ ] implement MCMC multiple chains starting & mean value of them for limit calculation
• [ ] fix segfault when multithreading
• [X] compute correct "depth" for raytracing focal spot
• [X] split window strongback from signal & position uncertainty Implemented the axion image w/o window & strongback separately. Still have to be implemented into limit calc.

• [X] implement MCMC multiple chains starting & mean value of them for limit calculation
• [ ] fix segfault when multithreading
• [X] implement strongback / signal split into limit calc
• [ ] Timepix3 background rate!

Questions for meeting with Klaus today:

1. Did you hear something from Igor? -> Nope he hasn't either. Apparently Igor is very busy currently. But Klaus doesn't think there will be any showstoppers regarding making the data available.
2. For reference distributions and logL morphing: We morph bin wise on pre-binned data. This leads to jumps in the logL cut value. Maybe a good idea after all not use a histogram, but a smooth KDE? Unbinned is not directly possible, because we don't have data to compute an unbinned distribution for everything outside main fluorescence lines! -> Klaus had a good idea here: We can estimate the systematic effect of our binning by moving the bin edges by half a bin width to the left / right and computing the expected limit based on these. If the effective limit changes, we know there is some systematic effect going on. More likely though, the expected limit remains unchanged (within variance) and therefore the systematic impact is smaller than the variance of the limit.

3. About septem veto and line veto: What to do with random coincidences? Is it honest to use those clusters? -> Klaus had an even better idea here: we can estimate the dead time by doing the following:

   1. read full septemboard data
   2. shuffle center + outer chip event data around such that we know the two are not correlated
   3. compute the efficiency of the septem veto.

   In theory 0% of all events should trigger either the septem or the line veto. The percentage that does anyway is our random coincidence!

All files that were in /tmp/playground (/t/playground) referenced here and in the meeting notes are backed up in ~/development_files/07_03_2023/playground (to not make sure we lose something / for reference to recreate some in development behavior etc.)

Just because I'll likely shut down the computer for the first time in 26 days soon and not sure if everything was backed up from there. I believe so, but who knows.

Let's rerun the likelihood after adding the tracking information back to the H5 files and fixing how the total duration is calculated from the data files.

Previously we used the total duration in every case, even when excluding tracking information and thus having less time in actuality. 'Fortunately', all background rate plots in the thesis as of today ran without any tracking info in the H5 files, meaning they include the solar tracking itself. Therefore the total duration is correct in those cases.

Run-2 testing (all vetoes):

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/playground/test_run2.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --lineveto --scintiveto --fadcveto --septemveto \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

Run-3 testing (all vetoes):

likelihood -f ~/CastData/data/DataRuns2018_Reco.h5 \
           --h5out /tmp/playground/test_run3.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --lineveto --scintiveto --fadcveto --septemveto \
           --calibFile ~/CastData/data/CalibrationRuns2018_Reco.h5

The likelihood outputs are here: ./resources/background_rate_test_correct_time_no_tracking/

Background:

plotBackgroundRate \
 /tmp/playground/test_run2.h5 \
 /tmp/playground/test_run3.h5 \
 --combName 2017/18 \
 --combYear 2017 \
 --centerChip 3 \
 --title "Background rate from CAST data, all vetoes" \
 --showNumClusters \
 --showTotalTime \
 --topMargin 1.5 \
 --energyDset energyFromCharge \
 --outfile background_rate_crGold_all_vetoes.pdf \
 --outpath /tmp/playground/ \
 --quiet

The number coming out (see title of the generated plot) is now 3158.57 h, which matches our (new :( ) expectation.

which is also in the same directory: ./resources/background_rate_test_correct_time_no_tracking/background_rate_crGold_all_vetoes.pdf

• [X] Rerun writeRunList and update the statusAndProgress and thesis tables about times!
• [ ]

  Update data in thesis! Run-2:

  ./writeRunList -b ~/CastData/data/DataRuns2017_Reco.h5 -c ~/CastData/data/CalibrationRuns2017_Reco.h5
  

  Type: rtBackground total duration: 14 weeks, 6 days, 11 hours, 25 minutes, 59 seconds, 97 milliseconds, 615 microseconds, and 921 nanoseconds In hours: 2507.433082670833 active duration: 2238.783333333333 trackingDuration: 4 days, 10 hours, and 20 seconds In hours: 106.0055555555556 active tracking duration: 94.12276972527778 nonTrackingDuration: 14 weeks, 2 days, 1 hour, 25 minutes, 39 seconds, 97 milliseconds, 615 microseconds, and 921 nanoseconds In hours: 2401.427527115278 active background duration: 2144.666241943055

  Solar tracking [h]	Background [h]	Active tracking [h]	Active background [h]	Total time [h]	Active time [h]
  106.006	2401.43	94.1228	2144.67	2507.43	2238.78

  Type: rtCalibration total duration: 4 days, 11 hours, 25 minutes, 20 seconds, 453 milliseconds, 596 microseconds, and 104 nanoseconds In hours: 107.4223482211111 active duration: 2.601388888888889 trackingDuration: 0 nanoseconds In hours: 0.0 active tracking duration: 0.0 nonTrackingDuration: 4 days, 11 hours, 25 minutes, 20 seconds, 453 milliseconds, 596 microseconds, and 104 nanoseconds In hours: 107.4223482211111 active background duration: 2.601391883888889

  Solar tracking [h]	Background [h]	Active tracking [h]	Active background [h]	Total time [h]	Active time [h]
  0	107.422	0	2.60139	107.422	2.60139

  Run-3:

  ./writeRunList -b ~/CastData/data/DataRuns2018_Reco.h5 -c ~/CastData/data/CalibrationRuns2018_Reco.h5
  

  Type: rtBackground total duration: 7 weeks, 23 hours, 13 minutes, 35 seconds, 698 milliseconds, 399 microseconds, and 775 nanoseconds In hours: 1199.226582888611 active duration: 1079.598333333333 trackingDuration: 3 days, 2 hours, 17 minutes, and 53 seconds In hours: 74.29805555555555 active tracking duration: 66.92306679361111 nonTrackingDuration: 6 weeks, 4 days, 20 hours, 55 minutes, 42 seconds, 698 milliseconds, 399 microseconds, and 775 nanoseconds In hours: 1124.928527333056 active background duration: 1012.677445774444

  Solar tracking [h]	Background [h]	Active tracking [h]	Active background [h]	Total time [h]	Active time [h]
  74.2981	1124.93	66.9231	1012.68	1199.23	1079.6

  Type: rtCalibration total duration: 3 days, 15 hours, 3 minutes, 47 seconds, 557 milliseconds, 131 microseconds, and 279 nanoseconds In hours: 87.06321031416667 active duration: 3.525555555555556 trackingDuration: 0 nanoseconds In hours: 0.0 active tracking duration: 0.0 nonTrackingDuration: 3 days, 15 hours, 3 minutes, 47 seconds, 557 milliseconds, 131 microseconds, and 279 nanoseconds In hours: 87.06321031416667 active background duration: 3.525561761944445

  Solar tracking [h]	Background [h]	Active tracking [h]	Active background [h]	Total time [h]	Active time [h]
  0	87.0632	0	3.52556	87.0632	3.52556

• [X]

  Rerun the createAllLikelihoodCombinations now that tracking information is there. -> Currently running (Update: <2023-03-10 Fri 21:43> We could now combine the below with the one further down that excludes the FADC!).

  ./createAllLikelihoodCombinations \
      --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
      --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
      --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --regions crGold --regions crAll \
      --vetoSets "{fkScinti, fkFadc, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
      --fadcVetoPercentiles 0.9 --fadcVetoPercentiles 0.95 --fadcVetoPercentiles 0.99 \
      --out /t/lhood_outputs_adaptive_fadc \
      --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --multiprocessing
  

  Found here: ./resources/lhood_limits_automation_correct_duration/

• [ ]

  Now generate the other likelihood outputs we need for more expected limit cases from sec. [BROKEN LINK: sec:meetings:10_03_23] in StatusAndProgress:

  • [ ] Calculate expected limits also for the following cases:
    • [X] Septem, line combinations without the FADC
    • [ ] Best case (lowest row of below) with lnL efficiencies of:
      • [ ] 0.7
      • [ ] 0.9

  The former (septem, line without FADC) will be done using createAllLikelihoodCombinations:

  ./createAllLikelihoodCombinations \
      --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
      --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
      --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --regions crGold --regions crAll \
      --vetoSets "{+fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
      --out /t/lhood_outputs_adaptive_fadc \
      --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --multiprocessing
  

  For simplicity, this will regenerate some of the files already generated (i.e. the no vetoes & the scinti case)

  These files are also found here: ./resources/lhood_limits_automation_correct_duration/ together with a rerun of the regular cases above.

Plot the background clusters to see if we indeed have less over the whole chip.

plotBackgroundClusters \
    /t/lhood_outputs_adaptive_duplicated_fadc_stuff/likelihood_cdl2018_Run2_crAll_scinti_vetoPercentile_0.99_fadc_vetoPercentile_0.99_septem_vetoPercentile_0.99_line_vetoPercentile_0.99.h5 \
    /t/lhood_outputs_adaptive_duplicated_fadc_stuff/likelihood_cdl2018_Run3_crAll_scinti_vetoPercentile_0.99_fadc_vetoPercentile_0.99_septem_vetoPercentile_0.99_line_vetoPercentile_0.99.h5 \
    --zMax 5 \
    --title "X-ray like clusters of CAST data after all vetoes" \
    --outpath /tmp/playground/ \
    --filterNoisyPixels \
    --filterEnergy 12.0 \
    --suffix "_all_vetoes"

Available here: resources/background_rate_test_correct_time_no_tracking/background_cluster_centers_all_vetoes.pdf

Where we can see that indeed we now have less than 10,000 clusters compared to the ~10,500 we had when using all data (including tracking).

Continue from yesterday:

• [X]

  Now generate the other likelihood outputs we need for more expected limit cases from sec. [BROKEN LINK: sec:meetings:10_03_23] in StatusAndProgress: -> All done, path to files below.

  • [X] Calculate expected limits also for the following cases:
    • [X] Septem, line combinations without the FADC (done yesterday)
    • [X] Best case (lowest row of below) with lnL efficiencies of:
      • [X] 0.7
      • [X] 0.9

  The latter has now also been implemented as functionality in likelihood and createAllLikelihoodCombinations (adjust signal efficiency from command line and add options to runner). Now run:

  ./createAllLikelihoodCombinations \
      --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
      --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
      --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --regions crGold --regions crAll \
      --signalEfficiency 0.7 --signalEfficiency 0.9 \
      --vetoSets "{+fkScinti, +fkFadc, +fkSeptem, fkLineVeto}" \
      --fadcVetoPercentile 0.9 \
      --out /t/lhood_outputs_adaptive_fadc \
      --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --multiprocessing \
      --dryRun
  

  to reproduce numbers for the best expected limit case together with a lnL signal efficiency of 70 and 90%.

  Finally, these files are also in: ./resources/lhood_limits_automation_correct_duration/ which means now all the setups we initially care about are there.

Let's look at the background rate we get from the 70% vs the 90% case:

plotBackgroundRate \
    likelihood_cdl2018_Run2_crGold_signalEff_0.7_scinti_fadc_septem_line_vetoPercentile_0.9.h5 \
    likelihood_cdl2018_Run3_crGold_signalEff_0.7_scinti_fadc_septem_line_vetoPercentile_0.9.h5 \
    likelihood_cdl2018_Run2_crGold_signalEff_0.9_scinti_fadc_septem_line_vetoPercentile_0.9.h5 \
    likelihood_cdl2018_Run3_crGold_signalEff_0.9_scinti_fadc_septem_line_vetoPercentile_0.9.h5 \
    --names "0.7" --names "0.7" \
    --names "0.9" --names "0.9" \
    --centerChip 3 \
    --title "Background rate from CAST data, incl. all vetoes, 70vs90" \
    --showNumClusters --showTotalTime \
    --topMargin 1.5 --energyDset energyFromCharge \
    --outfile background_rate_cast_all_vetoes_70p_90p.pdf \
    --outpath . \
    --quiet

[INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 4.1861e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 3.4884e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 9.3221e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 7.7684e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 8.6185e-06 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 4.3093e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 1.4775e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 7.3873e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.8116e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 4.0259e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 3.4650e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 7.7000e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 1.4423e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 5.7691e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 2.4273e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 9.7090e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 4.3972e-06 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 1.0993e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 1.1785e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 2.9461e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 2.7790e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 3.4738e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 5.3998e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 6.7497e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.7 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 1.3895e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 2.3159e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: 0.9 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 3.0956e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 5.1594e-06 keV⁻¹·cm⁻²·s⁻¹

yielded which shows quite the incredibly change especially in the 8 keV peak!

And in the 4 to 8 keV range we even almost got to the 1e-7 range for the 70% case (however note that in this case the total efficiency is only about 40% or so!).

• [X] Verify that the signal efficiency used is written to output logL file
  • [X] If not implement -> It was not, now implemented.
• [X] read signal efficiency from logL file in mcmc_limit and stop using the efficiency including ε in the context. Instead merge with the calculator for veto efficiencies. -> Implemented.

From meeting notes:

• [ ] Verify that those elements with lower efficiency indeed have \(R_T = 0\) at higher values! -> Just compute \(R_T = 0\) for all input files and output result, easiest.

Old model from March 2022 ./resources/mlp_trained_march2022.pt :

Test set: Average loss: 0.9876 | Accuracy: 0.988
Cut value: 2.483305978775025

Test set: Average loss: 0.9876 | Accuracy: 0.988
Total efficiency = 0.8999892098945267

Test set: Average loss: 0.9995 | Accuracy: 0.999
Target: Ag-Ag-6kV eff = 0.9778990694345026

Test set: Average loss: 0.9916 | Accuracy: 0.992
Target: Al-Al-4kV eff = 0.9226669690441093

Test set: Average loss: 0.9402 | Accuracy: 0.940
Target: C-EPIC-0.6kV eff = 0.6790938280413843

Test set: Average loss: 0.9941 | Accuracy: 0.994
Target: Cu-EPIC-0.9kV eff = 0.8284986713906112

Test set: Average loss: 0.9871 | Accuracy: 0.987
Target: Cu-EPIC-2kV eff = 0.8687534321801208

Test set: Average loss: 1.0000 | Accuracy: 1.000
Target: Cu-Ni-15kV eff = 0.9939449541284404

Test set: Average loss: 0.9999 | Accuracy: 1.000
Target: Mn-Cr-12kV eff = 0.9938112429087158

Test set: Average loss: 1.0000 | Accuracy: 1.000
Target: Ti-Ti-9kV eff = 0.9947166683932456

New model from yesterday ./resources/mlp_trained_bsz8192_hidden_5000.pt :

Test set: Average loss: 0.9714 | Accuracy: 0.971
Cut value: 1.847297704219818

Test set: Average loss: 0.9714 | Accuracy: 0.971
Total efficiency = 0.8999892098945267

Test set: Average loss: 0.9945 | Accuracy: 0.994
Target: Ag-Ag-6kV eff = 0.9525769506084467

Test set: Average loss: 0.9804 | Accuracy: 0.980
Target: Al-Al-4kV eff = 0.9097403333711305

Test set: Average loss: 0.8990 | Accuracy: 0.899
Target: C-EPIC-0.6kV eff = 0.7640920442383161

Test set: Average loss: 0.9584 | Accuracy: 0.958
Target: Cu-EPIC-0.9kV eff = 0.8211913197519929

Test set: Average loss: 0.9636 | Accuracy: 0.964
Target: Cu-EPIC-2kV eff = 0.8543657331136738

Test set: Average loss: 0.9980 | Accuracy: 0.998
Target: Cu-Ni-15kV eff = 0.981651376146789

Test set: Average loss: 0.9978 | Accuracy: 0.998
Target: Mn-Cr-12kV eff = 0.9807117070654977

Test set: Average loss: 0.9982 | Accuracy: 0.998
Target: Ti-Ti-9kV eff = 0.9776235367243344

Calculated with determineCdlEfficiency.

To get a background rate estimate we use the simple functionality in the NN training tool itself. Note that it needs the total time as an input to scale the data correctly (hardcoded has a value of Run-2, but for Run-3 need to use the totalTime argument!)

Background rate at 95% for Run-2 data:

cd ~/CastData/ExternCode/TimepixAnalysis/Tools/NN_playground
./train_ingrid -f ~/CastData/data/DataRuns2017_Reco.h5 --ε 0.95

which yields:

Background rate at 90% for Run-2 data:

cd ~/CastData/ExternCode/TimepixAnalysis/Tools/NN_playground
./train_ingrid -f ~/CastData/data/DataRuns2017_Reco.h5 --ε 0.9

Test set: Average loss: 0.9714 | Accuracy: 0.971 Cut value: 1.847297704219818

Test set: Average loss: 0.9714 | Accuracy: 0.971 Total efficiency = 0.8999892098945267

Test set: Average loss: 0.9945 | Accuracy: 0.994 Target: Ag-Ag-6kV eff = 0.9525769506084467

Test set: Average loss: 0.9804 | Accuracy: 0.980 Target: Al-Al-4kV eff = 0.9097403333711305

Test set: Average loss: 0.8990 | Accuracy: 0.899 Target: C-EPIC-0.6kV eff = 0.7640920442383161

Test set: Average loss: 0.9584 | Accuracy: 0.958 Target: Cu-EPIC-0.9kV eff = 0.8211913197519929

Test set: Average loss: 0.9636 | Accuracy: 0.964 Target: Cu-EPIC-2kV eff = 0.8543657331136738

Test set: Average loss: 0.9980 | Accuracy: 0.998 Target: Cu-Ni-15kV eff = 0.981651376146789

Test set: Average loss: 0.9978 | Accuracy: 0.998 Target: Mn-Cr-12kV eff = 0.9807117070654977

Test set: Average loss: 0.9982 | Accuracy: 0.998 Target: Ti-Ti-9kV eff = 0.9776235367243344

which yields:

Background rate at 80% for Run-2 data:

cd ~/CastData/ExternCode/TimepixAnalysis/Tools/NN_playground
./train_ingrid -f ~/CastData/data/DataRuns2017_Reco.h5 --ε 0.8

Test set: Average loss: 0.9714 | Accuracy: 0.971 Cut value: 3.556154251098633

Test set: Average loss: 0.9714 | Accuracy: 0.971 Total efficiency = 0.799991907420895

Test set: Average loss: 0.9945 | Accuracy: 0.994 Target: Ag-Ag-6kV eff = 0.856209735146743

Test set: Average loss: 0.9804 | Accuracy: 0.980 Target: Al-Al-4kV eff = 0.7955550515931511

Test set: Average loss: 0.8990 | Accuracy: 0.899 Target: C-EPIC-0.6kV eff = 0.6546557260078487

Test set: Average loss: 0.9584 | Accuracy: 0.958 Target: Cu-EPIC-0.9kV eff = 0.7030558015943313

Test set: Average loss: 0.9636 | Accuracy: 0.964 Target: Cu-EPIC-2kV eff = 0.7335529928610653

Test set: Average loss: 0.9980 | Accuracy: 0.998 Target: Cu-Ni-15kV eff = 0.9102752293577981

Test set: Average loss: 0.9978 | Accuracy: 0.998 Target: Mn-Cr-12kV eff = 0.9036616812790098

Test set: Average loss: 0.9982 | Accuracy: 0.998 Target: Ti-Ti-9kV eff = 0.8947477468144618

which yields:

Background rate at 70% for Run-2 data:

cd ~/CastData/ExternCode/TimepixAnalysis/Tools/NN_playground
./train_ingrid -f ~/CastData/data/DataRuns2017_Reco.h5 --ε 0.7

Test set: Average loss: 0.9714 | Accuracy: 0.971 Cut value: 5.098100709915161

Test set: Average loss: 0.9714 | Accuracy: 0.971 Total efficiency = 0.6999946049472634

Test set: Average loss: 0.9945 | Accuracy: 0.994 Target: Ag-Ag-6kV eff = 0.7380100214745884

Test set: Average loss: 0.9804 | Accuracy: 0.980 Target: Al-Al-4kV eff = 0.6931624900782402

Test set: Average loss: 0.8990 | Accuracy: 0.899 Target: C-EPIC-0.6kV eff = 0.5882090617195861

Test set: Average loss: 0.9584 | Accuracy: 0.958 Target: Cu-EPIC-0.9kV eff = 0.6312001771479185

Test set: Average loss: 0.9636 | Accuracy: 0.964 Target: Cu-EPIC-2kV eff = 0.6507413509060955

Test set: Average loss: 0.9980 | Accuracy: 0.998 Target: Cu-Ni-15kV eff = 0.7844036697247706

Test set: Average loss: 0.9978 | Accuracy: 0.998 Target: Mn-Cr-12kV eff = 0.7790613718411552

Test set: Average loss: 0.9982 | Accuracy: 0.998 Target: Ti-Ti-9kV eff = 0.7758209882937946

which yields:

• [X] Add background rates for new model using 95%, 90%, 80%
  • [X] And the outputs as above for the local efficiencies
• [ ] Implement selection of global vs. local target efficiency
• [ ] Check the efficiency we get when applying the model to the 55Fe calibration data (need same efficiency!). Cross check our helper program that does this for lnL method
• [ ] Check the background rate from the Run-3 data. Is it compatible? Or does it break down?

Practical:

• [ ] Move the model logic over to likelihood.nim as a replacement for the fkAggressive veto
  • [ ] including the selection of the target efficiency
• [X] Clean up veto system in likelihood for better insertion of NN
  • [X] add lnL as a form of veto
• [X] add vetoes for MLP and ConvNet (in principle)
• [ ] move NN code to main ingrid module
• [ ] make MLP / ConvNet types accessible in likelihood if compiled on cpp backend (and with CUDA?)
• [ ] add path to model file
• [ ] adjust CutValueInterpolator to make it work for both lnL as well as NN. Idea is the same!

Questions:

• [ ] Is there still a place for something like an equivalent for the likelihood morphing? In this case based on likely just interpolating the cut values?

TODOs from yesterday:

• [X] Add background rates for new model using 95%, 90%, 80%
  • [X] And the outputs as above for the local efficiencies
• [X] Implement selection of global vs. local target efficiency
• [X] Check the efficiency we get when applying the model to the 55Fe calibration data (need same efficiency!). Cross check our helper program that does this for lnL method -> Wrote effective_eff_55fe.nim in NN_playground -> Efficiency in 55Fe data is abysmal! ~40-55 % at 95% !
• [ ] Check the background rate from the Run-3 data. Is it compatible? Or does it break down?

Practical:

• [ ] Move the model logic over to likelihood.nim as a replacement for the fkAggressive veto
  • [ ] including the selection of the target efficiency
• [X] Clean up veto system in likelihood for better insertion of NN
  • [X] add lnL as a form of veto
• [X] add vetoes for MLP and ConvNet (in principle)
• [ ] move NN code to main ingrid module
• [X] make MLP / ConvNet types accessible in likelihood if compiled on cpp backend (and with CUDA?)
• [X] add path to model file
• [X] adjust CutValueInterpolator to make it work for both lnL as well as NN. Idea is the same!

Questions:

• [ ] Is there still a place for something like an equivalent for the likelihood morphing? In this case based on likely just interpolating the cut values?

With the refactor of likelihood we can now do things like disable the lnL cut itself and only use the vetoes.

NOTE: All outputs below that were placed in /t/testing can be found in ./resources/nn_testing_outputs/.

For example look at only using the FADC (with a much harsher cut than usual):

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_fadc.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --fadcveto \
           --vetoPercentile 0.75 \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5
plotBackgroundRate /t/testing/test_run2_only_fadc.h5 \
                   --combName "onlyFadc" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only FADC veto" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_fadc_veto.pdf \
                   --outpath /t/testing/ --quiet

The plot is here: The issue is there's a couple of runs that have no FADC / in which the FADC was extremely noisy, hence all we really see is the background distribution of those.

But for the more interesting stuff, let's try to create the background rate using the NN veto at 95% efficiency!:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_0.95.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.95 \
           --nnCutKind global \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_0.95.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ 95%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_0.95.pdf \
                   --outpath /t/testing/ --quiet

At 90% global efficiency:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_0.9.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.9 \
           --nnCutKind global \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_0.9.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ 90%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_0.9.pdf \
                   --outpath /t/testing/ --quiet

At 80% global efficiency:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_0.8.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.8 \
           --nnCutKind global \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_0.8.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ 90%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_0.8.pdf \
                   --outpath /t/testing/ --quiet

At 70% global efficiency:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_0.7.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.7 \
           --nnCutKind global \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_0.7.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ 90%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_0.7.pdf \
                   --outpath /t/testing/ --quiet

NOTE: Make sure to set neuralNetCutKind to local in the config file! And local 95%:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_local_0.95.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.95 \
           --nnCutKind local \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_local_0.95.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ local 95%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_local_0.95.pdf \
                   --outpath /t/testing/ --quiet

And local 80%:

likelihood -f ~/CastData/data/DataRuns2017_Reco.h5 \
           --h5out /tmp/testing/test_run2_only_mlp_local_0.8.h5 \
           --region crGold \
           --cdlYear 2018 \
           --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
           --mlp ~/org/resources/mlp_trained_bsz8192_hidden_5000.pt \
           --nnSignalEff 0.8 \
           --nnCutKind local \
           --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5

plotBackgroundRate /t/testing/test_run2_only_mlp_local_0.8.h5 \
                   --combName "onlyMLP" \
                   --combYear 2017 \
                   --centerChip 3 \
                   --title "Background rate from CAST data, only MLP @ local 80%" \
                   --showNumClusters \
                   --showTotalTime \
                   --topMargin 1.5 \
                   --energyDset energyFromCharge \
                   --outfile background_rate_only_mlp_local_0.8.pdf \
                   --outpath /t/testing/ --quiet

Continue on from yesterday:

• [ ] Implement 55Fe calibration data into the training process. E.g. add about 1000 events per calibration run to the training data as signal target to have a wider distribution of what real X-rays should look like. Hopefully that increases our efficiency!
• [ ]

  It seems like only very few events pass the cuts in readCalibData (e.g. what we use for effective efficiency check and in mixed data training). Why is that? Especially for escape peak often less than 300 events are valid! So little statistics, really? Looking at spectra, e.g. in ~/CastData/ExternCode/TimepixAnalysis/Analysis/ingrid/out/CalibrationRuns2018_Raw_2020-04-28_15-06-54 there is really this little statistics in the escape peak often. (peaks at less than per bin 50!) How do these spectra look without any cuts? Are our cuts rubbish? Quick look:

  plotData --h5file ~/CastData/data/CalibrationRuns2017_Reco.h5 \
           --runType rtCalibration \
           --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
           --region crSilver \
           --ingrid --separateRuns
  

  seems to support that there simply isn't much more statistics available!

• [X]

  First training with a mixed set of data, using, per run:

  • min(500, total) escape peak (after cuts)
  • min(500, total) photo peak (after cuts)
  • 6000 background
  • all CDL

  -> all of these are of course shuffled and then split into training and test datasets The resulting model is in: ./resources/nn_devel_mixing/trained_mlp_mixed_data.pt

  The generated plots are in: ./Figs/statusAndProgress/neuralNetworks/development/mixing_data/

  Looking at these figures we can see mainly that the ROC curve is extremely 'clean', but fitting for the separation seen in the training and validation output distributions.

  • [ ] effective efficiencies for 55Fe
  • [ ] efficiencies of CDL data!
• [X] make loss / accuracy curves log10
• [ ] Implement snapshots of the model during training whenever the training and test (or only test) accuracy improves

As discussed in the meeting today (sec. [BROKEN LINK: sec:meetings:17_03_23] in notes), let's rerun all expected limits and add the new two, namely:

• [ ] redo all expected limit calculations with the following new cases:
  • 0.9 lnL + scinti + FADC@0.98 + line
  • 0.8 lnL + scinti + FADC@0.98 + line ε_cut:
    • 1.0, 1.2, 1.4, 1.6

The standard cases (lnL 80 + all veto combinations with different FADC settings):

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --vetoSets "{fkLogL, fkScinti, fkFadc, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --fadcVetoPercentiles 0.9 --fadcVetoPercentiles 0.95 --fadcVetoPercentiles 0.99 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing

The no septem veto + different lnL efficiencies:

• [X] 0.9 lnL + scinti + FADC@0.98 + line

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --vetoSets "{+fkLogL, +fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --signalEfficiency 0.7 --signalEfficiency 0.9 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing

The older case of changing lnL efficiency with different lnL efficiency

• [X]

  add a case with less extreme FADC veto

  ./createAllLikelihoodCombinations \
      --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
      --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
      --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --regions crGold --regions crAll \
      --signalEfficiency 0.7 --signalEfficiency 0.9 \
      --vetoSets "{+fkLogL, +fkScinti, +fkFadc, +fkSeptem, fkLineVeto}" \
      --fadcVetoPercentile 0.95 --fadcVetoPercentile 0.99 \
      --out /t/lhood_outputs_adaptive_fadc \
      --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --multiprocessing
  

And finally different eccentricity cutoffs for the line veto:

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --signalEfficiency 0.7 --signalEfficiency 0.8 --signalEfficiency 0.9 \
    --vetoSets "{+fkLogL, +fkScinti, +fkFadc, fkLineVeto}" \
    --eccentricityCutoff 1.0 --eccentricityCutoff 1.2 --eccentricityCutoff 1.4 --eccentricityCutoff 1.6 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing

The output H5 files will be placed in: ./resources/lhood_limits_automation_with_nn_support

With the likelihood output files generated over night in resources/lhood_limits_automation_with_nn_support it's now time to let the limits run. I noticed something else was missing from these files: I forgot to re add the actual vetoes in use to the output (because those were written manually).

• [X] add flags to LikelihoodContext to auto serialize them
• [X] update mcmc limit code to use the new serialized data as veto efficiency and veto usage
• [X] rerun all limits with all the different setups.
• [X] update runLimits to be smarter about what has been done. In principle we can now quit the limit calculation and it should continue automatically on a restart (with the last file worked on!)

The script we actually ran today. This will be part of the thesis (or a variation thereof).

#!/usr/bin/zsh
cd ~/CastData/ExternCode/TimepixAnalysis/Analysis/
./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --vetoSets "{fkLogL, fkScinti, fkFadc, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --fadcVetoPercentiles 0.9 --fadcVetoPercentiles 0.95 --fadcVetoPercentiles 0.99 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing \
    --jobs 12
./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --vetoSets "{+fkLogL, +fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --signalEfficiency 0.7 --signalEfficiency 0.9 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing \
    --jobs 12
./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --signalEfficiency 0.7 --signalEfficiency 0.9 \
    --vetoSets "{+fkLogL, +fkScinti, +fkFadc, +fkSeptem, fkLineVeto}" \
    --fadcVetoPercentile 0.95 --fadcVetoPercentile 0.99 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing \
    --jobs 12
./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crGold --regions crAll \
    --signalEfficiency 0.7 --signalEfficiency 0.8 --signalEfficiency 0.9 \
    --vetoSets "{+fkLogL, +fkScinti, +fkFadc, fkLineVeto}" \
    --eccentricityCutoff 1.0 --eccentricityCutoff 1.2 --eccentricityCutoff 1.4 --eccentricityCutoff 1.6 \
    --out /t/lhood_outputs_adaptive_fadc \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing \
    --jobs 12

(all in ./resources/lhood_limits_automation_with_nn_support/)

And currently running:

./runLimits --path ~/org/resources/lhood_limits_automation_with_nn_support --nmc 1000

Train NN:

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 \
               ~/CastData/data/CalibrationRuns2018_Reco.h5 \
               --back ~/CastData/data/DataRuns2017_Reco.h5 \
               --back ~/CastData/data/DataRuns2018_Reco.h5 \
               --ε 0.95 \
               --modelOutpath /tmp/trained_mlp_mixed_data.pt

NOTE: The one thing I just realized is: the accuracy we print is of course related to the actual prediction of the network, i.e. which of the two output neurons is the maximum value. So maybe our approach of only looking at one and adjusting based on that is just dumb and the network actually is much better than we think?

The numbers we see as accuracy actually make sense. Consider predictBackground output:

Pred set: Average loss: 0.0169 | Accuracy: 0.9956
p inds len 1137 compared to all 260431

The 1137 clusters left after cuts correspond exactly to 99.56% (this is based on using the network's real prediction and not the output + cut value). The question here is: at what signal efficiency is this? From the CDL data it would seem to be at ~99%.

The network we trained today, including checkpoints is here: ./resources/nn_devel_mixing/18_03_23/

• [X] Check what efficiency we get for calibration data instead of background -> Yeah, it is also at over 99% efficiency. So we get a 1e-5 background rate at 99% efficiency. At least that's not too bad.

<2023-03-18 Sat 23:13>: The processed expected limit files are now here: ./resources/lhood_limits_automation_with_nn_support/limits/ with the logL output files in the lhood folder.

We'll continue on later with the processed.txt file as a guide there.

The expected limits for resources/lhood_limits_automation_with_nn_support/limits/ are still running because our processed continuation check was incorrect (looking at full path & not actually skipping files!).

Back to the NN: Let's look at the output of the network for both output neurons. Are they really essentially a mirror of one another?

predictAll in train_ingrid creates a plot of the different data kinds (55Fe, CDL, background) and each neurons output prediction. This yields the following plot by running:

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 \
               ~/CastData/data/CalibrationRuns2018_Reco.h5 \
               --back ~/CastData/data/DataRuns2017_Reco.h5 \
               --back ~/CastData/data/DataRuns2018_Reco.h5 \
               --ε 0.95 \
               --modelOutpath ~/org/resources/nn_devel_mixing/18_03_23/trained_mlp_mixed_datacheckpoint_epoch_95000_loss_0.0117_acc_0.9974.pt \
               --predict

where we can see the following points:

1. the two neurons are almost perfect mirrors, but not exactly
2. selecting the argmax of the two neurons gives us that neuron, which has a positive value almost certainly, due to the mirror nature around 0. It could be different (e.g. both neurons giving a positive or a negative value), but looking at the data this does not seem to happen (if then very rarely).
3. a cut value of 0 should reproduce pretty exactly the standard neural network prediction of picking the argmax

Question: Can the usage of both neurons be beneficial given the small but existing differences in the distributions? Not sure how if so.

An earlier checkpoint (the one before the extreme jump in the loss value based on the loss figure; need to regenerate it, but similar to except as log10) yields the following neuron output:

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 \
               ~/CastData/data/CalibrationRuns2018_Reco.h5 \
               --back ~/CastData/data/DataRuns2017_Reco.h5 \
               --back ~/CastData/data/DataRuns2018_Reco.h5 \
               --ε 0.95 \
               --modelOutpath ~/org/resources/nn_devel_mixing/18_03_23/trained_mlp_mixed_datacheckpoint_epoch_65000_loss_0.0103_acc_0.9977.pt \
               --predict

We can clearly see that at this stage in training the two types of signal data are predicted quite differently! In that sense the latest model is actually much more like what we want, i.e. same prediction for all different kinds of X-rays!

What does the case with the worst loss look like?

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 \
               ~/CastData/data/CalibrationRuns2018_Reco.h5 \
               --back ~/CastData/data/DataRuns2017_Reco.h5 \
               --back ~/CastData/data/DataRuns2018_Reco.h5 \
               --ε 0.95 \
               --modelOutpath ~/org/resources/nn_devel_mixing/18_03_23/trained_mlp_mixed_datacheckpoint_epoch_70000_loss_0.9683_acc_0.9977.pt \
               --predict

Interestingly essentially the same. But the accuracy is the same as before, only the loss is different. Not sure why that might be.

Training the network again after the charge bug was fixed:

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 \
               ~/CastData/data/CalibrationRuns2018_Reco.h5 \
               --back ~/CastData/data/DataRuns2017_Reco.h5 \
               --back ~/CastData/data/DataRuns2018_Reco.h5 \
               --ε 0.95 \
               --modelOutpath /t/nn_training/trained_model_charge_cut_bug_fixed.pt

which are stored here: ./resources/nn_devel_mixing/19_03_23_charge_bug_fixed/ with the generated plots: ./Figs/statusAndProgress/neuralNetworks/development/charge_cut_bug_fixed Looking at the loss plot, at around epoch 83000 the training data started to outpace the test data (test didn't get any worse though and test accuracy improved slightly). Also the ~all_prediction.pdf plot showing how CDL and 55Fe data is predicted is interesting. The CDL data is skewed significantly more to the right than 55Fe, explaining the again prevalent difference in 55Fe efficiency for a given CDL eff:

./effective_eff_55fe \
    -f ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/19_03_23_charge_bug_fixed/trained_model_charge_cut_bug_fixedcheckpoint_epoch_100000_loss_0.0102_acc_0.9978.pt \
    --ε 0.95

Run: 83 for target: signal Keeping : 759 of 916 = 0.8286026200873362

Run: 88 for target: signal Keeping : 763 of 911 = 0.8375411635565313

Run: 93 for target: signal Keeping : 640 of 787 = 0.8132147395171537

Run: 96 for target: signal Keeping : 4591 of 5635 = 0.8147293700088731

Run: 102 for target: signal Keeping : 1269 of 1588 = 0.7991183879093199

Run: 108 for target: signal Keeping : 2450 of 3055 = 0.8019639934533551

Run: 110 for target: signal Keeping : 1244 of 1554 = 0.8005148005148005

Run: 116 for target: signal Keeping : 1404 of 1717 = 0.8177052999417589

Run: 118 for target: signal Keeping : 1351 of 1651 = 0.8182919442761962

Run: 120 for target: signal Keeping : 2784 of 3413 = 0.8157046586580721

Run: 122 for target: signal Keeping : 4670 of 5640 = 0.8280141843971631

Run: 126 for target: signal Keeping : 2079 of 2596 = 0.8008474576271186

Run: 128 for target: signal Keeping : 6379 of 7899 = 0.8075705785542474

Run: 145 for target: signal Keeping : 2950 of 3646 = 0.8091058694459682

Run: 147 for target: signal Keeping : 1670 of 2107 = 0.7925961082107261

Run: 149 for target: signal Keeping : 1536 of 1936 = 0.7933884297520661

Run: 151 for target: signal Keeping : 1454 of 1839 = 0.790647090810223

Run: 153 for target: signal Keeping : 1515 of 1908 = 0.7940251572327044

Run: 155 for target: signal Keeping : 1386 of 1777 = 0.7799662352279122

Run: 157 for target: signal Keeping : 1395 of 1817 = 0.7677490368739681

Run: 159 for target: signal Keeping : 2805 of 3634 = 0.7718767198679142

Run: 161 for target: signal Keeping : 2825 of 3632 = 0.7778083700440529

Run: 163 for target: signal Keeping : 1437 of 1841 = 0.7805540467137425

Run: 165 for target: signal Keeping : 3071 of 3881 = 0.7912909044060809

Run: 167 for target: signal Keeping : 1557 of 2008 = 0.775398406374502

Run: 169 for target: signal Keeping : 4644 of 5828 = 0.7968428277282087

Run: 171 for target: signal Keeping : 1561 of 1956 = 0.7980572597137015

Run: 173 for target: signal Keeping : 1468 of 1820 = 0.8065934065934066

Run: 175 for target: signal Keeping : 1602 of 2015 = 0.7950372208436725

Run: 177 for target: signal Keeping : 1557 of 1955 = 0.7964194373401534

Run: 179 for target: signal Keeping : 1301 of 1671 = 0.7785757031717534

Run: 181 for target: signal Keeping : 2685 of 3426 = 0.7837127845884413

Run: 183 for target: signal Keeping : 2821 of 3550 = 0.7946478873239436

Run: 185 for target: signal Keeping : 3063 of 3856 = 0.7943464730290456

Run: 187 for target: signal Keeping : 2891 of 3616 = 0.7995022123893806

This is for a local efficiency. So once again 95% in CDL correspond to about 80% in 55Fe. Not ideal.

Let's try to train a network that also includes the total charge, so it has some idea of the gas gain in the events.

Otherwise we leave the settings as is:

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --ε 0.95 \
    --modelOutpath /t/nn_training/trained_model_incl_totalCharge.pt

Interestingly, when including the total charge, the loss of the test data remains lower than of the training set! Models: ./resources/nn_devel_mixing/19_03_23_with_total_charge/ and plots: ./Figs/statusAndProgress/neuralNetworks/development/with_total_charge/ Looking at the total charge, we see the same behavior of CDL and 55Fe data essentially. The background distribution has changed a bit.

We could attempt to change our definition of our loss function. Currently we in now way enforce that our result should be close to our target [1, 0] and [0, 1]. Using a MSE loss for example would make sure of that.

Now training with MSE loss. -> Couldn't get anything sensible out of MSE loss. Chatted with BingChat and it couldn't quite help me (different learning rate etc.), but it did suggest to try L1 loss (mean absolute error), which I am running with now.

L1 loss: ./resources/nn_devel_mixing/19_03_23_l1_loss/ ./Figs/statusAndProgress/neuralNetworks/development/l1_loss/ The all prediction plot is interesting. We see the same-ish behavior in this case as in the cross entropy loss. In the training dataset we can even more clearly see two distinct peaks. However, especially the effective efficiencies in the 55Fe data are all over the place:

./effective_eff_55fe \
    -f ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/19_03_23_l1_loss/trained_model_incl_totalCharge_l1_losscheckpoint_epoch_100000_loss_0.0157_acc_0.9971.pt \
    --ε 0.95

Run: 83 for target: signal Keeping : 781 of 916 = 0.8526200873362445

Run: 88 for target: signal Keeping : 827 of 911 = 0.9077936333699231

Run: 93 for target: signal Keeping : 612 of 787 = 0.7776365946632783

Run: 96 for target: signal Keeping : 4700 of 5635 = 0.8340727595385981

Run: 102 for target: signal Keeping : 1292 of 1588 = 0.8136020151133502

Run: 108 for target: signal Keeping : 2376 of 3055 = 0.7777414075286416

Run: 110 for target: signal Keeping : 1222 of 1554 = 0.7863577863577863

Run: 116 for target: signal Keeping : 1453 of 1717 = 0.8462434478741991

Run: 118 for target: signal Keeping : 1376 of 1651 = 0.8334342822531798

Run: 120 for target: signal Keeping : 2966 of 3413 = 0.8690301787283914

Run: 122 for target: signal Keeping : 5049 of 5640 = 0.8952127659574468

Run: 126 for target: signal Keeping : 2157 of 2596 = 0.8308936825885979

Run: 128 for target: signal Keeping : 6546 of 7899 = 0.8287124952525636

Run: 145 for target: signal Keeping : 2729 of 3646 = 0.7484914975315414

Run: 147 for target: signal Keeping : 1517 of 2107 = 0.7199810156620788

Run: 149 for target: signal Keeping : 1152 of 1936 = 0.5950413223140496

Run: 151 for target: signal Keeping : 1135 of 1839 = 0.6171832517672649

Run: 153 for target: signal Keeping : 1091 of 1908 = 0.5718029350104822

Run: 155 for target: signal Keeping : 974 of 1777 = 0.5481148002250985

Run: 157 for target: signal Keeping : 978 of 1817 = 0.5382498624105668

Run: 159 for target: signal Keeping : 2083 of 3634 = 0.5731975784259769

Run: 161 for target: signal Keeping : 2152 of 3632 = 0.5925110132158591

Run: 163 for target: signal Keeping : 1264 of 1841 = 0.6865833785985878

Run: 165 for target: signal Keeping : 2929 of 3881 = 0.7547023962896161

Run: 167 for target: signal Keeping : 1467 of 2008 = 0.7305776892430279

Run: 169 for target: signal Keeping : 4458 of 5828 = 0.7649279341111874

Run: 171 for target: signal Keeping : 1495 of 1956 = 0.7643149284253579

Run: 173 for target: signal Keeping : 1401 of 1820 = 0.7697802197802198

Run: 175 for target: signal Keeping : 1566 of 2015 = 0.7771712158808933

Run: 177 for target: signal Keeping : 1561 of 1955 = 0.7984654731457801

Run: 179 for target: signal Keeping : 1105 of 1671 = 0.6612806702573309

Run: 181 for target: signal Keeping : 2425 of 3426 = 0.7078225335668418

Run: 183 for target: signal Keeping : 2543 of 3550 = 0.716338028169014

Run: 185 for target: signal Keeping : 3033 of 3856 = 0.7865663900414938

Run: 187 for target: signal Keeping : 2712 of 3616 = 0.75

So definitely worse in that aspect.

Let's try cross entropy again, but with L1 or L2 regularization.

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --ε 0.95 \
    --modelOutpath ~/org/resources/nn_devel_mixing/19_03_23_l2_regularization/trained_model_incl_totalCharge_l2_regularization.pt

First attempt with:

        SGDOptions.init(0.005).momentum(0.2).weight_decay(0.01)

does not really converge. I guess that is too large.. :) Trying again with 0.001. This seems to work better. Oh, it broke between epoch 10000 and 15000, but better again at 20000 (but worse than before). It stayed on a plateau above the previous afterwards until the end. Also the distributions of the outputs are quite different now.

./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 ~/CastData/data/CalibrationRuns2018_Reco.h5 --back ~/CastData/data/DataRuns2017_Reco.h5 --back ~/CastData/data/DataRuns2018_Reco.h5 --ε 0.95 --modelOutpath ~/org/resources/nn_devel_mixing/19_03_23_l2_regularization/trained_model_incl_totalCharge_l2_regularizationcheckpoint_epoch_100000_loss_0.0261_acc_0.9963.pt --predict --plotPath ~/org/Figs/statusAndProgress/neuralNetworks/development/l2_regularization
./train_ingrid ~/CastData/data/CalibrationRuns2017_Reco.h5 ~/CastData/data/CalibrationRuns2018_Reco.h5 --back ~/CastData/data/DataRuns2017_Reco.h5 --back ~/CastData/data/DataRuns2018_Reco.h5 --ε 0.95 --modelOutpath ~/org/resources/nn_devel_mixing/19_03_23_l2_regularization/trained_model_incl_totalCharge_l2_regularizationcheckpoint_epoch_10000_loss_0.0156_acc_0.9964.pt --predict --plotPath ~/org/Figs/statusAndProgress/neuralNetworks/development/l2_regularization

./resources/nn_devel_mixing/19_03_23_l2_regularization/ ./Figs/statusAndProgress/neuralNetworks/development/l2_regularization/ Looking at the prediction of the final checkpoint (*_final_checkpoint.pdf) we see that we still have the same kind of shift in the data. However, after epoch 10000 () we see a much clearer overlap between the two (but likely also more background?). Still interesting, maybe L2 reg is useful if optimized to a good parameter. Let's look at the effective efficiencies of this particular checkpoint and comparing with the very last one: First the last:

./effective_eff_55fe \
    -f ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/19_03_23_l2_regularization/trained_model_incl_totalCharge_l2_regularizationcheckpoint_epoch_100000_loss_0.0261_acc_0.9963.pt \
    --ε 0.95

Run: 118 for target: signal Keeping : 1357 of 1651 = 0.8219261053906723

Run: 120 for target: signal Keeping : 2828 of 3413 = 0.8285965426311164

Run: 122 for target: signal Keeping : 4668 of 5640 = 0.8276595744680851

Run: 126 for target: signal Keeping : 2097 of 2596 = 0.8077812018489985

Run: 128 for target: signal Keeping : 6418 of 7899 = 0.8125079123939739

Run: 145 for target: signal Keeping : 2960 of 3646 = 0.811848601206802

Run: 147 for target: signal Keeping : 1731 of 2107 = 0.8215472235405791

Run: 149 for target: signal Keeping : 1588 of 1936 = 0.8202479338842975

Run: 151 for target: signal Keeping : 1482 of 1839 = 0.8058727569331158

Run: 153 for target: signal Keeping : 1565 of 1908 = 0.820230607966457

Run: 155 for target: signal Keeping : 1434 of 1777 = 0.806978052898143

Run: 157 for target: signal Keeping : 1457 of 1817 = 0.8018712162905889

Run: 159 for target: signal Keeping : 2914 of 3634 = 0.8018712162905889

Run: 161 for target: signal Keeping : 2929 of 3632 = 0.8064427312775331

Run: 163 for target: signal Keeping : 1474 of 1841 = 0.8006518196632265

Run: 165 for target: signal Keeping : 3134 of 3881 = 0.8075238340633857

Run: 167 for target: signal Keeping : 1609 of 2008 = 0.8012948207171314

Run: 169 for target: signal Keeping : 4738 of 5828 = 0.8129718599862732

Run: 171 for target: signal Keeping : 1591 of 1956 = 0.8133946830265849

Run: 173 for target: signal Keeping : 1465 of 1820 = 0.804945054945055

Run: 175 for target: signal Keeping : 1650 of 2015 = 0.8188585607940446

Run: 177 for target: signal Keeping : 1576 of 1955 = 0.8061381074168797

Run: 179 for target: signal Keeping : 1339 of 1671 = 0.8013165769000599

Run: 181 for target: signal Keeping : 2740 of 3426 = 0.7997664915353182

Run: 183 for target: signal Keeping : 2856 of 3550 = 0.8045070422535211

Run: 185 for target: signal Keeping : 3146 of 3856 = 0.8158713692946058

Run: 187 for target: signal Keeping : 2962 of 3616 = 0.8191371681415929

Once again in the ballpark of 80% while at 95% for CDL. And for epoch 10,000?

./effective_eff_55fe \
    -f ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/19_03_23_l2_regularization/trained_model_incl_totalCharge_l2_regularizationcheckpoint_epoch_10000_loss_0.0156_acc_0.9964.pt \
    --ε 0.95

Run: 118 for target: signal Keeping : 1357 of 1651 = 0.8219261053906723

Run: 120 for target: signal Keeping : 2790 of 3413 = 0.8174626428362145

Run: 122 for target: signal Keeping : 4626 of 5640 = 0.8202127659574469

Run: 126 for target: signal Keeping : 2092 of 2596 = 0.8058551617873652

Run: 128 for target: signal Keeping : 6377 of 7899 = 0.8073173819470819

Run: 145 for target: signal Keeping : 2974 of 3646 = 0.8156884256719693

Run: 147 for target: signal Keeping : 1735 of 2107 = 0.8234456573327005

Run: 149 for target: signal Keeping : 1606 of 1936 = 0.8295454545454546

Run: 151 for target: signal Keeping : 1485 of 1839 = 0.8075040783034257

Run: 153 for target: signal Keeping : 1575 of 1908 = 0.8254716981132075

Run: 155 for target: signal Keeping : 1444 of 1777 = 0.8126055149127743

Run: 157 for target: signal Keeping : 1478 of 1817 = 0.8134287286736379

Run: 159 for target: signal Keeping : 2932 of 3634 = 0.8068244358833242

Run: 161 for target: signal Keeping : 2942 of 3632 = 0.8100220264317181

Run: 163 for target: signal Keeping : 1484 of 1841 = 0.8060836501901141

Run: 165 for target: signal Keeping : 3134 of 3881 = 0.8075238340633857

Run: 167 for target: signal Keeping : 1612 of 2008 = 0.8027888446215139

Run: 169 for target: signal Keeping : 4700 of 5828 = 0.8064516129032258

Run: 171 for target: signal Keeping : 1582 of 1956 = 0.8087934560327198

Run: 173 for target: signal Keeping : 1469 of 1820 = 0.8071428571428572

Run: 175 for target: signal Keeping : 1630 of 2015 = 0.8089330024813896

Run: 177 for target: signal Keeping : 1571 of 1955 = 0.8035805626598466

Run: 179 for target: signal Keeping : 1366 of 1671 = 0.817474566128067

Run: 181 for target: signal Keeping : 2734 of 3426 = 0.7980151780502043

Run: 183 for target: signal Keeping : 2858 of 3550 = 0.8050704225352112

Run: 185 for target: signal Keeping : 3122 of 3856 = 0.8096473029045643

Run: 187 for target: signal Keeping : 2937 of 3616 = 0.8122234513274337

Interesting! Despite the much nicer overlap in the prediction at this checkpoint, the end result is not that different. Not quite sure what to make of that.

Next we try Adam, starting with this:

      var optimizer = Adam.init(
        model.parameters(),
        AdamOptions.init(0.005)
      )

./resources/nn_devel_mixing/19_03_23_adam_optim/ ./Figs/statusAndProgress/neuralNetworks/development/adam_optim/

-> Outputs are very funny. Extremely wide, need --clampOutput O(10000) or more. CDL and 55Fe are quite separated though! Enough for today.

• [X] Try L1 and L2 regularization of the network (weight decay parameter)
• [ ] Try L1 regularization
• [ ] Try Adam optimizer
• [ ] Try L2 with a value slightly larger and slightly smaller than 0.001

Continuing from yesterday, but before we do that, we need to generate the new expected limits table using the script in StatusAndProgress.org sec. [BROKEN LINK: sec:limit:expected_limits_different_setups_test].

• [X] Generate limits table
• [ ] Regenerate all limits once more to have them with the correct eccentricity cut off value in the files -> Should be done, but not priority right now. Our band aid fix relying on the filename is fine for now.
• [ ] continue NN training / investigation
• [ ] Update systematics due to determineEffectiveEfficiency using fixed code (correct energies & data frames) in thesis
  • [ ] fix that same code for fitByRun

From yesterday open TODOs:

• [ ] Regenerate all limits once more to have them with the correct eccentricity cut off value in the files -> Should be done, but not priority right now. Our band aid fix relying on the filename is fine for now.
• [ ] continue NN training / investigation
• [ ] Update systematics due to determineEffectiveEfficiency using fixed code (correct energies & data frames) in thesis
  • [ ] fix that same code for fitByRun

Additional:

• [X] look at prediction of our best trained network (and maybe the lnL variable one?) for all the different CDL datasets separately. Maybe a ridgeline plot of the different "sets", i.e. background, 55Fe photo, 55Fe escape, CDL sets
  • [ ] Do the same thing with the Run-2 and Run-3 calibration / background data split?
  • [ ] Do the same thing, but using the likelihood distributions for each instead of the NN predictions!
• [ ] Investigate whether effective efficiency (from tool) is correlated to mean gas gain of each calibration run. create a plot of the effective efficiency vs the mean gas gain of each run, per photo & escape type -> If this is strongly correlated it means we understand where the fluctuation comes from! If true, then can look at CDL data as well and check if this explains the variation.

Continue with the jobs from yesterday:

Additional:

• [X] look at prediction of our best trained network (and maybe the lnL variable one?) for all the different CDL datasets separately. Maybe a ridgeline plot of the different "sets", i.e. background, 55Fe photo, 55Fe escape, CDL sets
  • [ ] Do the same thing with the Run-2 and Run-3 calibration / background data split?
  • [X] Do the same thing, but using the likelihood distributions for each instead of the NN predictions!
• [X] Investigate whether effective efficiency (from tool) is correlated to mean gas gain of each calibration run. create a plot of the effective efficiency vs the mean gas gain of each run, per photo & escape type -> If this is strongly correlated it means we understand where the fluctuation comes from! If true, then can look at CDL data as well and check if this explains the variation.

That is: implement the lnL variant into the 'prediction' ridge line plots. And potentially look at the Run-2 vs Run-3 predictions.

• [ ] Look at the background rate of the 90% lnL cut variant. How much background do we have in that case? How does it compare to 99% accuracy MLP prediction?
• [ ] maybe try even larger MLP?

As a bonus:

• [ ] look at the hidden_2500neuron network for the background rate
• [ ] try to use the neural network for a limit calculation it its "natural" prediction. i.e close to 99% accuracy! That should give us quite the amazing signal (but of course decent background!). Still, as alternative combined with line veto and/or FADC could be very competitive!
• [X] Make notes about ROC curve plots
• [X] Next up: -> Look at effective efficiency again and how it varies -> Implement CAST ⇔ CDL transformation for cut values

From yesterday:

• [ ] Look at the background rate of the 90% lnL cut variant. How much background do we have in that case? How does it compare to 99% accuracy MLP prediction?
• [ ] maybe try even larger MLP?

As a bonus:

• [ ] look at the hidden_2500neuron network for the background rate
• [ ] try to use the neural network for a limit calculation it its "natural" prediction. i.e close to 99% accuracy! That should give us quite the amazing signal (but of course decent background!). Still, as alternative combined with line veto and/or FADC could be very competitive!

And in addition:

• [ ] Implement a fit that takes the effective efficiency and gas gain correlation into account and try it to correct for the efficiencies at CAST!
• [ ] Look at how the distributions change between different CDL runs with different gas gains.

First look into the energyFromCharge for the CDL data and see if it changes the cut values.

Important thought: ~3 keV escape events are not equivalent to 3 keV X-rays! escape events are effectively 5.9 keV photons that only deposit 3 keV! Real 3 keV X-rays have a much longer absorption length. That explains why 55Fe 3 keV data is shifted to a lower cut value, but CDL 3 keV data to a higher cut value when compared to the 5.9 keV data in each!

So our prediction of the cut value for the escape events via the slope of the 5.9 keV data is therefore too large, because the real events look "less" like X-rays to the network.

Generate two sets of fake events:

• [ ] Events of the same energy, but at an effectively different diffusion length, by taking transverse diffusion and the distance and 'pulling' all electrons of an event towards the center of the cluster -> e.g. generate 'real' 3 keV events from the escape peak 3 keV events
• [ ] Events of an artificially lower Timepix threshold of same energy. Look at how many electrons calibration set 1 has compared to 2. Then throw away this many electrons, biased by those of the lowest charges (or inversely fix a threshold in electrons and remove all pixels below and see where we end up). Problem is that total number of recorded electrons (i.e. ToT value) itself also depends on the threshold.
• [ ] (potential) a third set could just be looking at fake lower energy events generated in the way that we already do it in the lnL effective efficiency code!

These can then all be used to evaluate the MLP with.

• [ ] Investigate fake events!

Continuing from yesterday… Fake events and other stuff..

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/21_03_23_hidden_2500neurons/trained_model_hidden_2500neurons.pt \
    --ε 0.95 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit --plotDatasets

• [ ] Why is skewness of all CAST data (incl. fake data!) slightly positive? It should be centered around 0, no?
• [ ] rmsTransverse, length and width are as expected: fake data and escape peak data is larger than real 3 keV data! Due to different absorption lengths

And more continue, continue, continue…!

• [X] First start with reordering the ridgeline plots according to energy
• [ ] Then implement two different other fake data generation

Generation of data at different diffusion: Change the effective diffusion of the drawn event.

1. get the diffusion (transverse) coefficient σT = getDiffusion()

2. using the existing energy and target energies compute the required distance we 'move' the cluster from and to. That is: Assuming we have a diffusion equivalent to 3cm (conversion at cathode) and want to 'move' that to a diffusion of 2 cm (conversion 1cm away from cathode). We can compute the transverse diffusion by σ_T * √x cm (x ∈ [2, 3]) Each of the resulting numbers is the standard deviation of a normal distribution around the center position of the cluster!

   • [X] Verify how this relates to (see below) At 3cm the standard deviation is σ = √(6 D t) (3 dim)

   With the distributions we expect we now have a few options to generate new events

   1. simplest and deterministic: push all electrons to the equivalent value of the PDF (longer distance: shallower PDF. Find P(x_i) = P'(x_i') and move each x_i to x_i'.
   2. draw from the P' distribution for each pixel. The resulting x' is the location in distance from existing cluster center to place the pixel at.
   3. We could maybe somehow generate a 'local' PDF for each pixel (based on how far away each already is) and drawing from that. So a mix of 1 and 2?

   For now let's go with 2. Simpler to implement as we don't need to find an equivalent point on the PDF (which would be doable using lowerBound)

3. define the a gaussian with mean of resulting diffusion around that distance (what sigma does it have?) -> Or: define gaussian of the transverse diffusion coefficient and simply multiply!
4. for each pixel sample it and move each pixel the resulting distance towards the center / away from the center (depending)

And continue working on the fake event generation…!

• [X] Adjusting the diffusion down to low values (e.g. 400) does not move the fractionInTransverseRms to lower values!
• [X] Implement logL and tracking support into runAnalysisChain to really make it do (almost) everything

• [ ] Rerun the tracking log info!
• [ ] re evaluate the fake datasets & efficiencies in general using the non skewed data!

• [X] Understand weird filtering leaving some unwanted events in there -> Ohhh, the reason is that we filter on event numbers only, and not the actual properties of the clusters! Because one event can have multiple clusters. One of them will be passing, but the other likely not.

Ok, finally done with all the hick hack of coming up with working fake event generation etc etc.

As it seems right now:

Fake 5.9 keV data using correct diffusion gives 'identical' results in terms of efficiency than real CAST data. I.e. using the gain fit is correct and useful.

For the 3.0 keV case it is not as simple. The fake data is considered 'more' X-ray like (larger cut values), but quite clearly they don't fit onto the same slope as the 5.9 keV data!

What therefore might be a reasonable option:

• Generate X-ray data for all lines below 5.9 keV
• Use the generated 'runs' to fit the gas gain curve for each dataset
• Use that gas gain curve for each energy range. Lower energy ranges are likely to have somewhat shallower gas gain dependencies? Or maybe it's events with shorter absorption length. We'll have to test.
• [ ] Include one other energy, e.g. 930 eV, due to the very low absorption length. See how that behaves.

• [ ] Generalize the effective efficiency code to also include other CDL lines

• [ ] Rerun effective eff code and place efficiencies and plots somewhere
• [ ] Rerun old skewness model with --predict option in train_ingrid
• [X]

  Train a new MLP with same parameters as 2500 hidden neuron model, but using the corrected skewness data!

  ./train_ingrid \
      ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --back ~/CastData/data/DataRuns2017_Reco.h5 \
      --back ~/CastData/data/DataRuns2018_Reco.h5 \
      --modelOutpath ~/org/resources/nn_devel_mixing/28_03_23_hidden_2500_fixed_skew/trained_model_hidden_2500_fixed_skew.pt \
      --plotPath ~/org/Figs/statusAndProgress/neuralNetworks/development/hidden_2500neurons_fixed_skew/ \
      --numHidden 2500
  

  -> This was done, but we haven't looked at it yet!

• [ ] In https://www.youtube.com/watch?v=kCc8FmEb1nY Andrej mentions 3e-4 is a good learning rate for AdamW. We've only tried Adam. Let's try AdamW as well.

We'll continue from yesterday:

• [X] Rerun the tracking log info!
• [ ] re evaluate the fake datasets & efficiencies in general using the non skewed data!
• [X] Include one other energy, e.g. 930 eV, due to the very low absorption length. See how that behaves.
• [ ] Generalize the effective efficiency code to also include other CDL lines
• [ ] Rerun effective eff code and place efficiencies and plots somewhere
• [ ] Rerun old skewness model with --predict option in train_ingrid
• [X] In https://www.youtube.com/watch?v=kCc8FmEb1nY Andrej mentions 3e-4 is a good learning rate for AdamW. We've only tried Adam. Let's try AdamW as well.

New TODOs for today:

• [ ] Generate a plot of the cut values that contains all types of data in one with color being the data type

Can we make a fit to the rms transverse data of each 55Fe run, then use that upper limit of the transverse RMS to compute the diffusion and finally determine the cut value based on fake data with that diffusion and gas gain?

Goal: Determine NN cut value to use for a given cluster to achieve an efficiency of a desired value.

I have:

• CDL data that can be used to determine cut values for different energy. They have different gas gains and diffusion coefficients.
• Fake data for arbitrary energies, absorption lengths and diffusion coefficients
• real 5.9 keV data at different gas gains and diffusion coefficients.

What do I need: A relationship that maps a cut value from a CDL energy range to a cluster of different diffusion and gas gain.

How do I understand the dependence of diffusion and gas gain on NN cut value? Use upper percentile (e.g. 95) of rmsTransverse as an easy to use proxy for the diffusion in a dataset. Compute that value for each CDL run. Do same for every 55Fe run. Plot gas gain vs rmsTransverse for all this data.

Thoughts on diffusion & gas gain:

• [ ] INSERT PLOT OF RMS T VS GAS GAIN & CUT VAL

For hidden 2500 with skewed fixed:

• ./Figs/statusAndProgress/neuralNetworks/development/hidden_2500neurons_fixed_skew/rmsTransverse_vs_NN_cutVal.pdf
• ./Figs/statusAndProgress/neuralNetworks/development/hidden_2500neurons_fixed_skew/rmsTransverse_vs_gas_gain.pdf

For AdamW@3e-4 lr:

• ./../../../org/Figs/statusAndProgress/neuralNetworks/development/adamW_2500_3e-4/rmsTransverse_vs_NN_cutVal.pdf
• ./../../../org/Figs/statusAndProgress/neuralNetworks/development/adamW_2500_3e-4/rmsTransverse_vs_gas_gain.pdf

Higher gas gains are associated with longer transverse RMS values and lower NN cut values. Higher gas gains are associated with lower temperatures. Lower temperatures mean higher densities at the same pressure. Higher densities imply shorter absorption lengths. Shorter absorption lengths imply longer drift distances. Longer drift distances imply larger diffusion. Larger diffusion implies larger rms transverse.

So it may not actually be that the gas diffusion changes significantly (or only?), but that the change in density implies a change in average diffusion value.

Keep in mind that the crosses for 3 keV are the escape photons and not real 3 keV data!

The two are certainly related though.

Maybe rms transverse of CDL data is different due to different pressure?

• [ ] Make the same plot but instead of 3 keV escape photons generate 3 keV events at different diffusions

• [ ] Verify plotDatasets distributions of all fake data events! -> make this plot comparing CDL & fake data of same 'kind' by each set

Interesting observation:

Today we were trying to debug the reason why our sampled data seems to have a mean in the centerX and in particular centerY position that is not centered at ~7mm (128). Turns out, our original real data from which we sample has the same bias in the centerY position. Why is that? Don't we apply the same cuts in fake gen data reading as in the effective efficiency code?

-> The only difference between the two sets of cuts in the data reading is that in the fake data reading we do not apply an energy cut to the data. We only apply the xray cleaning cuts!

Finally implemented NN veto with run local cut values based on fake data in likelihood

• [ ]

  INSERT FIGURES OF EFFECTIVE EFFICIENCY USING FAKE CUT VALUES

  likelihood \
      -f ~/CastData/data/DataRuns2017_Reco.h5 \
      --h5out /tmp/testing/run2_only_mlp_local_0.95.h5 \
      --region crGold \
      --cdlYear 2018 \
      --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --mlp ~/org/resources/nn_devel_mixing/30_03_23_SGD_2500_3e-4_largerRmsT/trained_model_SGD_2500_3e-4_largerRmsT.pt \
      --nnSignalEff 0.95 \
      --nnCutKind runLocal \
      --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      --readonly
  

Ah, but we cannot directly sample fake events in our current approach for the background datasets, because we still rely on real events as a starting point for the sampling.

Therefore: develop sampling of full events using gas gain and polya sampling.

1. check what info is stored in H5 files regarding polya and gas gain -> nothing useful about threshold
2. SCurve should in theory tell us something about threshold
   • [ ] investigate if all polyas in a single file (Run-2 or Run-3) have the same threshold if determined from data, e.g. using quantile 1
3. sample from existing fixed parameter polya and check it looks reasonable
4. open H5 file run by run, look at minima and quantile 1 data for each gas gain slice. -> The minima is fixed for each set of chip calibrations! Run-2: Run-2: minimum: 1027.450870326596, quantile 1: 1247.965130336289 or 1282.375318596229 Run-3: minimum: 893.4812944899318, quantile 1: 1014.494678806292 or 1031.681385601604 -> this is useful! Means we can just sample with fixed cutoffs for each dataset and don't need raw charge data for anything for a run, but only gas gain slice fit parameters! -> Raises question though: How does then this relate to # of hits in photo peak? Lower gain, closer to cutoff, thus less pixels, but opposite? How do we get to too many pixels?
5. plot photo peak hit positions against gas gain
6. sample from a poyla using parameters as read from gas gain slices, plot against polya dataset -> Looks reasonable. Now need to include the cutoff. -> Looks even better. Note that the real data looks ugly due to the equal bin widths that are not realistic.

Ahh! Idea: we can reuse the gas gain vs energy calibration factor fit! We get the gas gain of a run for which to generate fake data. That gas gain is fed into that fit. The result is a calibration factor that tells us the charge corresponding to a given energy (or its inverse). Place our target energy into the function to get the desired charge. Then sample from a normal distribution around the target charge as the target for each event. Generate pixels until total charge close to target. -> Until total charge close or alternatively: Given gas gain calculate number based on target charge: e.g. 600,000 target and gain of 3,000) then target hits is target / per hit = 200 in this case. In this case we can get less due to threshold effects, but not more. So: better to do the former? Draw from target charge distribution and then accumulate pixels until the total drawn is matched?

Question: In the real data, are the total charge and the number of hits strongly correlated? It's important to understand how to go from a deposited energy to a number of electrons. There are multiple reasons to lose electrons and charge from a fixed input energy:

• amplifications below the threshold of a pixel
• higher / lower than average ionization events than W_i = 26 eV

We can model the former, but the latter is tricky.

Meeting with Klaus today:

Main take away: Train an MLP using purely generated data (for X-rays) with:

• an extra input neuron that corresponds to the transverse diffusion of the 'dataset' each.
• different MLPs, one for each input of the diffusion parameter

The former seems better to me. Just generate data with a uniform distribution in some range of diffusion parameters. Then each input actually has different values. When actually applying the network then we only have few distinct values (one for each run) of course. But that's fine.

Start a training run with the σT dataset!

With 100,000 generated events for each calibration file (why?). Still using the fixed cutoff of Run-3! This also uses a uniform distribution between G = 2400 .. 4500 σT = 550 .. 700 θ = 2.1 .. 2.4

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_data_diffusion/trained_mlp_sgd_sim_data_diffusion.pt \
    --plotPath ~/Sync/sgd_sim_data_diffusion/ \
    --numHidden 2500 \
    --datasets eccentricity \
    --datasets skewnessLongitudinal \
    --datasets skewnessTransverse \
    --datasets kurtosisLongitudinal \
    --datasets kurtosisTransverse \
    --datasets length \
    --datasets width \
    --datasets rmsLongitudinal \
    --datasets rmsTransverse \
    --datasets lengthDivRmsTrans \
    --datasets rotationAngle \
    --datasets fractionInTransverseRms \
    --datasets σT \
    --learningRate 3e-4 \
    --simulatedData 

At first glance it seems like the cut values determined from generating more fake data of the σT and gain of the real runs yields effective efficiencies that are all too high (95%, only some CDL runs approaching 80% target).

Now with normal distribution in G and σT:

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_gauss_diffusion/trained_mlp_sgd_sim_gauss_diffusion.pt \
    --plotPath ~/Sync/sgd_sim_gauss_diffusion/ \
    --numHidden 2500 \
    --datasets eccentricity \
    --datasets skewnessLongitudinal \
    --datasets skewnessTransverse \
    --datasets kurtosisLongitudinal \
    --datasets kurtosisTransverse \
    --datasets length \
    --datasets width \
    --datasets rmsLongitudinal \
    --datasets rmsTransverse \
    --datasets lengthDivRmsTrans \
    --datasets rotationAngle \
    --datasets fractionInTransverseRms \
    --datasets σT \
    --learningRate 3e-4 \
    --simulatedData 

-> Continue training this if this is useful. Still a lot of downward trend in loss visible!

In addition it may be a good idea to also try it with the gain as another input?

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_diffusion_gain/trained_mlp_sgd_sim_diffusion_gain.pt \
    --plotPath ~/Sync/sgd_sim_diffusion_gain/ \
    --numHidden 2500 \
    --datasets eccentricity \
    --datasets skewnessLongitudinal \
    --datasets skewnessTransverse \
    --datasets kurtosisLongitudinal \
    --datasets kurtosisTransverse \
    --datasets length \
    --datasets width \
    --datasets rmsLongitudinal \
    --datasets rmsTransverse \
    --datasets lengthDivRmsTrans \
    --datasets rotationAngle \
    --datasets fractionInTransverseRms \
    --datasets σT \
    --datasets gasGain \
    --learningRate 3e-4 \
    --simulatedData 

• [ ] Understand effective efficiency with gas gain parameter of FakeDesc. E.g. rerun the gauss diffusion network again after changing code etc
• [ ]

  verify the number of neighboring pixels that are actually active anywhere! -> We'll write a short script that extracts that information from a given file.

  • [ ] Correlate with gas gain and extracted diffusion!
  • [ ] Compare with generated fake data without artificial neighbor activation and with!
  • [ ] How do neighbors relate in their charge?

  In ./../CastData/ExternCode/TimepixAnalysis/Tools/countNeighborPixels:

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2017_Reco.h5

Average neighbors in run 110 = 55.09144736842105 Average neighbors in run 175 = 57.15461200585651 Average neighbors in run 122 = 47.4565240584658 Average neighbors in run 126 = 52.52208235545125 Average neighbors in run 183 = 59.16529930112428 Average neighbors in run 161 = 66.33755847291384 Average neighbors in run 116 = 51.58539765319426 Average neighbors in run 155 = 68.88583638583638 Average neighbors in run 173 = 57.67820710973725 Average neighbors in run 151 = 66.58298001211386 Average neighbors in run 153 = 68.5710128055879 Average neighbors in run 108 = 55.27580484566877 Average neighbors in run 93 = 55.01024327784891 Average neighbors in run 147 = 61.76867469879518 Average neighbors in run 179 = 63.32529743268628 Average neighbors in run 159 = 66.94385479157053 Average neighbors in run 163 = 63.48872858431019 Average neighbors in run 118 = 52.85084521047398 Average neighbors in run 102 = 54.66866666666667 Average neighbors in run 177 = 56.74977000919963 Average neighbors in run 181 = 60.51956253850894 Average neighbors in run 165 = 59.47996965098634 Average neighbors in run 167 = 61.3946587537092 Average neighbors in run 185 = 57.19969558599696 Average neighbors in run 149 = 67.73385167464114 Average neighbors in run 157 = 69.89698937426211 Average neighbors in run 187 = 58.77159520807061 Average neighbors in run 171 = 59.13408330799635 Average neighbors in run 128 = 52.05005608524958 Average neighbors in run 169 = 58.22026431718061 Average neighbors in run 145 = 61.53792611101196 Average neighbors in run 83 = 51.74536148432502 Average neighbors in run 88 = 46.48597521200261 Average neighbors in run 120 = 49.65804645033767 Average neighbors in run 96 = 51.49257633765991

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5

Average neighbors in run 304 = 83.96286764705883 Average neighbors in run 286 = 83.78564713256033 Average neighbors in run 294 = 84.31031159653068 Average neighbors in run 277 = 89.82936363636364 Average neighbors in run 241 = 81.1551888289432 Average neighbors in run 284 = 89.4854306756324 Average neighbors in run 260 = 85.9630966706779 Average neighbors in run 255 = 88.56640625 Average neighbors in run 292 = 85.56208945886769 Average neighbors in run 288 = 80.55608820709492 Average neighbors in run 247 = 77.27826358525921 Average neighbors in run 257 = 88.13119266055045 Average neighbors in run 239 = 73.20087064676616 Average neighbors in run 302 = 86.77864992150707 Average neighbors in run 249 = 79.53715365239294 Average neighbors in run 271 = 79.5505486808312 Average neighbors in run 300 = 90.83248730964468 Average neighbors in run 296 = 85.71424050632912 Average neighbors in run 243 = 80.90277344967279 Average neighbors in run 264 = 83.18354637823664 Average neighbors in run 280 = 88.5948709880428 Average neighbors in run 253 = 85.09293373659609 Average neighbors in run 251 = 77.57475909232204 Average neighbors in run 262 = 82.7622203811102 Average neighbors in run 290 = 79.6481004507405 Average neighbors in run 275 = 89.19542053956019 Average neighbors in run 269 = 82.60137931034483 Average neighbors in run 266 = 83.25234248788368 Average neighbors in run 273 = 81.36689741976086 Average neighbors in run 245 = 78.7516608668143 Average neighbors in run 259 = 86.49635416666666 Average neighbors in run 282 = 90.84639199809479

And using fake data, for the case of no simulated neighbors: Run-2

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2017_Reco.h5 --fake

Average neighbors in run 161 = 38.16306522609044 Average neighbors in run 128 = 29.87515006002401 Average neighbors in run 183 = 34.27851140456183 Average neighbors in run 185 = 33.6796718687475 Average neighbors in run 88 = 28.7222 Average neighbors in run 179 = 34.7138 Average neighbors in run 187 = 32.83693477390956 Average neighbors in run 155 = 39.0502 Average neighbors in run 163 = 35.952 Average neighbors in run 118 = 30.9526 Average neighbors in run 171 = 32.7596 Average neighbors in run 126 = 29.8646 Average neighbors in run 151 = 37.3554 Average neighbors in run 169 = 32.7694 Average neighbors in run 120 = 28.37575030012005 Average neighbors in run 102 = 30.4712 Average neighbors in run 159 = 37.93177270908363 Average neighbors in run 153 = 38.0882 Average neighbors in run 157 = 40.6742 Average neighbors in run 167 = 34.4742 Average neighbors in run 96 = 29.2946 Average neighbors in run 175 = 33.7786 Average neighbors in run 177 = 33.1542 Average neighbors in run 93 = 31.7502 Average neighbors in run 116 = 29.2908 Average neighbors in run 83 = 29.9446 Average neighbors in run 145 = 34.6312 Average neighbors in run 147 = 34.9114 Average neighbors in run 108 = 31.13945578231293 Average neighbors in run 122 = 27.9868 Average neighbors in run 181 = 34.3359343737495 Average neighbors in run 165 = 32.90716286514606 Average neighbors in run 173 = 32.5428 Average neighbors in run 149 = 38.0322 Average neighbors in run 110 = 31.2116

Run-3

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --fake

Average neighbors in run 284 = 41.75550220088035 Average neighbors in run 259 = 39.1072 Average neighbors in run 292 = 39.506 Average neighbors in run 286 = 39.2786 Average neighbors in run 239 = 37.4598 Average neighbors in run 288 = 37.3952 Average neighbors in run 251 = 36.9704 Average neighbors in run 255 = 40.5622 Average neighbors in run 262 = 39.1088 Average neighbors in run 260 = 40.5608 Average neighbors in run 294 = 38.376 Average neighbors in run 280 = 40.51660664265706 Average neighbors in run 271 = 37.5112 Average neighbors in run 296 = 39.0084 Average neighbors in run 275 = 40.3642 Average neighbors in run 269 = 38.63645458183273 Average neighbors in run 302 = 38.698 Average neighbors in run 304 = 38.2386 Average neighbors in run 266 = 37.8624 Average neighbors in run 243 = 37.6958 Average neighbors in run 264 = 39.4124 Average neighbors in run 257 = 40.6936 Average neighbors in run 282 = 41.2402 Average neighbors in run 290 = 37.7166 Average neighbors in run 277 = 40.49369369369369 Average neighbors in run 253 = 38.238 Average neighbors in run 249 = 38.3926 Average neighbors in run 273 = 38.073 Average neighbors in run 247 = 37.3596 Average neighbors in run 245 = 36.2114 Average neighbors in run 241 = 37.1948 Average neighbors in run 300 = 39.7052

So the number of neighbors:

Activating neighbor sharing we can get those numbers up, but that's for later.

Next: Plot charge of pixels with neighbors against something.

Note: For run 288 the following:

    if charge > 3500.0: # whatever
      # possibly activate a neighbor pixel!
      let activateNeighbor = rnd.rand(1.0) < 0.5
      if activateNeighbor:
        let neighbor = rand(3) # [up, down, right, left]
        let totNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
        case neighbor

yields very good agreement already, with the exception of the hits data if the gain.G * 0.85 is not used (but instead G itself).

In ./../CastData/ExternCode/TimepixAnalysis/Tools/countNeighborPixels we can now produce also histograms of the ToT values / charges recorded by the center chip (split by number of neighbors of a pixel). Three different versions, density of ToT, raw charge and density of charge.

The three are important because of the non-linearity of the ToT to charge conversion. The question is: which distribution should we really sample from to get the correct distribution of charges seen by the detector?

Running:

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2017_Reco.h5 --chargePlots

yields the plots in:

We can see that the density based version of the charge histograms looks very much not like a Polya. The histogram of ToT looks somewhat like it and the raw charge one looks closest I would say.

The other thing we see here is that the real data shows a shift to larger charges for the data with more neighbors. The effect is not extreme, but visible. Much more so in the Run-3 data than in the Run-2 data though, which matches our expectations (see table from yesterday).

For the fake data:

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake

yields: Figs/statusAndProgress/gasGainAndNeighbors/charges_neighbors_fake.pdf

At the very least the distributions currently generated do not match. Multiple reasons:

• our scaling of gas gain using G * 0.85 is bad (yes)
• currently we're sampling from a polya of the ToT values

We will now study the gas gains on a specific run, say 241. We'll try to imitate the look of the real 241 data in the fake data. Using:

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241

we'll make changes and try to get it to look better.

Starting point:

#let gInv = invert(gain.G * 0.85, calibInfo)
let gInv = invert(gain.G, calibInfo) # not using scaling
# ...
  let ToT = rnd.sample(psampler) # sampling from ToT
# and no neighbor activation logic

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_start.pdf

First step, what does our 0.85 scaling actually do?

let gInv = invert(gain.G * 0.85, calibInfo)
# ...
  let ToT = rnd.sample(psampler) # sampling from ToT
# and no neighbor activation logic

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_0.85_scaling.pdf

As expected, it makes the gains a bit smaller.

Notice how the ToT histogram of the real data has a much shorter tail than the fake data. At ToT = 150 it is effectively 0, but fake data still has good contribution there! Try scaling further down to see what that does.

First step, what does our 0.85 scaling actually do?

let gInv = invert(gain.G * 0.6, calibInfo)
# ...
  let ToT = rnd.sample(psampler) # sampling from ToT
# and no neighbor activation logic

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_0.6_scaling.pdf

The tail now looks more correct (almost zero at 150 ToT in fake), but:

• the peak is way too far left compared to real data
• [X]

  why does the real data have ToT values up to 0, but the generated data has a sharp cutoff at ~10 or so?

    charge 3516.613062845426 from ToT 64.00852228399629 for 1 is 893.4812944899318
    charge 1405.389895032497 from ToT 26.41913859339533 for 1 is 893.4812944899318
    charge 6083.063184722204 from ToT 88.29237432771937 for 1 is 893.4812944899318
  

  This matches my expectation, but not the data. Ah! It's because of our 1.15 scaling, no?

    charge 2295.67900197231 from ToT 47.40724170682395 for 1 is 893.4812944899318 for cutoff: 1027.503488663421 inverted cutoff: 9.0
    charge 4940.614026607494 from ToT 78.30781008751393 for 1 is 893.4812944899318 for cutoff: 1027.503488663421 inverted cutoff: 9.0
    charge 1318.261093845941 from ToT 23.29152025807923 for 1 is 893.4812944899318 for cutoff: 1027.503488663421 inverted cutoff: 9.0
    charge 2643.805171369868 from ToT 52.90409160403075 for 1 is 893.4812944899318 for cutoff: 1027.503488663421 inverted cutoff: 9.0
    charge 2102.98283203606 from ToT 43.92681742421028 for 1 is 893.4812944899318 for cutoff: 1027.503488663421 inverted cutoff: 9.0
  

  Exactly!

The ToT behavior makes me worried in one sense: I really feel like the issue is that certain pixels have different thresholds, which makes the distribution so ugly.

Next, go back to sampling from actual polya and see what that looks like (without any scaling of gas gain):

  let params = @[gain.N, gain.G, gain.theta]
  let psampler = initPolyaSampler(params, frm = 0.0, to = 20000.0) #invert(20000.0, calibInfo))
# ...
  let charge = rnd.sample(psampler) # sampling from ToT
  let ToT = invert(charge, calibInfo)

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_from_polya.pdf

This does look much more realistic! So sampling from the real polya seems more sensible after all I think. Maybe it's a bit too wide though?

Let's look at scaling the theta parameter down by 30% (50%)

  let params = @[gain.N, gain.G, gain.theta * 0.7] # 0.5
  let psampler = initPolyaSampler(params, frm = 0.0, to = 20000.0)
# ...
  let charge = rnd.sample(psampler) # sampling from ToT
  let ToT = invert(charge, calibInfo)

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_polya_theta_smaller30.pdf

In the 50% versin of the density charge plot a difference becomes visible. The behavior towards the peak on the left is different. The peak moves a bit to smaller values approaching the cutoff more. But the width towards the tail does not really change.

In case the individual pixel thresholds are very different we could maybe approximate that using an exponential distribution from which we draw, which implies 100% not activating a pixel at cutoff and exponentially more likely to activate approaching the peak.

We implemented this by using the function:

proc expShift(x: float, u, l, p: float): float =
  ## Shifted exponential distribution that satisfies:
  ## f(l) = p, f(u) = 1.0
  let λ = (u - l) / ln(p)
  result = exp(- (x - u) / λ)

# ...
  ## XXX: make dependent on Run-2 or Run-3 data!
  const cutoffs = @[1027.450870326596, # Run-2
                    893.4812944899318] # Run-3
  let cutoff = cutoffs[1] * 1.15

  let actSampler = (proc(rnd: var Rand, x: float): bool =
                      let activateThreshold = expShift(x, gain.G, cutoff, 1e-1)
                      result = x > cutoff and rnd.rand(1.0) < activateThreshold
                      echo "Charge ", x, " threshold: ", activateThreshold, " is ", result, " and cutoff ", cutoff
  )
    let activatePixel = actSampler(rnd, charge)
    if not activatePixel: #charge < cutoff:
      continue

This is way too extreme:

pdfunite *neighbor_run_charges_241* ~/org/Figs/statusAndProgress/gasGainAndNeighbors/run_241_comparison_polya_exp_activation.pdf

the latter with p = 0.3 instead of 0.1.

I think this is the wrong approach. It pushes us even stronger to lager values. In the real charge density plot we need to be lower rather than higher.

Not using an exponential cutoff, but modifying both the gas gain and theta parameters yields the best result. But it's still ugly. Fortunately our latest network does not rely on the total charge anymore.

Let's run the effective efficiency for the dataset plots of run 241 comparison:

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_diffusion_gain/trained_mlp_sgd_sim_diffusion_gain.pt \
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/sgd_sim_diffusion_gain/effective_eff \
    --run 241

Fake:

• too many hits
• energy mismatch
• eccentricity way too small!
• charge a bit too large

-> very bad match. Some of it is due to not using neighbors. But realistically at this gas gain we cannot justify more hits! So we need a higher gain to get less hits again.

Using

  let params = @[gain.N, gain.G * 0.75, gain.theta / 3.0]

yields a good agreement in the histograms. with showing similar issues as above.

But first let's try to look at neighbors again:

    if charge > 3500.0: # whatever
      # possibly activate a neighbor pixel!
      let activateNeighbor = rnd.rand(1.0) < 0.5
      if activateNeighbor:
        let neighbor = rand(3) # [right, left, up, down]
        let totNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
        case neighbor
        of 0: insert(xp + 1, yp,     totNeighbor)
        of 1: insert(xp - 1, yp,     totNeighbor)
        of 2: insert(xp,     yp + 1, totNeighbor)
        of 3: insert(xp,     yp - 1, totNeighbor)
        else: doAssert false
        totalCharge += calib(totNeighbor, calibInfo)
    insert(xp, yp, ToT)

same settings of gain etc as last above yields The issues are obvious. The neighbor histogram shows stark jumps (expected now that I think about it) and the eccentricity is still too low while energy a bit too high already.

Next: let's implement a strategy for smoother activation of neighbors. Let's go with a linear activation from 1000 to 10000 electrons, 0 to 1 chance and see what it looks like.

  let neighSampler = (proc(rnd: var Rand, x: float): bool =
                        let m = 1.0 / 9000.0
                        let activateThreshold = m * x - 1000 * m
                        result = rnd.rand(1.0) < activateThreshold
  )
# ...
      if neighSampler(rnd, charge): # charge > 3500.0: # whatever
      # possibly activate a neighbor pixel!
      #let activateNeighbor = rnd.rand(1.0) < 0.5
      let chargeNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
      let activateNeighbor = actSampler(rnd, chargeNeighbor)
      if true: # activateNeighbor:
        let neighbor = rand(3) # [right, left, up, down]
        #let chargeNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
        let totNeighbor = invert(chargeNeighbor, calibInfo)
        case neighbor
        of 0: insert(xp + 1, yp,     totNeighbor)
        of 1: insert(xp - 1, yp,     totNeighbor)
        of 2: insert(xp,     yp + 1, totNeighbor)
        of 3: insert(xp,     yp - 1, totNeighbor)
        else: doAssert false
        totalCharge += chargeNeighbor

This yields:

./countNeighborPixels -f ~/CastData/data/CalibrationRuns2018_Reco.h5 --chargePlots --fake --run 241
./effective_eff_55fe ~/CastData/data/CalibrationRuns2018_Reco.h5 --model ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_diffusion_gain/trained_mlp_sgd_sim_diffusion_gain.pt --ε 0.8 --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 --evaluateFit --plotDatasets --plotPath ~/Sync/sgd_sim_diffusion_gain/effective_eff --run 241

Looking better. The neighbor distribution is more sensible, in particular in the ToT plot we can see a similar drift for larger neighboring cases as in the real data. However, the energy is now even bigger and the number of hits way too large (peak 300 compared to 260).

The total charge and energy being too large is probably simply an issue with our gas gain. Given that we don't use the charge or energy (but the gain in one network!) it's not that big of an issue. But the number of hits better be believable.

Given that we still have probably less higher neighbor cases, too low eccentricity and too many hits, likely implies:

• add higher neighbor cases with lower chance
• scale down the target charge by our modification of the gain?

The former first implementation:

  let neighSampler = (proc(rnd: var Rand, x: float): int =
    ## Returns the number of neighbors to activate!
    let m = 1.0 / 9000.0
    let activateThreshold = m * x - 1000 * m
    let val = rnd.rand(1.0)
    if val * 4.0 < activateThreshold:
      result = 4
    elif val * 3.0 < activateThreshold:
      result = 3
    elif val * 2.0 < activateThreshold:
      result = 2
    elif val       < activateThreshold:
      result = 1
    else:
      result = 0
    #result = rnd.rand(1.0) < activateThreshold
  )
    # ...
    let numNeighbors = neighSampler(rnd, charge)
    if numNeighbors > 0: # charge > 3500.0: # whatever
      # possibly activate a neighbor pixel!
      #let activateNeighbor = rnd.rand(1.0) < 0.5
      var count = 0
      type Neighbor = enum Right, Left, Up, Down
      var seen: array[Neighbor, bool]
      while count < numNeighbors:
        let chargeNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
        #let activateNeighbor = actSampler(rnd, chargeNeighbor)
        let neighbor = block:
          var num = Neighbor(rnd.rand(3))
          while seen[num]:
            num = Neighbor(rnd.rand(3)) # [right, left, up, down]
          num
        seen[neighbor] = true
        #let chargeNeighbor = rnd.sample(psampler) / 2.0 # reduce amount
        let totNeighbor = invert(chargeNeighbor, calibInfo)
        case neighbor
        of Right: insert(xp + 1, yp,     totNeighbor)
        of Left:  insert(xp - 1, yp,     totNeighbor)
        of Up:    insert(xp,     yp + 1, totNeighbor)
        of Down:  insert(xp,     yp - 1, totNeighbor)
        totalCharge += chargeNeighbor
        insert(xp, yp, ToT)
        totalCharge += charge
        inc count

Yeah, a look at the count histogram shows this is clearly rubbish. At the very least we managed to produce something significantly more eccentric than the real data, yet has still significantly less hits! Sigh, bug in the insertion code…! The insert(xp, yp, ToT) shouldn't be in that part! This is looking half way reasonable for the neighbor histograms! Wow, this looks pretty convincing in the property histograms! Aside from having too many hits still, it looks almost perfect.

Let's be so crazy and run the effective efficiencies on all runs:

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_gauss_diffusion/trained_mlp_sgd_sim_gauss_diffusion.pt \
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/sgd_sim_gauss_diffusion/effective_eff
./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/12_04_23_sgd_sim_diffusion_gain/trained_mlp_sgd_sim_diffusion_gain.pt
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/sgd_sim_diffusion_gain/effective_eff

For the gauss diffusion network The results are not as good as I hoped. But there is a high chance the numbers are also bad because the training data was not the way it looks now? Didn't run with the network knowing the gain due to the problems of having the gain be equivalent between real and fake data.

Let's look at the worst performing run, 282 from this plot: In this case the eccentricity is still too small compared to the real case.

Let's first fix the number of hits we see.

  let gainToUse = gain.G * 0.75
  let calibFactor = linearFunc(@[calibInfo.bL, calibInfo.mL], gainToUse) * 1e-6

introduced this variable that's used everywhere instead of gain.G. There were still some places we only used gain.G! (Note: oh: we achieved an efficiency of 81% on this run with the gain diffusion network) Takeaway:

• energy way too low
• hits way too low (190 compared to 260)
• eccentricity slightly too large

(Note: weird: changing the gain the reconstruct fake event call to use the gainInfo.G * 0.75 we currently use, makes the effective efficiency be 18% only, but the energy histogram does not change?? Ahh, it's because of the gain being an input neuron! That of course means it breaks the expectation of the trained network. Implies same for network without gain should yield good result. -> Checked, it does. Bigger question: why does the energy reconstruction still yield such a low energy? Shouldn't that change if we change the apparent gain? -> OHH, it's because we use the hardcoded gain from the calibInfo argument in the computeEnergy function! Updating… Yup, energy is looking perfect now: but of course the effective efficiency is still completely off in this network!)

So, next target is getting the hits in order. The reason must be the fact that our neighbor tot/charge histograms have a longer tail than the real data. Hmm, but the tail looks almost identical now…

I tried different gain values after all again and played around with counting the added neighboring pixels in the total charge or not counting them. It seems to me like going back to the original gain is the better solution. We still produce too few hits (to be thought about more), but at least we natively get the correct energy & charge. The histograms of the gain curves also don't look so horrible with our theta/3 hack.

and the properties: Fraction in transverse RMS is slightly too low though. But the effective efficiencies look quite good:

So next points to do:

• [ ] Maybe number of neighbors should include diagonal elements after all? If UV photons are the cause the distance can be larger maybe. Currently our code says > 110 neighbors for the fake data but O(90) for the real data in Run-3!
• [ ] investigate more into the number of hits we see. Can we fix it?
• [ ] Look into how properties look now for CDL data. Maybe our slope etc of the neighbor logic needs to be adjusted for lower gains?
• [ ] train a new network that uses the correct neighboring charges as training data

First we changed the amount of charge for neighbors from div 2 to div 3:

        let chargeNeighbor = rnd.sample(psampler) / 3.0 # reduce amount

This already improved the numbers a bit, also lowering the number of neighbors to a mean of 100.

Then we lowered the gas gain from 1.0 to 0.9 but not for the target charge.

  let gainToUse = gain.G * 0.9

  let calibFactor = linearFunc(@[calibInfo.bL, calibInfo.mL], gain.G) * 1e-6

i.e. leaving the calibration factor unchanged. Why? I have no idea, but it gets the job done.

Running the effective efficiencies of all runs now: -> This makes the efficiencies worse! But that could be an effect of having a network trained on the wrong data.

So now let's try to train a new network using our now better fake data and trying again. This now uses all our new changes and:

• gain 0.9 for everything but target charge
• possibly up to 4 neighbors
• neighbors receive 1/3 of sampled charge
• linear chance to activate neighbors

Also: we increase the range of theta parameters in the generation of fake events from 2.1 .. 2.4 to 0.8 to 2.4.

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/16_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusion.pt \
    --plotPath ~/Sync/16_04_23_sgd_sim_diffusion_gain/ \
    --numHidden 2500 \
    --datasets eccentricity \
    --datasets skewnessLongitudinal \
    --datasets skewnessTransverse \
    --datasets kurtosisLongitudinal \
    --datasets kurtosisTransverse \
    --datasets length \
    --datasets width \
    --datasets rmsLongitudinal \
    --datasets rmsTransverse \
    --datasets lengthDivRmsTrans \
    --datasets rotationAngle \
    --datasets fractionInTransverseRms \
    --datasets σT \
    --learningRate 3e-4 \
    --simulatedData

(Note: we had a short run by accident with only 10000 photons and theta 2.1 to 2.4 that is stored in a directory with a 10000photons part in the name) <2023-04-16 Sun 13:46> There's some weirdness going on with the training. The outputs look fine and the accuracy is good, but the training data loss function is O(>100,000). Uhh.

In the meantime: looking at run 340 for energy 0.525 keV that has effective efficiency of ~90%. The biggest difference is in the length dataset! Our fake events are too long at the moment. Is the diffusion still bad for these datasets? YES. Hmmm.

• [ ] Why energy so wide in CDL data?
• [ ] why length too long?

SO to understand this:

• [ ] check if CDL data we use in plotting of all properties is representative. Shouldn't it have cuts to CDL charges?
• [ ] check rmsTransverse fit plots for this run. Is the fit good? Does it have still too much rubbish in there? Explanation of too much data for energy could be.

shows the energy seen in the CuEPIC0.9kV dataset (0.525 keV data) split by run based on CDL energy and charge energy. Run 340 has "all" energies, and is not correctly filtered?? What the heck. The code should be filtering both though. Is it our calibration of the energy being broken? I don't understand. I think it's the following:

  if not fileIsCdl:
    ## We don't need to write the calibrated dataset to the `calibration-cdl` file, because
    ## the data is copied over from the `CDL_Reco` file!
    let df = df.filter(f{idx("Cut?") == "Raw"})
      .clone() ## Clone to make sure we don't access full data underlying the DF!
    let dset = h5f.create_dataset(grpName / igEnergyFromCdlFit.toDset(),
                                  df.len,
                                  dtype = float,
                                  overwrite = true,
                                  filter = filter)
    dset.unsafeWrite(df["Energy", float].toUnsafeView(), df.len)

in cdl_spectrum_creation. Combining the toUnsafeView with the filter call is unsafe. Hence I added the clone line there now! See all plots with the bug present here: UPDATE: Upon further inspection after fixing the code in other ways again, it rather seems to be a different issue. The updated plots are here: (Note that the missing of the energy in Mn-Cr and Cu-Nu-15 in the former plots is a cut on < 3 keV data in the plotting code). With one exception (run 347) all runs that show behavior of very varying energies in the charge energy dataset are those from runs without the FADC! Why is that?

UPDATE <2023-04-16 Sun 18:43> I finally figured out the real problem. It was related to the above snippet, but not in the way I initially thought! The unsafe view was not the real culprit. The actual issue is that the calculated energies were not ordered correctly! We were just writing the data in the order in which it was in the data frame, which at this point in the code sometimes does not align correctly with the other datasets (due to a groups call over run numbers it did line up correctly for exactly one run I think). This is now fixed and the histograms of the energy (And yes, some of the runs really have that little statistics after the cuts! Check out: for an overview of the fits etc)

Let's look at the plot of the rmsTransverse from run 340 again: It has changed a bit. Can't say it's very obvious, but the shape is clearly less gaussian than before and the width to the right is a bit larger? What do the properties look like comparing run 340 now? IMPORTANT The size difference has actually become worse now! However, the effective efficiency still improves despite that. We probably should look at the determination of the diffusion again now that we actually look at the correct data for the CDL! Maybe the gauss fit is now actually viable.

Now that we actually look at the correct data for all the CDL runs, let's look at the effective efficiencies again using the gauss diffusion network without gain from <2023-04-12 Wed>: Aside from the 270 eV datasets all effective efficiencies have come much closer to the target efficiency! This is pretty great. See below for again looking at datasets comparing properties for the worst.

Let's also look at this using the new network that was trained on the new fake data!

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/16_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusion.pt \
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/16_04_23_sgd_sim_diffusion_gain/effective_eff

It looks surprisingly similar to the other network. It seems like regardless they learned similar things in spite of the other network having seen quite different data.

• [X] Look at rmsTransverse fit of run from earlier today again!
• [X] compare number of counts in the histograms before and after fix!

We'll continue the training of the model we created above

• [ ] LOOK At diffusion determination again after energy bug is fixed! Maybe the gauss fit method is now viable. In particular think about run 340 and how its length in the generated data is still too large! Also try to understand why runs for 250 eV are still so bad. ^– These are the properties of the worst run in terms of effective efficiency. By far the biggest issues are related to the size of the generated events and their eccentricity. We apparently generate way too many neighbors in this case? Investigate maybe.
• [ ] STILL have to adjust the gas gain and diffusion ranges for the fake data in training to be more sensible towards the real numbers!

Let's start by fixing the code to continue the training process. We've added the accuracy and loss fields to the MLPDesc, updated serialize.nim to deal with missing data better (just left empty) and will no rerun the training of the network from yesterday. We'll remove the full network and start from 0 to 100,000 epochs and then start it again to test (but this time containing the loss and accuracy values in MLPDesc H5 file).

So again (see MLP section below) While MLP is running back to the diffusion from the bad runs, 340 and also 342. 342 the sizes and diffusion is actually quite acceptable, but the eccentricity is pretty wrong, likely due to too many neighbors? Maybe need a lower threshold after all? What is the gas gain of 342 compared to other CDL runs? 2350 or so, compared to 2390 for 340. So comparable.

First 340: using gaussian fit with scale 1.65 again instead of limiting to 10% of height of peak. Better length & diffusion, but still too large. Let's check run 241 though before we lower further: -> This is already too small! So the fit itself is not going to be a good solution. Comparing the RMS plots: Given that drop off on RHS is MUCH stronger in 55Fe run, let's try a fixed value as peak + fixed. The run 241 data indicates something like 0.15 to the peak position (0.95 + 0.15 = 1.1). Run 241: And 340:

This matches well for run 241, but in case of 340 it is way too short a cutoff.

What defines the hardness of the cutoff? Number of pixels? Not quite, run 342 is narrower than 340. But 342 used FADC and 340 did not! Fixed cutoff is also a fail for 342 though! Could we determine it by simulation? Do a simple optimization that uses fit or something as start, rmsTransverse of events using diffusion and energy, then if small enough, cutoff? Then maybe use Kolmogorov-Smirnov to determine if agreement? Seems like maybe the most sane approach?

• [ ] Should we also store the training data of the net in MLPDesc?

We've implemented the ideas from yesterday to determine the diffusion based on simulating events itself using some optimization strategy (currently Newton combined with Cramér-von Mises).

The annoying part is we need to look at not only one axis, but actually both and then determine the long and short axes. Otherwise our estimate of the rms transverse ends up too large, because we look at the mean RMS instead of transverse RMS.

Having implemented it (with probably improvable parameters) and running it yields: which is actually quite good.

Let's look at run 333 (1.49 keV), one of the lowest ones in the efficiency. Was our value for the diffusion bad in that one? Hmm, not really. It looks quite reasonable. But, as expected: this run did not use the FADC! Probably NN filters out double photon events. Let's check if all "bad" ones are no FADC runs (y = FADC, n = no FADC): 0.93 keV runs:

• 335 (y), 336 (n), 337 (n) -> 335 is the best one, 336 and 337 are indeed the ones with low efficiency!

0.525 keV runs:

• 339 (y), 340 (n) -> 339 is the good one, 340 is the bad one!

So yes, this really seems to be the cause!

Let's check our new network from yesterday on the same:

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --ε 0.8 --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/17_04_23_sgd_sim_diffusion_gain/effective_eff/

As one could have hoped, the results look a bit better still.

Let's also look at the Run-2 data:

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --ε 0.8 --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/17_04_23_sgd_sim_diffusion_gain/effective_eff/

shows a plot of both run periods together.

We switched from Newton's method to gradient descent. It is much more stable in the approach to the minimum: n The outliers to very low values are again all CDL runs without the FADC! The 5.9 CAST data does indeed still have values quite a bit below target. But the spread is quite nice. The CDL data generally comes closer to the target, but overall has a larger spread (but notice not within a target! with the exception of the no FADC runs).

From here: look at some more distributions, but generally this is good to go in my opinion. Need to implement logic into application of the MLP method now that takes care of not only calculating the correct cut for each run & energy, but also writes the actual efficiency that the real data saw for the target to the output files, as well as the variance of all runs maybe.

Idea: Can we somehow find out which events are thrown out by the network that should be kept according to the simulated data cut value? More or less the 'difference' between the 80% cut based on the real data and the 80% cut based on the simulated data. The events in that region are thrown out. Look at some and generate distributions of those events to understand what they are? Maybe they are indeed "rubbish"? -> Events in the two 80% cuts of real and simulated. UPDATE: <2023-04-18 Tue 16:19> Hmm, at a first glance there does not seem to be much to see there. To a "human eye" they look like X-rays I would say. It's likely the correlation between different variables that is making them more background like than X-ray like. I implemented a short function that filters the data between the cut region of real and simulated data at target efficiency:

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/analyze_5.9/ \
    --run 241

which gives the following property plots:

Next: Check how many elements we really need to generate for our Newton optimization. And can we:

1. make the difference smaller again? -> Let's try with diffusion difference of 1 in numerical derivative (at 10k samples)
2. use Dual after all? -> No. Tried again, but Cramér-von Mises also destroys all derivative information.
3. only use about 1000 elements instead of 10000? -> at 1000 the spread becomes quite a bit bigger than at 10k!
4. [ ] Put 3.0 escape photon data for CAST back into plot!

So, from yesterday to summarize: Using gradient descent makes the determination of the diffusion from real data fast enough and stable.

First continue a short look at the 'intermediate' data between the cut values of simulated and real data. We will plot the intermediate against the data that passes both. Running the same command as yesterday for run 241 yields: In direct comparison we can see that the events that fail to pass the simulated cut are generally a bit more eccentric and longer / wider / slightly bigger RMS. Potentially events including a captured escape photon? Let's extract the event numbers of those intermediate events:

Event numbers that are intermediate: Column of type: int with length: 107
  contained Tensor: Tensor[system.int] of shape "[107]" on backend "Cpu"
    104     109     123     410     537     553     558    1042    1272    1346    1390    1447    1527    1583    1585    1594    1610    1720    1922    1965    2082    2155    2176    2198    2419    2512    2732    2800    2801    3038    3072    3095    3296    3310    3473    3621    3723    3820    4088    4145    4184    4220    4250    4308    4347    4353    4466    4558    4590    4725    4843    4988    5204    5234    5288    5497    5637    5648    5661    5792    5814    5848    5857    6090    6175    6187    6312    6328    6359    6622    6657    6698    6843    6944    6951    7121    7137    7162    7192    7350    7436    7472    7545    7633    7634    7730    7749    7788    7936    8003    8014    8075    8214    8364    8441    8538    8549    8618    8827    8969    8991    9008    9026    9082    9102    9260    9292

Now we'll use plotData to generate event displays for them:

plotData \
    --h5file ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --runType rtCalibration \
    --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
    --runs 241 \
    --eventDisplay \
    --septemboard \
    --events 104 --events 109 --events 123 --events 410 --events 537 --events 553 --events 558 --events 1042 --events 1272 --events 1346 --events 1390 --events 1447 --events 1527 --events 1583 --events 1585 --events 1594 --events 1610 --events 1720 --events 1922 --events 1965 --events 2082 --events 2155 --events 2176 --events 2198 --events 2419 --events 2512 --events 2732 --events 2800 --events 2801 --events 3038 --events 3072 --events 3095 --events 3296 --events 3310 --events 3473 --events 3621 --events 3723 --events 3820 --events 4088 --events 4145 --events 4184 --events 4220 --events 4250 --events 4308 --events 4347 --events 4353 --events 4466 --events 4558 --events 4590 --events 4725 --events 4843 --events 4988 --events 5204 --events 5234 --events 5288 --events 5497 --events 5637 --events 5648 --events 5661 --events 5792 --events 5814 --events 5848 --events 5857 --events 6090 --events 6175 --events 6187 --events 6312 --events 6328 --events 6359 --events 6622 --events 6657 --events 6698 --events 6843 --events 6944 --events 6951 --events 7121 --events 7137 --events 7162 --events 7192 --events 7350 --events 7436 --events 7472 --events 7545 --events 7633 --events 7634 --events 7730 --events 7749 --events 7788 --events 7936 --events 8003 --events 8014 --events 8075 --events 8214 --events 8364 --events 8441 --events 8538 --events 8549 --events
8618 --events 8827 --events 8969 --events 8991 --events 9008 --events 9026 --events 9082 --events 9102 --events 9260 --events 9292

It really seems like the main reason for them being a bit more eccentric is the presence of potentially higher number of single pixels that are outliers? I suppose it makes sense that these would be rejected with a larger likelihood. However, why would such events appear more likely in the 5.9 keV CAST data than in simulated data?

For now I'm not sure what to make of this. Of course one could attempt to do something like use a clustering algorithm that ignores outliers etc. But all that doesn't seem that interesting, at least not if we don't understand why we have this distinction from simulated to real data in the first place.

So for now let's just continue.

First: add 3.0 keV data with 5.9 keV absorption length into effective efficiency plot and see how the efficiency fares:

Run-2 @ 99%:

likelihood \
    -f ~/CastData/data/DataRuns2017_Reco.h5 \
    --h5out ~/Sync/run2_17_04_23_mlp_local_0.99.h5 \
    --region crGold \
    --cdlYear 2018 \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --mlp ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --nnSignalEff 0.99 \
    --nnCutKind runLocal \
    --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --readonly

Run-3 @ 99%:

likelihood \
    -f ~/CastData/data/DataRuns2018_Reco.h5 \
    --h5out ~/Sync/run3_17_04_23_mlp_local_0.99.h5 \
    --region crGold \
    --cdlYear 2018 \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --mlp ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --nnSignalEff 0.99 \
    --nnCutKind runLocal \
    --calibFile ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --readonly

Need to run effective efficiency at 99% as well to know the real efficiencies.

NOTE: The gradient descent code unfortunately is still rather slow. Should we aim to cache the results for each run? Generally I'd prefer no caching and just making it faster though. Maybe we can pick an example that is pretty hard to converge and use that as a reference to develop? Ah, and we can cache stuff within a single program run, so that we only need to compute the stuff once for each run?

• we know cache the results of the diffusion calculation, but still have to verify it works as intended
• for the background datasets the diffusion coefficient is now also ~correctly determined by sampling both the energy and position of the data, i.e.
  • uniform in drift distance (muons enter anywhere between cathode and anode)
  • energy exponential distribution that has ~20-30% of flux at 10 keV of the flux at ~0 eV
  • [ ] ADD PLOT
• verify the diffusion coefficients are correctly determined by calculating all of them for background and calibration and creating a plot showing the determined numbers -> maybe a good option to show the plots of the RMS data used and the generated data for each?
• [ ] implement calculation of variance and mean value for effective efficiencies for each run period (write to H5 out in likelihood)
• [X] TAKE NOTES of the diffusion vs run plot
• [ ] ADD EFFECTICE EFFICIENCY PLOT WITH 3 keV escape data

Yesterday we turned determineDiffusion into a usable standalone tool to compute the diffusion parameters of all runs. It generates a plot showing what the parameters are, colored by the 'loss' of the best estimate (via Cramér-von Mises).

./determineDiffusion \
    ~/CastData/data/DataRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/DataRuns2018_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5

yields the following plot of the diffusion values deduced using gradient descent:

We can clearly see two important aspects:

1. the background diffusion parameters are typically a good 20-30 μm/√cm larger than the 5.9 Fe points
2. there are some cases where:
   • the loss is very large (background, runs close to run 80)
   • the diffusion is exactly 660, our starting value

The latter implies that the loss is already smaller than 2 (our current stopping criterion). The former implies the distributions look very different than the simulated distribution. I've already seen one reason: in some runs there is an excessive contribution at rmsTransvers < 0.2 or so which has a huge peak. Some kind of noise signal? We'll check that.

For now here are all the plots of the rmsTransverse dataset for sim and real data that were generated from the best estimates:

The worst run is run 86 (see page 158 in the PDF). Let's investigate what kind of events contribute such data to the transverse RMS. My assumption would be noisy events at the top of the chip or something like that? First the ingrid distributions:

plotData \
    --h5file ~/CastData/data/DataRuns2017_Reco.h5 \
    --runType rtBackground \
    --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
    --runs 86 \
    --ingrid \
    --cuts '("rmsTransverse", 0.0, 0.2)' \
    --applyAllCuts \
    --chips 3

So clearly defined in centerX and centerY. The sparky point on the grid?

Yeah, that seems to be it. They are all pretty much below 20 hits anyway. So we'll just introduce an additional cut for non CDL data that clusters must be larger 20 hits. They match the pixels we consider noisy for the background cluster plot and in limit!

Added such a filter (hits > 20) and running again. The diffusion values are now these. The largest is on the order of 3 and they look more reasonable. Of course the distinction between background and calibration is still presetn.

Next we've removed the limit of ks > 2.0 as a stopping criterion and added logic to reset to best estimate if half number of bad steps already done.

Let's rerun again: Now it actually looks pretty nice!

The difference between 5.9 and background is now pretty much exactly 40. So we'll use that in getDiffusion as a constant for now.

Finally, let's try to run likelihood on one dataset (Run-2) at a fixed efficiency to see what is what (it should read all diffusion values from the cache).

The generated files are:

• /Sync/run2_17_04_23_mlp_local_0.80.h5
• /Sync/run2_17_04_23_mlp_local_0.90.h5
• /Sync/run3_17_04_23_mlp_local_0.80.h5
• /Sync/run3_17_04_23_mlp_local_0.90.h5

Time to plot a background rate from both, we compare it to the LnL only method at 80%:

plotBackgroundRate \
    ~/Sync/run2_17_04_23_mlp_local_0.90.h5 \
    ~/Sync/run3_17_04_23_mlp_local_0.90.h5 \
    ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run2_crGold.h5 \
    ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run3_crGold.h5 \
    --names "MLP@90" \
    --names "MLP@90" \
    --names "LnL@80" \
    --names "LnL@80" \
    --centerChip 3 \
    --title "Background rate from CAST data, LnL@80%, MLP@90%" \
    --showNumClusters \
    --showTotalTime \
    --topMargin 1.5 \
    --energyDset energyFromCharge \
    --outfile background_rate_only_lnL_0.8_mlp_0.9.pdf \
    --outpath ~/Sync/ \
    --quiet

[INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.3288e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.9406e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.7298e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 2.2748e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 6.2088e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 3.1044e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 3.4122e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 1.7061e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.1626e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.5836e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 9.6914e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.1537e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9703e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5881e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 1.0413e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 4.1650e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.5958e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 3.3067e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 8.2667e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.6551e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.0689e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.9383e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.4229e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.1788e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.3631e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@90 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 9.3572e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.5595e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:DataFrame with 7 columns and 116 rows: Idx Energy Rate totalTime RateErr Dataset yMin yMax dtype: float float constant float string float float 0 0.4 5.98 3159 1.026 LnL@80 4.955 7.006 1 0.6 6.332 3159 1.055 LnL@80 5.277 7.387 2 0.8 5.98 3159 1.026 LnL@80 4.955 7.006 3 1 5.98 3159 1.026 LnL@80 4.955 7.006 4 1.2 3.166 3159 0.7462 LnL@80 2.42 3.912 5 1.4 2.99 3159 0.7252 LnL@80 2.265 3.715 6 1.6 2.814 3159 0.7036 LnL@80 2.111 3.518 7 1.8 3.166 3159 0.7462 LnL@80 2.42 3.912 8 2 1.583 3159 0.5277 LnL@80 1.055 2.111 9 2.2 1.759 3159 0.5562 LnL@80 1.203 2.315 10 2.4 0.8794 3159 0.3933 LnL@80 0.4861 1.273 11 2.6 1.583 3159 0.5277 LnL@80 1.055 2.111 12 2.8 2.287 3159 0.6342 LnL@80 1.652 2.921 13 3 5.101 3159 0.9472 LnL@80 4.154 6.048 14 3.2 5.277 3159 0.9634 LnL@80 4.313 6.24 15 3.4 4.221 3159 0.8617 LnL@80 3.36 5.083 16 3.6 3.342 3159 0.7667 LnL@80 2.575 4.109 17 3.8 2.111 3159 0.6093 LnL@80 1.501 2.72 18 4 0.7036 3159 0.3518 LnL@80 0.3518 1.055 19 4.2 1.055 3159 0.4308 LnL@80 0.6245 1.486

[INFO]:INFO: storing plot in /home/basti/Sync/background_{rateonlylnL0.8mlp0.9.pdf}

As we can see the rates are quite comparable. What is nice to see is that the Argon peak is actually more visible in the MLP data than in the LnL cut method. Despite higher efficiency the network performs essentially the same though. That's very nice! The question is what is the effective efficiency based on 5.9 keV data using 90%?

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/17_04_23_sgd_gauss_diffusion/trained_mlp_sgd_gauss_diffusioncheckpoint_epoch_400000_loss_0.0115_acc_0.9975.pt \
    --ε 0.9 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/run2_run3_gradient_descent_eff_0.9/

So we're actually looking at about 85% realistically.

And background rate with MLP@85%:

plotBackgroundRate \
    ~/Sync/run2_17_04_23_mlp_local_0.85.h5 \
    ~/Sync/run3_17_04_23_mlp_local_0.85.h5 \
    ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run2_crGold.h5 \
    ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run3_crGold.h5 \
    --names "MLP@85" \
    --names "MLP@85" \
    --names "LnL@80" \
    --names "LnL@80" \
    --centerChip 3 \
    --title "Background rate from CAST data, LnL@80%, MLP@85%" \
    --showNumClusters \
    --showTotalTime \
    --topMargin 1.5 \
    --energyDset energyFromCharge \
    --outfile background_rate_only_lnL_0.8_mlp_0.85.pdf \
    --outpath ~/Sync/ \
    --quiet

[INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.3288e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.9406e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.4501e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 2.0418e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 6.2088e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 3.1044e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 1.3192e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.1626e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.5836e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 8.4250e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 1.8722e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9703e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5881e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9879e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5952e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.5958e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.9022e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 7.2554e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.6551e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.0689e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.7114e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.1392e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.1788e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.3631e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.5130e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.4188e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:DataFrame with 7 columns and 116 rows: Idx Energy Rate totalTime RateErr Dataset yMin yMax dtype: float float constant float string float float 0 0.4 5.98 3159 1.026 LnL@80 4.955 7.006 1 0.6 6.332 3159 1.055 LnL@80 5.277 7.387 2 0.8 5.98 3159 1.026 LnL@80 4.955 7.006 3 1 5.98 3159 1.026 LnL@80 4.955 7.006 4 1.2 3.166 3159 0.7462 LnL@80 2.42 3.912 5 1.4 2.99 3159 0.7252 LnL@80 2.265 3.715 6 1.6 2.814 3159 0.7036 LnL@80 2.111 3.518 7 1.8 3.166 3159 0.7462 LnL@80 2.42 3.912 8 2 1.583 3159 0.5277 LnL@80 1.055 2.111 9 2.2 1.759 3159 0.5562 LnL@80 1.203 2.315 10 2.4 0.8794 3159 0.3933 LnL@80 0.4861 1.273 11 2.6 1.583 3159 0.5277 LnL@80 1.055 2.111 12 2.8 2.287 3159 0.6342 LnL@80 1.652 2.921 13 3 5.101 3159 0.9472 LnL@80 4.154 6.048 14 3.2 5.277 3159 0.9634 LnL@80 4.313 6.24 15 3.4 4.221 3159 0.8617 LnL@80 3.36 5.083 16 3.6 3.342 3159 0.7667 LnL@80 2.575 4.109 17 3.8 2.111 3159 0.6093 LnL@80 1.501 2.72 18 4 0.7036 3159 0.3518 LnL@80 0.3518 1.055 19 4.2 1.055 3159 0.4308 LnL@80 0.6245 1.486

[INFO]:INFO: storing plot in /home/basti/Sync/background_{rateonlylnL0.8mlp0.85.pdf} [WARNING]: Printing total background time currently only supported for single datasets.

So at this point the effective efficiency is pretty much 80%. This means the network is significantly better at lower energies (aside from first bin), but essentially the same everywhere else. Pretty interesting. Question is how it fares at higher signal efficiencies!

Wrote a small helper studyMLP that runs:

• effective efficiency
• likelihood for 2017 and 2018 data

for a list of target efficiencies for a given MLP model.

./studyMLP 0.85 0.95 0.99 \
           --model ~/org/resources/nn_devel_mixing/22_04_23_sgd_gauss_diffusion_2hidden/trained_mlp_sgd_gauss_diffusion_2hiddencheckpoint_epoch_400000_loss_0.0067_acc_0.9984.pt \
           --plotPath ~/Sync/22_04_23_sgd_sim_diffusion_gain_2hidden/ \
           --h5out ~/Sync/22_04_23_tanh_hidden2_100.h5

Let's see!

Hmm, we have some issue that the likelihood output files are on the order of 1.5 - 2.5 GB! Investigating… Ah: I didn't recompile likelihood, so it still used the single layer MLP instead of the two layer. Therefore the weights obviously didn't properly match.

The above ran properly now and produced the following effective efficiency plots:

Background rate tanh 2 layer 100 neurons 85%:

plotBackgroundRate ~/Sync/22_04_23_tanh_hidden2_100eff_0.85_run2.h5 ~/Sync/22_04_23_tanh_hidden2_100eff_0.85_run3.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run2_crGold.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run3_crGold.h5 --names "MLP@85" --names "MLP@85" --names "LnL@80" --names "LnL@80" --centerChip 3 --title "Background rate from CAST data, LnL@80%, SGD tanh100 MLP@85%" --showNumClusters --showTotalTime --topMargin 1.5 --energyDset energyFromCharge --outfile background_rate_only_lnL_0.8_sgd_tanh100_22_04_23_mlp_0.85.pdf --outpath ~/Sync/ --quiet

[INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.3288e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.9406e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 1.9594e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.6328e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 6.2088e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 3.1044e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 2.0579e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 1.0289e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.1626e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.5836e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 7.4577e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 1.6573e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9703e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5881e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 5.5405e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 2.2162e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.5958e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.5328e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.3320e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.6551e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.0689e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.3033e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 1.6292e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.1788e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.3631e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@85 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 7.7918e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.2986e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:DataFrame with 7 columns and 116 rows: Idx Energy Rate totalTime RateErr Dataset yMin yMax dtype: float float constant float string float float 0 0.4 5.98 3159 1.026 LnL@80 4.955 7.006 1 0.6 6.332 3159 1.055 LnL@80 5.277 7.387 2 0.8 5.98 3159 1.026 LnL@80 4.955 7.006 3 1 5.98 3159 1.026 LnL@80 4.955 7.006 4 1.2 3.166 3159 0.7462 LnL@80 2.42 3.912 5 1.4 2.99 3159 0.7252 LnL@80 2.265 3.715 6 1.6 2.814 3159 0.7036 LnL@80 2.111 3.518 7 1.8 3.166 3159 0.7462 LnL@80 2.42 3.912 8 2 1.583 3159 0.5277 LnL@80 1.055 2.111 9 2.2 1.759 3159 0.5562 LnL@80 1.203 2.315 10 2.4 0.8794 3159 0.3933 LnL@80 0.4861 1.273 11 2.6 1.583 3159 0.5277 LnL@80 1.055 2.111 12 2.8 2.287 3159 0.6342 LnL@80 1.652 2.921 13 3 5.101 3159 0.9472 LnL@80 4.154 6.048 14 3.2 5.277 3159 0.9634 LnL@80 4.313 6.24 15 3.4 4.221 3159 0.8617 LnL@80 3.36 5.083 16 3.6 3.342 3159 0.7667 LnL@80 2.575 4.109 17 3.8 2.111 3159 0.6093 LnL@80 1.501 2.72 18 4 0.7036 3159 0.3518 LnL@80 0.3518 1.055 19 4.2 1.055 3159 0.4308 LnL@80 0.6245 1.486

[INFO]:INFO: storing plot in /home/basti/Sync/background_{rateonlylnL0.8sgdtanh100220423mlp0.85.pdf} [WARNING]: Printing total background time currently only supported for single datasets.

plotBackgroundRate ~/Sync/22_04_23_tanh_hidden2_100eff_0.95_run2.h5 ~/Sync/22_04_23_tanh_hidden2_100eff_0.95_run3.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run2_crGold.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run3_crGold.h5 --names "MLP@95" --names "MLP@95" --names "LnL@80" --names "LnL@80" --centerChip 3 --title "Background rate from CAST data, LnL@80%, SGD tanh100 MLP@95%" --showNumClusters --showTotalTime --topMargin 1.5 --energyDset energyFromCharge --outfile background_rate_only_lnL_0.8_sgd_tanh100_22_04_23_mlp_0.95.pdf --outpath ~/Sync/ --quiet

[INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.3288e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.9406e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.6278e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 2.1898e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 6.2088e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 3.1044e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 3.3771e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 1.6885e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.1626e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.5836e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 9.8673e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.1927e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9703e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5881e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 9.1462e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.6585e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.5958e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 3.2715e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 8.1788e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.6551e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.0689e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.8275e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.2843e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.1788e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.3631e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@95 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 9.5155e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.5859e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:DataFrame with 7 columns and 116 rows: Idx Energy Rate totalTime RateErr Dataset yMin yMax dtype: float float constant float string float float 0 0.4 5.98 3159 1.026 LnL@80 4.955 7.006 1 0.6 6.332 3159 1.055 LnL@80 5.277 7.387 2 0.8 5.98 3159 1.026 LnL@80 4.955 7.006 3 1 5.98 3159 1.026 LnL@80 4.955 7.006 4 1.2 3.166 3159 0.7462 LnL@80 2.42 3.912 5 1.4 2.99 3159 0.7252 LnL@80 2.265 3.715 6 1.6 2.814 3159 0.7036 LnL@80 2.111 3.518 7 1.8 3.166 3159 0.7462 LnL@80 2.42 3.912 8 2 1.583 3159 0.5277 LnL@80 1.055 2.111 9 2.2 1.759 3159 0.5562 LnL@80 1.203 2.315 10 2.4 0.8794 3159 0.3933 LnL@80 0.4861 1.273 11 2.6 1.583 3159 0.5277 LnL@80 1.055 2.111 12 2.8 2.287 3159 0.6342 LnL@80 1.652 2.921 13 3 5.101 3159 0.9472 LnL@80 4.154 6.048 14 3.2 5.277 3159 0.9634 LnL@80 4.313 6.24 15 3.4 4.221 3159 0.8617 LnL@80 3.36 5.083 16 3.6 3.342 3159 0.7667 LnL@80 2.575 4.109 17 3.8 2.111 3159 0.6093 LnL@80 1.501 2.72 18 4 0.7036 3159 0.3518 LnL@80 0.3518 1.055 19 4.2 1.055 3159 0.4308 LnL@80 0.6245 1.486

[INFO]:INFO: storing plot in /home/basti/Sync/background_{rateonlylnL0.8sgdtanh100220423mlp0.95.pdf} [WARNING]: Printing total background time currently only supported for single datasets.

plotBackgroundRate ~/Sync/22_04_23_tanh_hidden2_100eff_0.99_run2.h5 ~/Sync/22_04_23_tanh_hidden2_100eff_0.99_run3.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run2_crGold.h5 ~/org/resources/lhood_limits_automation_correct_duration/likelihood_cdl2018_Run3_crGold.h5 --names "MLP@99" --names "MLP@99" --names "LnL@80" --names "LnL@80" --centerChip 3 --title "Background rate from CAST data, LnL@80%, SGD tanh100 MLP@99%" --showNumClusters --showTotalTime --topMargin 1.5 --energyDset energyFromCharge --outfile background_rate_only_lnL_0.8_sgd_tanh100_22_04_23_mlp_0.99.pdf --outpath ~/Sync/ --quiet

[INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 2.3288e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 1.9406e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 0.0 .. 12.0: 3.3911e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 12.0: 2.8259e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 6.2088e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 3.1044e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 0.5 .. 2.5: 5.9450e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 2.5: 2.9725e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.1626e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.5836e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 0.5 .. 5.0: 1.3192e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.5 .. 5.0: 2.9315e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 8.9703e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 3.5881e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 0.0 .. 2.5: 1.4599e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 2.5: 5.8395e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 2.6383e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 6.5958e-06 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 4.0 .. 8.0: 4.0630e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 4.0 .. 8.0: 1.0158e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 1.6551e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 2.0689e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 0.0 .. 8.0: 2.5082e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 0.0 .. 8.0: 3.1352e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: LnL@80 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 8.1788e-05 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.3631e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:Dataset: MLP@99 [INFO]: Integrated background rate in range: 2.0 .. 8.0: 1.0923e-04 cm⁻²·s⁻¹ [INFO]: Integrated background rate/keV in range: 2.0 .. 8.0: 1.8204e-05 keV⁻¹·cm⁻²·s⁻¹ [INFO]:DataFrame with 7 columns and 116 rows: Idx Energy Rate totalTime RateErr Dataset yMin yMax dtype: float float constant float string float float 0 0.4 5.98 3159 1.026 LnL@80 4.955 7.006 1 0.6 6.332 3159 1.055 LnL@80 5.277 7.387 2 0.8 5.98 3159 1.026 LnL@80 4.955 7.006 3 1 5.98 3159 1.026 LnL@80 4.955 7.006 4 1.2 3.166 3159 0.7462 LnL@80 2.42 3.912 5 1.4 2.99 3159 0.7252 LnL@80 2.265 3.715 6 1.6 2.814 3159 0.7036 LnL@80 2.111 3.518 7 1.8 3.166 3159 0.7462 LnL@80 2.42 3.912 8 2 1.583 3159 0.5277 LnL@80 1.055 2.111 9 2.2 1.759 3159 0.5562 LnL@80 1.203 2.315 10 2.4 0.8794 3159 0.3933 LnL@80 0.4861 1.273 11 2.6 1.583 3159 0.5277 LnL@80 1.055 2.111 12 2.8 2.287 3159 0.6342 LnL@80 1.652 2.921 13 3 5.101 3159 0.9472 LnL@80 4.154 6.048 14 3.2 5.277 3159 0.9634 LnL@80 4.313 6.24 15 3.4 4.221 3159 0.8617 LnL@80 3.36 5.083 16 3.6 3.342 3159 0.7667 LnL@80 2.575 4.109 17 3.8 2.111 3159 0.6093 LnL@80 1.501 2.72 18 4 0.7036 3159 0.3518 LnL@80 0.3518 1.055 19 4.2 1.055 3159 0.4308 LnL@80 0.6245 1.486

[INFO]:INFO: storing plot in /home/basti/Sync/background_{rateonlylnL0.8sgdtanh100220423mlp0.99.pdf} [WARNING]: Printing total background time currently only supported for single datasets.

All together:

So clearly we also mainly gain at low energies here. Even at 99% we still are effectively not seeing any degradation in most regions!

IMPORTANT:

• See if adjusting the rms transverse values down to the old ones (~1.0 to 1.1 instead of 1.2 - 1.3) gets us closer to the 80% in case of the e.g. 80% network efficiency!!!

Finished the implementation of mapping all parameters of interest for the model layout and training to the MLPDesc and loading / applying the correct thing at runtime.

In order to update any of the existing mlp_desc.h5 files, we need to provide the settings used for each network! At least the new settings that is.

For example let's update the last tanh network we trained with a sigmoid output layer. We do it by also continuing training on it for another 100k epochs:

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
    --plotPath ~/Sync/23_04_23_sgd_sim_diffusion_gain_tanh300_mse_loss/ \
    --numHidden 300 \
    --numHidden 300 \
    --activationFunction tanh \
    --outputActivation sigmoid \
    --lossFunction MSE \
    --optimizer SGD \
    --learningRate 7e-4 \
    --simulatedData \
    --continueAfterEpoch 500000

As we have the existing H5 file for the MLPDesc we don't need to supply the datasets. But the additional new fields are required to update them in the file.

We implemented that the code raises an exception if only an old MLPDesc file is found in effective_eff_fe55. Also the file name now contains the version number itself for convenience. Also the version is another field now that is serialized.

Let's check the following TODO from yesterday:

IMPORTANT:

• See if adjusting the rms transverse values down to the old ones (~1.0 to 1.1 instead of 1.2 - 1.3) gets us closer to the 80% in case of the e.g. 80% network efficiency!!!

-> We'll reset the cuts to their old values and then run the above network on effective_eff_fe55

./effective_eff_55fe \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --model ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
    --ε 0.8 \
    --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
    --evaluateFit \
    --plotDatasets \
    --plotPath ~/Sync/run2_run3_23_04_23_tanh300_rmsT_old/

-> The plot generally doesn't look too bad, but it's hard to read because there are 2 escape peak values with negative efficiency! Need to fix that.

We'll look at only run 88, which should show that behavior.

Fake data for 3.0 at run 88 and energy 2.98 keV target Ag-Ag-6kV has a cut value 0.9976992011070251 and effective eff -0.01075268817204301

That's fishy. Computed from (kept vs total in that dataset):

Number kept: -1 vs 93

Only 93 escape events? And -1 kept? Ohhh! It's because of this line:

 if pred.len < 100: return (-1, @[])

in predictCut!

Removed that line as it's not very useful nowadays. We still care about the number.

Rerunning and storing in: This looks pretty much like the other case with wider rmsT, no? Let's make sure and change again and run again.

-> Ok, it does improve the situation. The mean is now rather at 75-76 percent compared to 73% or so.

What to make of this? The more rms transverse we allow, the more ugly non X-ray events we will have in our data? Or just that we reproduce the ugly events with our fake data?

Next steps:

1. implement computing the effective efficiencies for:

   • the given network
   • the desired efficiency

   These will also be stored in some form of a cache (given that the calculation takes O(15 min) -> Implemented including caching. Just needs to be called from likelihood now. VetoSettings has the fields required. Should we compute it first or last? -> Compute in initLikelihoodContext. -> Not in the init, due to problems with circular imports. The effective eff fields are now filled in likelihood itself. They are mainly needed for that anyway, so that should be fine.

2. fix the application of the septem veto etc when using an MLP. Use that to decide if a cluster passes instead of logL -> up next. -> Implemented now!
3. automate MLP for likelihood in createAllLikelihoodCombinations
   • this should be done via a --mlp option which takes seq[string] which are paths to the model file. For each we add a fkMLP flag that will perform a call to likelihood as a replacement fkLogL
4. [X] IMPORTANT: We were still using the old getDiffusion working with the fit for any code using readValidDsets and readCdlDset! Means both the diffusion values used in training for the background data as well as for anything with prediction based on real data likely had a wrong diffusion value! FIX and check the impact.

• Fixed the usage of the old getDiffusion in io_helpers. Now using the correct value from the cache table. Let's see the effect of that on the effective efficiency!

  ./effective_eff_55fe \
      ~/CastData/data/CalibrationRuns2017_Reco.h5 \
      ~/CastData/data/CalibrationRuns2018_Reco.h5 \
      --model ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
      --ε 0.8 \
      --cdlFile ~/CastData/data/CDL_2019/CDL_2019_Reco.h5 \
      --evaluateFit \
      --plotDatasets \
      --plotPath ~/Sync/run2_run3_23_04_23_tanh300_correct_diffusion
  

  The result is pretty much unchanged, if not possibly actually a bit worse than before. Yeah, comparing with the old version it's a bit worse even. :/

We've implemented MLP support into createAllLikelihoodCombinations now. The result is not super pretty, but it should work. The filename of the MLP is added to the output file name, leading to very long filenames. Ideally we'd have a better solution there.

Let's give it a test run, but before we do that, let's run a likelihood combining the MLP with FADC, septem and line veto.

Run2:

likelihood \
    -f ~/CastData/data/DataRuns2017_Reco.h5 \
    --h5out /t/run_2_mlp_all_vetoes.h5 \
    --region crGold \
    --cdlYear 2018 \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --mlp ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/23_04_23_tanh300_mse.pt \
    --nnSignalEff 0.99 \
    --vetoPercentile 0.99 \
    --lineVeto \
    --septemVeto \
    --calibFile ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --readonly

likelihood \
    -f ~/CastData/data/DataRuns2018_Reco.h5 \
    --h5out /t/run_3_mlp_all_vetoes.h5 \
    --region crGold \
    --cdlYear 2018 \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --mlp ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/23_04_23_tanh300_mse.pt \
    --nnSignalEff 0.99 \
    --vetoPercentile 0.99 \
    --lineVeto \
    --septemVeto \
    --calibFile ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --readonly

• [ ] Well, issue:
  1. need to find a way to better deal with the cache table for effective efficiencies! It uses strings in the keys, which we currently don't correctly support. Either change them to fixed length or fix the code for variable length strings as part compound types.
  2. when running likelihood our center event index is not assigned correctly for the septem veto! Something is up there with the predicted values or the cut values, such that we never enter the branch the sets centerEvIdx.

• [X] fixed the centerEvIdx bug: the problem was our comparison of the NN prediction value to the cut value interpolator. I missed that the LnL and NN values are inverted!

Running the program to create all likelihood combinations:

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crAll --regions crGold \
    --signalEfficiency 0.8 --signalEfficiency 0.85 --signalEfficiency 0.9 --signalEfficiency 0.95 \
    --fadcVetoPercentile 0.99 \
    --vetoSets "{fkMLP, fkFadc, fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --mlpPath ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
    --out ~/org/resources/lhood_limits_23_04_23_mlp/ \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing

So we output the files to: ~/org/resources/lhood_{limits230423mlp}

Running the code last night more or less worked, but in some cases there was a segfault running the code.

I spent most of today debugging this issue. gdb while useful was still pretty confusing. valgrind is still running…

It turns out that

1. confusion and bad stack traces are due to injected destructor calls. those break the line information
2. the issue is a crash in the destructor of a H5Group
3. only triggered with our try/except code. We've changed that to if/else now. Saner anyway.

With this fixed, we can rerun the likelihood combinations. Both using the sigmoid output layer MLP from <2023-04-23 Sun> and if possible one of the linear output ones. For now start with one of them.

• [ ] Note: when running with another MLP still need to regenerate the MLPDesc H5 file!

We'll rerun the command from yesterday. Still, can we speed up the process somehow?

I think we should introduce caching also for the CutValueInterpolator.

Finally fixed all the HDF5 data writing stuff.. Will write more about that tomorrow.

Will test likelihood producing a sensible file now: First a run to generate the file. Will take a short look, then read it in a separate script to see what it contains. Then finally rerun the same command to see if the fake event generation is then skipped.

likelihood -f /home/basti/CastData/data/DataRuns2018_Reco.h5 --h5out /t/blabla_broken.h5 --region=crGold --cdlYear=2018 --scintiveto --fadcveto --septemveto --lineveto --mlp /home/basti/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt --cdlFile=/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 --readOnly --calibFile=/home/basti/CastData/data/CalibrationRuns2018_Reco.h5 --vetoPercentile=0.99 --nnSignalEff=0.95

• [X] segfault due to string
• [ ]

    HDF5-DIAG: Error detected in HDF5 (1.10.5) thread 0:
    #000: H5Tcompound.c line 354 in H5Tinsert(): unable to insert member
      major: Datatype
      minor: Unable to insert object
    #001: H5Tcompound.c line 446 in H5T__insert(): member extends past end of compound type
      major: Datatype
      minor: Unable to insert object
  

  -> partially fixed it but:

-> We have alignment issues. The H5 library also seems to align data in some cases when necessary. Our code currently assumes there is no such thing.

So I guess we need to go back to an approach that does actually generate some helper type or what?

-> I think I got it all working now. Merged the two approaches of copying data to a Buffer and

Ok, the caching seems to work I think. And the generated file looks good.

Let's read the data and see what it contains:

import nimhdf5, tables

const CacheTabFile = "/dev/shm/cacheTab_runLocalCutVals.h5"
type
  TabKey = (int, string, float)
  #         ^-- run number
  #              ^-- sha1 hash of the NN model `.pt` file
  #                      ^-- target signal efficiency
  TabVal = seq[(string, float)]
  #             ^-- CDL target
  #                     ^-- MLP cut value
  CacheTabTyp = Table[TabKey, TabVal]
var tab = deserializeH5[CacheTabTyp](CacheTabFile)
for k, v in tab:
  echo "Key: ", k, " = ", v

Key: (295, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.974539938569069), ("Al-Al-4kV", 0.9415906369686127), ("C-EPIC-0.6kV", 0.7813522785902023), ("Cu-EPIC-0.9kV", 0.8288368076086045), ("Cu-EPIC-2kV", 0.8751996099948883), ("Cu-Ni-15kV", 0.9722824782133103), ("Mn-Cr-12kV", 0.9686738938093186), ("Ti-Ti-9kV", 0.9641254603862762)] Key: (270, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9771067887544632), ("Al-Al-4kV", 0.9460293263196945), ("C-EPIC-0.6kV", 0.7850535094738007), ("Cu-EPIC-0.9kV", 0.8130916118621826), ("Cu-EPIC-2kV", 0.8937846541404724), ("Cu-Ni-15kV", 0.9715777307748794), ("Mn-Cr-12kV", 0.9703928083181381), ("Ti-Ti-9kV", 0.9622162997722625)] Key: (285, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9717331320047379), ("Al-Al-4kV", 0.9252435654401779), ("C-EPIC-0.6kV", 0.7428321331739426), ("Cu-EPIC-0.9kV", 0.7788086831569672), ("Cu-EPIC-2kV", 0.8636017471551896), ("Cu-Ni-15kV", 0.9687924206256866), ("Mn-Cr-12kV", 0.9663917392492294), ("Ti-Ti-9kV", 0.9504937410354615)] Key: (267, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.975317457318306), ("Al-Al-4kV", 0.9462290376424789), ("C-EPIC-0.6kV", 0.776226544380188), ("Cu-EPIC-0.9kV", 0.834472405910492), ("Cu-EPIC-2kV", 0.8766408234834671), ("Cu-Ni-15kV", 0.9714102745056152), ("Mn-Cr-12kV", 0.9717013716697693), ("Ti-Ti-9kV", 0.9653892040252685)] Key: (240, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9765782624483108), ("Al-Al-4kV", 0.9469905078411103), ("C-EPIC-0.6kV", 0.8053286731243133), ("Cu-EPIC-0.9kV", 0.8526969790458679), ("Cu-EPIC-2kV", 0.8962102770805359), ("Cu-Ni-15kV", 0.9759411454200745), ("Mn-Cr-12kV", 0.9700609385967255), ("Ti-Ti-9kV", 0.9632380992174149)] Key: (244, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9737197607755661), ("Al-Al-4kV", 0.936311411857605), ("C-EPIC-0.6kV", 0.7668724238872529), ("Cu-EPIC-0.9kV", 0.8007214874029159), ("Cu-EPIC-2kV", 0.8769152045249939), ("Cu-Ni-15kV", 0.9682861983776092), ("Mn-Cr-12kV", 0.9659816890954971), ("Ti-Ti-9kV", 0.9604099780321121)] Key: (287, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9702932089567184), ("Al-Al-4kV", 0.9345487594604492), ("C-EPIC-0.6kV", 0.7561751186847686), ("Cu-EPIC-0.9kV", 0.7987294286489487), ("Cu-EPIC-2kV", 0.8735345602035522), ("Cu-Ni-15kV", 0.9701228439807892), ("Mn-Cr-12kV", 0.9675930917263031), ("Ti-Ti-9kV", 0.9575580269098282)] Key: (278, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9720843970775604), ("Al-Al-4kV", 0.9305856913328171), ("C-EPIC-0.6kV", 0.741540789604187), ("Cu-EPIC-0.9kV", 0.7682088732719422), ("Cu-EPIC-2kV", 0.864213228225708), ("Cu-Ni-15kV", 0.9648653626441955), ("Mn-Cr-12kV", 0.9665578484535218), ("Ti-Ti-9kV", 0.9552424371242523)] Key: (250, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9766805320978165), ("Al-Al-4kV", 0.9418972045183182), ("C-EPIC-0.6kV", 0.76584292948246), ("Cu-EPIC-0.9kV", 0.8238443195819855), ("Cu-EPIC-2kV", 0.8775055557489395), ("Cu-Ni-15kV", 0.9738869041204452), ("Mn-Cr-12kV", 0.9712955445051193), ("Ti-Ti-9kV", 0.9667880535125732)] Key: (283, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9675714403390885), ("Al-Al-4kV", 0.9285876244306565), ("C-EPIC-0.6kV", 0.71473089158535), ("Cu-EPIC-0.9kV", 0.7696082562208175), ("Cu-EPIC-2kV", 0.8377785980701447), ("Cu-Ni-15kV", 0.9635965615510941), ("Mn-Cr-12kV", 0.9641305923461914), ("Ti-Ti-9kV", 0.9575453609228134)] Key: (301, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9705011487007141), ("Al-Al-4kV", 0.9279218137264251), ("C-EPIC-0.6kV", 0.7229891419410706), ("Cu-EPIC-0.9kV", 0.793598598241806), ("Cu-EPIC-2kV", 0.8582320868968963), ("Cu-Ni-15kV", 0.9688272416591645), ("Mn-Cr-12kV", 0.9668716788291931), ("Ti-Ti-9kV", 0.9546888172626495)] Key: (274, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.974557614326477), ("Al-Al-4kV", 0.9371200978755951), ("C-EPIC-0.6kV", 0.7444291532039642), ("Cu-EPIC-0.9kV", 0.7895265400409699), ("Cu-EPIC-2kV", 0.8598116517066956), ("Cu-Ni-15kV", 0.9712087035179138), ("Mn-Cr-12kV", 0.9688791006803512), ("Ti-Ti-9kV", 0.9589674681425094)] Key: (242, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9761280864477158), ("Al-Al-4kV", 0.9446888864040375), ("C-EPIC-0.6kV", 0.766591414809227), ("Cu-EPIC-0.9kV", 0.8117899149656296), ("Cu-EPIC-2kV", 0.8900630325078964), ("Cu-Ni-15kV", 0.971353754401207), ("Mn-Cr-12kV", 0.9718088060617447), ("Ti-Ti-9kV", 0.9645727574825287)] Key: (306, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9756848812103271), ("Al-Al-4kV", 0.9455605119466781), ("C-EPIC-0.6kV", 0.7938183635473252), ("Cu-EPIC-0.9kV", 0.8287457168102265), ("Cu-EPIC-2kV", 0.8792453199625015), ("Cu-Ni-15kV", 0.9696165889501571), ("Mn-Cr-12kV", 0.972235518693924), ("Ti-Ti-9kV", 0.9627663731575012)] Key: (303, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9724529892206192), ("Al-Al-4kV", 0.9305596590042114), ("C-EPIC-0.6kV", 0.7230549484491349), ("Cu-EPIC-0.9kV", 0.7953008472919464), ("Cu-EPIC-2kV", 0.8613291561603547), ("Cu-Ni-15kV", 0.9646958172321319), ("Mn-Cr-12kV", 0.9663623839616775), ("Ti-Ti-9kV", 0.9565762877464294)] Key: (291, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9746044754981995), ("Al-Al-4kV", 0.945956015586853), ("C-EPIC-0.6kV", 0.7661843031644822), ("Cu-EPIC-0.9kV", 0.8199316382408142), ("Cu-EPIC-2kV", 0.8820369154214859), ("Cu-Ni-15kV", 0.972326734662056), ("Mn-Cr-12kV", 0.9720319092273713), ("Ti-Ti-9kV", 0.9641989678144455)] Key: (281, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9678589403629303), ("Al-Al-4kV", 0.9328784346580505), ("C-EPIC-0.6kV", 0.7547005414962769), ("Cu-EPIC-0.9kV", 0.7789339125156403), ("Cu-EPIC-2kV", 0.8745017945766449), ("Cu-Ni-15kV", 0.9656407535076141), ("Mn-Cr-12kV", 0.9641033113002777), ("Ti-Ti-9kV", 0.9552679359912872)] Key: (276, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9702300161123276), ("Al-Al-4kV", 0.9325366169214249), ("C-EPIC-0.6kV", 0.7419947564601898), ("Cu-EPIC-0.9kV", 0.7747547417879105), ("Cu-EPIC-2kV", 0.8381759762763977), ("Cu-Ni-15kV", 0.9623401463031769), ("Mn-Cr-12kV", 0.9613375902175904), ("Ti-Ti-9kV", 0.9547658741474152)] Key: (279, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9683325439691544), ("Al-Al-4kV", 0.9248780459165573), ("C-EPIC-0.6kV", 0.7393383264541626), ("Cu-EPIC-0.9kV", 0.7804951041936874), ("Cu-EPIC-2kV", 0.8707629531621933), ("Cu-Ni-15kV", 0.9666939914226532), ("Mn-Cr-12kV", 0.9646541595458984), ("Ti-Ti-9kV", 0.9589365422725677)] Key: (268, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9719461470842361), ("Al-Al-4kV", 0.9368112564086915), ("C-EPIC-0.6kV", 0.7573535948991775), ("Cu-EPIC-0.9kV", 0.8275747984647751), ("Cu-EPIC-2kV", 0.8748720288276672), ("Cu-Ni-15kV", 0.9710170537233352), ("Mn-Cr-12kV", 0.9695031344890594), ("Ti-Ti-9kV", 0.9605758100748062)] Key: (258, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9745412111282349), ("Al-Al-4kV", 0.9339351028203964), ("C-EPIC-0.6kV", 0.7520094603300095), ("Cu-EPIC-0.9kV", 0.8004794657230377), ("Cu-EPIC-2kV", 0.8735634952783584), ("Cu-Ni-15kV", 0.9667491674423218), ("Mn-Cr-12kV", 0.9659819602966309), ("Ti-Ti-9kV", 0.9595000624656678)] Key: (248, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.977547001838684), ("Al-Al-4kV", 0.943942391872406), ("C-EPIC-0.6kV", 0.7833251833915711), ("Cu-EPIC-0.9kV", 0.8251269578933715), ("Cu-EPIC-2kV", 0.8876170873641968), ("Cu-Ni-15kV", 0.9738454699516297), ("Mn-Cr-12kV", 0.9718320488929748), ("Ti-Ti-9kV", 0.9669605851173401)] Key: (246, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9743694633245468), ("Al-Al-4kV", 0.9456877648830414), ("C-EPIC-0.6kV", 0.7678399056196212), ("Cu-EPIC-0.9kV", 0.8119639933109284), ("Cu-EPIC-2kV", 0.8899865686893463), ("Cu-Ni-15kV", 0.9739743441343307), ("Mn-Cr-12kV", 0.9730724036693573), ("Ti-Ti-9kV", 0.966443908214569)] Key: (254, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9705976217985153), ("Al-Al-4kV", 0.9223994761705399), ("C-EPIC-0.6kV", 0.7517837852239608), ("Cu-EPIC-0.9kV", 0.7894137501716614), ("Cu-EPIC-2kV", 0.8740812391042709), ("Cu-Ni-15kV", 0.9649749875068665), ("Mn-Cr-12kV", 0.9622029483318328), ("Ti-Ti-9kV", 0.9560768663883209)] Key: (299, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9735205113887787), ("Al-Al-4kV", 0.9309952080249786), ("C-EPIC-0.6kV", 0.7400766223669052), ("Cu-EPIC-0.9kV", 0.7843391090631485), ("Cu-EPIC-2kV", 0.8648605406284332), ("Cu-Ni-15kV", 0.967012819647789), ("Mn-Cr-12kV", 0.9662521809339524), ("Ti-Ti-9kV", 0.9518016576766968)] Key: (256, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9687688469886779), ("Al-Al-4kV", 0.9329577207565307), ("C-EPIC-0.6kV", 0.7124274164438248), ("Cu-EPIC-0.9kV", 0.7790258109569549), ("Cu-EPIC-2kV", 0.8657959043979645), ("Cu-Ni-15kV", 0.9670290321111679), ("Mn-Cr-12kV", 0.9633038669824601), ("Ti-Ti-9kV", 0.9567578852176666)] Key: (293, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9724262475967407), ("Al-Al-4kV", 0.9334424823522568), ("C-EPIC-0.6kV", 0.7237139195203781), ("Cu-EPIC-0.9kV", 0.8003135979175567), ("Cu-EPIC-2kV", 0.8604774057865143), ("Cu-Ni-15kV", 0.9647800117731095), ("Mn-Cr-12kV", 0.9671240687370301), ("Ti-Ti-9kV", 0.9593446969985961)] Key: (261, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.971569812297821), ("Al-Al-4kV", 0.9378552496433258), ("C-EPIC-0.6kV", 0.7479894399642945), ("Cu-EPIC-0.9kV", 0.8020979523658752), ("Cu-EPIC-2kV", 0.8605112731456757), ("Cu-Ni-15kV", 0.9698976039886474), ("Mn-Cr-12kV", 0.9648890495300293), ("Ti-Ti-9kV", 0.9587083220481872)] Key: (289, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9720789223909378), ("Al-Al-4kV", 0.9372004061937332), ("C-EPIC-0.6kV", 0.7466429948806763), ("Cu-EPIC-0.9kV", 0.7778348356485367), ("Cu-EPIC-2kV", 0.8603413850069046), ("Cu-Ni-15kV", 0.967164334654808), ("Mn-Cr-12kV", 0.9689937591552734), ("Ti-Ti-9kV", 0.9565973937511444)] Key: (298, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9661869823932647), ("Al-Al-4kV", 0.9258407413959503), ("C-EPIC-0.6kV", 0.7028941214084625), ("Cu-EPIC-0.9kV", 0.7626788705587387), ("Cu-EPIC-2kV", 0.8486940711736679), ("Cu-Ni-15kV", 0.9655997604131699), ("Mn-Cr-12kV", 0.9632629603147507), ("Ti-Ti-9kV", 0.9513149082660675)] Key: (265, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9717926532030106), ("Al-Al-4kV", 0.9275272369384766), ("C-EPIC-0.6kV", 0.7334349215030671), ("Cu-EPIC-0.9kV", 0.7848327666521072), ("Cu-EPIC-2kV", 0.879998528957367), ("Cu-Ni-15kV", 0.9686383992433548), ("Mn-Cr-12kV", 0.9665170550346375), ("Ti-Ti-9kV", 0.9604606479406357)] Key: (272, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9744620800018311), ("Al-Al-4kV", 0.9406907320022583), ("C-EPIC-0.6kV", 0.7583581209182739), ("Cu-EPIC-0.9kV", 0.803549861907959), ("Cu-EPIC-2kV", 0.8741284489631653), ("Cu-Ni-15kV", 0.9701647937297821), ("Mn-Cr-12kV", 0.9710357129573822), ("Ti-Ti-9kV", 0.9592038929462433)] Key: (297, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.970538979768753), ("Al-Al-4kV", 0.9283351331949234), ("C-EPIC-0.6kV", 0.7336910545825959), ("Cu-EPIC-0.9kV", 0.7939026236534119), ("Cu-EPIC-2kV", 0.8525126844644546), ("Cu-Ni-15kV", 0.9671412736177445), ("Mn-Cr-12kV", 0.9655984342098236), ("Ti-Ti-9kV", 0.9555086642503738)] Key: (263, "D7DBC196401F3CAC564FBE899306343BEB5022BA", 0.95) = @[("Ag-Ag-6kV", 0.9716018617153168), ("Al-Al-4kV", 0.9364291220903397), ("C-EPIC-0.6kV", 0.7137390047311782), ("Cu-EPIC-0.9kV", 0.7985885977745056), ("Cu-EPIC-2kV", 0.8676451265811921), ("Cu-Ni-15kV", 0.9684803396463394), ("Mn-Cr-12kV", 0.9639465093612671), ("Ti-Ti-9kV", 0.9574855715036392)]

And finally rerun all likelihood combinations:

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crAll --regions crGold \
    --signalEfficiency 0.8 --signalEfficiency 0.85 --signalEfficiency 0.9 --signalEfficiency 0.95 \
    --fadcVetoPercentile 0.99 \
    --vetoSets "{fkMLP, fkFadc, fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --mlpPath ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
    --out ~/org/resources/lhood_limits_23_04_23_mlp/ \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing

Running the command was somewhat of a failure:

• the default of 8 jobs per process is not enough. We use way more memory than expected. Even with only 4 jobs, there is a risk of one or more being killed!

• sometimes (at least in crAll combinations) we still get the "broken event" quit error. In this case however we can spy the following:

  Cluster: 1 of chip 3 has val : 0.9999758005142212 copmare: 0.9467523097991943
  Cluster: 0 of chip 5 has val : 0.02295641601085663 copmare: 0.9783756822347641
  Cluster: 1 of chip 3 has val : 0.9999998807907104 copmare: 0.9467523097991943
  Cluster: 0 of chip 4 has val : 5.220065759203862e-06 copmare: 0.99133480489254
  Cluster: 1 of chip 3 has val : 0.9999876022338867 copmare: 0.9467523097991943
  Cluster: 0 of chip 3 has val : nan copmare: 0.9545457273721695
  Broken event! DataFrame with 3 columns and 1 rows:
       Idx     eventIndex    eventNumber     chipNumber
    dtype:            int            int          float
         0          10750          35811              3
  

  -> Note the nan value for cluster 0 in the last line before the DF print! This means the MLP predicted a value of nan for a cluster!

-> We need to debug the likelihood call for one of the cases and isolate the events for that.

• [X]

  get a combination that causes this:

    lhood_c18_R2_crAll_sEff_0.85_scinti_fadc_septem_line_mlp_trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss_vQ_0.99.log
  

• [X] extract which run is at fault -> Run 92 in this case
• [X] Does nan appear more often? -> Yes! O(5) times in the same log file
• [ ]

  Rerun the equivalent command for testing and try to debug the cause. -> Reconstruct the command:

  likelihood \
      -f /home/basti/CastData/data/DataRuns2017_Reco.h5 \
      --h5out /t/debug_broken_event.h5 \
      --region=crAll \
      --cdlYear=2018 \
      --scintiveto \
      --fadcveto \
      --septemveto \
      --lineveto \
      --mlp /home/basti/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
      --cdlFile=/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
      --readOnly \
      --calibFile=/home/basti/CastData/data/CalibrationRuns2017_Reco.h5 \
      --vetoPercentile=0.99 \
      --nnSignalEff=0.85
  

• [X] Modify createAllLikelihoodCombinations such that each job does not stop on a failure, but retries failed jobs -> Done
• [ ] I'll restart the create likelihood combinations command now that we've implemented restart on failure & disabled failing on NaN septem events (i.e. "no cluster found") -> FIX ME

• [X] We've also implemented a version of toH5 and deserializeH5 that

retries the writing / reading if the file is locked. This should make it safer to run multiple instances in parallel which might try to access the same file.

• [X] Look into the origin of events with a NaN value! -> Likely they just pass by accident due to having NaN? That would be "good"
• [X] look at what such clusters look like

Having inserted printing the event number for run 92 with the NaN event, the output is:

Broken event! DataFrame with 3 columns and 1 rows:
     Idx     eventIndex    eventNumber     chipNumber
  dtype:            int            int          float
       0          10750          35811              3

The event is: @[("eventNumber", (kind: VInt, num: 35811))]

Let's plot that event:

plotData \
    --h5file ~/CastData/data/DataRuns2017_Reco.h5 \
    --runType rtBackground \
    --config ~/CastData/ExternCode/TimepixAnalysis/Plotting/karaPlot/config.toml \
    --runs 92 \
    --eventDisplay \
    --septemboard \
    --events 35811

The event is:

Well, even some of the properties are NaN! No wonder it yields a NaN result for the NN prediction. We have to be careful not to accidentally consider them "passing" though! (Which seems to be the case currently)

We've changed the code to always veto NaN events. And as a result the broken event should never happen again, which is why we reintroduced the quit condition. Let's see if the command from yesterday now runs correctly.

Well, great. The process was killed due to memory usage of running 5 jobs concurrently. Guess I'll just wait for all the other jobs to finish now. We're ~half way done or so.

Small problem: In the limit calculation given that our effective efficiency is Run-2/3 specific we need to think about how to handle that. Can we use different efficiencies for different parts of the data? Certainly, but makes everything more complicated. -> Better to compute an average of the two different runs.

Sigh, it seems like even 2 jobs at the same time can use too much memory!

(oh, what a dummy: I didn't even add the 99% case to the createLogL call!)

Therefore:

For now let's try to see what we can get with exactly one setup:

• 99% MLP
• all vetoes except septem veto

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crAll --regions crGold \
    --signalEfficiency 0.99 \
    --fadcVetoPercentile 0.99 \
    --vetoSets "{fkMLP, +fkFadc, +fkScinti, fkLineVeto}" \
    --mlpPath ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/trained_mlp_sgd_gauss_diffusion_tanh300_mse_loss.pt \
    --out ~/org/resources/lhood_limits_23_04_23_mlp_sEff_0.99/ \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    #--multiprocessing

We're running without multiprocessing now. :( At least for the crAll cases. We'll do crGold with 2 jobs first.

Gold is finished and now running crAll with a single job. :/

So: Things we want for the meeting:

1. background rate in gold region
2. background clusters of MLP@99
3. background clusters of MLP@99 + vetoes
4. expected limits for MLP@99 & MLP@99+vetoes

UPDATE: <2023-10-30 Mon 11:14> We changed the train_ingrid code slightly today, so that our previous hardcoded change to replace subsetPerRun by 6 * subsetPerRun for background data was taken out. That's why I modified the command below to accommodate for that by adding the --subsetPerRun 6000 argument!

• the fact that mcmc_limit gets stuck on toH5 is because of the KDtree, which is massive with O(350 k) entries!
• [X] add filter option (energyMin and energyMax) -> shows that in Run-2 alone without any energy filtering & without vetoes there are almost 200k cluster! -> cutting 0.2 < E < 12 and using vetoes cuts it in half O(40k) compared to O(75k) when 0 < E < 12.
• [X] Check data selection for background training sample -> OUCH: We still filter to THE GOLD REGION in the background sample. This likely explains why the background is so good there, but horrible towards the edges!
• [X] write a custom serializer for KDTree to avoid extremely nested H5 data structure -> We now only serialize the actual data of the tree. The tree can be rebuilt from that after all.
• [X] implement filtering on energy < 200 eV for mcmc_limit
  • [X] also check the impact of that on the limit! -> running right now -> very slight improvement to g_ae² = 6.28628480082639e-21 from g_ae² = 6.333984435685045e-21
• [ ] analyze memory usage of likelihood when using NN veto
• [X] Train a new MLP of same architecture as currently, but using background data over the whole chip! -> [X] CenterX or Y is not part of inputs, right? -> nope, it's not

Here we go:

./train_ingrid \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --back ~/CastData/data/DataRuns2017_Reco.h5 \
    --back ~/CastData/data/DataRuns2018_Reco.h5 \
    --modelOutpath ~/org/resources/nn_devel_mixing/10_05_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_tanh300_mse.pt \
    --plotPath ~/Sync/10_05_23_sgd_tanh300_mse/ \
    --datasets eccentricity \
    --datasets skewnessLongitudinal \
    --datasets skewnessTransverse \
    --datasets kurtosisLongitudinal \
    --datasets kurtosisTransverse \
    --datasets length \
    --datasets width \
    --datasets rmsLongitudinal \
    --datasets rmsTransverse \
    --datasets lengthDivRmsTrans \
    --datasets rotationAngle \
    --datasets fractionInTransverseRms \
    --datasets totalCharge \
    --datasets σT \
    --numHidden 300 \
    --numHidden 300 \
    --activation tanh \
    --outputActivation sigmoid \
    --lossFunction MSE \
    --optimizer SGD \
    --learningRate 7e-4 \
    --simulatedData \
    --backgroundRegion crAll \
    --nFake 250_000 \
    --subsetPerRun 6000

-> Trained up to 500k on <2023-05-11 Thu 10:21>.

Continue from yesterday:

• [ ] analyze memory usage of likelihood when using NN veto -> try with nimprof -> useless

Regenerating the cache tables:

./determineDiffusion \
    ~/CastData/data/DataRuns2017_Reco.h5 \
    ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    ~/CastData/data/DataRuns2018_Reco.h5 \
    ~/CastData/data/CalibrationRuns2018_Reco.h5

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crAll --regions crGold \
    --signalEfficiency 0.99 \
    --fadcVetoPercentile 0.99 \
    --vetoSets "{fkMLP, +fkFadc, +fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --mlpPath ~/org/resources/nn_devel_mixing/10_05_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.pt \
    --out ~/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/ \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5

With these plots generated and the corresponding sections written, we are almost done with that part of the thesis.

• [ ] Mini section about background rate with MLP
• [ ] Extend MLP w/ all vetoes to include MLP

• [X]

  Let's see if running likelihood with lnL instead of the MLP also eats more and more memory!

  likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out /t/dm_lnl_slowit.h5 --region=crAll --cdlYear=2018 --scintiveto --fadcveto --septemveto --lineveto --lnL --cdlFile=/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 --readOnly --calibFile=/home/basti/CastData/data/CalibrationRuns2017_Reco.h5 --vetoPercentile=0.99 --signalEff=0.8
  

  -> Not really. It stays nice and low.

Running the following script:

import nimhdf5, os, seqmath, sequtils, datamancer

import ingrid / ingrid_types

proc readAllChipData(h5f: H5File, group: H5Group, numChips: int): AllChipData =
  ## Read all data for all chips of this run that we need for the septem veto
  let vlenXY = special_type(uint8)
  let vlenCh = special_type(float64)
  result = AllChipData(x: newSeq[seq[seq[uint8]]](numChips),
                       y: newSeq[seq[seq[uint8]]](numChips),
                       ToT: newSeq[seq[seq[uint16]]](numChips),
                       charge: newSeq[seq[seq[float]]](numChips))
  for i in 0 ..< numChips:
    result.x[i] =      h5f[group.name / "chip_" & $i / "x", vlenXY, uint8]
    result.y[i] =      h5f[group.name / "chip_" & $i / "y", vlenXY, uint8]
    result.ToT[i] =    h5f[group.name / "chip_" & $i / "ToT", vlenCh, uint16]
    result.charge[i] = h5f[group.name / "chip_" & $i / "charge", vlenCh, float]


var h5f = H5open("~/CastData/data/DataRuns2017_Reco.h5", "r")
let grp = h5f["/reconstruction/run_186/".grp_str]
echo grp
var df = newDataFrame()
for i in 0 ..< 10:
  let data = readAllChipData(h5f, grp, 7)
  df = toDf({"x" : data.x.flatten.mapIt(it.float)})
  echo "Read: ", i, " = ", getOccupiedMem().float / 1e6, " MB", " df len: ", df.len
discard h5f.close()

under valgrind right now. -> In this setup valgrind does not see any leaks.

Let's also run it under heaptrack to see where the 8GB of memory come from! -> Running without -d:useMalloc yields pretty much nothing (as it seems to intercept malloc calls) -> With -d:useMalloc it doesn't look unusual at all -> the 8GB seen when running without -d:useMalloc seem to be the standard Nim allocator doing its thing -> So this code snippet seems fine.

• [ ] let's try to trim down likelihood to still reproduce the memleak issue

  1. see the memory usage of a cpp MLP run without vetoes -> so far seems to run without growing endlessely. -> it grows slightly, approaching 9GB Maybe we can run this under heaptrack? -> Yeah, no problems without septem&line veto -> Now trying with only septem veto -> This already crashes heaptrack! -> Let's try likelihood with lnL

          heaptrack likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out /t/more_lnl_slowit.h5 --region=crAll --cdlYear=2018 --scintiveto --fadcveto --septemveto --lnL --cdlFile=/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 --readOnly --calibFile=/home/basti/CastData/data/CalibrationRuns2017_Reco.h5 --vetoPercentile=0.99 --signalEff=0.80
     

     • [X] Even with lnL and septem veto heaptrack crashes! Let's see if we understand where.
• [ ]

  Debug the crashing of heaptrack -> cause is in septem veto -> in geometry -> in DBScan reconstruction -> in nearestNeighbor call

  • [X] serialized a bunch of pixels that cause the crash in /tmp/serial_pixels.h5

  -> let's only run DBSCAN on these pixels!

  import arraymancer
  
  
  proc callDb(data: Tensor[float]) =
    echo dbscan(data, 65.0, 3)
  
  import nimhdf5
  type
    Foo = object
      shape: seq[int]
      data: seq[float]
  
  let foo = deserializeH5[Foo]("/tmp/serial_pixels.h5")
  var pT = foo.data.toTensor.reshape(foo.shape)
  
  for _ in 0 ..< 1000:
    callDb(pT)
  

  -> cannot reproduce the problem this way

• [X] Ok, given that heaptrack always seems to crash when using DBSCAN in that context, I just ran it using the default cluster algo for the septem veto logic. This ran correctly and showed no real memory leak. And a peak heap mem size of ~6 GB.
• [X] Running again now using default clusterer but using MLP! -> Reproduces the problem and heaptrack tracks the issue. We're leaking in the H5 library due to identifiers not being closed. I've changed the code now to automatically close identifiers by attaching them to =destroy calls. -> Also changed the logic to only read the MLP H5 file a single time for the septem veto, because this was the main origin (reading the file for each single cluster!)

UPDATE: -> We've replaced the distinct hid_t logic in nimhdf5 by an approach that wraps the identifiers in a ref object to make sure we destroy every single identifiers when it goes out of scope. This fixed the memory leak that could finally be tracked properly with the following snippet:

import nimhdf5

type
  MLPDesc* = object
    version*: int # Version of this MLPDesc object
    path*: string # model path to the checkpoint files including the default model name!
    modelDir*: string # the parent directory of `path`
    plotPath*: string # path in which plots are placed
    calibFiles*: seq[string] ## Path to the calibration files
    backFiles*: seq[string] ## Path to the background data files
    simulatedData*: bool
    numInputs*: int
    numHidden*: seq[int]
    numLayers*: int
    learningRate*: float
    datasets*: seq[string] # Not `InGridDsetKind` to support arbitrary new columns
    subsetPerRun*: int
    rngSeed*: int
    backgroundRegion*: string
    nFake*: int # number of fake events per run period to generate
    #
    activationFunction*: string
    outputActivation*: string
    lossFunction*: string
    optimizer*: string
    # fields that store training information
    epochs*: seq[int] ## epochs at which plots and checkpoints are generated
    accuracies*: seq[float]
    testAccuracies*: seq[float]
    losses*: seq[float]
    testLosses*: seq[float]


proc getNumber(file: string): int =
  let desc = deserializeH5[MLPDesc](file)
  result = desc.numInputs

proc main(fname: string) =
  for i in 0 ..< 50_000:
    echo "Number: ", i, " gives ", getNumber(fname)

when isMainModule:
  import cligen
  dispatch main

running via:

heaptrack ./testmem3 -f ~/org/resources/nn_devel_mixing/23_04_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_desc.h5

The issue really was the many deserializeH5 calls, which suddenly highlighted the memory leak!

In the meantime let's look at background clusters of the likelihood outputs from last night:

• [ ]

  Try again to rerun heaptrack on the DBSCAN code:

  heaptrack likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out /t/mlp_isit_faster_noleak_dbscan.h5 --region=crAll --cdlYear=2018 --scintiveto --fadcveto --septemveto --mlp ~/org/resources/nn_devel_mixing/10_05_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.pt --cdlFile=/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 --readOnly --calibFile=/home/basti/CastData/data/CalibrationRuns2017_Reco.h5 --vetoPercentile=0.99 --nnSignalEff=0.99
  

  -> Still the same problem! So either a bug in heaptrack or a bug in our code still. :/

• [ ]

  Run perf on the code using default clustering algo. Then use hotspot to check the report:

  perf record -o /t/log_noleaks_slow.data --call-graph dwarf -- likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out /t/log_noleaks_debug_slowit.
  h5 --region=crAll --cdlYear=2018 --scintiveto --fadcveto --septemveto --lineveto --mlp ~/org/resources/nn_devel_mixing/10_05_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.pt --cdlFile
  =/home/basti/CastData/data/CDL_2019/calibration-cdl-2018.h5 --readOnly --calibFile=/home/basti/CastData/data/CalibrationRuns2017_Reco.h5 --vetoPercentile=0.99 --nnSignalEff=0.99
  

  Killed it after about 20 runs processed as the file already grew to over 15GB. Should be enough statistics… :) -> Ok, the perf data shows that time is spent in a variety of different places. mergeChain of the cluster algo is a bigger one, so is forward of the MLP and generally a whole lot of copying data. Nothing in particular jumps out though.

  I guess performance for this is acceptable for now.

  Let's also look at DBSCAN (same command as above): -> Yeah, as expected. The queryImpl call is by far the dominating part of the performance report when using DBSCAN. The heapqueues used make up >50% of the time with the distance and ~index_select and toTensorTuple the rest. In toTensorTuple the dominating factor is pop.

  -> We've replaced the HeapQueue by a SortedSeq now for testing. With it the pop procedure is much faster (but insert is a bit slower). We've finally added a --run option to likelihood to now time the performance of the heapqueue vs the sorted seq. We'll run on run 168 as it is one of the longest runs. Sorted seq:

    likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out           314.52s user 3.16s system 100%
  

  Heap queue:

    likelihood -f /home/basti/CastData/data/DataRuns2017_Reco.h5 --h5out           324.78s user 3.12s system 100%
  

  -> Well. Guess it is a bit faster than the heap queue approach, but it's a small improvement. I guess we mostly traded building the sorted seq by popping from it. NOTE: There was a segfault after the run was done?

• [X] just ran heaptrack again on the default clustering case for the full case, because I didn't remember if we did that with the ref object IDs. All looking good now. Peak memory at 9.5GB. High, but explained by the 3.7GB overhead of the CUDA (that is "leaked"). Oh, I just googled for the CUDA "leak": https://discuss.pytorch.org/t/memory-leak-in-libtorch-extremely-simple-code/38149/3 -> It's because we don't use "no grad" mode! Hmm, but we *are* using our ~no_grad_mode template in nn_types and effectively in nn_predict. -> I introduced the "NoGradGuard" into Flamebau and used it in places in addition to the no_grad_mode and run heaptrack again. Let's see. -> Didn't change anything!

So I guess that means we continue with our actual work. Performance is deemed acceptable now.

Let's go with a test run of 5 different jobs in createAllLikelihoodCombinations:

./createAllLikelihoodCombinations \
    --f2017 ~/CastData/data/DataRuns2017_Reco.h5 \
    --f2018 ~/CastData/data/DataRuns2018_Reco.h5 \
    --c2017 ~/CastData/data/CalibrationRuns2017_Reco.h5 \
    --c2018 ~/CastData/data/CalibrationRuns2018_Reco.h5 \
    --regions crAll --regions crGold \
    --signalEfficiency 0.95 --signalEfficiency 0.90 \
    --fadcVetoPercentile 0.99 \
    --vetoSets "{fkMLP, +fkFadc, +fkScinti, fkSeptem, fkLineVeto, fkExclusiveLineVeto}" \
    --mlpPath ~/org/resources/nn_devel_mixing/10_05_23_sgd_gauss_diffusion_tanh300_mse_loss/mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.pt \
    --out ~/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/ \
    --cdlFile ~/CastData/data/CDL_2019/calibration-cdl-2018.h5 \
    --multiprocessing \
    --jobs 5

• [X] We get errors for dataspaces. Apparently our nimhdf5 code is problematic in some cases. Let's try to fix that up, run all its tests. -> tread_write reproduces it -> fixed all the issues by another rewrite of the ID logic

NOTE: I'm stopping the following limit calculation at <2023-05-16 Tue 08:10>:

shellCmd: mcmc_limit_calculation limit -f /home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/lhood_c18_R2_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5 -f /home/basti/org/resources/lhood_li
mits_10_05_23_mlp_sEff_0.99/lhood_c18_R3_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5 --years 2017 --years 2018 --σ_p 0.05 --limitKind lkMCMC --nmc 1000 --suffix=_sEff_0.99_mlp_mlp_tanh300_msecheckpo
int_epoch_485000_loss_0.0055_acc_0.9933 --path "" --outpath /home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/limits/
shell 5127> files @["/home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/lhood_c18_R2_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5", "/home/basti/org/resources/lhood_limits_10_05_23_mlp_sEf
f_0.99/lhood_c18_R3_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5"]
shell 5127>
shell 5127> @["/home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/lhood_c18_R2_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5", "/home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99
/lhood_c18_R3_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5"]
shell 5127> [INFO]: Read a total of 340911 input clusters.

because each MCMC takes too long to build:

shell 5127> MC index 12
shell 5127> Building chain of 150000 elements took 126.9490270614624 s
shell 5127> Acceptance rate: 0.2949866666666667 with last two states of chain: @[@[1.250804325635904e-21, 0.02908566421900792, -0.004614450435544698, 0.03164360408307887, 0.03619049049344585], @[1.250804325635904e-21, 0.02908566421900792,
-0.004614450435544698, 0.03164360408307887, 0.03619049049344585]]
shell 5127> Limit at 7.666033305250432e-21
shell 5127> Number of candidates: 16133
shell 5127> INFO: The integer column `Hist` has been automatically determined to be continuous. To overwrite this behavior add a `+ scale_x/y_discrete()` call to the plotting chain. Choose `x` or `y` depending on which axis this column refers to. Or apply a `factor` to the column name in the `aes` call, i.e. `aes(..., factor("Hist"), ...)`.
shell 5127> MC index 12
shell 5127> Building chain of 150000 elements took 133.913836479187 s
shell 5127> Acceptance rate: 0.3029 with last two states of chain: @[@[5.394302558724233e-21, 0.005117702932866319, -0.003854151251413666, 0.1140589851066757, -0.01650063836525805], @[5.394302558724233e-21, 0.005117702932866319, -0.003854151251413666, 0.1140589851066757, -0.01650063836525805]]
shell 5127> Limit at 8.815294304890919e-21
shell 5127> Number of candidates: 16518
shell 5127> INFO: The integer column `Hist` has been automatically determined to be continuous. To overwrite this behavior add a `+ scale_x/y_discrete()` call to the plotting chain. Choose `x` or `y` depending on which axis this column refers to. Or apply a `factor` to the column name in the `aes` call, i.e. `aes(..., factor("Hist"), ...)`.
shell 5127> MC index 12

(due to the ~16k candidates)

I'll add the file to processed.txt and start the rest now.

shell 5127>             Initial chain state: @[3.325031213438127e-21, -0.005975150812670194, 0.2566543411242529, 0.1308918272537833, 0.3838098582402962]
^C
Interrupted while running processing of file: /home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/lhood_c18_R2_crAll_sEff_0.99_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5
Check the `processed.txt` file in the output path to see which files were processed successfully!

Added that file manually to the processed.txt file now. Restarting:

./runLimits --path ~/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/ --outpath ~/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/limits/ --prefix lhood_c18_R2_crAll --nmc 1000

Doing the same with 90% MLP:

^C
Interrupted while running processing of file: /home/basti/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/lhood_c18_R2_crAll_sEff_0.9_mlp_mlp_tanh300_msecheckpoint_epoch_485000_loss_0.0055_acc_0.9933.h5
Check the `processed.txt` file in the output path to see which files were processed successfully!

Restarted with the veto based ones for 90% left.

Running the expected limits table generator:

• [X] I updated it to include the MLP efficiencies. Still have to change it so that it prints whether MLP or LnL was in use!

./generateExpectedLimitsTable \
    -p ~/org/resources/lhood_limits_10_05_23_mlp_sEff_0.99/limits/ \
    --prefix mc_limit_lkMCMC_skInterpBackground_nmc_1000

yields:

Continue on from meeting with Klaus on <2023-05-16 Tue 13:00>:

• [ ]

  Start by fixing the systematics -> start by extracting the systematic values used to compute the current number. The default values are:

      σ_sig = 0.04692492913207222,
      σ_back = 0.002821014576353691,
      σ_p = 0.05,
  

  computed from section sec:systematics:combined_uncertainties in StatusAndProgress via the following code (stripped down to signal):

  import math, sequtils
  let ss = [3.3456,
            0.5807,
            1.0,
            0.2159,
            2.32558,
            0.18521]
            #2.0] #1.727]
  ##          ^-- This is the effective efficiency for 55Fe apparently.
  proc total(vals: openArray[float]): float =
    for x in vals:
      result += x * x
    result = sqrt(result)
  
  let ss0 = total(ss)
  let ss17 = total(concat(@ss, @[1.727]))
  let ss2 = total(concat(@ss, @[2.0]))
    
  echo "Combined uncertainty signal (Δ software eff = 0%): ",     ss0 / 100.0  
  echo "Combined uncertainty signal (Δ software eff = 1.727%): ", ss17 / 100.0
  echo "Combined uncertainty signal (Δ software eff = 2%): ",     ss2 / 100.0
  

  • [X] There is one mystery here: The value that comes out of that

  calculation is 0.458 instead of the ~0.469 used in the code. I don't know why that is exactly. -> *SOLVED*: The reason is the ~0.469 is from assuming 2%, but the ~0.458 is from using the value ~1.727!

  So now all we need to do is to combine the value without the software efficiency numbers with the effective efficiency from the MLP! (or using the 2% for LnL).

  We do this by:

  import math
  let ss = [3.3456,
            0.5807,
            1.0,
            0.2159,
            2.32558,
            0.18521]
  proc total(vals: openArray[float]): float =
    for xf in vals:
      let x = xf / 100.0
      result += x * x
    result = sqrt(result)
  
  let sStart = total(ss) # / 100.0
  doAssert abs(sStart - 0.04244936953654317) < 1e-5
  
  # add new value:
  let seff = 1.727 / 100.0
  let s17 = sqrt( (sStart)^2 + seff^2 )
  echo "Uncertainty including Δseff = 1.727% after the fact: ", s17
  

  • [X] Implemented

• [ ]

  IMPORTANT: When using MLP in classification, how do we deal with inputs that have NaN values??? Do we filter them out, i.e. reject them?

    let data = h5f.readValidDsets(grp, filterNan = false)
  

  -> How do we deal with them? Ah, from likelihood:

               (classify(nnPred[ind]) == fcNaN or # some clusters are NaN due to certain bad geometry, kick those out!
                                                  # -> clusters of sparks on edge of chip
  

  However, is this actually enough? Maybe there are cases where the result is not NaN for some bizarre reason? Shouldn't be allowed to happen, because NaN "infects" everything.

  Still: check the mapping of input NaN values to output NaN values! Any weirdness? What do outputs of all those many clusters in the corners look like? Maybe there is something to learn there? Same with very low energy behavior.

I noticed that there seems to be a bug in the sorted_seq based k-d tree implementation of arraymancer! Some test cases do not pass!

• [ ] INVESTIGATE

This might be essential for the limit calculation!

Also: Yesterday evening in discussion with Cris she told me that their group thinks the focal spot of the telescope is actually in the center of the detector chamber and not the focal plane as we always assumed!

• [ ] FIND OUT -> I wrote a message to Johanna to ask if she knows anything more up to date

• [ ] recompile MCMC limit program with different leaf sizes of the k-d tree (try 64)
• [ ] Redo limit calculations for the best 3 cases of MLP & LnL with 10000 toys
• [ ] compute the axion image again in the actual focal spot, then use that input to generate the limit for the best case of the above!
• [X] I moved the data from ./../CastData/data/ to ./../../../mnt/1TB/Uni/DataFromHomeCastData/ to make more space on the home partition!

Yesterday I sent a mail to multiple people from CAST who might know about the beamline design behind the LLNL telescope to find out the intended idea for the focal spot. The question is where in the detector the focal spot was intended to be.

I was always under the impression the focal spot was on the readout plane, but in a discussion with Cristina she mentioned that she thinks it's in the center of the gas volume.

See the mail titled "CAST detectors behind LLNL telescope" (sent from uni bonn address via fastmail).

Turns out from Igor's and Juan's answer the idea was indeed to place the focal spot in the center of the volume!

This is massive and means we need to recompute the raytracing image!

• [X] recompute the axion image! -> DONE!
• [ ] Rerun the correct expected limit calculations incl vetoes! What we ran yesterday didn't include vetoes, hence so slow!
• [ ] Rerun the same case (best case) with the correct axion image!
• [ ] Think about whether we ideally really should have a systematic for the z position of the detector. I.e. varying it change the size of the axion image