search

LDMX-Software / ldmx-sw

The Light Dark Matter eXperiment simulation and reconstruction framework.
https://ldmx-software.github.io
GNU General Public License v3.0
20 stars 15 forks source link

Large positions and energy in Hcal reconstruction #1348

Open EinarElen opened 2 years ago

EinarElen commented 2 years ago

Strange Hcal Reconstruction output

Describe the bug

There seems to be two strange parts of the Hcal reconstruction output (at least with the prototype geometry).

  1. Some hits have positions that are wildly outside of the geometry

    This seems to happen when we have one of the TOA values (TOA_posend or TOA_negend) be negative and the other positive. I suspect that negative TOA values aren’t good in the first place, but I’m not sure.

  2. The distribution of the reconstructed energy is significantly different from the distribution of energy per bar from simulated hits. Specifically, the total reconstructed energy is significantly larger than the total simulated energy

    I haven’t looked into this much but wanted to get the bug report up so I’ll update the issue here. I’m keeping them together because I suspect that they might be related. If they turn out not to be, ill make a separate issue for that.

To reproduce

I don’t think that this should be a problem that is specific to the prototype geometry, unless we have some part of the digitization code that relies on particularities of the regular detector. I haven’t tested with the full detector yet but in principle, you’d just have to change the corresponding line in the configuration file below.

Fastest way to see it is to just run a simulation with the configuration file below, and then just inspecting it either manually or with a script like below.

ldmx fire iss21.py
ldmx python3 checkIss21.py

import ROOT
from ROOT import gSystem
from ROOT import TFile
gSystem.Load('libFramework.so')
file = TFile('iss21.root')
tree = file.Get('LDMX_Events')
# Position tests
for i,event in enumerate(tree):
    recHits = event.HcalRecHits_iss21
    simHits = event.HcalSimHits_iss21
    for j,hit in enumerate(recHits):
        x = hit.getXPos()
        y = hit.getYPos()
        t = hit.getTime()
        if (x > 1000):
            print("X position {:.4f} out of bounds, event {}, hit {}".format(x,i+1,j))
        if (y > 1000):
            print("Y position {:.4f} out of bounds, event {}, hit {}".format(y,i+1,j))
# Energy tests
for i,event in enumerate(tree):
    recHits = event.HcalRecHits_iss21
    simHits = event.HcalSimHits_iss21
    totalSimEnergy = 0.
    totalRecEnergy = 0.
    for hit in simHits:
        totalSimEnergy += hit.getEdep()
    for hit in recHits:
        totalRecEnergy += hit.getEnergy()
    deltaEnergy = totalRecEnergy - totalSimEnergy
    print("Event: {}, E_{{sim}}: {:.2f}, E_{{rec}}: {:.2f}".format(i+1, totalSimEnergy, totalRecEnergy))
    print("E_{{rec}} - E_{{sim}}: {:.2f}".format(deltaEnergy))
    # Don't need all the output
    if (i+1 >= 5):
        break

this should output

X position 8037.7466 out of bounds, event 3, hit 12
Y position 8069.7676 out of bounds, event 53, hit 31
X position 7619.1831 out of bounds, event 106, hit 49
Y position 7347.0024 out of bounds, event 136, hit 5
X position 8186.4165 out of bounds, event 169, hit 38
Y position 7930.2466 out of bounds, event 193, hit 0
Event: 1, E_{sim}: 113.80, E_{rec}: 187.65
E_{rec} - E_{sim}: 73.85
Event: 2, E_{sim}: 540.15, E_{rec}: 620.73
E_{rec} - E_{sim}: 80.58
Event: 3, E_{sim}: 454.35, E_{rec}: 579.76
E_{rec} - E_{sim}: 125.41
Event: 4, E_{sim}: 50.94, E_{rec}: 95.35
E_{rec} - E_{sim}: 44.42
Event: 5, E_{sim}: 95.47, E_{rec}: 176.14
E_{rec} - E_{sim}: 80.67

With GDB inside of the container you can explore it as follows, which lets us see that the large positions comes from the TOA values. I’m setting three breakpoints, you only need the first one. The second and third will break whenever one (but not two) of the TOA values is negative.

ldmx use dev v3.2
ldmx gdb fire
# Note if it is a fresh session, you'll have to respond "yes" to creating a breakpoint despite the shared library not being loaded yet
break HcalRecProducer.cxx:238 if position_bar > 1000

break HcalRecProducer.cxx:217 if TOA_posend < 0 && TOA_negend > 0

break HcalRecProducer.cxx:221 if TOA_negend < 0 && TOA_posend > 0

run iss21.py

# When one of the breakpoints triggers, the following commands can be used to check out the variables

