Skip to content

Commit

Permalink
Implement module-theta readout (with per-layer merging) (#56)
Browse files Browse the repository at this point in the history 
* implement module-theta (with per-layer merging) readout

* implement module-theta (with per-layer merging) readout

* cleanup

* cleanup

* cleanup

* performance improvements

* little fix for granularity

* read number of modules from xml

* add protection to check that calculated numPlanes is consistent with constant in xml

* small updates after review of PR

* add a comment to the code

* add number of layers info (currently only for printout of merging config but could be used to enforce that length of vector passed via xml to configure merging per layer is correct

---------

Co-authored-by: Giovanni Marchiori <[email protected]>
  • Loading branch information
@giovannimarchiori
giovannimarchiori and Giovanni Marchiori authored Sep 26, 2023
1 parent 70a989a commit 75c3e84
Show file tree
Hide file tree
Showing 10 changed files with 639 additions and 5 deletions.
143 changes: 143 additions & 0 deletions Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
<?xml version="1.0" encoding="UTF-8"?>
<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">

<info name="FCCee_ECalBarrel"
title="Settings for FCCee Inclined ECal Barrel Calorimeter"
author="M.Aleksa,J.Faltova,A.Zaborowska, V. Volkl"
url="no"
status="development"
version="1.0">
<comment>
Settings for the inclined EM calorimeter.
The barrel is filled with liquid argon. Passive material includes lead in the middle and steal on the outside, glued together.
Passive plates are inclined by a certain angle from the radial direction.
In between of two passive plates there is a readout.
Space between the plate and readout is of trapezoidal shape and filled with liquid argon.
Definition of sizes, visualization settings, readout and longitudinal segmentation are specified.
</comment>
</info>

<define>
<!-- Inclination angle of the lead plates -->
<constant name="InclinationAngle" value="50*degree"/>
<!-- thickness of active volume between two absorber plates at barrel Rmin, measured perpendicular to the readout plate -->
<constant name="LArGapThickness" value="1.239749*2*mm"/>

<!-- Air margin, thicknesses of cryostat and LAr bath -->
<constant name="AirMarginThickness" value="54*mm"/> <!-- Space holder for air gap between cryostat vessels -->

<constant name="CryoBarrelFrontWarm" value="10*mm"/> <!-- Al solid corresponding to 0.11 X0 -->
<constant name="CryoBarrelFrontCold" value="3.8*mm"/> <!-- Al solid equivalent of 0.043 X0 sandwich CFRP -->
<constant name="CryoBarrelFront" value="CryoBarrelFrontWarm+CryoBarrelFrontCold"/>

<constant name="CryoBarrelBackCold" value="30*mm"/> <!-- Al solid corresponding to 0.34 X0 -->
<constant name="CryoBarrelBackWarm" value="2.7*mm"/> <!-- Al solid equivalent of 0.03 X0 sandwich CFRP -->
<constant name="SolenoidBarrel" value="70*mm"/> <!-- Al solenoid with thickness of 0.8 X0 -->
<constant name="CryoBarrelBack" value="CryoBarrelBackWarm+SolenoidBarrel+CryoBarrelBackCold"/>

<constant name="CryoBarrelSideWarm" value="30*mm"/>
<constant name="CryoBarrelSideCold" value="3.8*mm"/>
<constant name="CryoBarrelSide" value="CryoBarrelSideWarm+CryoBarrelSideCold"/>

<constant name="LArBathThicknessFront" value="5*mm"/>
<constant name="LArBathThicknessBack" value="40*mm"/>

<!-- air margin around calorimeter -->
<constant name="BarCryoECal_rmin" value="BarECal_rmin+AirMarginThickness"/>
<constant name="BarCryoECal_rmax" value="BarECal_rmax-AirMarginThickness"/>
<constant name="BarCryoECal_dz" value="BarECal_dz"/>
<!-- calorimeter active volume -->
<constant name="EMBarrel_rmin" value="BarCryoECal_rmin+CryoBarrelFront+LArBathThicknessFront"/>
<constant name="EMBarrel_rmax" value="BarCryoECal_rmax-CryoBarrelBack-LArBathThicknessBack"/>
<constant name="EMBarrel_dz" value="BarECal_dz-CryoBarrelSide"/>
<!-- thickness of active volume between two absorber plates at EMBarrel_rmin, measuring perpendicular to the readout plate -->
<constant name="LAr_thickness" value="LArGapThickness"/>
<!-- passive layer consists of lead in the middle and steel on the outside, glued -->
<!-- When employing trapezoidal planes Pb_thickness corresponds to the minimum thickness, i.e at the front of the calo -->
<constant name="Pb_thickness" value="1.80*mm"/>
<constant name="planeLength" value="-EMBarrel_rmin*cos(InclinationAngle) + sqrt(EMBarrel_rmax*EMBarrel_rmax - EMBarrel_rmin*EMBarrel_rmin*sin(InclinationAngle)*sin(InclinationAngle))"/>
<constant name="ECalBarrelNumPlanes" value="1545"/>
<constant name="ECalBarrelNumLayers" value="12"/>
<constant name="phi" value="asin(planeLength / EMBarrel_rmax * sin(InclinationAngle))"/>
<!-- use a different value for Pb_thickness_max when employing trapezoidal planes -->
<!-- approximate constant sampling fraction: make the absorber grow linearly with the radius,
taking into account the angular projection effect -->
<!-- <constant name="Pb_thickness_max" value="1.3 * Pb_thickness * EMBarrel_rmax/EMBarrel_rmin *
cos(InclinationAngle - phi) / cos(InclinationAngle)" />-->
<constant name="Pb_thickness_max" value="Pb_thickness" />
<!-- total amount of steel in one passive plate: it is divided for the outside layer on top and bottom -->
<constant name="Steel_thickness" value="0.1*mm"/>
<!-- total amount of glue in one passive plate: it is divided for the outside layer on top and bottom -->
<constant name="Glue_thickness" value="0.1*mm"/>
<!-- readout in between two absorber plates -->
<constant name="readout_thickness" value="1.2*mm"/>
</define>

<display>
<vis name="ecal_envelope" r="0.1" g="0.2" b="0.6" alpha="1" showDaughers="false" visible="true" />
</display>

<readouts>
<!-- readout for the simulation, with the baseline merging: 2x along the module direction in each layer; 4x along theta in each layer except layer 1 -->
<!-- the lists mergedCells_Theta and mergedModules define the number of cells to group together in the theta and module direction as a function of the layer -->
<readout name="ECalBarrelModuleThetaMerged">
<segmentation type="FCCSWGridModuleThetaMerged" nModules="1545" mergedCells_Theta="4 1 4 4 4 4 4 4 4 4 4 4" mergedModules="2 2 2 2 2 2 2 2 2 2 2 2" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>
<id>system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10</id>
</readout>

<!-- example of adding a second readout for the reconstruction, to compare the two -->
<readout name="ECalBarrelModuleThetaMerged2">
<segmentation type="FCCSWGridModuleThetaMerged" nModules="1545" mergedCells_Theta="2 4 2 1 2 1 2 2 1 1 1 2" mergedModules="2 1 1 2 2 1 1 1 2 2 1 1" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>
<id>system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10</id>
</readout>
</readouts>

<detectors>
<detector id="BarECal_id" name="ECalBarrel" type="EmCaloBarrelInclined" readout="ECalBarrelModuleThetaMerged">
<type_flags type=" DetType_CALORIMETER + DetType_ELECTROMAGNETIC + DetType_BARREL"/>
<sensitive type="SimpleCalorimeterSD"/>
<dimensions rmin="BarCryoECal_rmin" rmax="BarCryoECal_rmax" dz="BarCryoECal_dz" vis="ecal_envelope"/>
<cryostat name="ECAL_Cryo">
<material name="Aluminum"/>
<dimensions rmin1="BarCryoECal_rmin" rmin2="BarCryoECal_rmin+CryoBarrelFront" rmax1="BarCryoECal_rmax-CryoBarrelBack" rmax2="BarCryoECal_rmax" dz="BarCryoECal_dz"/>
<front sensitive="false"/> <!-- inner wall of the cryostat -->
<side sensitive="false"/> <!-- both sides of the cryostat -->
<back sensitive="false"/> <!-- outer wall of the cryostat -->
</cryostat>
<calorimeter name="EM_barrel">
<!-- 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"/>
<active thickness="LAr_thickness">
<material name="LAr"/>
<!-- overlap offset is a specific feature of the construction; do not change! -->
<!-- one volume for a gap on both side of the readout) -->
<overlap offset="0.5"/>
</active>
<passive>
<rotation angle="InclinationAngle"/> <!-- inclination angle -->
<inner thickness="Pb_thickness" sensitive="false">
<material name="Lead"/>
</inner>
<innerMax thickness="Pb_thickness_max" sensitive="false">
<material name="Lead"/>
</innerMax>
<glue thickness="Glue_thickness" sensitive="false">
<material name="lArCaloGlue"/>
</glue>
<outer thickness="Steel_thickness" sensitive="false">
<material name="lArCaloSteel"/>
</outer>
</passive>
<readout thickness="readout_thickness" sensitive="false">
<material name="PCB"/>
</readout>
<layers> <!-- pcb electrode segmentation in the radial direction -->
<layer thickness="1.5*cm" repeat="1"/>
<layer thickness="3.5*cm" repeat="11"/>
</layers>
</calorimeter>
</detector>
</detectors>
</lccdd>
45 changes: 45 additions & 0 deletions Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
<?xml version="1.0" encoding="UTF-8"?>
<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">

<info name="FCCee-IDEA-LAr Master"
title="FCCee-IDEA-LAr Master"
author="Valentin Volkl"
url="no"
status="development"
version="1.0">
<comment>
Master compact file describing the latest developments of the FCCee IDEA detector concept with a LAr calorimeter.
</comment>
</info>

<include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml" />

<includes>
<gdmlFile ref="../../DetCommon/compact/elements.xml"/>
<gdmlFile ref="../../DetCommon/compact/materials.xml"/>
</includes>

<define>
<constant name="world_size" value="25*m"/>
<constant name="world_x" value="world_size"/>
<constant name="world_y" value="world_size"/>
<constant name="world_z" value="world_size"/>
</define>

<include ref="./FCCee_DectDimensions.xml" />

<include ref="../../DetFCCeeIDEA/compact/Beampipe.xml"/>
<include ref="../../DetFCCeeIDEA/compact/BeamInstrumentation.xml"/>
<include ref="../../DetFCCeeIDEA/compact/LumiCal.xml"/>
<include ref="../../DetFCCeeIDEA/compact/HOMAbsorber.xml"/>
<include ref="../../DetFCCeeCLD/compact/FCCee_o2_v02/Vertex.xml"/>
<include ref="../../DetFCCeeIDEA/compact/SimplifiedDriftChamber.xml"/>
<include ref="../../DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml"/>
<include ref="../../DetFCCeeHCalTile/compact/FCCee_HCalBarrel_TileCal.xml"/>
<include ref="../../DetFCCeeCalDiscs/compact/FCCee_EcalEndcaps_coneCryo.xml"/>
<include ref="../../DetFCCeeHCalTile/compact/FCCee_HCalEndcaps_ThreeParts_TileCal.xml"/>
<include ref="MuonTagger.xml"/>

</lccdd>
42 changes: 42 additions & 0 deletions Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,25 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
layersTotalHeight += layer.repeat() * layer.thickness();
}
lLog << MSG::DEBUG << "Number of layers: " << numLayers << " total thickness " << layersTotalHeight << endmsg;
// The following code checks if the xml geometry file contains a constant defining
// the number of layers the barrel. In that case, it makes the program abort
// if the number of planes in the xml is different from the one calculated from
// the geometry. This is because the number of layers is needed
// in other parts of the code (the readout for the FCC-ee ECAL with
// inclined modules).
int nLayers = -1;
try {
nLayers = aLcdd.constant<int>("ECalBarrelNumLayers");
}
catch(...) {
;
}
if (nLayers > 0 && nLayers != numLayers) {
lLog << MSG::ERROR << "Incorrect number of layers (ECalBarrelNumLayers) in xml file!" << endmsg;
// todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure"));
// make the code crash (incidentSvc does not work)
assert(nLayers == numLayers);
}

