LayeredCalorimeterData::Layer::distance set wrong for ALLEGRO ECal?

[scott-snyder]

Check duplicate issues.

• Checked for duplicates

Goal

hi -

The per-layer distance field of LayeredCalorimeterData is documented as
(https://github.com/AIDASoft/DD4hep/blob/0aafe1c25ab1c834482f109f08a91441510af66a/DDRec/include/DDRec/DetectorData.h#L433)

	/// distance from Origin (or the z-axis) to the inner-most face of the layer
	double distance;

However, createECalBarrelInclined (

k4geo/detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v04_geo.cpp

Line 771 in 211e0de

double rad_first = Rmin;

) seems to fill this in with the distance
to the outer face, not the inner face:

  double rad_first = Rmin;
  double rad_last = 0;
  ...
  for (size_t il = 0; il < layerHeight.size(); il++) {
    ...
    rad_last = rad_first + (layerHeight[il] * scale_fact);
    ...
    rad_first = rad_last;
    ...
    caloLayer.distance = rad_first;

Is this intentional, or should it be fixed?

Operating System and Version

n/a

compiler

n/a

The version of the key4hep stack

n/a

Package Version

d5a7656

Reproducer

After

. /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh

This python script

import DDG4
import ROOT
import os

k4geo = os.environ['K4GEO']
kernel = DDG4.Kernel()
geo_file = "file:" + k4geo + "/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml"
kernel.loadGeometry(geo_file)

dd = kernel.detectorDescription()
print ('rmin', dd.constantAsDouble('EMBarrel_rmin'),
       'rmax', dd.constantAsDouble('EMBarrel_rmax'))

code = """
#include "DDRec/DetectorData.h"
#include "DD4hep/DetElement.h"
const dd4hep::rec::LayeredCalorimeterStruct* getLayerInfo (const dd4hep::DetElement& de)
{
  return de.extension<dd4hep::rec::LayeredCalorimeterData>();
}
"""
ROOT.gInterpreter.Declare(code)

emb = dd.detector('ECalBarrel')
li = ROOT.getLayerInfo (emb)
print ('first distance', li.layers[0].distance,
       'last distance', li.layers[-1].distance)

prints

rmin 217.28 rmax 257.83000000000004
first distance 218.9304102909949 last distance 257.83000347886093

(the crash at the end is an unrelated deletion order problem).

Additional context

No response

[scott-snyder]

Further, after failing to reproduce it, i don't think the calculation
of the layer radii is quite right.

If the distance along the signal board is l and the cylindrical radius is r,
then we have

l = sqrt(r^2 - (rmin*sin(angle))^2) - rmin*cos(angle)

For r=rmax, this agrees with the calculation of the signal board length:

  // Electrode length, given inclination angle and min/max radius of the active calorimeter volume
  double planeLength = -Rmin * cos(angle) + sqrt(pow(Rmax, 2) - pow(Rmin * sin(angle), 2));

call this L.

The calculation in the loop above is then finding the cylindrical radius
corresponding to a length l along the signal board as the linear scaling

r = l/L * (rmax-rmin)  +  rmin

But l is not a linear function of r, so this is not correct.
The correct expression is

r = sqrt(rmin**2 + l**2 + 2*rmin*l*cos(angle))

Is there a reason why this linear approximation was used here?
I don't really see a downside to doing this correctly.

[atolosadelgado]

Hi @scott-snyder @SwathiSasikumar

A sketch from both of you would help clarify the issue (at least for me).

From Scott algebra, I inferred this sketch:

However, I’m unsure if angle in the code refers to something else.

Thanks for your time!

[scott-snyder]

hi -

I set this up as in this (hopefully legible) sketch:

2227_001.pdf

As far as i can make out, an inclination angle of 0 corresponds to the boards being oriented radially. In that case, all the expressions discussed here reduce to r = rmin+l.

The functions r(l) and l(r) don't depend on rmax, only rmin and the angle (alpha). The expressions in the code do depend on rmax because they treat r(l) and l(r) as straight lines between r=rmin, l=0 and r=rmax, l=L.

[giovannimarchiori]

The correct calculations are given in slide 5 of this old talk https://indico.cern.ch/event/985994/contributions/4153095/attachments/2162942/3653083/Exploring%20FCC%20EMB%20Geometry.pdf
The numbers have changed a bit but the formula is correct i.e. r^2 = r_in^2 + L^2 + 2Lr_in cos(alpha).
I think the fact that the wrong formula is used is a (small) mistake that we should fix.

Cheers,
Giovanni

[SwathiSasikumar]

Dear @scott-snyder

Thanks for pointing this out. You are right that the linear scaling is a mistake in this case. I will correct this.

Cheers,
Swathi

[scott-snyder]

hi -

Thanks for the pointer to that talk. Something like this should
really be referenced from the code.

A few more comments/questions. These are all small things, but it would
be good to get things nailed down precisely. Some of this is a bit
orthogonal to above, so if desired, i could move to a separate issue.

I've never seen it mentioned anywhere, but curves of constant theta are in
fact not in general straight lines on the readout boards. Rather, they
are hyperbolas --- the intersection of the cone defined by theta=constant
and the plane formed by the readout board. It's small, but if i plot it out,
i can perceive the curvature for the curves close to the edge of the board.
This does not seem to have been taken into account in, eg, the figure on p19
of that talk. If instead of heaving perfectly projective cells,
the lines are instead meant to be straight on the readout board, then
in principle the segmentation code should be adjusted.

Second, i don't understand where the exact numbers used in the segmentation
come from:

      <segmentation type="FCCSWGridModuleThetaMerged_k4geo" nModules="ECalBarrelNumPlanes" mergedCells_Theta="2 4 2 1 2 1 2 2 1 1 1" mergedModules="2 1 1 2 2 1 1 1 2 2 1" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>

The theta offset is given as 0.5902785.
But the corner of the calorimeter seems to be at rho=217.28 and z=306.62,
which corresponds to theta=atan2(217.28,306.62)=0.6164941753384752.
If i instead use the numbers given on p. 3 of that talk, i get
theta=atan2(216, 300) = 0.6240230529767569. Indeed, when i plot out
the cells, the first theta bins completely miss the calorimeter.
The bin width is given as 0.009817477 (for all except the strip layer).
This seems to match the value of 0.5625deg = 0.009817477042468103
given on p. 4 of the talk. It seems to have been chosen to give
100 coarse bins in each half, and does seem to work out
(once one realizes that the offset is to the center of the first bin):

(pi/2 - 0.5902785 + 0.009817477/8) / 0.009817477 = 99.99972614347827

I guess in the end it doesn't really matter if there are unused theta
bins at the beginning and end, but it seems a bit funny then to set
the number of bins, inclusive of the unused ones, to a round number.

Finally, it is a bit surprising that module/phi index of 0 is not
at (nor even close to) phi = 0.

>>> volid = 4  # mod=0, layer=0, theta=0
>>> detectorDescription = kernel.detectorDescription()
>>> vman=DDG4.VolumeManager.getVolumeManager(detectorDescription)
>>> vc = vman.lookupContext(volid)
>>> vc.localToWorld([0,0,0]).phi()
-0.8716917291774724

This seems to come from this line in xml geometry:

        <!-- offset defines the numbering of the modules: module==0 for phi=0 direction -->
        <dimensions rmin="EMBarrel_rmin" rmax="EMBarrel_rmax" dz="EMBarrel_dz" offset="-InclinationAngle"/>

but the comment here does not seem correct.
If i change the offset to 0, i get instead

>>> vc.localToWorld([0,0,0]).phi()
0.004114489473282187

which is just about the angular width of one unit.

[giovannimarchiori]

I've never seen it mentioned anywhere, but curves of constant theta are in
fact not in general straight lines on the readout boards. Rather, they
are hyperbolas --- the intersection of the cone defined by theta=constant
and the plane formed by the readout board. It's small, but if i plot it out,
i can perceive the curvature for the curves close to the edge of the board.
This does not seem to have been taken into account in, eg, the figure on p19
of that talk. If instead of heaving perfectly projective cells,
the lines are instead meant to be straight on the readout board, then
in principle the segmentation code should be adjusted.

good point, the lines of the readout board are straight lines, but I wonder how large the correction would be and if somebody already thought about this in the past (the concept has been there for quite some time) - maybe @gartrog or @BrieucF ?

Second, i don't understand where the exact numbers used in the segmentation
come from:

  <segmentation type="FCCSWGridModuleThetaMerged_k4geo" nModules="ECalBarrelNumPlanes" mergedCells_Theta="2 4 2 1 2 1 2 2 1 1 1" mergedModules="2 1 1 2 2 1 1 1 2 2 1" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>

The theta offset is given as 0.5902785.
But the corner of the calorimeter seems to be at rho=217.28 and z=306.62,
which corresponds to theta=atan2(217.28,306.62)=0.6164941753384752.
If i instead use the numbers given on p. 3 of that talk, I get
theta=atan2(216, 300) = 0.6240230529767569. Indeed, when i plot out
the cells, the first theta bins completely miss the calorimeter.
The bin width is given as 0.009817477 (for all except the strip layer).
This seems to match the value of 0.5625deg = 0.009817477042468103
given on p. 4 of the talk. It seems to have been chosen to give
100 coarse bins in each half, and does seem to work out
(once one realizes that the offset is to the center of the first bin):

(pi/2 - 0.5902785 + 0.009817477/8) / 0.009817477 = 99.99972614347827
I guess in the end it doesn't really matter if there are unused theta
bins at the beginning and end, but it seems a bit funny then to set
the number of bins, inclusive of the unused ones, to a round number.

the dimensions of the calorimeter have changed a bit from time to time as some details of the calo or other parts of the detector changed, we've kept in the simulation the original theta granularity to avoid changing these small details for every change of parameters, in the end it's not a problem if there are a few unused cells at the beginning, it would be worse if there were parts of the acceptance not covered by the readout. Of course in the future when dimensions are stable we could refine these. One more thing to consider is that we want the readout to cover the full envelope of the calorimeter so that we can have hits also in the inactive parts when we make them active in the simulation to study sampling fractions and leakage up/downstream

Finally, it is a bit surprising that module/phi index of 0 is not
at (nor even close to) phi = 0.

volid = 4 # mod=0, layer=0, theta=0
detectorDescription = kernel.detectorDescription()
vman=DDG4.VolumeManager.getVolumeManager(detectorDescription)
vc = vman.lookupContext(volid)
vc.localToWorld([0,0,0]).phi()
-0.8716917291774724
This seems to come from this line in xml geometry:

    <!-- offset defines the numbering of the modules: module==0 for phi=0 direction -->
    <dimensions rmin="EMBarrel_rmin" rmax="EMBarrel_rmax" dz="EMBarrel_dz" offset="-InclinationAngle"/>

but the comment here does not seem correct.
If i change the offset to 0, i get instead

vc.localToWorld([0,0,0]).phi()
0.004114489473282187
which is just about the angular width of one unit.

This comment in the xml has been there since a long time
HEP-FCC/FCCDetectors@3e03649
maybe @faltovaj can comment

My understanding is that module==0 is oriented along phi=0 direction (i.e. is parallel to the x axis).

[scott-snyder]

hi Giovanni -

Thanks for your reply.

good point, the lines of the readout board are straight lines, but I wonder how large the correction would be and if somebody already thought about this in the past (the concept has been there for quite some time) - maybe @gartrog or @BrieucF ?

The largest distance between the constant-theta curves and the straight-line approximation on the surface of the
readout board is ~ 1cm. (Largest distance is 0.78cm, largest displacement in the z-direction is 1.02 cm.) Here's a diagram; black curves show constant theta, and the red curves are straight lines:

pcbdraw.pdf

the dimensions of the calorimeter have changed a bit from time to time as some details of the calo or other parts of the detector changed, we've kept in the simulation the original theta granularity to avoid changing these small details for every change of parameters, in the end it's not a problem if there are a few unused cells at the beginning, it would be worse if there were parts of the acceptance not covered by the readout. Of course in the future when dimensions are stable we could refine these. One more thing to consider is that we want the readout to cover the full envelope of the calorimeter so that we can have hits also in the inactive parts when we make them active in the simulation to study sampling fractions and leakage up/downstream

Fair enough, thanks.

This comment in the xml has been there since a long time HEP-FCC/FCCDetectors@3e03649 maybe @faltovaj can comment

My understanding is that module==0 is oriented along phi=0 direction (i.e. is parallel to the x axis).

I agree that seems to be the intention, but it does not seem to be what is actually done.

[giovannimarchiori]

Hi Scott,

This comment in the xml has been there since a long time HEP-FCC/FCCDetectors@3e03649 maybe @faltovaj can comment
My understanding is that module==0 is oriented along phi=0 direction (i.e. is parallel to the x axis).

I agree that seems to be the intention, but it does not seem to be what is actually done.

What I meant is that module = 0 is horizontal so it has 0 azimuth - if you print the coordinates of all layers with module = 0 they all have the same y = -166.88408469920424

[scott-snyder]

Ah, thank you, it makes sense now. The language was a bit ambiguous, and i was picturing it wrongly.

[giovannimarchiori]

Hi @scott-snyder
this PR #433 should fix the radius calculation for the layered calorimeter data for the ALLEGRO ecal in v03/04