lorenzetti-hep / lorenzetti

Lorenzetti: Empowering Physics Performance and Analysis with Low-level Calorimetry Data
GNU General Public License v3.0
8 stars 13 forks source link

Implements Optimal Filter #2

Closed ingoncalves closed 1 year ago

ingoncalves commented 2 years ago

Until now, Optimal Filter's weights are being designed by an external tool and hardcoded into Lorenzetti's code. In this PR, we fully implement the Optimal Filter estimator, including some useful features:

Each of these features is detailed below.

We had to refactor and implement some tools. Pulse Shape code was extracted from PulseGenerator into a new tool called CaloPulseShapeMaker. This allows us to reuse the pulse shape in the Optimal Filter design.

Optimal Filter and Pulse Shape

Optimal Filter now uses the pulse shape associated with each cell through CaloPulseShapeMaker. It is initialized in the geometry/DetectorATLASModel/python/CaloCellBuilder.py files as:

pulse_shape = CaloPulseShapeMaker("CaloPulseShapeMaker",
                                  ShaperFile  = seg.ShaperFile,
                                  OutputLevel = self.__outputLevel
                                 )

Control the number of constraints (OF1 and OF2)

The number of constraints used in the Optimal Filter optimization process can be controlled. In the geometry/DetectorATLASModel/python/DetectorLayers.py file, we useUseOF2 = True to configure the Optimal Filter with two restrictions (phase and pedestal immunity). Setting this flag to False gives us the OF1 implementation, which uses only the phase immunity restriction. See an example below.

# using OF1
psb = CaloSamplingMaker( "PSB", "Collection_PSB", cells["PSB"],
                    ShaperFile      = basepath + "/pulseLar.dat",
                    BunchIdStart    = -21, # -525ns
                    BunchIdEnd      = 3, # +75ns
                    StartSamplingBC = -2, 
                    NSamples        = 5,
                    EletronicNoise  = 90, # MeV
                    UseOF2          = False,
                    )

# using OF2
tilecal1 = CaloSamplingMaker( "TileCal1", "Collection_TileCal1", cells['TileCal1'],
                    ShaperFile      = basepath+"/pulseTile.dat",
                    BunchIdStart    = -6, # -150ns
                    BunchIdEnd      = 4, # +100ns
                    StartSamplingBC = -3, 
                    NSamples        = 7,
                    EletronicNoise  = 20, # MeV
                    UseOF2          = True,
                    )

Calibration with Minimum Bias datasets

We can also calibrate the Optimal Filter by using the covariance matrix of the noise (noise whitening). To do so, we need to provide a noise dataset file, which is a text file composed of a matrix (each cell separated by whitespace, and one row per line).

For instance, if a cell uses a pulse with 7 samples, the noise dataset should be like

# file geometry/DetectorATLASModel/data/noise-dataset-tile.dat
30.6789 30.4616 30.5791 47.2937 100.1157 121.6271 78.5617
47.1437 35.8586 34.7989 44.4713 148.5912 256.1586 187.6970
95.3952 72.3061 49.1841 38.0381 29.2746 29.9340 31.8294
30.9098 31.2579 31.4419 58.2861 89.3281 68.1710 70.5107
...

with as many rows as you need.

In possession of these datasets, we can enable the calibration by using the OFCalibrationFile property as follows:

# file geometry/DetectorATLASModel/python/DetectorLayers.py
tilecal1 = CaloSamplingMaker( "TileCal1", "Collection_TileCal1", cells['TileCal1'],
                    ShaperFile        = basepath+"/pulseTile.dat",
                    BunchIdStart      = -6, # -150ns
                    BunchIdEnd        = 4, # +100ns
                    StartSamplingBC   = -3, 
                    NSamples          = 7,
                    EletronicNoise    = 20, # MeV
                    UseOF2            = True,
                    OFCalibrationFile = basepath+"/noise-dataset-tile.dat",
                    )

Speedup with static coefficients

As the Optimal Filter coefficients are now calculated in runtime, it may slow down the simulation process. If this is an issue for you, you can use pre-calculated static weights. This skips the entire weight calculation step, speeding up the simulation. To do so, use the OFWeights property as follows:

# file geometry/DetectorATLASModel/python/DetectorLayers.py
tilecal1 = CaloSamplingMaker( "TileCal1", "Collection_TileCal1", cells['TileCal1'],
                    ShaperFile      = basepath+"/pulseTile.dat",
                    BunchIdStart    = -6, # -150ns
                    BunchIdEnd      = 4, # +100ns
                    StartSamplingBC = -3, 
                    NSamples        = 7,
                    EletronicNoise  = 20, # MeV
                    OFWeights       =  [-0.3781, -0.3572, 0.1808, 0.8125, 0.2767, -0.2056, -0.3292]
                    )
micaelverissimo commented 2 years ago

@jodafons, @jlieberm, @edmaregidio for me it's clear!

ingoncalves commented 2 years ago

Thanks @micaelverissimo! Shall we merge the PR then?

jodafons commented 2 years ago

Super clear description. only one comment to change. About the changes, did you check with the output (e.g Zee et) before the changes is the same after you introduce this modification?

Unfortunately, we don't have any CI to compare with some reference file. I would like to merge this after this comparison. Example: Zee et energy from the cluster vs after, for this test you just only need to run HIT->ESD->AOD (don't need to rerun the Geant4 step, only the reconstruction).

ingoncalves commented 2 years ago

Hi, @jodafons. Thanks for your feedback!

I totally agree with you. I'll do this test and as soon as I get the result I'll let you know.

micaelverissimo commented 2 years ago

@ingoncalves, could please check the conflicts? Or can I close this?

micaelverissimo commented 1 year ago

There are several conflicts with this pull request. Do we need this for now @jodafons?

jodafons commented 1 year ago

Not sure, @jlieberm could you check if this PR will be include into the release (today)?

micaelverissimo commented 1 year ago

Closing this since there is no activity.