dd4hep::xml::DetElement readout = calo.child(_Unicode(readout));
std::string readoutMaterial = readout.materialStr();
Expand Down Expand Up @@ -185,6 +204,29 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
<< " rotation angle = " << angle << endmsg;
uint numPlanes =
round(M_PI / asin((passiveThickness + activeThickness + readoutThickness) / (2. * caloDim.rmin() * cos(angle))));

// The following code checks if the xml geometry file contains a constant defining
// the number of planes in the barrel. In that case, it makes the program abort
// if the number of planes in the xml is different from the one calculated from
// the geometry. This is because the number of plane information (retrieved from the
// xml) is used in other parts of the code (the readout for the FCC-ee ECAL with
// inclined modules). In principle the code above should be refactored so that the number
// of planes is one of the inputs of the calculation and other geometrical parameters
// are adjusted accordingly. This is left for the future, and we use the workaround
// below to enforce for the time being that the number of planes is "correct"
int nModules = -1;
try {
nModules = aLcdd.constant<int>("ECalBarrelNumPlanes");
}
catch(...) {
;
}
if (nModules > 0 && nModules != numPlanes) {
lLog << MSG::ERROR << "Incorrect number of planes (ECalBarrelNumPlanes) in xml file!" << endmsg;
// todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure"));
// make the code crash (incidentSvc does not work)
assert(nModules == numPlanes);
}
double dPhi = 2. * M_PI / numPlanes;
lLog << MSG::INFO << "number of passive plates = " << numPlanes << " azim. angle difference = " << dPhi << endmsg;
lLog << MSG::INFO << " distance at inner radius (cm) = " << 2. * M_PI * caloDim.rmin() / numPlanes << "\n"
Expand Down
117 changes: 117 additions & 0 deletions Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
#ifndef DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H
#define DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H

