Commit 3732550a authored by Tadej Novak's avatar Tadej Novak
Browse files

Merge branch 'hgtd_fix_sihit_coords_21p9' into '21.9'

Fixing SiHit xEta and xPhi for HGTD and corresponding HGTD_ModuleDesign setup

See merge request !47155
parents 43ec3b4a e3aede33
......@@ -996,14 +996,12 @@ InDetDD::HGTD_ModuleDesign* HGTD_DetectorFactory::createHgtdDesign( double thick
int diodeRowsPerCircuit = cellRowsPerCircuit;
InDetDD::PixelDiodeMatrix* normalCell = new InDetDD::PixelDiodeMatrix(phiPitch, etaPitch);
InDetDD::PixelDiodeMatrix* singleRow = new InDetDD::PixelDiodeMatrix(InDetDD::PixelDiodeMatrix::etaDir, 0,
InDetDD::PixelDiodeMatrix* singleRow = new InDetDD::PixelDiodeMatrix(InDetDD::PixelDiodeMatrix::phiDir, 0,
normalCell, diodeColumnsPerCircuit, 0);
InDetDD::PixelDiodeMatrix* fullMatrix = new InDetDD::PixelDiodeMatrix(InDetDD::PixelDiodeMatrix::phiDir, 0,
InDetDD::PixelDiodeMatrix* fullMatrix = new InDetDD::PixelDiodeMatrix(InDetDD::PixelDiodeMatrix::etaDir, 0,
singleRow, 2*diodeRowsPerCircuit, 0); // note 30 = 2*15 rows adopted
DetectorDesign::Axis yDirection = InDetDD::DetectorDesign::xAxis;
if (m_geomVersion == 0 )
yDirection = InDetDD::DetectorDesign::yAxis;
DetectorDesign::Axis yDirection = InDetDD::DetectorDesign::yAxis;
InDetDD::HGTD_ModuleDesign* design = new InDetDD::HGTD_ModuleDesign(thickness,
circuitsPerColumn, circuitsPerRow,
......
......@@ -22,8 +22,8 @@ HGTD_ModuleDesign::HGTD_ModuleDesign(const double thickness,
int readoutSide,
DetectorDesign::Axis yDirection,
DetectorDesign::Axis depthDirection):
DetectorDesign(thickness,
true, true, true, // phi,eta,depth axes symmetric
DetectorDesign(thickness,
false, false, true, // phi,eta,depth axes symmetric
carrierType,
readoutSide,
yDirection,
......@@ -40,14 +40,14 @@ HGTD_ModuleDesign::~HGTD_ModuleDesign() {
delete m_bounds;
}
// Returns distance to nearest detector edge
// Returns distance to nearest detector edge
// +ve = inside
// -ve = outside
void
HGTD_ModuleDesign::distanceToDetectorEdge(const SiLocalPosition & localPosition,
double & etaDist, double & phiDist) const
{
// This assume element is centered at 0,0
{
// This assume element is centered at 0,0
// As the calculation is symmetric around 0,0 we only have to test it for one side.
double xEta = abs(localPosition.xEta());
double xPhi = abs(localPosition.xPhi());
......@@ -57,7 +57,7 @@ HGTD_ModuleDesign::distanceToDetectorEdge(const SiLocalPosition & localPosition,
// Distance to top/bottom
etaDist = xEtaEdge - xEta;
// Distance to right/left edge
phiDist = xPhiEdge - xPhi;
}
......@@ -112,7 +112,7 @@ double HGTD_ModuleDesign::widthFromRowRange(const int rowMin, const int rowMax)
double
HGTD_ModuleDesign::phiPitch() const
{
// Average pitch.
// Average pitch.
return width() / rows();
}
......@@ -132,7 +132,7 @@ HGTD_ModuleDesign::etaPitch() const
return length() / columns();
}
const Trk::SurfaceBounds &
const Trk::SurfaceBounds &
HGTD_ModuleDesign::bounds() const
{
// We create on demand as width and length are 0 when HGTD_ModuleDesign first gets
......@@ -148,7 +148,7 @@ SiDiodesParameters HGTD_ModuleDesign::parameters(const SiCellId & cellId) const
SiLocalPosition HGTD_ModuleDesign::localPositionOfCell(const SiCellId & cellId) const
{
return m_diodeMap.parameters(cellId).centre();
return m_diodeMap.parameters(cellId).centre();
}
int HGTD_ModuleDesign::numberOfConnectedCells(const SiReadoutCellId & readoutId) const
......
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
from AthenaCommon import CfgMgr
# The earliest bunch crossing time for which interactions will be sent
# to the HGTD Digitization code.
def HGTD_FirstXing():
return -50
# The latest bunch crossing time for which interactions will be sent
# to the HGTD Digitization code.
def HGTD_LastXing():
return 25
# NOTE: related to 3BC mode?
######################################################################################
def getHGTD_SurfaceChargesGenerator(name="HGTD_SurfaceChargesGenerator", **kwargs):
from HGTD_Digitization.HGTD_DigitizationConf import HGTD_SurfaceChargesGenerator
from AthenaCommon import CfgMgr
return CfgMgr.HGTD_SurfaceChargesGenerator(name, **kwargs)
######################################################################################
def getHGTD_FrontEndTool(name="HGTD_FrontEndTool", **kwargs):
from Digitization.DigitizationFlags import digitizationFlags
from HGTD_Digitization.HGTD_DigitizationConf import HGTD_FrontEndTool
from AthenaCommon import CfgMgr
return CfgMgr.HGTD_FrontEndTool(name, **kwargs)
######################################################################################
def commonHGTD_DigitizationConfig(name,**kwargs):
from Digitization.DigitizationFlags import digitizationFlags
kwargs.setdefault("InputCollecionName", "HGTD_Hits")
kwargs.setdefault("RndmSvc", digitizationFlags.rndmSvc() )
streamName = "HGTD_Digitization" # FIXME ideally random stream name would be configurable
if not digitizationFlags.rndmSeedList.checkForExistingSeed(streamName):
digitizationFlags.rndmSeedList.addSeed(streamName, 49261510, 105132394 )
if digitizationFlags.doXingByXingPileUp():
kwargs.setdefault("FirstXing", HGTD_FirstXing())
kwargs.setdefault("LastXing", HGTD_LastXing() )
from AthenaCommon import CfgMgr
return CfgMgr.HGTD_DigitizationTool(name,**kwargs)
######################################################################################
def getHGTD_DigitizationTool(name="HGTD_DigitizationTool", **kwargs):
kwargs.setdefault("OutoutRDOContName", "HGTD_RDOs")
kwargs.setdefault("OutputSDOName", "HGTD_SDO_Map")
return commonHGTD_DigitizationConfig(name,**kwargs)
######################################################################################
## FIXME needs to go back in later
# def getHGTD_Amp(name="HGTD_Amp", **kwargs):
# kwargs.setdefault("CrossFactor2sides", 0.1)
# kwargs.setdefault("CrossFactorBack", 0.07)
# kwargs.setdefault("PeakTime", 21)
# kwargs.setdefault("deltaT", 1.0)
# kwargs.setdefault("Tmin", -25.0)
# kwargs.setdefault("Tmax", 150.0)
# kwargs.setdefault("NbAverage", 0)
# from HGTD_Digitization.HGTD_DigitizationConf import HGTD_Amp
# return HGTD_Amp(name, **kwargs)
######################################################################################
def getHGTD_FrontEnd(name="HGTD_FrontEnd", **kwargs):
from Digitization.DigitizationFlags import digitizationFlags
from HGTD_Digitization.HGTD_DigitizationConf import HGTD_FrontEnd
return HGTD_FrontEnd(name, **kwargs)
......@@ -71,12 +71,16 @@ void HGTD_SurfaceChargesGenerator::createSurfaceChargesFromHit(
double sensor_thickness = element->design().thickness();
int readout_side = element->design().readoutSide();
double pixel_size_x = element->design().phiPitch();
double pixel_size_y = element->design().etaPitch();
float pixel_size_xphi = element->design().phiPitch();
float pixel_size_xeta = element->design().etaPitch();
const CLHEP::Hep3Vector start_pos(hit.localStartPosition());
const CLHEP::Hep3Vector end_pos(hit.localEndPosition());
ATH_MSG_DEBUG("start_pos xEta=" << start_pos[SiHit::xEta]
<< ", xPhi=" << start_pos[SiHit::xPhi]
<< ", xDep=" << start_pos[SiHit::xDep]);
CLHEP::Hep3Vector direction = end_pos - start_pos;
double deposit_length = direction.mag();
int n_steps = deposit_length / m_small_step_length + 1;
......@@ -94,8 +98,11 @@ void HGTD_SurfaceChargesGenerator::createSurfaceChargesFromHit(
// FIXME needed to check for deposits in guardrings. This should be taken over
// by the module design class and not hardcoded here!
float x_offset = 9.75;
float y_offset = 19.5;
float xphi_offset = 9.75;
float xeta_offset = 19.5;
// FIXME this should be handled by the module design class in the future
float interpad = 50 * CLHEP::micrometer;
for (int i_step = 0; i_step < n_steps; i_step++) {
CLHEP::Hep3Vector surface_pos = start_pos + i_step * direction;
......@@ -103,7 +110,7 @@ void HGTD_SurfaceChargesGenerator::createSurfaceChargesFromHit(
// Distance between charge and readout side. p_design->readoutSide() is
// +1 if readout side is in +ve depth axis direction and visa-versa.
// FIXME ask Noemi about what happens here
double spess =
float spess =
0.5 * sensor_thickness - readout_side * surface_pos[SiHit::xDep];
if (spess < 0) {
spess = 0; // FIXME this means I am on the surface already?
......@@ -114,42 +121,39 @@ void HGTD_SurfaceChargesGenerator::createSurfaceChargesFromHit(
// position at the surface, adding smearing
// FIXME currently no Lorentz angle considered, can be studied in the future
double surf_pos_x = surface_pos[SiHit::xEta] +
rdif * CLHEP::RandGaussZiggurat::shoot(rndm_engine);
double surf_pos_y = surface_pos[SiHit::xPhi] +
rdif * CLHEP::RandGaussZiggurat::shoot(rndm_engine);
float surf_pos_xphi = surface_pos[SiHit::xPhi] +
rdif * CLHEP::RandGaussZiggurat::shoot(rndm_engine);
float surf_pos_xeta = surface_pos[SiHit::xEta] +
rdif * CLHEP::RandGaussZiggurat::shoot(rndm_engine);
// if the deposit is outside the guard ring, don't consider it
if (fabs(surf_pos_x) > x_offset or fabs(surf_pos_y) > y_offset) {
if (fabs(surf_pos_xphi) > xphi_offset or
fabs(surf_pos_xeta) > xeta_offset) {
ATH_MSG_DEBUG("Hit in guard ring");
continue;
}
// FIXME this should be handled by the module design class in the future
float interpad = 50 * CLHEP::micrometer;
int bin_xphi = floor(fabs(surf_pos_xphi + xphi_offset) / pixel_size_xphi);
int bin_xeta = floor(fabs(surf_pos_xeta + xeta_offset) / pixel_size_xeta);
int bin_x = floor(fabs(surf_pos_x + x_offset) / pixel_size_x);
int bin_y = floor(fabs(surf_pos_y + y_offset) / pixel_size_y);
float pos_xphi_inpixel =
fabs(surf_pos_xphi + xphi_offset) - float(bin_xphi) * pixel_size_xphi;
float pos_xeta_inpixel =
fabs(surf_pos_xeta + xeta_offset) - float(bin_xeta) * pixel_size_xeta;
float pos_x_inpixel =
fabs(surf_pos_x + x_offset) - float(bin_x) * pixel_size_x;
float pos_y_inpixel =
fabs(surf_pos_y + y_offset) - float(bin_y) * pixel_size_y;
bool is_interpad_x =
(pos_x_inpixel < interpad or pos_x_inpixel > (pixel_size_x - interpad));
bool is_interpad_y =
(pos_y_inpixel < interpad or pos_y_inpixel > (pixel_size_y - interpad));
bool is_interpad_xphi = (pos_xphi_inpixel < interpad or
pos_xphi_inpixel > (pixel_size_xphi - interpad));
bool is_interpad_xeta = (pos_xeta_inpixel < interpad or
pos_xeta_inpixel > (pixel_size_xeta - interpad));
// check if the charge is sitting in the interpad region
if (is_interpad_x or is_interpad_y) {
ATH_MSG_DEBUG("INTERPAD DETECTED!!!");
ATH_MSG_DEBUG("surf_pos_x=" << surf_pos_x
<< ", surf_pos_y=" << surf_pos_y);
if (is_interpad_xphi or is_interpad_xeta) {
ATH_MSG_DEBUG("Hit in interpad region");
continue;
}
// charges deposited within the active sensor get added
const InDetDD::SiLocalPosition position(
element->hitLocalToLocal(surf_pos_x, surf_pos_y));
element->hitLocalToLocal(surf_pos_xeta, surf_pos_xphi));
SiSurfaceCharge surface_charge(
position, SiCharge(charge_per_step, time_of_flight, hitproc,
......@@ -157,7 +161,7 @@ void HGTD_SurfaceChargesGenerator::createSurfaceChargesFromHit(
InDetDD::SiCellId cell_id =
element->cellIdOfPosition(surface_charge.position());
ATH_MSG_DEBUG("cell_id x=" << cell_id);
if (cell_id.isValid()) {
// add this charge to the collection (or merge in existing charged
// diode)
......
......@@ -273,41 +273,37 @@ StatusCode HGTD_SmearedDigitizationTool::digitize() {
HepGeom::Point3D<double> hit_loc_start_pos = hit->localStartPosition();
HepGeom::Point3D<double> hit_loc_end_pos = hit->localEndPosition();
HepGeom::Point3D<double> hit_glob_start_pos =
curr_det_element->hitLocalToLocal3D(hit_loc_start_pos);
HepGeom::Point3D<double> hit_glob_end_pos =
curr_det_element->hitLocalToLocal3D(hit_loc_end_pos);
double globa_entry_x = hit_glob_start_pos.x();
double globa_entry_y = hit_glob_start_pos.y();
ATH_MSG_DEBUG("hit meanTime=" << hit->meanTime());
ATH_MSG_DEBUG("hit_loc_start_pos xEta="
<< hit_loc_start_pos[SiHit::xEta]
<< ", xPhi=" << hit_loc_start_pos[SiHit::xPhi]
<< ", xDep=" << hit_loc_start_pos[SiHit::xDep]);
double globa_entry_x = hit_loc_start_pos[SiHit::xPhi];
double globa_entry_y = hit_loc_start_pos[SiHit::xEta];
m_x_entry_hit = globa_entry_x;
m_y_entry_hit = globa_entry_y;
m_z_entry_hit = hit_glob_start_pos.z();
m_z_entry_hit = hit_loc_start_pos[SiHit::xDep];
double globa_exit_x = hit_glob_end_pos.x();
double globa_exit_y = hit_glob_end_pos.y();
double globa_exit_x = hit_loc_end_pos[SiHit::xPhi];
double globa_exit_y = hit_loc_end_pos[SiHit::xEta];
m_x_exit_hit = globa_exit_x;
m_y_exit_hit = globa_exit_y;
m_z_exit_hit = hit_glob_end_pos.z();
m_z_exit_hit = hit_loc_end_pos[SiHit::xDep];
double dist_x = std::abs(std::abs(globa_exit_x) - std::abs(globa_entry_x));
double dist_y = std::abs(std::abs(globa_exit_y) - std::abs(globa_entry_y));
Amg::Vector2D local_entry(globa_entry_x, globa_entry_y);
Amg::Vector2D local_exit(globa_exit_x, globa_exit_y);
// get the identifier of the entry and the exit
Identifier entry_id = curr_det_element->identifierOfPosition(local_entry);
Identifier exit_id = curr_det_element->identifierOfPosition(local_exit);
const InDetDD::SiLocalPosition position_entry(
curr_det_element->hitLocalToLocal(globa_entry_y, globa_entry_x));
const InDetDD::SiLocalPosition position_exit(
curr_det_element->hitLocalToLocal(globa_exit_y, globa_exit_x));
// now get the cellIds and check whether they're valid
InDetDD::SiCellId entry_cell_id =
curr_det_element->cellIdFromIdentifier(entry_id);
curr_det_element->cellIdOfPosition(position_entry);
InDetDD::SiCellId exit_cell_id =
curr_det_element->cellIdFromIdentifier(exit_id);
curr_det_element->cellIdOfPosition(position_exit);
ATH_MSG_DEBUG("--- HGTD_SmearedDigitizationTool: entryId "
<< entry_id << " --- exitId " << exit_id);
ATH_MSG_DEBUG("--- HGTD_SmearedDigitizationTool: entryCellId "
<< entry_cell_id << " --- exitCellId " << exit_cell_id);
......
......@@ -30,7 +30,7 @@
#include "CxxUtils/make_unique.h"
HGTDSensorSD::HGTDSensorSD(const std::string& name)
: G4VSensitiveDetector( name ),
: G4VSensitiveDetector( name ),
m_HitColl( "HGTD_Hits" )
{
......@@ -109,15 +109,15 @@ G4bool HGTDSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
G4cout << "DEBUG HGTDG4SD : LOCAL: x: " << pos_current.x()*CLHEP::mm << ", y: " << pos_current.y()*CLHEP::mm << ", z: " << pos_current.z()*CLHEP::mm << G4endl;
if ( i >= 1) {
const G4AffineTransform transformation2 = myTouch->GetHistory()->GetTransform( i-1 ); // Transformation from global to 1 up
const G4AffineTransform transformation2 = myTouch->GetHistory()->GetTransform( i-1 ); // Transformation from global to 1 up
G4AffineTransform transformation_up; // Transformation from current to 1 up
transformation_up.Product( transformationInverse, transformation2 );
G4ThreeVector pos_up = transformation_up.TransformPoint( pos_current );
G4RotationMatrix rotmat = transformation_up.NetRotation(); // https://www-zeuthen.desy.de/geant4/clhep-2.0.4.3/classCLHEP_1_1HepRotation.html
G4ThreeVector translation = transformation_up.NetTranslation(); // https://www-zeuthen.desy.de/ILC/geant4/clhep-2.0.4.3/classCLHEP_1_1Hep3Vector.html
G4ThreeVector translation = transformation_up.NetTranslation(); // https://www-zeuthen.desy.de/ILC/geant4/clhep-2.0.4.3/classCLHEP_1_1Hep3Vector.html
G4cout << "DEBUG HGTDG4SD : Rotation:"
<< "| xx:" << rotmat.xx() << ", xy: " << rotmat.xy() << ", xz: " << rotmat.xz()
<< "| yx:" << rotmat.yx() << ", yy: " << rotmat.yy() << ", yz: " << rotmat.yz()
<< "| xx:" << rotmat.xx() << ", xy: " << rotmat.xy() << ", xz: " << rotmat.xz()
<< "| yx:" << rotmat.yx() << ", yy: " << rotmat.yy() << ", yz: " << rotmat.yz()
<< "| zx:" << rotmat.zx() << ", zy: " << rotmat.zy() << ", zz: " << rotmat.zz() << " | " << G4endl;
G4cout << "DEBUG HGTDG4SD : Translation: x: " << translation.x() << ", y:" << translation.y() << ", z:" << translation.z() << G4endl;
G4cout << "DEBUG HGTDG4SD : TRANSFORMED: x:" << pos_up.x()*CLHEP::mm << ", y:" << pos_up.y()*CLHEP::mm << ", z:" << pos_up.z()*CLHEP::mm << G4endl;
......@@ -141,27 +141,27 @@ G4bool HGTDSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
G4ThreeVector localPosition1 = transformation.TransformPoint(startCoord);
G4ThreeVector localPosition2 = transformation.TransformPoint(endCoord);
HepGeom::Point3D<double> lP1,lP2;
lP1[SiHit::xEta] = localPosition1[2]*CLHEP::mm;
lP1[SiHit::xPhi] = localPosition1[1]*CLHEP::mm;
lP1[SiHit::xDep] = localPosition1[0]*CLHEP::mm;
lP2[SiHit::xEta] = localPosition2[2]*CLHEP::mm;
lP2[SiHit::xPhi] = localPosition2[1]*CLHEP::mm;
lP2[SiHit::xDep] = localPosition2[0]*CLHEP::mm;
lP1[SiHit::xEta] = localPosition1[1]*CLHEP::mm; //long edge of the module
lP1[SiHit::xPhi] = localPosition1[0]*CLHEP::mm; //short edge of the module
lP1[SiHit::xDep] = localPosition1[2]*CLHEP::mm; //depth (z)
lP2[SiHit::xEta] = localPosition2[1]*CLHEP::mm;
lP2[SiHit::xPhi] = localPosition2[0]*CLHEP::mm;
lP2[SiHit::xDep] = localPosition2[2]*CLHEP::mm;
std::string module_indices = myTouch->GetVolume(1)->GetLogicalVolume()->GetName();
std::size_t found = module_indices.find_last_of("_");
// get indices from the volume name
// nomenclature is expected to be e.g. "HGTDModule0_layer_0_1_2"
// for layer=0, phi=1, eta=2 (defined from HGTD_DetectorFactory)
int eta = atoi((module_indices.substr(found+1)).c_str());
module_indices.erase(found);
module_indices.erase(found);
found = module_indices.find_last_of("_");
int phi = atoi((module_indices.substr(found+1)).c_str());
module_indices.erase(found);
module_indices.erase(found);
found = module_indices.find_last_of("_");
int layer = atoi((module_indices.substr(found+1)).c_str());
......
......@@ -45,6 +45,11 @@ StatusCode PadClusterizationAlg::execute() {
const HGTD::HGTD_RDOContainer* rdo_container = rdo_container_handle.cptr();
if (not rdo_container) {
//assume this is fast digi input and ignore
return StatusCode::SUCCESS;
}
// register the PRD container in storegate
SG::WriteHandle<HGTD::HGTD_ClusterContainer> prd_container_handle(
m_prd_wh_key);
......
......@@ -98,9 +98,8 @@ StatusCode TrackTimeExtensionAlg::execute() {
m_sdo_coll_rh_key);
const InDetSimDataCollection* sdo_collection = sdo_collection_handle.cptr();
if (not sdo_collection) {
ATH_MSG_ERROR(
"[TrackTimeExtensionAlg] SDO Collection not found, aborting execute!");
return StatusCode::FAILURE;
ATH_MSG_WARNING("[TrackTimeExtensionAlg] SDO Collection not found, no "
"truth info available!");
}
SG::ReadHandle<McEventCollection> mc_collection_handle(m_mc_coll_rh_key);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment