-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathlikelihood.nim
1953 lines (1794 loc) · 91.6 KB
/
likelihood.nim
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import std / [os, tables, strutils, strformat, algorithm, sets,
stats, sequtils, typetraits, random]
from std/sugar import dup
# import docopt except Value
import cligen / macUt # for `docCommentAdd`
import nimhdf5, tos_helpers, seqmath, arraymancer
import helpers/utils
import ingrid / [ingrid_types, calibration]
import ingrid/calibration/fit_functions
from ingrid / private / geometry import recoEvent
import ingrid / private / [likelihood_utils, fadc_utils]
import ingridDatabase / [databaseRead, databaseDefinitions]
import parsetoml except `{}`
from std / envvars import getEnv # for some additional config, plotSeptem cutoff ...
when defined(cpp):
## What do we need? Flambeau ? NN submodule?
import ../../Tools/NN_playground/nn_predict
from flambeau/flambeau_nn import Device, DeviceKind, Torch, cuda_is_available, init, to, load, NoGradGuard
var Model: MLP
var Desc: MLPDesc
when defined(linux):
const commitHash = staticExec("git rev-parse --short HEAD")
else:
const commitHash = ""
# get date using `CompileDate` magic
const compileDate = CompileDate & " at " & CompileTime
const versionStr = "Version: $# built on: $#" % [commitHash, compileDate]
const IgnoreCdlFile {.booldefine.} = false
const h5cdl_file = currentSourcePath() / "../../../resources/calibration-cdl.h5"
const cdlExists = fileExists(h5cdl_file)
when not cdlExists and not IgnoreCdlFile:
{.fatal: "CAST CDL reference file `calibration-cdl.h5` does not exist at: " &
$h5cdl_file.}
# cut performed regardless of logL value on the data, since transverse
# rms > 1.5 cannot be a physical photon, due to diffusion in 3cm drift
# distance
const RmsCleaningCut = 1.5
# the filter we use globally in this file
let filter = H5Filter(kind: fkZlib, zlibLevel: 4)
when not defined(useMalloc):
{.fatal: "Please compile with `-d:useMalloc` to reduce the amount of memory kept by the " &
"program. This allows to use more jobs in `createAllLikelihoodCombinations`.".}
template withConfig(body: untyped): untyped =
const sourceDir = currentSourcePath().parentDir
let config {.inject.} = parseToml.parseFile(sourceDir / "config.toml")
body
proc readMorphKind(): MorphingKind =
## reads the `morphingKind` field from the TOML file
withConfig:
result = parseEnum[MorphingKind](config["Likelihood"]["morphingKind"].getStr)
proc readSearchRadius(): int =
## reads the cluster finding `searchRadius` field from the TOML file
withConfig:
result = config["Likelihood"]["searchRadius"].getInt
proc readClusterAlgo(): ClusteringAlgorithm =
## Reads the clustering algorithm to use for the septem veto clustering
withConfig:
result = parseEnum[ClusteringAlgorithm](config["Likelihood"]["clusterAlgo"].getStr)
proc readDbscanEpsilon(): float =
## Reads the `ε` to use for the DBSCAN algorithm septem veto clustering
withConfig:
result = config["Likelihood"]["epsilon"].getFloat
proc writeCdlAttributes[T: H5DataSet | H5Group](dset: var T,
cdlFile: string,
year: YearKind) =
## writes information about what datasets were used to calculate the likelihood
## dataset
dset.attrs["year of CDL data"] = $year
dset.attrs["calibration CDL file"] = cdlFile
proc calcLogLikelihood*(h5f: var H5File,
ctx: LikelihoodContext) =
## - read all data of single run
## - get energy dataset
## - create energy bins
## - use digitize to get the correct energy bins for each cluster
## - create individual seq's for each energy bin (containing all
## clusters and the corresponding properties)
## - create histogram for each of these energy binned properties
## in XrayLikelihoodProcessor we have
## - 1 histogram for each property and each energy bin
## - fill them with the log likelihood of all events falling into that histogram
## - logLikelihood:
## get number of elements in histogram at the bin for the element for which we get
## the logL
## # elements / # total in histogram = likelihood. Take log
for (num, group) in runs(h5f):
echo &"Start logL calc of run {group}"
var logL_chips = newSeq[seq[float64]](ctx.vetoCfg.numChips)
# iterate over all chips and perform logL calcs
for (_, chipNumber, grp) in chipGroups(h5f, group):
let logL = calcLikelihoodDataset(h5f, grp, ctx)
# after walking over all events for this chip, add to correct
# index for logL
logL_chips[chipNumber] = logL
# after we added all logL data to the seqs, write it to the file
var logL_dsets = mapIt(toSeq(0 ..< ctx.vetoCfg.numChips),
h5f.create_dataset((group / &"chip_{it}/likelihood"),
logL_chips[it].len,
float64,
overwrite = true,
filter = filter))
# write the data to the file
echo &"Writing data of run {group} to file {h5f.name}"
for tup in zip(logL_dsets, logL_chips):
var (dset, logL) = tup
dset[dset.all] = logL
dset.writeCdlAttributes(ctx.cdlFile, ctx.year)
import datamancer / serialize
proc writeInfos(h5f: H5File, grp: H5Group, ctx: LikelihoodContext,
fadcVetoCount, scintiVetoCount: int,
flags: set[LogLFlagKind]) =
## writes information about used vetoes and the number of events removed by
## the vetos
# write the CDL file info
var mgrp = grp
mgrp.writeCdlAttributes(ctx.cdlFile, ctx.year)
# now serialize the veto settings we used
h5f.toH5(ctx, "logLCtx", grp.name, exclude = @["refSetTuple", "refDf", "refDfEnergy"])
if fadcVetoCount >= 0:
mgrp.attrs["# removed by FADC veto"] = fadcVetoCount
if scintiVetoCount >= 0:
mgrp.attrs["# removed by scinti veto"] = scintiVetoCount
proc writeLikelihoodData(h5f: var H5File,
h5fout: var H5File,
group: var H5Group,
chipNumber: int,
cutTab: CutValueInterpolator,
nnCutTab: CutValueInterpolator,
passedInds: OrderedSet[int],
fadcVetoCount, scintiVetoCount: int, flags: set[LogLFlagKind],
ctx: LikelihoodContext
) =
#durations: (float64, float64)) =
## writes all relevant cluster data of events corresponding to `passedInds` in
## the group given by `group` for chip `chipNumber` to the output file `h5f`
when false:
let (totalDuration, totalDurationPassed) = durations
# TODO: add copy of attributes from energyFromPixel dataset, which contains
# the calibration factor!
# read all float datasets, which we want to write to the output file
var dsetNames = concat(@(getFloatDsetNames()), @(getIntDsetNames()))
# add the final two datasets, which we'd like to write
dsetNames.add @["likelihood", "energyFromPixel", "energyFromCharge", "x", "y", "ToT"]
if ctx.timepix == Timepix3:
dsetNames.add getFloatToANames()
dsetNames.add @["ToA", "ToACombined", "toaMin"]
let chpGrpName = group.name / &"chip_{chipNumber}"
# get mutable group for this chip to copy attributes
var chpGrpIn = h5f[chpGrpName.grp_str]
# now get the event numbers to compare against the indices
let evNumbers = h5f[(chpGrpName / "eventNumber"), int64]
# then write to the new likelihood group for this chip
# create the groups for the likelihood file, run group and chip group
let
runGrpName = group.name.replace("reconstruction", "likelihood")
# group name of this chip's group
logLgroup = runGrpName / &"chip_{chip_number}"
var
runGrp = h5fout.create_group(runGrpName)
chpGrpOut = h5fout.create_group(logLgroup)
## walk all datasets and read data, filter to indices and write to output
for dsetName in dsetNames: # all datasets including vlen
let path = chpGrpName / dsetName
if path in h5f: ## check if it exists so that we can try to write both energy dsets
let h5dset = h5f[(path).dset_str]
withDset(h5dset): # reads full dataset into injected `dset` variable
var data = newSeqOfCap[elementType(dset)](passedInds.card)
for idx in passedInds:
data.add dset[idx]
when not (elementType(dset) is string or elementType(dset) is seq[string]):
var outDset = h5fout.create_dataset((logLgroup / dsetName),
passedInds.len,
elementType(dset),
filter = filter)
outDset[outDset.all] = data
# get all event numbers from hash set by using elements as indices for event numbers
let evNumsPassed = mapIt(passedInds, evNumbers[it]).sorted
# create dataset for allowed indices
var evDset = h5fout.create_dataset((logLgroup / "eventNumber"),
evNumsPassed.len,
int,
filter = filter)
# write event numbers
evDset[evDset.all] = evNumsPassed
# now write the regular datasets, need to match up the event numbers!
if chipNumber == 3: # centerChip # XXX: DO THIS
let passedEvSet = evNumsPassed.toSet
proc copyOver(h5f, h5fout: H5File, grp: H5Group, baseName: string) =
let grpEvNums = h5f[grp.name / "eventNumber", int]
# get the indices matching the event numbers that passed on center chip!
let grpIdxs = toSeq(0 ..< grpEvNums.len).filterIt(grpEvNums[it] in passedEvSet)
for dataset in items(grp): # copy one level (nested is not copied)
if dataset.name.extractFilename() == "pedestalRun": continue # do not copy it a) not useful b) different size
echo "Copying dataset ", dataset
withDset(dataset): # reads full dataset into injected `dset` variable
var data: Tensor[elementType(dset)]
type T = elementType(dset)
when T is SomeNumber: # do not want bool, char, string etc! Also because they break `arraymancer`...
if dataset.shape.len == 1:
data = newTensor[T]([grpIdxs.len]) #newSeqOfCap[elementType(dset)](passedInds.card)
for i, idx in grpIdxs:
data[i] = dset[idx]
else:
doAssert dataset.shape.len == 2, "Datasets with shape " & $dataset.shape & " not yet supported!"
data = newTensor[T]([grpIdxs.len, dataset.shape[1]]) #newSeqOfCap[elementType(dset)](passedInds.card)
let dsetShaped = dset.toTensor.reshape(dataset.shape)
for i, idx in grpIdxs:
data[i, _] = dsetShaped[idx, _]
if data.size > 0: # if there's no data to write, don't create dataset
var outDset = h5fout.create_dataset((baseName / dataset.name.extractFilename),
data.shape.toSeq,
elementType(dset),
filter = filter)
outDset.unsafeWrite(data.toUnsafeView(), data.size)
h5f.copyOver(h5fout, group, runGrpName) # copy over the common datasets
if group.name / "fadc" in h5f:
let fadcGroup = h5f[(group.name / "fadc").grp_str]
h5f.copyOver(h5fout, fadcGroup, runGrpName / "fadc") # copy over the FADC datasets
# finally write all interesting attributes
chpGrpOut.attrs["MorphingKind"] = $cutTab.morphKind
if cutTab.morphKind == mkNone:
for key, val in pairs(cutTab.lnLcutTab):
chpGrpOut.attrs[&"logL cut value: {key}"] = val
chpGrpOut.attrs["SpectrumType"] = "background"
else:
## TODO: add dataset of energies & cut values for cut
echo "WARN: Add information about used energies and cut values for mkLinear!"
discard
# write total number of clusters
chpGrpOut.attrs["Total number of cluster"] = evNumbers.len
when false:
# still missing per chip information
runGrp.attrs["totalDurationRun"] = totalDuration
runGrp.attrs["totalDurationPassed"] = totalDurationPassed
# copy attributes over from the input file
runGrp.copy_attributes(group.attrs)
chpGrpOut.copy_attributes(chpGrpIn.attrs)
h5fout.writeInfos(chpGrpOut, ctx, fadcVetoCount, scintiVetoCount, flags)
runGrp.writeCdlAttributes(ctx.cdlFile, ctx.year)
func isVetoedByFadc(ctx: LikelihoodContext, run, eventNumber: int, fadcTrigger, fadcEvNum: seq[int64],
fadcRise, fadcFall: seq[uint16], fadcSkew: seq[float]): bool =
## returns `true` if the event of `ind` is vetoed by the FADC based on cuts on
## the rise and fall time. Vetoed means the event must be thrown out
## because it does ``not`` conform to the X-ray hypothesis.
## ------------ NOTE --------------
## looking at ~/org/Figs/SPSC_Jan_2019/test/Run3/{rise,fall}Time_normalizedPDF_run3.pdf
## makes it seem like anything above ~130 is probably not an X-ray. Take that for
## now.
## NOTE: these values do not match the typical values seen in the 2017 data taking!
## (which to be fair had different FADC settings etc!)
## NOTE: Numbers adjusted quickly based on the 1, 5, 95, 99-th percentiles of the
## rise / fall time FADC data (see `statusAndProgress` notes)
##
## For a motivation see sec. [[#sec:fadc:noisy_events_and_fadc_veto]] in ~statusAndProgress.org~
##
## UPDATE: The values nowadays are chosen based on a custom percentile of the rise/fall time
## distributions as seen in the 55Fe calibration data (which is handed as the `--calibFile`
## argument).
# get correct veto cut values for the current FADC settings (of this run)
let cut = ctx.vetoCfg.fadcVetoes[ run.toFadcSetting() ]
doAssert cut.active, "The cuts for FADC setting " & $run.toFadcSetting() & " is not active, i.e. " &
"we did not read corresponding run data in the given calibration file: " & $ctx.vetoCfg.calibFile & "!"
result = false
let fIdx = fadcEvNum.lowerBound(eventNumber)
if fIdx < fadcEvNum.high and
fadcEvNum[fIdx] == eventNumber: # and
#fadcTrigger[fIdx] == 1: # cannot use trigger, as it's a global dataset and `fIdx` does not match! Anyhow not needed.
if fadcRise[fIdx].float >= cut.riseLow and
fadcRise[fIdx].float <= cut.riseHigh and
fadcFall[fIdx].float >= cut.fallLow and
fadcFall[fIdx].float <= cut.fallHigh and
fadcSkew[fIdx] <= cut.skewness:
result = false
else:
# outside either of the cuts, throw it out
result = true
func isVetoedByScintis(eventNumber: int,
scintEvNum: seq[int64],
scinti1, scinti2: seq[int64]): bool =
## returns `true` if the event of `ind` is vetoed by the scintillators based
## on the fact that one of the two scintillators had a non trivial scintillator
## count value ( > 0 and < 4095; in practice ~< 400).
## Restrict to 150 to avoid the 255 scintillator 'bug' in the large veto paddle.
## Vetoed means the event must be thrown out, because the event was most
## likely induced by a muon
result = false
# throw out any event with a non trivial (> 0 and < 4095)
# scintillator trigger
const low = 0
const high = 150 # be pessimistic about high
let sIdx = scintEvNum.lowerBound(eventNumber)
if sIdx < scintEvNum.high and
scintEvNum[sIdx] == eventNumber and
((scinti1[sIdx] > low and scinti1[sIdx] < high) or
(scinti2[sIdx] > low and scinti2[sIdx] < high)):
# had a non trivial trigger, throw out
result = true
proc printRow(df: DataFrame) =
doAssert df.len == 1, "not 1 element"
for k in getKeys(df).sorted:
echo k, " = ", df[k, Value][0]
proc evaluateCluster(clTup: (int, ClusterObject[PixInt]),
septemFrame: var SeptemFrame,
septemGeometry: var SeptemEventGeometry,
centerData: CenterClusterData,
gainVals: seq[float],
calibTuple: tuple[b, m: float], ## contains the parameters required to perform energy calibration
σT: float, # diffusion of the run this cluster is part of
ctx: LikelihoodContext,
flags: set[LogLFlagKind],
eventNumber: int
): tuple[logL, nnPred, energy: float, lineVetoPassed: bool] =
## XXX: add these to config.toml and as a cmdline argument in addition
# total charge for this cluster
when defined(cpp): # initialize the Device
var device_type: DeviceKind
if Torch.cuda_is_available():
#echo "CUDA available! Training on GPU."
device_type = kCuda
else:
#echo "Training on CPU."
device_type = kCPU
let Dev = Device.init(device_type)
let clusterId = clTup[0]
let cl = clTup[1]
let clData = cl.data
# reset running stat for each cluster
var rs: RunningStat
var pixIdx = septemFrame.numRecoPixels
var totCharge = 0.0
for pix in clData:
# take each pixel tuple and reconvert it to chip based coordinates
# first determine chip it corresponds to
let pixChip = if ctx.vetoCfg.useRealLayout: determineRealChip(pix)
else: determineChip(pix)
# for this pixel and chip get the correct gas gain and push it to the `RunningStat`
rs.push gainVals[pixChip]
# taken the chip of the pixel, reconvert that to a local coordinate system
# given charge of this pixel, assign it to some intermediate storage
totCharge += septemFrame.charge[pix.y, pix.x]
# overwrite the `septemFrame` by `pix` and cluster id
septemFrame.pixels[pixIdx] = (x: pix.x, y: pix.y, ch: clusterId)
inc pixIdx
septemFrame.numRecoPixels = pixIdx
doAssert totCharge > 0.0, "Total charge of the cluster is 0. This should not happen. Lookup of charge buggy?"
# determine parameters of the lines through the cluster centers
# invert the slope
let slope = tan(Pi - cl.geometry.rotationAngle)
let cX = if ctx.vetoCfg.useRealLayout: toRealXPix(cl.centerX) else: toXPix(cl.centerX)
let cY = if ctx.vetoCfg.useRealLayout: toRealYPix(cl.centerY) else: toYPix(cl.centerY)
let intercept = cY.float - slope * cX.float
septemGeometry.centers.add (x: cX.float, y: cY.float)
septemGeometry.lines.add (m: slope, b: intercept)
# using total charge and `RunningStat` calculate energy from charge
## TODO: need to look *only* at gains from the corresponding chips each
## and compute the energy of the part of each chip indidually and add
let (b, m) = calibTuple
let energy = totCharge * linearFunc(@[b, m], rs.mean) * 1e-6
let logL = ctx.calcLikelihoodForEvent(energy, # <- TODO: hack for now!
cl.geometry.eccentricity,
cl.geometry.lengthDivRmsTrans,
cl.geometry.fractionInTransverseRms)
## XXX: Potentially we need to load the model only once. Might be costly?
var nnPred: float
when defined(cpp):
if ctx.vetoCfg.useNeuralNetworkCut:
# build single row DF from cluster
let clusterDf = clusterToDf(cl, logL, energy, totCharge, σT)
doAssert clusterDf.len == 1, "More than 1 row in cluster DF. How?"
# forward it through the network
nnPred = Model.predict(Dev, Desc, clusterDf)[0] #ctx.vetoCfg.nnModelPath.predict(clusterDf)[0]
## Check if the current cluster is in input chip region. If it is, either it is part of something
## super big that makes the center still fall into the gold region or it remains unchanged.
## In the unchanged case, let's compare the energy and cluster pixels
let inRegionOfInterest = inRegion(cl.centerX - TimepixSize, cl.centerY - TimepixSize, ctx.region)
var lineVetoPassed = true #
## As long as this cluster does not contain the original cluster data (depending on the line veto kind!)
## we perform the line veto
let containsOriginal = clTup[1].containsOriginal(septemFrame)
let isOriginal = clTup[1].isOriginal(septemFrame)
let ofInterest =
if isOriginal: false # never care about *only* the original cluster!
elif containsOriginal and ctx.vetoCfg.lineVetoKind in {lvRegular, lvCheckHLC}: true # contained, but we want it
elif containsOriginal and ctx.vetoCfg.lineVetoKind == lvRegularNoHLC: false # contained, but we *dont* want it
else: true # if it neither contains the original, nor is it, it is of interest
if fkLineVeto in flags and ofInterest and
cl.geometry.eccentricity >= ctx.vetoCfg.eccLineVetoCut:
# if this is not in region of interest, check its eccentricity and compute if line points to original
# center cluster
# compute line orthogonal to this cluster's line
let orthSlope = -1.0 / slope
# determine y axis intersection
## Center data must be converted as it doesn't go through reco!
let centerNew = if ctx.vetoCfg.useRealLayout: tightToReal((x: centerData.cX, y: centerData.cY))
else: (x: centerData.cX, y: centerData.cY)
septemGeometry.xCenter = if ctx.vetoCfg.useRealLayout: toRealXPix(centerNew.x)
else: toXPix(centerNew.x)
septemGeometry.yCenter = if ctx.vetoCfg.useRealLayout: toRealYPix(centerNew.y)
else: toYPix(centerNew.y)
let centerInter = septemGeometry.yCenter.float - orthSlope * septemGeometry.xCenter.float
# use orthogonal line to determine crossing point between it and this cluster's line
let xCut = (intercept - centerInter) / (orthSlope - slope)
let yCut = orthSlope * xCut + centerInter
septemGeometry.centers.add (x: xCut, y: yCut)
# compute distance between two points
let dist = sqrt( (septemGeometry.xCenter.float - xCut)^2 + (septemGeometry.yCenter.float - yCut)^2 )
# if distance is smaller than radius
septemGeometry.centerRadius = ((centerData.rmsT +
centerData.rmsL) /
2.0 * 3.0) /
TimepixSize * 256.0
if dist < septemGeometry.centerRadius:
lineVetoPassed = false
septemGeometry.lines.add (m: orthSlope, b: centerInter)
result = (logL: logL,
nnPred: nnPred,
energy: energy,
lineVetoPassed: lineVetoPassed) ## No line veto is equivalent to passing the line veto.
## Init to `true`, just use variable
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]
proc readCenterChipData(h5f: H5File, group: H5Group, ctx: LikelihoodContext): CenterChipData =
## Read all data of this run for the center chip that we need in the rest of
## the septem veto logic
let chpGrp = group.name / "chip_3"
result.logL = h5f[chpGrp / "likelihood", float]
result.energies = h5f[chpGrp / ctx.energyDset.toDset, float]
result.cX = h5f[chpGrp / "centerX", float]
result.cY = h5f[chpGrp / "centerY", float]
result.hits = h5f[chpGrp / "hits", int]
result.rmsT = h5f[chpGrp / "rmsTransverse", float]
result.rmsL = h5f[chpGrp / "rmsLongitudinal", float]
when defined(cpp):
if ctx.vetoCfg.useNeuralNetworkCut: # get NN prediction for center chip if needed
result.nnPred = ctx.predict(h5f, chpGrp)
proc buildSeptemEvent(evDf: DataFrame,
valToCut, energies: seq[float],
cutTab: CutValueInterpolator,
allChipData: AllChipData,
centerChip: int,
useRealLayout: bool,
isNNcuts: bool): SeptemFrame =
## Given a sub DF containing the indices for the clusters of a given event
## number assembles the full septemboard frame using `allChipData`.
##
## `valToCut` is either the `logL` data or the NN prediction for all clusters on
## the center chip. `cutTab` is the corresponding helper containing the cut values.
result = SeptemFrame(centerEvIdx: -1, numRecoPixels: 0)
# assign later to avoid indirection for each pixel
## XXX: only xsize pix, y size pix when dealing with real layout!!
let chargeSize = if useRealLayout: [YSizePix, XSizePix]
else: [256 * 3, 256 * 3]
var chargeTensor = zeros[float](chargeSize)
var pixels: PixelsInt = newSeqOfCap[PixInt](5000) # 5000 = 5k * 26eV = 130keV. Plenty to avoid reallocations
let chipData = evDf["chipNumber", int] # get data as tensors and access to avoid `Value` wrapping
let idxData = evDf["eventIndex", int]
for i in 0 ..< evDf.len:
# get the chip number and event index, dump corresponding event pixel data
# onto the "SeptemFrame"
let chip = chipData[i]
let idx = idxData[i]
# convert to septem coordinate and add to frame
let pixData = getPixels(allChipData, chip, idx, chargeTensor, useRealLayout)
pixels.add pixData
if chip == centerChip:
# NOTE: Index `idx` is `*only valid*` for `valToCut` when we look at the center chip!
let passCut = if isNNCuts: valToCut[idx.int] > cutTab[energies[idx.int]]
else: valToCut[idx.int] < cutTab[energies[idx.int]]
if passCut:
# assign center event index if this is the cluster that passes logL cut
result.centerEvIdx = idx.int
# in this case assign `pixData` to the result as a reference for data of original cluster
if useRealLayout:
result.centerCluster = newSeq[PixInt](pixData.len)
for i in 0 ..< pixData.len:
result.centerCluster[i] = septemPixToRealPix(pixData[i])
else:
result.centerCluster = pixData
result.charge = chargeTensor
result.pixels = pixels
proc buildEventIdxMap(septemDf: DataFrame): Table[int, seq[(int, int)]] =
## Creates a table which maps event numbers to a sequence of pairs (chip, eventIndex)
## such that we can later look up the indices given an event number.
result = initTable[int, seq[(int, int)]]()
for (tup, subDf) in groups(septemDf.group_by("eventNumber")):
let evNum = tup[0][1].toInt
result[evNum] = newSeq[(int, int)]()
let chips = subDf["chipNumber", int]
let idxs = subDf["eventIndex", int]
for i in 0 ..< chips.len:
result[evNum].add (chips[i], idxs[i]) # add chip number & event index for this event number
proc bootstrapFakeEvents(septemDf, centerDf: DataFrame,
passedEvs: var OrderedSet[int],
passedInds: var OrderedSet[int],
fixedCluster: bool): DataFrame =
## Bootstrap fake events to estimate the random coincidence of the septem and line veto
## by rewriting the `septemDf`, which is an index of indices in each chip group to the
## corresponding event number.
##
## We rewrite it by randomly drawing from center clusters and then picking a different
## event from which to assign the outer ring.
##
## Environment variables `NUM_SEPTEM_FAKE` and `SEPTEM_FAKE_FIXED_CLUSTER` can be used to
## adjust the number of bootstrapped events and the cluster index in case we run in
## "fixed center cluster" mode.
let numBootstrap = getEnv("NUM_SEPTEM_FAKE", "2000").parseInt
let fixedClusterIdx = getEnv("SEPTEM_FAKE_FIXED_CLUSTER", "5").parseInt
let idxs = passedInds.toSeq
# from the list of passed indices, we now need to rebuild the `septemDf`.
# We can't just draw a random sample of indices as we have to take care of
# the following things:
# - do not mix the outer chips. All events of these must remain together
# - do not accidentally combine correct events
# how do we get a) the indices for our new `idx`?
# best if we turn the `septemDf` into a `Table[int, seq[int]]` of sorts where the key
# is the event number and the argument corresponds to the *indices* of that event.
let evIdxMap = buildEventIdxMap(septemDf)
let maxEvent = septemDf["eventNumber", int].max # get max event to know up to where to draw
let validEvs = toSeq(keys(evIdxMap))
# We _can_ start by drawing a `numBootstrap` elements from the *passed indices*
# These then define the *center clusters* that we will use as a basis.
# Then: draw an outer ring that is *not* the same event number.
var evOuter = -1
# `i` will just be the new eventNumber for the fake events
result = newDataFrame()
passedInds.clear()
for i in 0 ..< numBootstrap:
# draw an event from which to build a fake event
let idxCenter = idxs[rand(0 .. idxs.high)]
# get event number for this index
let evCenter = centerDf["eventNumber", int][idxCenter]
# So: draw another event number, which is the outer ring we want to map to this
evOuter = rand(0 ..< maxEvent)
while evOuter == evCenter: # if same, draw again until different
evOuter = rand(0 ..< maxEvent)
# get indices for outer chips from map
let outerIdxs = if evOuter in evIdxMap: evIdxMap[evOuter]
else: newSeq[(int, int)]()
# now add this to the resulting DF
var indices = if fixedCluster: @[ idxs[fixedClusterIdx] ]
else: @[ idxCenter ]
var chips = @[3] ## XXX: replace by using center number
for (chip, index) in outerIdxs:
if chip == 3: continue # do not add center chip
indices.add index
chips.add chip
result.add toDf({"eventNumber" : i, "eventIndex" : indices, "chipNumber" : chips })
# finally reset the `passedEvs` to all events we just faked
passedEvs = toSeq(0 ..< numBootstrap).toOrderedSet
from ../../Tools/determineDiffusion/determineDiffusion import getDiffusionForRun
import std / cpuinfo
import pkg / weave
proc getBatchedIndices(septemDf: DataFrame, batchSize: int): seq[(int, int)] =
## Returns a sequence of indices that the `septemDf` can be sliced at
## to keep `batchSize` events together
##
## The slice is inclusive, to be used as `septemDf[slice[0] .. slice[1]]`
# # batch around events
var count = 0
var idx = 0
var lastIdx = 0
for (tup, subDf) in groups(septemDf.group_by("eventNumber")):
# increase counters first!
inc count
inc idx, subDf.len # increase by length
if count == batchSize:
result.add (lastIdx, idx - 1) # inclusive, so -1
lastIdx = idx
count = 0
if count > 0:
result.add (lastIdx, septemDf.high)
const sanity = false
if sanity:
var lastEv = -1
for slice in result:
let subDf = septemDf[slice[0] .. slice[1]]
#echo "Head::: ", subDf.head(10)
#echo "Tail::: ", subDf.tail(10)
#echo "slice", slice
let ev = subDf["eventNumber", int][0]
if lastEv > 0:
doAssert lastEv != ev, "Slice from " & $slice & " evnets " & $ev
lastEv = subDf["eventNumber", int][subDf.high]
proc buildSeptemEvents(ctx: LikelihoodContext, septemDf: DataFrame,
cutTab: CutValueInterpolator,
allChipData: AllChipData,
centerData: CenterChipData,
passedInds, passedEvs: OrderedSet[int]
): (seq[SeptemFrame], seq[RecoInputEvent[PixInt]]) =
result[0] = newSeqOfCap[SeptemFrame](passedInds.len)
result[1] = newSeqOfCap[RecoInputEvent[PixInt]](passedInds.len)
let septemGrouped = septemDf.group_by("eventNumber")
for (pair, evGroup) in groups(septemGrouped):
let evNum = pair[0][1].toInt
#if evNum != 98659: continue
if evNum in passedEvs:
# then grab all chips for this event
let isNNcut = ctx.vetoCfg.useNeuralNetworkCut
let valsToCut = if isNNcut: centerData.nnPred else: centerData.logL
var septemFrame = buildSeptemEvent(evGroup, valsToCut, centerData.energies,
cutTab, allChipData, ctx.vetoCfg.centerChip, ctx.vetoCfg.useRealLayout,
isNNcut)
if septemFrame.centerEvIdx == -1:
echo "Broken event! ", evGroup.pretty(-1)
echo "The event is: ", pair
#continue # skip "broken" events for now (i.e. NaN)
quit(1)
# given the full frame run through the full reconstruction for this cluster
# here we give chip number as -1, indicating "Septem"
## XXX: for full ToA support in Timepix3, need to also read the `toa` data and insert
## it here!
result[0].add septemFrame
result[1].add (pixels: septemFrame.pixels, eventNumber: evNum.int,
toa: newSeq[uint16](), toaCombined: newSeq[uint64](), ftoa: newSeq[uint8]())
proc reconstructSeptemEvents(ctx: LikelihoodContext, septemEvents: seq[RecoInputEvent[PixInt]],
runNumber: int): seq[RecoEvent[PixInt]] =
result = newSeq[RecoEvent[PixInt]](septemEvents.len)
var resBuf = cast[ptr UncheckedArray[RecoEvent[PixInt]]](result[0].addr)
var dataBuf = cast[ptr UncheckedArray[RecoInputEvent[PixInt]]](septemEvents[0].addr)
## XXX: I think this cannot work as `RecoEvent` contains a `seq[ClusterObject]`. Once we return
## from this scope they will be freed as they were created on a different thread.
## Update: or rather the issue is that a seq / tensor of `RecoEvent` is not a flat memory structure.
## Therefore it's not memcopyable and we cannot get a ptr to it in a sane way?
## Update 2: above 21 threads the code results in a segfault. This _seems_ reproducible.
putEnv("WEAVE_NUM_THREADS", $min(countProcessors(), 20))
init(Weave)
let
searchRadius = ctx.vetoCfg.searchRadius
dbscanEpsilon = ctx.vetoCfg.dbscanEpsilon
clusterAlgo = ctx.vetoCfg.clusterAlgo
useRealLayout = ctx.vetoCfg.useRealLayout
parallelFor event in 0 ..< septemEvents.len:
captures: {resBuf, dataBuf, runNumber, searchRadius, dbscanEpsilon, clusterAlgo, useRealLayout}
resBuf[event] = recoEvent(dataBuf[event], -1,
runNumber, searchRadius,
dbscanEpsilon = dbscanEpsilon,
clusterAlgo = clusterAlgo,
useRealLayout = useRealLayout)
echoCounted(event, 5000, msg = " clusters reconstructed")
syncRoot(Weave)
exit(Weave)
delEnv("WEAVE_NUM_THREADS")
# `result` modified via `resBuf`
var eventCounter = 0
proc performSeptemVeto(ctx: LikelihoodContext,
runNumber: int,
recoEvents: seq[RecoEvent[PixInt]],
septemFrames: seq[SeptemFrame],
cutTab: CutValueInterpolator,
allChipData: AllChipData,
centerData: CenterChipData,
passedInds: var OrderedSet[int],
passedEvs: var OrderedSet[int],
gains: seq[seq[GasGainIntervalResult]],
calibTuple: tuple[b, m: float], ## contains the parameters required to perform energy calibration
σT: float,
flags: set[LogLFlagKind]
) =
let PlotCutEnergyLow = getEnv("PLOT_SEPTEM_E_CUTOFF_LOW", "0.0").parseFloat
let PlotCutEnergyHigh = getEnv("PLOT_SEPTEM_E_CUTOFF", "5.0").parseFloat
let PlotEventNumber = getEnv("PLOT_SEPTEM_EVENT_NUMBER", "-1").parseInt # can be used to plot only a single event
let useLineVeto = fkLineVeto in flags
let useSeptemVeto = fkSeptem in flags
for i in 0 ..< recoEvents.len:
inc eventCounter
let recoEv = recoEvents[i]
var septemFrame = septemFrames[i] # <- will be modified
let evNum = recoEv.eventNumber
let centerClusterData = getCenterClusterData(septemFrame, centerData, recoEv, ctx.vetoCfg.lineVetoKind)
# extract the correct gas gain slices for this event
var gainVals: seq[float]
for chp in 0 ..< ctx.vetoCfg.numChips:
let gainSlices = gains[chp]
let gainEvs = gainSlices.mapIt(it.sliceStartEvNum)
let sliceIdx = gainEvs.lowerBound(evNum)
gainVals.add gainSlices[min(gainSlices.high, sliceIdx)].G # add the correct gas gain
# calculate log likelihood of all reconstructed clusters
# Note: `septem` and `line` vetoes are "inverse" in their logic. Only a *single* line needed
# to veto center cluster, but *any* cluster passing septem (hence passed vs rejected)
var septemVetoed = true
var lineVetoed = false
var septemGeometry: SeptemEventGeometry # no need for constructor. `default` is fine
if fkAggressive in flags:
# if there's more than 1 cluster, remove
if recoEv.cluster.len > 1: # <- this is a very bad idea knowing something about random coincidence rates now!
# (anyhow this is not really ever used)
passedInds.excl septemFrame.centerEvIdx
passedEvs.excl evNum
continue # skip to next iteration
for clusterTup in pairs(recoEv.cluster):
let (logL, nnPred, energy, lineVetoPassed) = evaluateCluster(
clusterTup, septemFrame, septemGeometry,
centerClusterData,
gainVals,
calibTuple,
σT,
ctx,
flags,
evNum)
let cX = toXPix(clusterTup[1].centerX)
let cY = toYPix(clusterTup[1].centerY)
var chipClusterCenter: int
if ctx.vetoCfg.useRealLayout:
chipClusterCenter = (x: cX, y: cY, ch: 0).determineRealChip(allowOutsideChip = true)
else:
chipClusterCenter = (x: cX, y: cY, ch: 0).determineChip(allowOutsideChip = true)
var
lnLVeto = false
nnVeto = false
## IMPORTANT: the following code relies on the fact that only *one* branch (lnL or NN)
## is being used. The `cutTab` corresponds to the CutValueInterpolator for *that* branch.
## NN cut
when defined(cpp):
if ctx.vetoCfg.useNeuralNetworkCut and
(classify(nnPred) == fcNaN or # some clusters are NaN due to certain bad geometry, kick those out!
# -> clusters of sparks on edge of chip
nnPred < cutTab[energy]):
# more background like if smaller than cut value, -> veto it
nnVeto = true
## LnL cut
if ctx.vetoCfg.useLnLCut and logL > cutTab[energy]:
# more background like if larger than cut value, -> veto it
lnLVeto = true
# needs to be on center chip & pass both cuts. Then septem veto does *not* remove it
#if evNum == 21860:
# echo "Chip center cluster: ", chipClusterCenter, " nnVeto? ", nnVeto
# echo "== center? ", chipClusterCenter == ctx.vetoCfg.centerChip, " for center: ", ctx.vetoCfg.centerChip
# echo "nnVeto : ", nnVeto, " for ", nnPred, " < ", cutTab[energy], " for energy: ", energy
# echo "Index was: ", septemFrame.centerEvIdx
if chipClusterCenter == ctx.vetoCfg.centerChip and # this cluster's center is on center chip
not lnLVeto and # passes lnL veto cut
not nnVeto: # passes NN veto cut
#if evNum == 98659:
# echo "SEPTEM VETOED, FALSE"
septemVetoed = false # <-- not vetoed!
if not lineVetoed and not lineVetoPassed:
lineVetoed = true
#if evNum == 98659:
# echo "EVENT:::::: lnL ", lnLVeto, " nn ", nnVeto, " septem: ", septemVetoed
if (useSeptemVeto and septemVetoed) or (useLineVeto and lineVetoed):
## If `passed` is still false, it means *no* reconstructed cluster passed the logL now. Given that
## the original cluster which *did* pass logL is part of the septem event, the conclusion is that
## it was now part of a bigger cluster that did *not* pass anymore.
passedInds.excl septemFrame.centerEvIdx
passedEvs.excl evNum
if fkPlotSeptem in flags:
if septemFrame.centerEvIdx < 0:
doAssert false, "this cannot happen. it implies no cluster found in the given event"
if centerData.energies[septemFrame.centerEvIdx] < PlotCutEnergyHigh and # only plot if below energy cut
centerData.energies[septemFrame.centerEvIdx] > PlotCutEnergyLow and # only plot if above lower energy cut
(PlotEventNumber < 0 or PlotEventNumber == evNum): # and given event matches target event (or all)
# shorten to actual number of stored pixels. Otherwise elements with ToT / charge values will remain
# in the `septemFrame`
septemFrame.pixels.setLen(septemFrame.numRecoPixels)
if not useLineVeto: # if no line veto, don't draw lines
septemGeometry.lines = newSeq[tuple[m, b: float]]()
## XXX: it might be a good idea to extend the plotting to include cluster information
## that affects the line veto. We could hand eccentricity & the eccentricity cut to make
## it clearer why a cluster was (not) cut (add the ecc. as text next to cluster center?)
## Ideally for that we would change the code to hand some object instead of centers, lines, ...
## and all that jazz. Instead have one object per cluster.
plotSeptemEvent(septemFrame.pixels, runNumber, evNum,
lines = septemGeometry.lines,
centers = septemGeometry.centers,
septemVetoed = septemVetoed,
lineVetoed = lineVetoed,
xCenter = septemGeometry.xCenter,
yCenter = septemGeometry.yCenter,
radius = septemGeometry.centerRadius,
energyCenter = centerData.energies[septemFrame.centerEvIdx],
useTeX = ctx.useTeX,
plotPath = ctx.plotPath)
#if evNum == 21860:
# echo "Done"
# quit()
proc applySeptemVeto(h5f: var H5File,
runNumber: int,
passedInds: var OrderedSet[int],
cutTab: CutValueInterpolator, ## Important: `cutTab` can either be for LogL or NN!
ctx: LikelihoodContext,
flags: set[LogLFlagKind]) =
## Applies the septem board veto to the given `passedInds` in `runNumber` of `h5f`.
## If an event does not pass the septem veto cut, it is excluded from the `passedInds`
## set.
##
## The environment variable `PLOT_SEPTEM_E_CUTOFF` can be used to adjust the energy
## cutoff for which events to plot when running with `--plotSeptem`. In addition the
## variable `USE_TEX` can be adjusted to generate TikZ TeX plots.
## XXX: add these to config.toml and as a cmdline argument in addition
echo "Passed indices before septem veto ", passedInds.card
let group = h5f[(recoBase() & $runNumber).grp_str]
var septemDf = h5f.getSeptemEventDF(runNumber)
# now filter events for `centerChip` from and compare with `passedInds`
let centerDf = septemDf.filter(f{int: `chipNumber` == ctx.vetoCfg.centerChip})
# The `passedEvs` simply maps the *indices* to the *event numbers* that pass *on the
# center chip*.
var passedEvs = passedInds.mapIt(centerDf["eventNumber", int][it]).sorted.toOrderedSet
## Read all the pixel data for all chips
let allChipData = readAllChipData(h5f, group, ctx.vetoCfg.numChips)
var centerData = readCenterChipData(h5f, group, ctx)
let estimateRandomCoinc = fkEstRandomCoinc in flags
if estimateRandomCoinc:
septemDf = bootstrapFakeEvents(septemDf, centerDf, passedEvs, passedInds, fkEstRandomFixedEvent in flags)
#echo septemDf.pretty(1000)
var fout = open(ctx.septemLineVetoEfficiencyFile, fmAppend)
let useLineVeto = fkLineVeto in flags
let useSeptemVeto = fkSeptem in flags
fout.write(&"Septem events before: {passedEvs.len} (S,L,F) = ({$useSeptemVeto}, {$useLineVeto}, {estimateRandomCoinc})\n")
# use diffusion cache to get diffusion for this run
let runType = parseEnum[RunTypeKind](group.attrs["runType", string])
let σT = if ctx.vetoCfg.useNeuralNetworkCut: getDiffusionForRun(runNumber, isBackground = (runType == rtBackground))
else: -1.0 # irrelevant without NN!
let chips = toSeq(0 ..< ctx.vetoCfg.numChips)
let gains = chips.mapIt(h5f[(group.name / "chip_" & $it / "gasGainSlices"), GasGainIntervalResult])
let septemHChips = chips.mapIt(getSeptemHChip(it))
let ccGrp = h5f[recoPath(runNumber, ctx.vetoCfg.centerChip)] # cc = center chip
let calibTuple = ccGrp.getCalibVsGasGainFactors(septemHChips[ctx.vetoCfg.centerChip], runNumber)
septemDf = septemDf.filter(f{int: `eventNumber` in passedEvs})
when false:
let batchedIdx = getBatchedIndices(septemDf, 500)
for slice in batchedIdx:
echo "Slicing to ", slice
let subDf = septemDf[slice[0] .. slice[1]] # slice is inclusive!
#echo "Sub : ", subDf, " of total length : ", septemDf.len
#echo "Tail of DF: ", subDf.tail(20)
# 1. extract all events as `'SeptemEvents'` from the DF
let (septemFrames, septemEvents) = ctx.buildSeptemEvents(subDf, cutTab, allChipData, centerData,
passedInds, passedEvs)
#echo "Built events"
# 2. reconstruct all events in parallel
let recoEvents = ctx.reconstructSeptemEvents(septemEvents, runNumber)
echo "Reconstructed events"
# 3. perform the actual septem veto
ctx.performSeptemVeto(runNumber,
recoEvents, septemFrames, cutTab, allChipData, centerData,
passedInds, passedEvs,
gains, calibTuple,
σT,
flags)
when false: # no batching
echo "Start building events"
# 1. extract all events as `'SeptemEvents'` from the DF
let (septemFrames, septemEvents) = ctx.buildSeptemEvents(septemDf, cutTab, allChipData, centerData,
passedInds, passedEvs)
echo "Built events"
# 2. reconstruct all events in parallel
let recoEvents = ctx.reconstructSeptemEvents(septemEvents, runNumber)
echo "Reconstructed events"
# 3. perform the actual septem veto
ctx.performSeptemVeto(runNumber,
recoEvents, septemFrames, cutTab, allChipData, centerData,
passedInds, passedEvs,
gains, calibTuple,
σT,
flags)
#echo "EVENT COUNTER : ", eventCounter
when false:
echo "Passed indices after septem veto ", passedEvs.card
fout.write("Septem events after fake cut: " & $passedEvs.len & "\n")
fout.close()
when true:
let PlotCutEnergyLow = getEnv("PLOT_SEPTEM_E_CUTOFF_LOW", "0.0").parseFloat
let PlotCutEnergyHigh = getEnv("PLOT_SEPTEM_E_CUTOFF", "5.0").parseFloat
let PlotEventNumber = getEnv("PLOT_SEPTEM_EVENT_NUMBER", "-1").parseInt # can be used to plot only a single event
let septemGrouped = septemDf.group_by("eventNumber")
for (pair, evGroup) in groups(septemGrouped):
let evNum = pair[0][1].toInt
#if evNum != 98659: continue
if evNum in passedEvs:
# then grab all chips for this event
let isNNcut = ctx.vetoCfg.useNeuralNetworkCut
let valsToCut = if isNNcut: centerData.nnPred else: centerData.logL
var septemFrame = buildSeptemEvent(evGroup, valsToCut, centerData.energies,
cutTab, allChipData, ctx.vetoCfg.centerChip, ctx.vetoCfg.useRealLayout,
isNNcut)
if septemFrame.centerEvIdx == -1:
echo "Broken event! ", evGroup.pretty(-1)
echo "The event is: ", pair
#continue # skip "broken" events for now (i.e. NaN)
quit(1)
# given the full frame run through the full reconstruction for this cluster
# here we give chip number as -1, indicating "Septem"
## XXX: for full ToA support in Timepix3, need to also read the `toa` data and insert
## it here!
let inp = (pixels: septemFrame.pixels, eventNumber: evNum.int,
toa: newSeq[uint16](), toaCombined: newSeq[uint64](), ftoa: newSeq[uint8]())
let recoEv = recoEvent(inp, -1,
runNumber, searchRadius = ctx.vetoCfg.searchRadius,
dbscanEpsilon = ctx.vetoCfg.dbscanEpsilon,
clusterAlgo = ctx.vetoCfg.clusterAlgo,
useRealLayout = ctx.vetoCfg.useRealLayout)
let centerClusterData = getCenterClusterData(septemFrame, centerData, recoEv, ctx.vetoCfg.lineVetoKind)
# extract the correct gas gain slices for this event
var gainVals: seq[float]
for chp in 0 ..< ctx.vetoCfg.numChips:
let gainSlices = gains[chp]
let gainEvs = gainSlices.mapIt(it.sliceStartEvNum)
let sliceIdx = gainEvs.lowerBound(evNum)
gainVals.add gainSlices[min(gainSlices.high, sliceIdx)].G # add the correct gas gain
# calculate log likelihood of all reconstructed clusters
# Note: `septem` and `line` vetoes are "inverse" in their logic. Only a *single* line needed
# to veto center cluster, but *any* cluster passing septem (hence passed vs rejected)
var septemVetoed = true
var lineVetoed = false
var septemGeometry: SeptemEventGeometry # no need for constructor. `default` is fine
if fkAggressive in flags:
# if there's more than 1 cluster, remove
if recoEv.cluster.len > 1: # <- this is a very bad idea knowing something about random coincidence rates now!
# (anyhow this is not really ever used)