print TOA_posend
print TOA_negend
# You might have to press next a couple of times to get to where this is calculated
print position_bar

# First break has negative TOA_negend, but only ~-9 so position is within bar (258)
continue
# Second break is similar
continue
# Third triggers the position > 1000 trigger
# Should output 8037.7464360643971
# Should output -53.27899194934173
# Should output 32.516454878479138

Configuration file

from LDMX.Framework import ldmxcfg
from LDMX.SimCore import generators
from LDMX.SimCore import simulator

process = ldmxcfg.Process("iss21")
process.logFrequency = 1
process.termLogLevel = 1
process.maxEvents = 200
process.run = 0
simulation = simulator.simulator("iss21")
# simulation.setDetector('ldmx-det-v13')
simulation.setDetector('ldmx-hcal-prototype-v1.0')
gun = generators.gun('particle_gun')

gun.particle  = 'proton'
gun.energy    = 8.
gun.position  = [0., 0., 0.]
gun.direction = [0., 0., 1.]

simulation.generators=[gun]

import LDMX.Hcal.HcalGeometry
import LDMX.Ecal.EcalGeometry

from LDMX.Hcal import digi
from LDMX.Hcal import hcal

digis = digi.HcalDigiProducer()
reco = digi.HcalRecProducer()

from LDMX.Hcal import hcal_hardcoded_conditions
process.sequence=[simulation,
                  digis,
                  reco,
                  ]
process.outputFiles = ['iss21.root']
EinarElen commented 2 years ago

I've been looking around a bit for these things, @cmantill is a negative TOA value ever reasonable? If not, I might have an idea of where to look into

cmantill commented 2 years ago

It shouldn't be, but is possible that the TOA correction is actually causing this? Can you check if it is positive before this line? https://github.com/LDMX-Software/Hcal/blob/trunk/src/Hcal/HcalRecProducer.cxx#L103

EinarElen commented 2 years ago

I'll mess around and see what I can find!

EinarElen commented 2 years ago

I went through one call to getTOA that yields a negative TOA. The loop variables unfold as seen in the end of this comment (the counter is the sample number).

This then gives the TOA as

double toa = (maxSample - toaSample) * clock_cycle_ - toaRelStartBX;
// toa = (3 - 3) * 25 - 8.935 = -8.935 
// time w.r.t to the SOI
toa += ((int)iSOI - maxSample) * clock_cycle_;
// toa = -8.935 + (3 - 3) * 25 = -8.935
EinarElen commented 2 years ago

It shouldn't have been possible that it would have been caused by this update 0d1e4b437909f3c1bfee220b52f065bf1f8466da but I checked to confirm. Issue is still there.

EinarElen commented 2 years ago

Ok the energy is just a factor 2 larger than before. I'm guessing this is because the voltage is defined as 

      // set voltage as the sum of both bars
      voltage = (voltage_posend + voltage_negend);

and the energy is given by 

    double num_mips_equivalent = voltage / voltage_per_mip_;
    double energy_deposited = num_mips_equivalent * mip_energy_ ;

Is this double counting intended? If so, I can just deal with it in my analysis

cmantill commented 2 years ago

hi @EinarElen, looking at this on more detail, we made this change of a factor of 2 because of issue 16. I think the num_mips_equivalent should stay as is since we want the total number of PEs saved as the sum, while we could introduce the facto of 2 in the energy calculation: 

 double energy_deposited = num_mips_equivalent * mip_energy_  / 2;

Or, we can leave it for the analysis. Do you have a strong opinion?

EinarElen commented 2 years ago

That would work for me, but wouldn't that create issues for single ended readout which would then get half the energy?

EinarElen commented 2 years ago

hitsOutsideDetector

Here's a plot showing the position issue. It shows results for 10k 2 GeV pi- events with the prototype geometry. The histogram contains all hits that have a position along the bar outside of the bounds of the detector

EinarElen commented 2 years ago

Ok my comments before were a bit confused (or maybe I'm confused now), but I guess that getTOA returns a relative time so that time being negative in and of itself wouldn't be an issue necessarily. However, the problem still remains. And it isn't caused (at least note solely) by the correction step.

For example, I have an event where fabs(TOA_posend - TOA_negend) = 80 (with correction) fabs(TOA_posend_copy - TOA_negend_copy) = 78 (without)

When multiplied with v, this gives you a position around 8000. Here, TOA_posend = -53, TOA_negend = 27, TOA_posend_copy = 75, and TOA_negend_copy = -3. So the correction makes a difference but its not obvious here that it is working correctly with bad input.

I'm not sure how to move forward with this on my own, @cmantill do you have any suggestions?