// FCCSW
#include "DetSegmentation/GridTheta.h"
#include "DD4hep/VolumeManager.h"

/** FCCSWGridModuleThetaMerged Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged.h FCCSWGridModuleThetaMerged.h
*
* Segmentation in theta and module.
* Based on GridTheta, merges modules and theta cells based on layer number
*
*/

namespace dd4hep {
namespace DDSegmentation {
class FCCSWGridModuleThetaMerged : public GridTheta {
public:
/// default constructor using an arbitrary type
FCCSWGridModuleThetaMerged(const std::string& aCellEncoding);
/// Default constructor used by derived classes passing an existing decoder
FCCSWGridModuleThetaMerged(const BitFieldCoder* decoder);

/// destructor
virtual ~FCCSWGridModuleThetaMerged() = default;

/// read n(modules) from detector metadata
void GetNModulesFromGeom();

/// read n(layers) from detector metadata
void GetNLayersFromGeom();

/** Determine the local position based on the cell ID.
* @param[in] aCellId ID of a cell.
* return Position (relative to R, phi of Geant4 volume it belongs to, scaled for R=1).
*/
virtual Vector3D position(const CellID& aCellID) const;
/** Determine the cell ID based on the position.
* @param[in] aLocalPosition (not used).
* @param[in] aGlobalPosition
* @param[in] aVolumeId ID of the Geant4 volume
* return Cell ID.
*/
virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition,
const VolumeID& aVolumeID) const;
/** Determine the azimuthal angle (relative to the G4 volume) based on the cell ID.
* @param[in] aCellId ID of a cell.
* return Phi.
*/
double phi(const CellID& aCellID) const;
/** Determine the polar angle (relative to the G4 volume) based on the cell ID.
* @param[in] aCellId ID of a cell.
* return Theta.
*/
double theta(const CellID& aCellID) const;
/** Determine the radius based on the cell ID.
* @param[in] aCellId ID of a cell.
* return Radius.
*/
// double radius(const CellID& aCellID) const;
/** Get the number of merged cells in theta for given layer
* @param[in] layer
* return The number of merged cells in theta
*/
inline int mergedThetaCells(const int layer) const {
if (layer<m_mergedCellsTheta.size())
return m_mergedCellsTheta[layer];
else
return 1;
}
/** Get the number of merged modules (inclined in phi)
* @param[in] layer
* return The number of merged modules
*/
inline int mergedModules(const int layer) const {
if (layer<m_mergedModules.size())
return m_mergedModules[layer];
else
return 1;
}
/** Get the total number of modules of detector
* return The number of modules (as it was set by the user in the xml file..)
*/
inline int nModules() const { return m_nModules; }
/** Get the number of layers
* return The number of layers
*/
inline int nLayers() const { return m_nLayers; }
/** Get the field name used for the layer
* return The field name for the layer.
*/
inline const std::string& fieldNameLayer() const { return m_layerID; }
/** Get the field name used for the module number
* return The field name for the module.
*/
inline const std::string& fieldNameModule() const { return m_moduleID; }

protected:
/// the field name used for layer
std::string m_layerID;
/// the field name used for the read-out module (can differ from module due to merging)
std::string m_moduleID;
/// vector of number of cells to be merged along theta for each layer
std::vector<int> m_mergedCellsTheta;
/// vector of number of modules to be merged for each layer
std::vector<int> m_mergedModules;

/// number of modules (or, equivalently, the deltaPhi between adjacent modules)
int m_nModules;

/// number of layers (from the geometry)
int m_nLayers;

};
}
}
#endif /* DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H */
Loading

0 comments on commit 75c3e84

Please sign in to comment.