From d8866f95015265f7247127de1b7cd5ab12a2393e Mon Sep 17 00:00:00 2001 From: Anthony Denis Frawley Date: Sat, 20 Jun 2026 01:15:46 -0400 Subject: [PATCH 1/7] Implements the ability to tilt the TPC envelope in sPHENIX, and attempts to handle the fallout in the simulation and reconstruction code. --- .../TrackingDiagnostics/TrackResiduals.cc | 4 +- offline/packages/tpc/TpcClusterMover.cc | 19 +- offline/packages/tpc/TpcClusterMover.h | 4 +- offline/packages/tpc/TpcClusterizer.cc | 42 +- offline/packages/trackbase/ActsGeometry.cc | 104 ++- offline/packages/trackbase/ActsGeometry.h | 6 + .../trackbase/AlignmentTransformation.cc | 41 +- .../trackbase_historic/TrackAnalysisUtils.cc | 7 +- .../packages/trackreco/MakeActsGeometry.cc | 82 ++- offline/packages/trackreco/MakeActsGeometry.h | 4 + offline/packages/trackreco/MakeSourceLinks.cc | 6 +- offline/packages/trackreco/MakeSourceLinks.h | 2 +- offline/packages/trackreco/PHActsTrkFitter.cc | 2 +- .../packages/trackreco/PHCosmicsTrkFitter.cc | 4 +- .../g4simulation/g4detectors/Makefile.am | 3 + .../g4simulation/g4detectors/PHG4TpcGeom.h | 53 +- .../g4simulation/g4detectors/PHG4TpcGeomv2.cc | 593 ++++++++++++++++++ .../g4simulation/g4detectors/PHG4TpcGeomv2.h | 158 +++++ .../g4detectors/PHG4TpcGeomv2LinkDef.h | 5 + .../g4simulation/g4eval/SvtxTruthEval.cc | 44 +- .../g4simulation/g4tpc/PHG4TpcDetector.cc | 39 +- .../g4tpc/PHG4TpcElectronDrift.cc | 25 +- .../g4simulation/g4tpc/PHG4TpcElectronDrift.h | 1 + .../g4simulation/g4tpc/TpcClusterBuilder.cc | 4 +- 24 files changed, 1130 insertions(+), 122 deletions(-) create mode 100644 simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc create mode 100644 simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h create mode 100644 simulation/g4simulation/g4detectors/PHG4TpcGeomv2LinkDef.h diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.cc b/offline/packages/TrackingDiagnostics/TrackResiduals.cc index 1a5a69fae0..9c9a09f651 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.cc +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.cc @@ -115,7 +115,9 @@ int TrackResiduals::InitRun(PHCompositeNode* topNode) m_globalPositionWrapper.set_suppressCrossing(m_convertSeeds); // clusterMover needs the correct radii of the TPC layers auto *tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); - m_clusterMover.initialize_geometry(tpccellgeo); + + auto *geometry = findNode::getClass(topNode, "ActsGeometry"); + m_clusterMover.initialize_geometry(tpccellgeo, geometry); m_clusterMover.set_verbosity(0); auto *se = Fun4AllServer::instance(); diff --git a/offline/packages/tpc/TpcClusterMover.cc b/offline/packages/tpc/TpcClusterMover.cc index 39c5e156d4..d1f391b9ea 100644 --- a/offline/packages/tpc/TpcClusterMover.cc +++ b/offline/packages/tpc/TpcClusterMover.cc @@ -45,13 +45,16 @@ TpcClusterMover::TpcClusterMover() } } -void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo) +void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo, ActsGeometry *tGeometry) { if (_verbosity > 0) { std::cout << "TpcClusterMover: Initializing layer radii for Tpc from cell geometry object" << std::endl; } + + _tGeometry = tGeometry; + int layer = 0; PHG4TpcGeomContainer::ConstRange layerrange = cellgeo->get_begin_end(); for (PHG4TpcGeomContainer::ConstIterator layeriter = layerrange.first; @@ -67,8 +70,9 @@ void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo) std::vector> TpcClusterMover::processTrack(const std::vector>& global_in) { // Get the global positions of the TPC clusters for this track, already corrected for distortions, and move them to the surfaces - // The input object contains all clusters for the track - + // The input object contains all clusters for the track in world coordinates + // The surface radii are in envelope coordinates, we transform the positions to envelope coordinates + std::vector> global_moved; std::vector tpc_global_vec; @@ -80,7 +84,8 @@ std::vector> TpcClusterMover::proces if (trkrid == TrkrDefs::tpcId) { tpc_cluskey_vec.push_back(ckey); - tpc_global_vec.push_back(global); + Acts::Vector3 env_global = _tGeometry->transformTpcWorldToEnvelope(global); + tpc_global_vec.push_back(env_global); } else { @@ -140,9 +145,9 @@ std::vector> TpcClusterMover::proces // now move the cluster to the surface radius // we keep the cluster key fixed, change the surface if necessary - Acts::Vector3 global_new(xnew, ynew, znew); - - // add the new position and surface to the return object + Acts::Vector3 env_global_new(xnew, ynew, znew); + // now we transform back to global coordinates and add the new position and surface to the return object + Acts::Vector3 global_new = _tGeometry->transformTpcEnvelopeToWorld(env_global_new); global_moved.emplace_back(cluskey, global_new); if (_verbosity > 2) diff --git a/offline/packages/tpc/TpcClusterMover.h b/offline/packages/tpc/TpcClusterMover.h index 92a053e990..7eaa0b3495 100644 --- a/offline/packages/tpc/TpcClusterMover.h +++ b/offline/packages/tpc/TpcClusterMover.h @@ -24,7 +24,7 @@ class TpcClusterMover //! Updates the assumed default geometry below to that contained in the //! cell geo - void initialize_geometry(PHG4TpcGeomContainer *cellgeo); + void initialize_geometry(PHG4TpcGeomContainer *cellgeo, ActsGeometry *tGeometry); private: int get_circle_circle_intersection(double target_radius, double R, double X0, double Y0, double xclus, double yclus, double &x, double &y) const; @@ -48,6 +48,8 @@ class TpcClusterMover double outer_tpc_spacing = 0.0; int _verbosity = 0; + + ActsGeometry *_tGeometry = nullptr; }; #endif diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index 687d5b58f9..a9ec66a219 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -29,7 +29,8 @@ #include #include // for SubsysReco -#include +//#include +#include #include #include @@ -657,30 +658,29 @@ namespace } } - // This is the global position + // This is the phi position in the tpc_envelope double clusiphi = iphi_sum / adc_sum; - double clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side); + double env_clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side); - double clusx = radius * cos(clusphi); - double clusy = radius * sin(clusphi); + // these positions are in the tpc_envelope + double env_clusx = radius * cos(env_clusphi); + double env_clusy = radius * sin(env_clusphi); double clust = t_sum / adc_sum; // needed for surface identification double zdriftlength = clust * my_data.tGeometry->get_drift_velocity(); // convert z drift length to z position in the TPC - double clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength; + double env_clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength; if (my_data.side == 0) { - clusz = -clusz; + env_clusz = -env_clusz; } - // std::cout << " side " << my_data.side << " clusz " << clusz << " clust " << clust << " driftmax " << my_data.m_tdriftmax << std::endl; const double phi_cov = (iphi2_sum / adc_sum - square(clusiphi)) * pow(my_data.layergeom->get_phistep(), 2); const double t_cov = t2_sum / adc_sum - square(clust); - // Get the surface key to find the surface from the - // TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(my_data.layer, my_data.sector, my_data.side); - Acts::Vector3 global(clusx, clusy, clusz); + // get_tpc_surface_from_coords expects a world global position + Acts::Vector3 env_global(env_clusx, env_clusy, env_clusz); + Acts::Vector3 global = my_data.tGeometry->transformTpcEnvelopeToWorld(env_global); TrkrDefs::subsurfkey subsurfkey = 0; - Surface surface = my_data.tGeometry->get_tpc_surface_from_coords( tpcHitSetKey, global, @@ -693,6 +693,8 @@ namespace hitkeyvec.clear(); return; } + // Acts::Vector3 surfcent = surface->center(my_data.tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm; + // std::cout << " surf center = " << surfcent.x() << " " << surfcent.y() << " " << surfcent.z() << std::endl; // Estimate the errors // Blow up error on single pixel clusters by a factor 3 to compensate for threshold effects @@ -744,6 +746,7 @@ namespace b_made_cluster = true; } + // This code needs to be reviewed in case of a non-zero TPC tilt - ADF 6/16/26 if (use_nn && clus_base && training_hits) { try @@ -1298,16 +1301,11 @@ int TpcClusterizer::InitRun(PHCompositeNode *topNode) AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep(); - std::cout << "FirstLayerCellGeomv1 streamer: " << std::endl; - auto *g1 = static_cast (geom->GetFirstLayerCellGeom()); // cast because << not in the base class - std::cout << *g1 << std::endl; - std::cout << "LayerCellGeomv1 streamer for layer 24: " << std::endl; - auto *g2 = static_cast (geom->GetLayerCellGeom(24)); // cast because << not in the base class - std::cout << *g2 << std::endl; - std::cout << "LayerCellGeomv1 streamer for layer 40: " << std::endl; - auto *g3 = static_cast (geom->GetLayerCellGeom(40)); // cast because << not in the base class - std::cout << *g3 << std::endl; - + // the identify now contains all information from the streamer for v2 + geom->GetFirstLayerCellGeom()->identify(); + geom->GetLayerCellGeom(24)->identify(); + geom->GetLayerCellGeom(40)->identify(); + if (m_maskDeadChannels) { m_deadChannelMap.clear(); diff --git a/offline/packages/trackbase/ActsGeometry.cc b/offline/packages/trackbase/ActsGeometry.cc index 8096be1896..bac8fdf558 100644 --- a/offline/packages/trackbase/ActsGeometry.cc +++ b/offline/packages/trackbase/ActsGeometry.cc @@ -7,6 +7,10 @@ #include +#include +#include +#include + namespace { /// square @@ -149,9 +153,16 @@ Surface ActsGeometry::get_tpc_surface_from_coords( Acts::Vector3 world, TrkrDefs::subsurfkey& subsurfkey) const { + // Assume that the world coordinates are in the sPHENIX frame, where the TPC is tilted + // We convert the position to tpc envelope coordinates, where we know where everything is + Acts::Vector3 world_envelope = transformTpcWorldToEnvelope(world); + double world_phi = atan2(world_envelope[1], world_envelope[0]); + unsigned int layer = TrkrDefs::getLayer(hitsetkey); unsigned int side = TpcDefs::getSide(hitsetkey); - + unsigned int sector = TpcDefs::getSectorId(hitsetkey); + + // returns an iterator to all of the surfaces for this layer auto mapIter = m_surfMaps.m_tpcSurfaceMap.find(layer); if (mapIter == m_surfMaps.m_tpcSurfaceMap.end()) @@ -160,11 +171,42 @@ Surface ActsGeometry::get_tpc_surface_from_coords( << hitsetkey << std::endl; return nullptr; } - double world_phi = atan2(world[1], world[0]); const auto& surf_vec = mapIter->second; unsigned int surf_index = 999; + // Apparently, tilting the TPC leads to the surfaces not being sorted in phi in the outer layers + // just test all surfaces in each layer for now + for(unsigned int isurf = 0; isurf < surf_vec.size(); ++isurf) + { + Surface this_surf = surf_vec[isurf]; + auto surf_center = this_surf->center(m_tGeometry.getGeoContext()); + surf_center /= 10.0; // convert from mm to cm + //this surface center includes the TPC tilt used in PHG4TpcDetector construction, transform it to tpc envelope coordinates + Acts::Vector3 surf_center_envelope = transformTpcWorldToEnvelope(surf_center); + double surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]); + double surfStepPhi = m_tGeometry.tpcSurfStepPhi; + + if ((world_phi > surf_phi - surfStepPhi / 2.0) && (world_phi < surf_phi + surfStepPhi / 2.0)) + { + if(surf_center.z() < 0 && side != 0) { continue; } + if(surf_center.z() > 0 && side != 1) { continue; } + surf_index = isurf; + subsurfkey = isurf; + break; + } + } + + if(surf_index == 999) + { + std::cout << "Error: surface not found in ActsGeometry::get_tpc_surface_from_coords " + << " layer " << layer << " side " << side << " sector " << sector << " world_phi " << world_phi << " world[0] " << world[0] << " world[1] " << world[1] << " hitsetkey " << hitsetkey << std::endl; + return nullptr; + } + + return surf_vec[surf_index]; + + /* // Predict which surface index this phi and side will correspond to // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi // we use TPC side from the hitsetkey, since z can be either sign in north and south, depending on crossing @@ -172,33 +214,38 @@ Surface ActsGeometry::get_tpc_surface_from_coords( double rounded_nsurf = std::round((double) (surf_vec.size() / 2) * fraction - 0.5); // NOLINT unsigned int nsurfm = (unsigned int) rounded_nsurf; - + std::cout << " surf_vec.size " << surf_vec.size() << " rounded_nsurf " << rounded_nsurf << " initial nsurfm " << nsurfm << std::endl; + if (side == 0) { nsurfm += surf_vec.size() / 2; } unsigned int nsurf = nsurfm % surf_vec.size(); Surface this_surf = surf_vec[nsurf]; - //std::cout << " world_phi " << world_phi << " fraction " << fraction << " rounded_nsurf " << rounded_nsurf << " nsurfm " << nsurfm << " nsurf " << nsurf << std::endl; - - auto vec3d = this_surf->center(m_tGeometry.getGeoContext()); - std::vector surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm - double surf_phi = atan2(surf_center[1], surf_center[0]); + auto surf_center = this_surf->center(m_tGeometry.getGeoContext()); + surf_center /= 10.0; // convert from mm to cm + //this surface center is from the default geometry, which includes the TPC tilt used in PHG4TpcDetector construction + // transform it to tpc envelope coordinates + Acts::Vector3 surf_center_envelope = m_tpc_world_envelope_transform * surf_center; + + double surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]); double surfStepPhi = m_tGeometry.tpcSurfStepPhi; - // std::cout << " surf_phi " << surf_phi << " surfStepPhi " << surfStepPhi << " nsurf " << nsurf << std::endl; - if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0)) + if ((world_phi > surf_phi - surfStepPhi / 2.0) && (world_phi < surf_phi + surfStepPhi / 2.0)) { surf_index = nsurf; subsurfkey = nsurf; + std::cout << "success, found nsurf = " << nsurf << std::endl; } else { // check for the periodic boundary condition auto firstsurf = *surf_vec.begin(); - auto firstsurfcenter = firstsurf->center(geometry().getGeoContext()); - float firstsurf_phi = atan2(firstsurfcenter[1], firstsurfcenter[0]); - if (world_phi < firstsurf_phi - surfStepPhi / 2.0) + auto firstsurfcenter = firstsurf->center(m_tGeometry.getGeoContext()); + firstsurfcenter /= 10.0; + auto firstsurfcenter_envelope = m_tpc_world_envelope_transform * firstsurfcenter; + double firstsurf_phi = atan2(firstsurfcenter_envelope[1], firstsurfcenter_envelope[0]); + if (world_phi < -M_PI) { world_phi += 2.0 * M_PI; } @@ -211,12 +258,16 @@ Surface ActsGeometry::get_tpc_surface_from_coords( } unsigned int new_nsurf = (nsurf+i) % surf_vec.size(); this_surf = surf_vec[new_nsurf]; - vec3d = this_surf->center(geometry().getGeoContext()); - surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm - surf_phi = atan2(surf_center[1], surf_center[0]); - //std::cout << " new world_phi " << world_phi << " new surf_phi " << surf_phi << " new_nsurf " << new_nsurf << std::endl; - if ((world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0)) + surf_center = this_surf->center(m_tGeometry.getGeoContext()); + surf_center /= 10.0; + surf_center_envelope = m_tpc_world_envelope_transform * surf_center; + surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]); + double this_philow = surf_phi - surfStepPhi / 2.0; + double this_phihigh = surf_phi + surfStepPhi / 2.0; + + if ((world_phi > this_philow) && (world_phi < this_phihigh)) { + std::cout << "success, found nsurf = " << new_nsurf << std::endl; surf_index = new_nsurf; subsurfkey = new_nsurf; return surf_vec[surf_index]; @@ -224,8 +275,9 @@ Surface ActsGeometry::get_tpc_surface_from_coords( } return nullptr; } + */ + - return surf_vec[surf_index]; } //________________________________________________________________________________________________ @@ -291,3 +343,17 @@ Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* c return local; } + + Acts::Vector3 ActsGeometry::transformTpcWorldToEnvelope(Acts::Vector3 world) const + { + Acts::Vector3 envelope = m_tpc_world_envelope_transform * world; + + return envelope; + } + + Acts::Vector3 ActsGeometry::transformTpcEnvelopeToWorld(Acts::Vector3 envelope) const + { + Acts::Vector3 world = m_tpc_world_envelope_transform.inverse() * envelope; + + return world; + } diff --git a/offline/packages/trackbase/ActsGeometry.h b/offline/packages/trackbase/ActsGeometry.h index df391d98bb..d92bf4af38 100644 --- a/offline/packages/trackbase/ActsGeometry.h +++ b/offline/packages/trackbase/ActsGeometry.h @@ -51,6 +51,7 @@ class ActsGeometry void set_CM_halfwidth(double val) { _CM_halfwidth = val; } void set_tpc_tzero(double tz) { _tpc_tzero = tz; } void set_sampa_tzero_bias(double tzb) { _sampa_tzero_bias = tzb; } + void set_tpc_world_envelope_transform(Acts::Transform3 transf) { m_tpc_world_envelope_transform = transf; } double get_tpc_tzero() const { return _tpc_tzero; } double get_sampa_tzero_bias() const { return _sampa_tzero_bias; } @@ -77,12 +78,17 @@ class ActsGeometry Acts::Transform3 makeAffineTransform(Acts::Vector3 rotation, Acts::Vector3 translation) const; + Acts::Vector3 transformTpcWorldToEnvelope(Acts::Vector3 vin) const ; + Acts::Vector3 transformTpcEnvelopeToWorld(Acts::Vector3 vin) const ; + Acts::Vector2 getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster) const; Acts::Vector2 getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster, short int crossing) const; private: ActsTrackingGeometry m_tGeometry; ActsSurfaceMaps m_surfMaps; + Acts::Transform3 m_tpc_world_envelope_transform; + Acts::Transform3 m_tpc_envelope_world_transform; double _drift_velocity = 8.0e-3; // cm/ns double _max_driftlength = 102.235; // cm double _CM_halfwidth = 0.28; // cm diff --git a/offline/packages/trackbase/AlignmentTransformation.cc b/offline/packages/trackbase/AlignmentTransformation.cc index 171d9684ac..78078b0fbc 100644 --- a/offline/packages/trackbase/AlignmentTransformation.cc +++ b/offline/packages/trackbase/AlignmentTransformation.cc @@ -254,14 +254,14 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) unsigned int side = TpcDefs::getSide(hitsetkey); unsigned int sector = TpcDefs::getSectorId(hitsetkey); - // std::cout << "New module hitsetkey " << hitsetkey << "test_layer " << test_layer << " side " << side << " sector " << sector << " nlayers " << nlayers << " layer_begin " << layer_begin << std::endl; + // std::cout << "New module hitsetkey " << hitsetkey << "test_layer " << test_layer << " side " << side << " sector " << sector << " nlayers " << nlayers << " layer_begin " << layer_begin << std::endl; // loop over layers in module for (unsigned int this_layer = layer_begin; this_layer < layer_begin + nlayers; ++this_layer) { TrkrDefs::hitsetkey this_hitsetkey = TpcDefs::genHitSetKey(this_layer, sector, side); - // std::cout << " *** module hitsetkey " << hitsetkey << " this_hitsetkey " << this_hitsetkey << " this layer " << this_layer << " side " << side << " sector " << sector << std::endl; + // std::cout << " *** module hitsetkey " << hitsetkey << " this_hitsetkey " << this_hitsetkey << " this layer " << this_layer << " side " << side << " sector " << sector << std::endl; // is this correct?????? int subsurfkey_min = (1 - side) * 144 + (144 - sector * 12) - 12 - 6; @@ -643,26 +643,24 @@ void AlignmentTransformation::extractModuleCenterPositions() for (int isector = 0; isector < 12; ++isector) { - double sectorphi = sectorPhi[iside][iregion]; + double sectorphi = sectorPhi[iside][isector]; TrkrDefs::hitsetkey hitsetkey_in = TpcDefs::genHitSetKey(lin, isector, iside); - if (localVerbosity) - { - std::cout << " hitsetkey_in " << hitsetkey_in << " lin " << lin << " sector " << isector << " side " << iside << " region " << iregion << std::endl; - } - double surf_rad_in = extractModuleCenter(hitsetkey_in, sectorphi); + double surf_rad_in = extractModuleCenter(hitsetkey_in, sectorphi); TrkrDefs::hitsetkey hitsetkey_out = TpcDefs::genHitSetKey(lout, isector, iside); - double surf_rad_out = extractModuleCenter(hitsetkey_out, sectorphi); + double surf_rad_out = extractModuleCenter(hitsetkey_out, sectorphi); + double mod_radius = (surf_rad_in + surf_rad_out) / 2.0; - TpcModuleRadii[iside][isector][iregion] = mod_radius; - if (localVerbosity) - { - std::cout << " hitsetkey_out " << hitsetkey_out << " lout " << lout << " sector " << isector << " side " << iside - << " region " << iregion << " module radius " << mod_radius << std::endl; - } + if (localVerbosity) + { + std::cout << " hitsetkey_in " << hitsetkey_in << " lin " << lin << " sector " << isector << " side " << iside << " region " << iregion << std::endl; + std::cout << " hitsetkey_out " << hitsetkey_out << " lout " << lout << " sector " << isector << " side " << iside << " region " << iregion << std::endl; + std::cout << " module radius " << mod_radius << std::endl; + } + } } } @@ -670,7 +668,7 @@ void AlignmentTransformation::extractModuleCenterPositions() double AlignmentTransformation::extractModuleCenter(TrkrDefs::hitsetkey hitsetkey, double sectorphi) { - // We want the module center position from the ideal geometry + // We want the module center position from the ideal geometry in the tpc envelope frame // the radius and z are not used, only the phi value double x = std::cos(sectorphi + 0.01) * 10.0; @@ -680,7 +678,12 @@ double AlignmentTransformation::extractModuleCenter(TrkrDefs::hitsetkey hitsetke Acts::Vector3 world(x, y, z); TrkrDefs::subsurfkey subsurfkey = 0; - Surface surface = m_tGeometry->get_tpc_surface_from_coords(hitsetkey, world, subsurfkey); + // std::cout << "extractModuleCenter: sectorphi " << sectorphi << " world " << world(0) << " " << world(1) << " " << world(2) << std::endl; + + // Note: the "world" position here is in pre-tilt tpc envelope coordinates, not global coordinates + // But, get_tpc_surface_from_coords() expects a global position as input, so we convert to world coordinates + Acts::Vector3 world_envelope = m_tGeometry->transformTpcEnvelopeToWorld(world); + Surface surface = m_tGeometry->get_tpc_surface_from_coords(hitsetkey, world_envelope, subsurfkey); if (!surface) { std::cout << PHWHERE << "Failed to find surface, quit " << std::endl; @@ -689,7 +692,9 @@ double AlignmentTransformation::extractModuleCenter(TrkrDefs::hitsetkey hitsetke Eigen::Vector3d surf_center = surface->center(m_tGeometry->geometry().getGeoContext()); surf_center /= 10.0; // convert from mm to cm - double surf_radius = std::sqrt(surf_center[0] * surf_center[0] + surf_center[1] * surf_center[1]); + // convert to tpc envelope coords + Acts::Vector3 surf_center_envelope = m_tGeometry->transformTpcWorldToEnvelope(surf_center); + double surf_radius = std::sqrt(surf_center_envelope[0] * surf_center_envelope[0] + surf_center_envelope[1] * surf_center_envelope[1]); return surf_radius; } diff --git a/offline/packages/trackbase_historic/TrackAnalysisUtils.cc b/offline/packages/trackbase_historic/TrackAnalysisUtils.cc index f7dbbb2b80..49bbd418bb 100644 --- a/offline/packages/trackbase_historic/TrackAnalysisUtils.cc +++ b/offline/packages/trackbase_historic/TrackAnalysisUtils.cc @@ -308,13 +308,14 @@ namespace TrackAnalysisUtils TpcGlobalPositionWrapper globalWrapper; globalWrapper.loadNodes(topNode); globalWrapper.set_suppressCrossing(true); + + auto* geometry = findNode::getClass(topNode, "ActsGeometry"); + TpcClusterMover mover; auto* tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); - mover.initialize_geometry(tpccellgeo); + mover.initialize_geometry(tpccellgeo, geometry); mover.set_verbosity(0); - auto* geometry = findNode::getClass(topNode, "ActsGeometry"); - std::vector> global_raw; for (const auto& key : get_cluster_keys(track)) diff --git a/offline/packages/trackreco/MakeActsGeometry.cc b/offline/packages/trackreco/MakeActsGeometry.cc index 0afea11b6c..7e3976aa68 100644 --- a/offline/packages/trackreco/MakeActsGeometry.cc +++ b/offline/packages/trackreco/MakeActsGeometry.cc @@ -180,10 +180,47 @@ int MakeActsGeometry::InitRun(PHCompositeNode *topNode) PHG4TpcGeom *layergeom = m_geomContainerTpc->GetLayerCellGeom(20); // z geometry is the same for all layers m_max_driftlength = layergeom->get_max_driftlength(); - m_CM_halfwidth = layergeom->get_CM_halfwidth(); - + m_CM_halfwidth = layergeom->get_CM_halfwidth(); m_maxSurfZ = m_max_driftlength - 0.0001; // add clearance from physical TPC gas volume length to avoid overlaps - + + // This transform will eventually be built using the tilt and placement variables that will be in layergeom + // TPC envelope to global transformation + + double rot_x = layergeom->get_rot_x(); + double rot_y = layergeom->get_rot_y(); + double rot_z = layergeom->get_rot_z(); + double place_x = layergeom->get_place_x(); + double place_y = layergeom->get_place_y(); + double place_z = layergeom->get_place_z(); + Eigen::Vector3d rot(rot_x, rot_y, rot_z); + Eigen::Vector3d trans(place_x, place_y, place_z); + + Eigen::AngleAxisd alpha(rot(0), Eigen::Vector3d::UnitX()); + Eigen::AngleAxisd beta(rot(1), Eigen::Vector3d::UnitY()); + Eigen::AngleAxisd gamma(rot(2), Eigen::Vector3d::UnitZ()); + Eigen::Quaternion q = gamma * beta * alpha; + m_tpc_envelope_world_transform.linear() = q.matrix(); + m_tpc_envelope_world_transform.translation() = trans; + // and the inverse + m_tpc_world_envelope_transform = m_tpc_envelope_world_transform.inverse(); + + // test + Acts::Vector3 test_env(10.0, 40.0, 80.0); + std::cout << "MakeActsGeometry::InitRun transform tests north" << std::endl; + std::cout << " test_env " << test_env.x() << " " << test_env.y() << " " << test_env.z() << std::endl; + Acts::Vector3 test_glob = m_tpc_envelope_world_transform * test_env; + std::cout << " test_glob " << test_glob.x() << " " << test_glob.y() << " " << test_glob.z() << std::endl; + Acts::Vector3 test_env_check = m_tpc_world_envelope_transform * test_glob; + std::cout << " test_env_check " << test_env_check.x() << " " << test_env_check.y() << " " << test_env_check.z() << std::endl; + + Acts::Vector3 test_envs(10.0, 40.0, -80.0); + std::cout << "MakeActsGeometry::InitRun transform tests south" << std::endl; + std::cout << " test_env " << test_envs.x() << " " << test_envs.y() << " " << test_envs.z() << std::endl; + Acts::Vector3 test_globs = m_tpc_envelope_world_transform * test_envs; + std::cout << " test_glob " << test_globs.x() << " " << test_globs.y() << " " << test_globs.z() << std::endl; + Acts::Vector3 test_env_checks = m_tpc_world_envelope_transform * test_globs; + std::cout << " test_env_check " << test_env_checks.x() << " " << test_env_checks.y() << " " << test_env_checks.z() << std::endl; + // Alignment Transformation declaration of instance - must be here to set initial alignment flag AlignmentTransformation alignment_transformation; alignment_transformation.createAlignmentTransformContainer(topNode); @@ -298,6 +335,7 @@ int MakeActsGeometry::InitRun(PHCompositeNode *topNode) m_actsGeometry->set_CM_halfwidth(m_CM_halfwidth); m_actsGeometry->set_tpc_tzero(m_tpc_tzero); m_actsGeometry->set_sampa_tzero_bias(m_sampa_tzero_bias); + m_actsGeometry->set_tpc_world_envelope_transform(m_tpc_world_envelope_transform); // transform world position to TPC envelope position // alignment_transformation.useInttSurveyGeometry(m_inttSurvey); if (Verbosity() > 1) @@ -846,11 +884,12 @@ void MakeActsGeometry::makeTpcMapPairs(TrackingVolumePtr &tpcVolume) { auto surf = j->getSharedPtr(); auto vec3d = surf->center(m_geoCtxt); - - // convert to cm - std::vector world_center = {vec3d(0) / 10.0, - vec3d(1) / 10.0, - vec3d(2) / 10.0}; + vec3d /= 10.0; + auto vec3d_envelope = m_tpc_world_envelope_transform * vec3d; // needs to be in TPC envelope coordinates due to tilt in sims + + std::vector world_center = {vec3d_envelope(0), + vec3d_envelope(1), + vec3d_envelope(2)}; TrkrDefs::hitsetkey hitsetkey = getTpcHitSetKeyFromCoords(world_center); unsigned int layer = TrkrDefs::getLayer(hitsetkey); @@ -1210,7 +1249,8 @@ void MakeActsGeometry::makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume) TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector &world) { - // Look up TPC surface index values from world position of surface center + // the input position is assumed to be in tpc envelope coords - i.e. tilt removed + // Look up TPC surface index values from tpc envelope position of surface center // layer unsigned int layer = 999; double layer_rad = sqrt(pow(world[0], 2) + pow(world[1], 2)); @@ -1220,6 +1260,7 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector= tpc_ref_radius_low && layer_rad < tpc_ref_radius_high) { layer = ilayer; @@ -1248,7 +1289,7 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector 3 && layer == 7) + { + std::cout << " layer_rad " << layer_rad << " m_layerRadius[layer] " << m_layerRadius[layer-7] << " found layer " << layer << " side " << side << " world " << world[0] << " " << world[1] << " " << world[2] << " phi_world " << phi_world << " readout_mod " << readout_mod << std::endl; + } + if (readout_mod >= m_nTpcModulesPerLayer) { std::cout << PHWHERE @@ -1271,16 +1318,13 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector 3) + if (Verbosity() > 3 && layer == 7) { - if (layer == 30) - { - std::cout << " world = " << world[0] << " " << world[1] - << " " << world[2] << " phi_world " - << phi_world * 180 / M_PI << " layer " << layer - << " readout_mod " << readout_mod << " side " << side - << " hitsetkey " << hitset_key << std::endl; - } + std::cout << " world = " << world[0] << " " << world[1] + << " " << world[2] << " phi_world " + << phi_world * 180 / M_PI << " layer " << layer + << " readout_mod " << readout_mod << " side " << side + << " hitsetkey " << hitset_key << std::endl; } return hitset_key; diff --git a/offline/packages/trackreco/MakeActsGeometry.h b/offline/packages/trackreco/MakeActsGeometry.h index 272182b713..4180db74f3 100644 --- a/offline/packages/trackreco/MakeActsGeometry.h +++ b/offline/packages/trackreco/MakeActsGeometry.h @@ -303,6 +303,10 @@ class MakeActsGeometry : public SubsysReco bool m_use_module_tilt_always = false; bool m_use_new_silicon_rotation_order = false; + + Acts::Transform3 m_tpc_world_envelope_transform; + Acts::Transform3 m_tpc_envelope_world_transform; + }; #endif diff --git a/offline/packages/trackreco/MakeSourceLinks.cc b/offline/packages/trackreco/MakeSourceLinks.cc index f2975ac952..74c9300d97 100644 --- a/offline/packages/trackreco/MakeSourceLinks.cc +++ b/offline/packages/trackreco/MakeSourceLinks.cc @@ -50,12 +50,12 @@ namespace } // namespace -void MakeSourceLinks::initialize(PHG4TpcGeomContainer* cellgeo) +void MakeSourceLinks::initialize(PHG4TpcGeomContainer* cellgeo, ActsGeometry *tGeometry) { // get the TPC layer radii from the geometry object - if (cellgeo) + if (cellgeo && tGeometry) { - _clusterMover.initialize_geometry(cellgeo); + _clusterMover.initialize_geometry(cellgeo, tGeometry); } } diff --git a/offline/packages/trackreco/MakeSourceLinks.h b/offline/packages/trackreco/MakeSourceLinks.h index 7cabf268dc..62b5031624 100644 --- a/offline/packages/trackreco/MakeSourceLinks.h +++ b/offline/packages/trackreco/MakeSourceLinks.h @@ -42,7 +42,7 @@ class MakeSourceLinks public: MakeSourceLinks() = default; - void initialize(PHG4TpcGeomContainer* cellgeo); + void initialize(PHG4TpcGeomContainer* cellgeo, ActsGeometry *tGeometry); void setVerbosity(int verbosity) { m_verbosity = verbosity; } diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index 1cf72bdfdc..f7a58fb813 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -436,7 +436,7 @@ void PHActsTrkFitter::loopTracks(Acts::Logging::Level logLevel) SourceLinkVec sourceLinks; MakeSourceLinks makeSourceLinks; - makeSourceLinks.initialize(_tpccellgeo); + makeSourceLinks.initialize(_tpccellgeo, m_tGeometry); makeSourceLinks.setVerbosity(Verbosity()); makeSourceLinks.set_pp_mode(m_pp_mode); makeSourceLinks.set_cluster_edge_rejection(m_cluster_edge_rejection); diff --git a/offline/packages/trackreco/PHCosmicsTrkFitter.cc b/offline/packages/trackreco/PHCosmicsTrkFitter.cc index fc6de17328..2c173d2ed2 100644 --- a/offline/packages/trackreco/PHCosmicsTrkFitter.cc +++ b/offline/packages/trackreco/PHCosmicsTrkFitter.cc @@ -297,7 +297,7 @@ void PHCosmicsTrkFitter::loopTracks(Acts::Logging::Level logLevel) SourceLinkVec sourceLinks; MakeSourceLinks makeSourceLinks; - makeSourceLinks.initialize(_tpccellgeo); + makeSourceLinks.initialize(_tpccellgeo, m_tGeometry); makeSourceLinks.setVerbosity(Verbosity()); makeSourceLinks.set_pp_mode(false); @@ -1116,4 +1116,4 @@ Acts::Vector3 PHCosmicsTrkFitter::calculateMomentum(TrackSeed* tpcseed, const st momentum.z() = pz; return momentum; -} \ No newline at end of file +} diff --git a/simulation/g4simulation/g4detectors/Makefile.am b/simulation/g4simulation/g4detectors/Makefile.am index 9e219fe9e0..acff3e0c83 100644 --- a/simulation/g4simulation/g4detectors/Makefile.am +++ b/simulation/g4simulation/g4detectors/Makefile.am @@ -100,6 +100,7 @@ pkginclude_HEADERS = \ PHG4TpcCylinderGeomContainer.h \ PHG4TpcGeom.h \ PHG4TpcGeomv1.h \ + PHG4TpcGeomv2.h \ PHG4TpcGeomContainer.h \ PHG4ZDCDefs.h \ PHG4ZDCSubsystem.h @@ -137,6 +138,7 @@ ROOTDICTS = \ PHG4TpcCylinderGeomContainer_Dict.cc \ PHG4TpcGeom_Dict.cc \ PHG4TpcGeomv1_Dict.cc \ + PHG4TpcGeomv2_Dict.cc \ PHG4TpcGeomContainer_Dict.cc pcmdir = $(libdir) @@ -177,6 +179,7 @@ libg4detectors_io_la_SOURCES = \ PHG4TpcCylinderGeomContainer.cc \ PHG4TpcGeom.cc \ PHG4TpcGeomv1.cc \ + PHG4TpcGeomv2.cc \ PHG4TpcGeomContainer.cc libg4detectors_la_SOURCES = \ diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeom.h b/simulation/g4simulation/g4detectors/PHG4TpcGeom.h index bdb3bfbae7..521cf2fa9c 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeom.h +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeom.h @@ -192,6 +192,38 @@ class PHG4TpcGeom : public PHObject return -99999; } + virtual double get_rot_x() const + { + PHOOL_VIRTUAL_WARN("get_rot_x()"); + return std::numeric_limits::quiet_NaN(); + } + virtual double get_rot_y() const + { + PHOOL_VIRTUAL_WARN("get_rot_y()"); + return std::numeric_limits::quiet_NaN(); + } + virtual double get_rot_z() const + { + PHOOL_VIRTUAL_WARN("get_rot_z()"); + return std::numeric_limits::quiet_NaN(); + } + virtual double get_place_x() const + { + PHOOL_VIRTUAL_WARN("get_place_x()"); + return std::numeric_limits::quiet_NaN(); + } + virtual double get_place_y() const + { + PHOOL_VIRTUAL_WARN("get_place_y()"); + return std::numeric_limits::quiet_NaN(); + } + virtual double get_place_z() const + { + PHOOL_VIRTUAL_WARN("get_place_z()"); + return std::numeric_limits::quiet_NaN(); + } + + virtual const std::array, 2> &get_sector_min_phi(); virtual const std::array, 2> &get_sector_max_phi(); @@ -211,6 +243,17 @@ class PHG4TpcGeom : public PHObject { PHOOL_VIRTUAL_WARN("set_phi_bias(const std::array, 2>&)"); } + + + + /* + double get_rot_x() const override { return rot_x; } + double get_rot_y() const override { return rot_y; } + double get_rot_z() const override { return rot_z; } + double get_place_x() const override { return place_x; } + double get_place_y() const override { return place_y; } + double get_place_z() const override { return place_z; } + */ virtual void set_layer(const int) { PHOOL_VIRTUAL_WARN("set_layer(const int)"); } virtual void set_radius(const double) { PHOOL_VIRTUAL_WARN("set_radius(const double)"); } @@ -237,7 +280,15 @@ class PHG4TpcGeom : public PHObject virtual void set_adc_clock(const double) { PHOOL_VIRTUAL_WARN("set_adc_clock(const double)"); } virtual void set_extended_readout_time(const double) { PHOOL_VIRTUAL_WARN("set_extended_readout_time(const double)"); } virtual void set_drift_velocity_sim(const double) { PHOOL_VIRTUAL_WARN("set_drift_velocity_sim(const double)"); } - + + virtual void set_rot_x(const double) { PHOOL_VIRTUAL_WARN("set_rot_x(const double)"); } + virtual void set_rot_y(const double) { PHOOL_VIRTUAL_WARN("set_rot_y(const double)"); } + virtual void set_rot_z(const double) { PHOOL_VIRTUAL_WARN("set_rot_z(const double)"); } + + virtual void set_place_x(const double) { PHOOL_VIRTUAL_WARN("set_place_x(const double)"); } + virtual void set_place_y(const double) { PHOOL_VIRTUAL_WARN("set_place_y(const double)"); } + virtual void set_place_z(const double) { PHOOL_VIRTUAL_WARN("set_place_z(const double)"); } + //! load parameters from PHParameters, which interface to Database/XML/ROOT files virtual void ImportParameters(const PHParameters & /*param*/) { return; } diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc new file mode 100644 index 0000000000..050baf86a9 --- /dev/null +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc @@ -0,0 +1,593 @@ +#include "PHG4TpcGeomv2.h" +#include "PHG4CylinderCellDefs.h" + +#include + +#include + +namespace +{ + // streamer for internal 2dimensional arrays + using array_t = std::array, PHG4TpcGeomv2::NSides>; + std::ostream& operator<<(std::ostream& out, const array_t& array) + { + out << "{ "; + for (const auto& iside : array) + { + out << "{"; + bool first = true; + for (const auto& value : iside) + { + if (!first) + { + out << ", "; + } + first = false; + out << value; + } + out << "} "; + } + out << " }"; + return out; + } +} // namespace + +std::ostream& operator<<(std::ostream& out, const PHG4TpcGeomv2& geom) +{ + out << "PHG4TpcGeomv2 - layer: " << geom.layer << std::endl; + out + << " binnig: " << geom.binning + << ", radius: " << geom.radius + << ", nzbins: " << geom.nzbins + << ", zmin: " << geom.zmin + << ", zstep: " << geom.zstep + << ", nphibins: " << geom.nphibins + << ", phimin: " << geom.phimin + << ", phistep: " << geom.phistep + << ", thickness: " << geom.thickness + << std::endl; + + out << " sector_R_bias: " << geom.sector_R_bias << std::endl; + out << " sector_Phi_bias: " << geom.sector_Phi_bias << std::endl; + out << " sector_min_Phi: " << geom.sector_min_Phi << std::endl; + out << " sector_max_Phi: " << geom.sector_max_Phi << std::endl; + + return out; +} + +void PHG4TpcGeomv2::set_zbins(const int i) +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + nzbins = i; +} + +void PHG4TpcGeomv2::set_zmin(const double z) +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + zmin = z; +} + +int PHG4TpcGeomv2::get_zbins() const +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return nzbins; +} + +double +PHG4TpcGeomv2::get_zmin() const +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return zmin; +} + +double +PHG4TpcGeomv2::get_zstep() const +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return zstep; +} + +void PHG4TpcGeomv2::set_zstep(const double z) +{ + check_binning_method(PHG4CylinderCellDefs::sizebinning); + zstep = z; +} + +int PHG4TpcGeomv2::get_phibins() const +{ + check_binning_method_phi("PHG4TpcGeomv2::get_phibins"); + return nphibins; +} + +double +PHG4TpcGeomv2::get_phistep() const +{ + check_binning_method_phi("PHG4TpcGeomv2::get_phistep"); + return phistep; +} + +double +PHG4TpcGeomv2::get_phimin() const +{ + check_binning_method_phi("PHG4TpcGeomv2::get_phimin"); + return phimin; +} + +void PHG4TpcGeomv2::set_phibins(const int i) +{ + check_binning_method_phi("PHG4TpcGeomv2::set_phibins"); + nphibins = i; +} + +void PHG4TpcGeomv2::set_phistep(const double phi) +{ + check_binning_method_phi("PHG4TpcGeomv2::set_phistep"); + phistep = phi; +} + +void PHG4TpcGeomv2::set_phimin(const double phi) +{ + check_binning_method_phi("PHG4TpcGeomv2::set_phimin"); + phimin = phi; +} + +int PHG4TpcGeomv2::get_etabins() const +{ + check_binning_method_eta("PHG4TpcGeomv2::get_etabins"); + return nzbins; +} + +double +PHG4TpcGeomv2::get_etastep() const +{ + check_binning_method_eta("PHG4TpcGeomv2::get_etastep"); + return zstep; +} +double +PHG4TpcGeomv2::get_etamin() const +{ + check_binning_method_eta("PHG4TpcGeomv2::get_etamin"); + return zmin; +} + +void PHG4TpcGeomv2::set_etamin(const double z) +{ + check_binning_method_eta("PHG4TpcGeomv2::set_etamin"); + zmin = z; +} + +void PHG4TpcGeomv2::set_etastep(const double z) +{ + check_binning_method_eta("PHG4TpcGeomv2::set_etastep"); + zstep = z; +} + +void PHG4TpcGeomv2::set_etabins(const int i) +{ + check_binning_method_eta("PHG4TpcGeomv2::set_etabins"); + nzbins = i; +} + +void PHG4TpcGeomv2::identify(std::ostream& os) const +{ + os << "PHG4TpcGeomv2::identify - layer: " << layer << std::endl; + + os + << " binning: " << binning + << ", radius: " << radius + << ", nzbins: " << nzbins + << ", zmin: " << zmin + << ", zstep: " << zstep + << ", nphibins: " << nphibins + << ", phimin: " << phimin + << ", phistep: " << phistep + << ", thickness: " << thickness + << std::endl; + + os << " sector_R_bias: " << sector_R_bias << std::endl; + os << " sector_Phi_bias: " << sector_Phi_bias << std::endl; + os << " sector_min_Phi: " << sector_min_Phi << std::endl; + os << " sector_max_Phi: " << sector_max_Phi << std::endl; +} + +std::pair +PHG4TpcGeomv2::get_zbounds(const int ibin) const +{ + if (ibin < 0 || ibin > nzbins) + { + std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl; + exit(1); + } + check_binning_method(PHG4CylinderCellDefs::sizebinning); + double zlow = zmin + ibin * zstep; + double zhigh = zlow + zstep; + return std::make_pair(zlow, zhigh); +} + +std::pair +PHG4TpcGeomv2::get_etabounds(const int ibin) const +{ + if (ibin < 0 || ibin > nzbins) + { + std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl; + exit(1); + } + check_binning_method_eta("PHG4TpcGeomv2::get_etabounds"); + // check_binning_method(PHG4CylinderCellDefs::etaphibinning); + double zlow = zmin + ibin * zstep; + double zhigh = zlow + zstep; + return std::make_pair(zlow, zhigh); +} + +std::pair +PHG4TpcGeomv2::get_phibounds(const int ibin) const +{ + if (ibin < 0 || ibin > nphibins) + { + std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; + exit(1); + } + + double philow = phimin + ibin * phistep; + double phihigh = philow + phistep; + return std::make_pair(philow, phihigh); +} + +int PHG4TpcGeomv2::get_zbin(const double z) const +{ + if (z < zmin || z > (zmin + nzbins * zstep)) + { + // cout << PHWHERE << "Asking for bin for z outside of z range: " << z << endl; + return -1; + } + + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return floor((z - zmin) / zstep); +} + +int PHG4TpcGeomv2::get_etabin(const double eta) const +{ + if (eta < zmin || eta > (zmin + nzbins * zstep)) + { + // cout << "Asking for bin for eta outside of eta range: " << eta << endl; + return -1; + } + check_binning_method_eta(); + return floor((eta - zmin) / zstep); +} + +int PHG4TpcGeomv2::get_phibin_new(const double phi) const +{ + double norm_phi = phi; + if (phi < phimin || phi > (phimin + nphibins * phistep)) + { + int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); + norm_phi += 2 * M_PI * nwraparound; + } + check_binning_method_phi(); + return floor((norm_phi - phimin) / phistep); +} + +int PHG4TpcGeomv2::find_phibin(const double phi, int side) const +{ + double norm_phi = phi; + if (phi < phimin || phi > (phimin + nphibins * phistep)) + { + int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); + norm_phi += 2 * M_PI * nwraparound; + } + // if (phi > M_PI){ + // norm_phi = phi - 2* M_PI; + // } + // if (phi < phimin){ + // norm_phi = phi + 2* M_PI; + // } + //side = 0; + + int phi_bin = -1; + + for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++) + { + if (norm_phi < sector_max_Phi[side][s] && norm_phi > sector_min_Phi[side][s]) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = (floor(std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s); + break; + } + if (s == 11) + { + if (norm_phi < sector_max_Phi[side][s] && norm_phi >= -M_PI) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = floor(std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s; + break; + } + if (norm_phi > sector_min_Phi[side][s] + 2 * M_PI) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = floor(std::abs(sector_max_Phi[side][s] - (norm_phi - 2 * M_PI)) / phistep) + nphibins / 12 * s; + break; + } + } + } + return phi_bin; +} + +float PHG4TpcGeomv2::get_pad_float(const double phi, int side) const +{ + double norm_phi = phi; + if (phi < phimin || phi > (phimin + nphibins * phistep)) + { + int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); + norm_phi += 2 * M_PI * nwraparound; + } + // if (phi > M_PI){ + // norm_phi = phi - 2* M_PI; + // } + // if (phi < phimin){ + // norm_phi = phi + 2* M_PI; + // } + //side = 0; + + float phi_bin = -1; + + for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++) + { + if (norm_phi < sector_max_Phi[side][s] && norm_phi > sector_min_Phi[side][s]) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = (std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s; + break; + } + if (s == 11) + { + if (norm_phi < sector_max_Phi[side][s] && norm_phi >= -M_PI) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = (std::abs(sector_max_Phi[side][s] - norm_phi) / phistep) + nphibins / 12 * s; + break; + } + if (norm_phi > sector_min_Phi[side][s] + 2 * M_PI) + { + // NOLINTNEXTLINE(bugprone-integer-division) + phi_bin = (std::abs(sector_max_Phi[side][s] - (norm_phi - 2 * M_PI)) / phistep) + nphibins / 12 * s; + break; + } + } + } + return phi_bin - 0.5; +} + +float PHG4TpcGeomv2::get_tbin_float(const double z) const +{ + if (z < zmin || z > (zmin + nzbins * zstep)) + { + // cout << PHWHERE << "Asking for bin for z outside of z range: " << z << endl; + return -1; + } + + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return ((z - zmin) / zstep) - 0.5; +} + +int PHG4TpcGeomv2::get_phibin(const double phi, int side) const +{ + double new_phi = phi; + if (phi > M_PI) + { + new_phi = phi - 2 * M_PI; + } + if (phi < phimin) + { + new_phi = phi + 2 * M_PI; + } + // Get phi-bin number + int phi_bin = find_phibin(new_phi, side); + + //side = 0; + // If phi-bin is not defined, check that it is in the dead area and put it to the edge of sector + if (phi_bin < 0) + { + // + for (std::size_t s = 0; s < sector_max_Phi[side].size(); s++) + { + double daPhi = 0; + if (s == 0) + { + daPhi = fabs(sector_min_Phi[side][11] + 2 * M_PI - sector_max_Phi[side][s]); + } + else + { + daPhi = fabs(sector_min_Phi[side][s - 1] - sector_max_Phi[side][s]); + } + + double min_phi = sector_max_Phi[side][s]; + double max_phi = sector_max_Phi[side][s] + daPhi; + if (new_phi <= max_phi && new_phi >= min_phi) + { + if (fabs(max_phi - new_phi) > fabs(new_phi - min_phi)) + { + new_phi = min_phi - phistep / 5; + } + else + { + new_phi = max_phi + phistep / 5; + } + } + } + // exit(1); + + phi_bin = find_phibin(new_phi, side); + if (phi_bin < 0) + { + std::cout << PHWHERE << "Asking for bin for phi outside of phi range: " << phi << std::endl; + exit(1); + // phi_bin=0; + } + } + return phi_bin; +} + +double +PHG4TpcGeomv2::get_zcenter(const int ibin) const +{ + if (ibin < 0 || ibin > nzbins) + { + std::cout << PHWHERE << "Asking for invalid bin in z: " << ibin << std::endl; + exit(1); + } + check_binning_method(PHG4CylinderCellDefs::sizebinning); + return zmin + (ibin + 0.5) * zstep; +} + +double +PHG4TpcGeomv2::get_etacenter(const int ibin) const +{ + if (ibin < 0 || ibin > nzbins) + { + std::cout << PHWHERE << "Asking for invalid bin in eta: " << ibin << std::endl; + std::cout << "minbin: 0, maxbin " << nzbins << std::endl; + exit(1); + } + check_binning_method_eta(); + return zmin + (ibin + 0.5) * zstep; +} + +double +PHG4TpcGeomv2::get_phicenter_new(const int ibin) const +{ + if (ibin < 0 || ibin > nphibins) + { + std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; + exit(1); + } + + check_binning_method_phi(); + + return (phimin + (ibin + 0.5) * phistep); +} + +double +PHG4TpcGeomv2::get_phicenter(const int ibin, const int side) const +{ + // double phi_center = -999; + if (ibin < 0 || ibin > nphibins) + { + std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; + exit(1); + } + + check_binning_method_phi(); + + //const int side = 0; + unsigned int pads_per_sector = nphibins / 12; + unsigned int sector = ibin / pads_per_sector; + double phi_center = (sector_max_Phi[side][sector] - (ibin + 0.5 - sector * pads_per_sector) * phistep); + if (phi_center <= -M_PI) + { + phi_center += 2 * M_PI; + } + return phi_center; +} + +double +PHG4TpcGeomv2::get_phi(const float ibin, const int side) const +{ + // double phi_center = -999; + if (ibin < 0 || ibin > nphibins) + { + std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; + exit(1); + } + + check_binning_method_phi(); + + //const int side = 0; + unsigned int pads_per_sector = nphibins / 12; + unsigned int sector = ibin / pads_per_sector; + double phi = (sector_max_Phi[side][sector] - (ibin + 0.5 - sector * pads_per_sector) * phistep); + if (phi <= -M_PI) + { + phi += 2 * M_PI; + } + return phi; +} + +std::string +PHG4TpcGeomv2::methodname(const int i) const +{ + switch (i) + { + case PHG4CylinderCellDefs::sizebinning: + return "Bins in cm"; + break; + case PHG4CylinderCellDefs::etaphibinning: + return "Eta/Phi bins"; + break; + case PHG4CylinderCellDefs::etaslatbinning: + return "Eta/numslat bins"; + break; + case PHG4CylinderCellDefs::spacalbinning: + return "SPACAL Tower bins"; + break; + default: + break; + } + return "Unknown"; +} + +void PHG4TpcGeomv2::check_binning_method(const int i) const +{ + if (binning != i) + { + std::cout << "different binning method used " << methodname(binning) + << ", not : " << methodname(i) + << std::endl; + exit(1); + } + return; +} + +void PHG4TpcGeomv2::check_binning_method_eta(const std::string& src) const +{ + if (binning != PHG4CylinderCellDefs::etaphibinning && + binning != PHG4CylinderCellDefs::etaslatbinning && + binning != PHG4CylinderCellDefs::spacalbinning) + { + if (!src.empty()) + { + std::cout << src << " : "; + } + + std::cout << "different binning method used " << methodname(binning) + << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning) + << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning) + << " or " << methodname(PHG4CylinderCellDefs::spacalbinning) + << std::endl; + exit(1); + } + return; +} + +void PHG4TpcGeomv2::check_binning_method_phi(const std::string& src) const +{ + if (binning != PHG4CylinderCellDefs::etaphibinning && + binning != PHG4CylinderCellDefs::sizebinning && + binning != PHG4CylinderCellDefs::etaslatbinning && + binning != PHG4CylinderCellDefs::spacalbinning) + { + if (!src.empty()) + { + std::cout << src << " : "; + } + + std::cout << "different binning method used " << methodname(binning) + << ", not : " << methodname(PHG4CylinderCellDefs::etaphibinning) + << " or " << methodname(PHG4CylinderCellDefs::sizebinning) + << " or " << methodname(PHG4CylinderCellDefs::etaslatbinning) + << " or " << methodname(PHG4CylinderCellDefs::spacalbinning) + << std::endl; + exit(1); + } + return; +} diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h new file mode 100644 index 0000000000..235256b67c --- /dev/null +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h @@ -0,0 +1,158 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef G4DETECTORS_PHG4TPCGEOMV1_H +#define G4DETECTORS_PHG4TPCGEOMV1_H + +#include "PHG4TpcGeom.h" + +#include +#include +#include // for cout, ostream +#include +#include // for pair + +class PHG4TpcGeomv2 : public PHG4TpcGeom +{ + public: + PHG4TpcGeomv2() = default; + + ~PHG4TpcGeomv2() override = default; + + // from PHObject + void identify(std::ostream& os = std::cout) const override; + + int get_layer() const override { return layer; } + double get_radius() const override { return radius; } + double get_thickness() const override { return thickness; } + int get_binning() const override { return binning; } + int get_zbins() const override; + int get_phibins() const override; + double get_zmin() const override; + double get_phistep() const override; + double get_phimin() const override; + double get_zstep() const override; + int get_etabins() const override; + double get_etastep() const override; + double get_etamin() const override; + + double get_max_driftlength() const override { return max_driftlength; } + double get_CM_halfwidth() const override { return CM_halfwidth; } + double get_adc_clock() const override { return adc_clock; } // default sim value + double get_extended_readout_time() const override { return extended_readout_time; } + double get_drift_velocity_sim() const override { return drift_velocity_sim; } + + double get_rot_x() const override { return rot_x; } + double get_rot_y() const override { return rot_y; } + double get_rot_z() const override { return rot_z; } + double get_place_x() const override { return place_x; } + double get_place_y() const override { return place_y; } + double get_place_z() const override { return place_z; } + + std::pair get_zbounds(const int ibin) const override; + std::pair get_phibounds(const int ibin) const override; + std::pair get_etabounds(const int ibin) const override; + double get_etacenter(const int ibin) const override; + double get_zcenter(const int ibin) const override; + double get_phicenter(const int ibin, const int side = 0) const override; + double get_phicenter_new(const int ibin) const override; + double get_phi(const float ibin, const int side = 0) const override; + + int get_etabin(const double eta) const override; + int get_zbin(const double z) const override; + int get_phibin(const double phi, int side = 0) const override; + int get_phibin_new(const double phi) const override; + + float get_pad_float(const double phi, int side = 0) const override; + float get_tbin_float(const double z) const override; + int find_phibin(const double phi, int side = 0) const override; + + void set_layer(const int i) override { layer = i; } + void set_binning(const int i) override { binning = i; } + void set_radius(const double r) override { radius = r; } + void set_thickness(const double t) override { thickness = t; } + void set_zbins(const int i) override; + void set_zmin(const double z) override; + void set_zstep(const double z) override; + void set_phibins(const int i) override; + void set_phistep(const double phi) override; + void set_phimin(const double phi) override; + void set_etabins(const int i) override; + void set_etamin(const double z) override; + void set_etastep(const double z) override; + // capture the z geometry related setup parameters + void set_max_driftlength(const double val) override { max_driftlength = val; } + void set_CM_halfwidth(const double val) override { CM_halfwidth = val; } + void set_adc_clock(const double val) override { adc_clock = val; } + void set_extended_readout_time(const double val) override { extended_readout_time = val; } + void set_drift_velocity_sim(const double val) override { drift_velocity_sim = val; } + void set_rot_x(const double val) override { rot_x = val; } + void set_rot_y(const double val) override { rot_y = val; } + void set_rot_z(const double val) override { rot_z = val; } + void set_place_x(const double val) override { place_x = val; } + void set_place_y(const double val) override { place_y = val; } + void set_place_z(const double val) override { place_z = val; } + + static const int NSides = 2; + + void set_r_bias(const std::array, NSides> &dr) override { sector_R_bias = dr; } + void set_phi_bias(const std::array, NSides> &dphi) override { sector_Phi_bias = dphi; } + + void set_sector_min_phi(const std::array, NSides> &s_min_phi) override + { + sector_min_Phi = s_min_phi; + } + void set_sector_max_phi(const std::array, NSides> &s_max_phi) override + { + sector_max_Phi = s_max_phi; + } + + const std::array, NSides> &get_sector_min_phi() override + { + return sector_min_Phi; + } + const std::array, NSides> &get_sector_max_phi() override + { + return sector_max_Phi; + } + + protected: + void check_binning_method(const int i) const; + void check_binning_method_eta(const std::string& src = "") const; + void check_binning_method_phi(const std::string& src = "") const; + std::string methodname(const int i) const; + int layer{-999}; + int binning{0}; + double radius{std::numeric_limits::quiet_NaN()}; + int nzbins{-1}; + double zmin{std::numeric_limits::quiet_NaN()}; + double zstep{std::numeric_limits::quiet_NaN()}; + int nphibins{-1}; + double phimin{-M_PI}; + double phistep{std::numeric_limits::quiet_NaN()}; + double thickness{std::numeric_limits::quiet_NaN()}; + + double max_driftlength{std::numeric_limits::quiet_NaN()}; + double CM_halfwidth{std::numeric_limits::quiet_NaN()}; + double adc_clock{std::numeric_limits::quiet_NaN()}; + double extended_readout_time{std::numeric_limits::quiet_NaN()}; + double drift_velocity_sim{std::numeric_limits::quiet_NaN()}; + + double rot_x{std::numeric_limits::quiet_NaN()}; + double rot_y{std::numeric_limits::quiet_NaN()}; + double rot_z{std::numeric_limits::quiet_NaN()}; + double place_x{std::numeric_limits::quiet_NaN()}; + double place_y{std::numeric_limits::quiet_NaN()}; + double place_z{std::numeric_limits::quiet_NaN()}; + + std::array, NSides> sector_R_bias; + std::array, NSides> sector_Phi_bias; + std::array, NSides> sector_min_Phi; + std::array, NSides> sector_max_Phi; + + // streamer + friend std::ostream& operator<<(std::ostream&, const PHG4TpcGeomv2&); + + ClassDefOverride(PHG4TpcGeomv2, 1) +}; + +#endif diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2LinkDef.h b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2LinkDef.h new file mode 100644 index 0000000000..e25664077d --- /dev/null +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2LinkDef.h @@ -0,0 +1,5 @@ +#ifdef __CINT__ + +#pragma link C++ class PHG4TpcGeomv2 + ; + +#endif /* __CINT__ */ diff --git a/simulation/g4simulation/g4eval/SvtxTruthEval.cc b/simulation/g4simulation/g4eval/SvtxTruthEval.cc index 2ecaccf41c..f2fcd220a7 100644 --- a/simulation/g4simulation/g4eval/SvtxTruthEval.cc +++ b/simulation/g4simulation/g4eval/SvtxTruthEval.cc @@ -457,8 +457,16 @@ void SvtxTruthEval::LayerClusterG4Hits(const std::set& truth_hits, std // we do not assume that the truth hits know what layer they are in for (auto *this_g4hit : truth_hits) { - float rbegin = std::sqrt(this_g4hit->get_x(0) * this_g4hit->get_x(0) + this_g4hit->get_y(0) * this_g4hit->get_y(0)); - float rend = std::sqrt(this_g4hit->get_x(1) * this_g4hit->get_x(1) + this_g4hit->get_y(1) * this_g4hit->get_y(1)); + // The truth hits are in world coordinates + // They have to be transformed to envelope coords to find what layer they are in + // Then the cluster positions have to be transformed back to world coordinates + Acts::Vector3 world0(this_g4hit->get_x(0), this_g4hit->get_y(0), this_g4hit->get_z(0)); + Acts::Vector3 env0 = _tgeometry->transformTpcWorldToEnvelope(world0); + Acts::Vector3 world1(this_g4hit->get_x(1), this_g4hit->get_y(1), this_g4hit->get_z(1)); + Acts::Vector3 env1 = _tgeometry->transformTpcWorldToEnvelope(world1); + + float rbegin = std::sqrt(env0.x() * env0.x() + env0.y() * env0.y()); + float rend = std::sqrt(env1.x() * env1.x() + env1.y() * env1.y()); // std::cout << " Eval: g4hit " << this_g4hit->get_hit_id() << " layer " << layer << " rbegin " << rbegin << " rend " << rend << std::endl; // make sure the entry point is at lower radius @@ -468,21 +476,21 @@ void SvtxTruthEval::LayerClusterG4Hits(const std::set& truth_hits, std if (rbegin < rend) { - xl[0] = this_g4hit->get_x(0); - yl[0] = this_g4hit->get_y(0); - zl[0] = this_g4hit->get_z(0); - xl[1] = this_g4hit->get_x(1); - yl[1] = this_g4hit->get_y(1); - zl[1] = this_g4hit->get_z(1); + xl[0] = env0.x(); + yl[0] = env0.y(); + zl[0] = env0.z(); + xl[1] = env1.x(); + yl[1] = env1.y(); + zl[1] = env1.z(); } else { - xl[0] = this_g4hit->get_x(1); - yl[0] = this_g4hit->get_y(1); - zl[0] = this_g4hit->get_z(1); - xl[1] = this_g4hit->get_x(0); - yl[1] = this_g4hit->get_y(0); - zl[1] = this_g4hit->get_z(0); + xl[0] = env1.x(); + yl[0] = env1.y(); + zl[0] = env1.z(); + xl[1] = env0.x(); + yl[1] = env0.y(); + zl[1] = env0.z(); std::swap(rbegin, rend); // std::cout << "swapped in and out " << std::endl; } @@ -659,6 +667,14 @@ void SvtxTruthEval::LayerClusterG4Hits(const std::set& truth_hits, std } } + // convert cluster position back to world coordinates + Acts::Vector3 clus_env(gx,gy,gz); + Acts::Vector3 clus_world = _tgeometry->transformTpcEnvelopeToWorld(clus_env); + gx = clus_world.x(); + gy = clus_world.y(); + gz = clus_world.z(); + // what is gr used for? + } // if TPC else { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc index 5d0eb1e66c..6f736ccc39 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc @@ -3,7 +3,7 @@ #include "PHG4TpcDisplayAction.h" #include -#include +#include #include #include @@ -118,10 +118,26 @@ void PHG4TpcDetector::ConstructMe(G4LogicalVolume *logicWorld) ConstructTpcCageVolume(tpc_envelope_logic); ConstructTpcGasVolume(tpc_envelope_logic); - new G4PVPlacement(nullptr, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm), + G4RotationMatrix rot; + rot.rotateX(m_Params->get_double_param("rot_x")*rad); + rot.rotateY(m_Params->get_double_param("rot_y")*rad); + rot.rotateZ(m_Params->get_double_param("rot_z")*rad); + + G4ThreeVector trans(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm); + + new G4PVPlacement( + G4Transform3D(rot, trans), tpc_envelope_logic, "tpc_envelope", - logicWorld, false, false, OverlapCheck()); - + logicWorld, + false, false, OverlapCheck()); + + /* + G4ThreeVector test_env(10.0, 40.0, 80.0); + std::cout << " test_env " << test_env.x() << " " << test_env.y() << " " << test_env.z() << std::endl; + G4ThreeVector test_glob = test_env.transform(rot); + std::cout << " test_glob " << test_glob.x() << " " << test_glob.y() << " " << test_glob.z() << std::endl; + */ + // geometry node add_geometry_node(); } @@ -538,7 +554,12 @@ void PHG4TpcDetector::add_geometry_node() auto *newNode = new PHIODataNode(geonode, geonode_name, "PHObject"); geomNode->addNode(newNode); } - + else + { + std::cout << "PHG4TpcGeomContainer already exists with name " << geonode_name << " and it should not! " << std::endl; + geonode->identify(); + } + m_cdb = CDBInterface::instance(); std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP"); @@ -695,7 +716,7 @@ void PHG4TpcDetector::add_geometry_node() << " phibins " << NPhiBins[iregion] << " phistep " << phi_bin_width_cdb[layer] << std::endl; } - auto *layerseggeo = new PHG4TpcGeomv1; + auto *layerseggeo = new PHG4TpcGeomv2; layerseggeo->set_layer(layer); double r_length = Thickness[iregion]; @@ -731,6 +752,12 @@ void PHG4TpcDetector::add_geometry_node() layerseggeo->set_adc_clock(m_Params->get_double_param("tpc_adc_clock")); layerseggeo->set_extended_readout_time(m_Params->get_double_param("extended_readout_time")); layerseggeo->set_drift_velocity_sim(m_Params->get_double_param("drift_velocity_sim")); + layerseggeo->set_rot_x(m_Params->get_double_param("rot_x")); + layerseggeo->set_rot_y(m_Params->get_double_param("rot_y")); + layerseggeo->set_rot_z(m_Params->get_double_param("rot_z")); + layerseggeo->set_place_x(m_Params->get_double_param("place_x")); + layerseggeo->set_place_y(m_Params->get_double_param("place_y")); + layerseggeo->set_place_z(m_Params->get_double_param("place_z")); } // Chris Pinkenburg: greater causes huge memory growth which causes problems diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index 9de4cf6639..24d753cb66 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -390,6 +390,14 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } } + /* + Acts::Vector3 test_env(10.0, 40.0, 80.0); + std::cout << " test_env " << test_env.x() << " " << test_env.y() << " " << test_env.z() << std::endl; + Acts::Vector3 test_glob = m_tGeometry-> transformTpcEnvelopeToWorld(test_env); + std::cout << " test_glob " << test_glob.x() << " " << test_glob.y() << " " << test_glob.z() << std::endl; + Acts::Vector3 test_env_check = m_tGeometry-> transformTpcWorldToEnvelope(test_glob); + std::cout << " test_env_check " << test_env_check.x() << " " << test_env_check.y() << " " << test_env_check.z() << std::endl; + */ return Fun4AllReturnCodes::EVENT_OK; } @@ -556,11 +564,22 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // values between 0 and 1 const double f = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0); - const double x_start = hiter->second->get_x(0) + f * (hiter->second->get_x(1) - hiter->second->get_x(0)); - const double y_start = hiter->second->get_y(0) + f * (hiter->second->get_y(1) - hiter->second->get_y(0)); - const double z_start = hiter->second->get_z(0) + f * (hiter->second->get_z(1) - hiter->second->get_z(0)); + const double x_start_glob = hiter->second->get_x(0) + f * (hiter->second->get_x(1) - hiter->second->get_x(0)); + const double y_start_glob = hiter->second->get_y(0) + f * (hiter->second->get_y(1) - hiter->second->get_y(0)); + const double z_start_glob = hiter->second->get_z(0) + f * (hiter->second->get_z(1) - hiter->second->get_z(0)); const double t_start = hiter->second->get_t(0) + f * (hiter->second->get_t(1) - hiter->second->get_t(0)); + Acts::Vector3 start_glob(x_start_glob, y_start_glob, z_start_glob); + Acts::Vector3 start = m_tGeometry->transformTpcWorldToEnvelope(start_glob); // we drift in tpc envelope coords, where E is in the z direction + + const double x_start = start.x(); + const double y_start = start.y(); + const double z_start = start.z(); + /* + std::cout << " xg " << x_start_glob << " x " << x_start + <<" yg " << y_start_glob << " y " << y_start + <<" zg " << z_start_glob << " z " << z_start << std::endl; + */ unsigned int side = 0; if (z_start > 0) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h index de2810e918..f290d7d92f 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h @@ -38,6 +38,7 @@ class DistortedTrackContainer; class TpcClusterBuilder; class PHG4TpcGeomContainer; class ClusHitsVerbose; +class ActsGeometry; class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface { diff --git a/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc b/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc index eb39ed8033..5fb3ee433b 100644 --- a/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc +++ b/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc @@ -330,9 +330,11 @@ void TpcClusterBuilder::cluster_hits(TrkrTruthTrack* track) } // end debug printing // get the global vector3 to then get the surface local phi and z - Acts::Vector3 global(clusx, clusy, clusz); + Acts::Vector3 global_env(clusx, clusy, clusz); TrkrDefs::subsurfkey subsurfkey = 0; + // get_tpc_surface_from coords and the Acts transform both expect coordinates in world (i.e. tilted TPC) coordinates + Acts::Vector3 global = m_tGeometry->transformTpcEnvelopeToWorld(global_env); Surface surface = m_tGeometry->get_tpc_surface_from_coords( hitsetkey, global, subsurfkey); From e08f803676c1f09f9d7eeccd35a72c5ef9841a94 Mon Sep 17 00:00:00 2001 From: Joe Osborn Date: Mon, 22 Jun 2026 11:52:22 -0400 Subject: [PATCH 2/7] fix compilation --- offline/QA/Tracking/TpcSeedsQA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/offline/QA/Tracking/TpcSeedsQA.cc b/offline/QA/Tracking/TpcSeedsQA.cc index dc78d7d51c..3a23895c18 100644 --- a/offline/QA/Tracking/TpcSeedsQA.cc +++ b/offline/QA/Tracking/TpcSeedsQA.cc @@ -73,7 +73,7 @@ int TpcSeedsQA::InitRun(PHCompositeNode *topNode) // global position wrapper m_globalPositionWrapper.loadNodes(topNode); - m_clusterMover.initialize_geometry(g4geom); + m_clusterMover.initialize_geometry(g4geom, actsgeom); m_clusterMover.set_verbosity(0); auto *hm = QAHistManagerDef::getHistoManager(); From a9a5359fd0203aa501a4926b1dac6b585f513017 Mon Sep 17 00:00:00 2001 From: Anthony Denis Frawley Date: Thu, 25 Jun 2026 12:00:41 -0400 Subject: [PATCH 3/7] Implemented coderabbit suggestions. TpcClusterMover initialization, add ActsGeom for envelope transform. TpcClusterizer, use envelope coords where appropriate. AlignmentTransformation, redo sskey finding. Use envelope center coords in module tilt section. Add PHG4TpcGeom, add g4detectors to makefile. PHG4TpcEndCapDetector, change rot_x etc from deg to rad. TpcClusterBuilder, add check for adc_sum = 0. --- offline/packages/tpc/TpcClusterMover.cc | 7 +- offline/packages/tpc/TpcClusterizer.cc | 6 +- offline/packages/trackbase/ActsGeometry.cc | 7 +- .../trackbase/AlignmentTransformation.cc | 148 +++++++++++------- .../trackbase/AlignmentTransformation.h | 4 +- offline/packages/trackbase/Makefile.am | 1 + .../trackbase_historic/TrackAnalysisUtils.cc | 4 +- .../packages/trackreco/MakeActsGeometry.cc | 26 +-- offline/packages/trackreco/PHActsTrkFitter.cc | 1 + offline/packages/trackreco/PHCASeeding.cc | 14 ++ .../g4simulation/g4detectors/PHG4TpcGeom.h | 12 +- .../g4detectors/PHG4TpcGeomContainer.h | 2 +- .../g4simulation/g4detectors/PHG4TpcGeomv2.cc | 58 +++++-- .../g4simulation/g4detectors/PHG4TpcGeomv2.h | 4 +- .../g4simulation/g4eval/SvtxTruthEval.cc | 3 + .../g4simulation/g4tpc/PHG4TpcDetector.cc | 4 +- .../g4tpc/PHG4TpcEndCapDetector.cc | 14 +- .../g4tpc/PHG4TpcEndCapSubsystem.cc | 3 +- .../g4simulation/g4tpc/PHG4TpcSubsystem.cc | 2 +- .../g4simulation/g4tpc/TpcClusterBuilder.cc | 6 + 20 files changed, 219 insertions(+), 107 deletions(-) diff --git a/offline/packages/tpc/TpcClusterMover.cc b/offline/packages/tpc/TpcClusterMover.cc index d1f391b9ea..a7663349d5 100644 --- a/offline/packages/tpc/TpcClusterMover.cc +++ b/offline/packages/tpc/TpcClusterMover.cc @@ -49,9 +49,14 @@ void TpcClusterMover::initialize_geometry(PHG4TpcGeomContainer* cellgeo, ActsGeo { if (_verbosity > 0) { - std::cout << "TpcClusterMover: Initializing layer radii for Tpc from cell geometry object" << std::endl; + std::cout << "TpcClusterMover: Getting ActsGeometry, and getting layer radii for Tpc from cell geometry object" << std::endl; } + if(!tGeometry || !cellgeo) + { + std::cout << PHWHERE << " Failed to get ActsGeometry or TPC cell geometry, cannot continue - quit!" << std::endl; + exit(1); + } _tGeometry = tGeometry; diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index a9ec66a219..ba62fda0f8 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -764,9 +764,11 @@ namespace double nn_z = training_hits->z + std::clamp(ten_pos[0][1][0].item(), -(double) nd, (double) nd) * training_hits->zstep; double nn_x = radius * std::cos(nn_phi); double nn_y = radius * std::sin(nn_phi); - Acts::Vector3 nn_global(nn_x, nn_y, nn_z); + + Acts::Vector3 nn_env_global(nn_x, nn_y, nn_z); + Acts::Vector3 nn_global = my_data.tGeometry->transformTpcEnvelopeToWorld(nn_env_global); nn_global *= Acts::UnitConstants::cm; - Acts::Vector3 nn_local = surface->localToGlobalTransform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global; + Acts::Vector3 nn_local = surface->localToGlobalTransform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global; nn_local /= Acts::UnitConstants::cm; double nn_t = my_data.m_tdriftmax - std::fabs(nn_z) / my_data.tGeometry->get_drift_velocity(); clus_base->setLocalX(nn_local(0)); diff --git a/offline/packages/trackbase/ActsGeometry.cc b/offline/packages/trackbase/ActsGeometry.cc index bac8fdf558..b57a8f1b79 100644 --- a/offline/packages/trackbase/ActsGeometry.cc +++ b/offline/packages/trackbase/ActsGeometry.cc @@ -187,10 +187,11 @@ Surface ActsGeometry::get_tpc_surface_from_coords( double surf_phi = atan2(surf_center_envelope[1], surf_center_envelope[0]); double surfStepPhi = m_tGeometry.tpcSurfStepPhi; - if ((world_phi > surf_phi - surfStepPhi / 2.0) && (world_phi < surf_phi + surfStepPhi / 2.0)) + const double dphi = std::atan2(std::sin(world_phi - surf_phi), std::cos(world_phi - surf_phi)); + if (std::abs(dphi) < surfStepPhi / 2.0) { - if(surf_center.z() < 0 && side != 0) { continue; } - if(surf_center.z() > 0 && side != 1) { continue; } + if(surf_center_envelope.z() < 0 && side != 0) { continue; } + if(surf_center_envelope.z() > 0 && side != 1) { continue; } surf_index = isurf; subsurfkey = isurf; break; diff --git a/offline/packages/trackbase/AlignmentTransformation.cc b/offline/packages/trackbase/AlignmentTransformation.cc index 78078b0fbc..8cfaa2bf88 100644 --- a/offline/packages/trackbase/AlignmentTransformation.cc +++ b/offline/packages/trackbase/AlignmentTransformation.cc @@ -21,6 +21,9 @@ #include #include +#include +#include + #include #include #include @@ -194,7 +197,7 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) Acts::GeometryIdentifier id = surf->geometryId(); - if (localVerbosity) + if (localVerbosity > 1) { std::cout << " Add transform for MVTX with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl; std::cout << " final mvtx transform:" << std::endl @@ -223,7 +226,7 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, localFrameTranslation, sensorAnglesGlobal, trkrId, use_intt_survey_geometry); Acts::GeometryIdentifier id = surf->geometryId(); - if (localVerbosity) + if (localVerbosity > 1) { std::cout << " Add transform for INTT with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl; } @@ -263,58 +266,82 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) // std::cout << " *** module hitsetkey " << hitsetkey << " this_hitsetkey " << this_hitsetkey << " this layer " << this_layer << " side " << side << " sector " << sector << std::endl; - // is this correct?????? - int subsurfkey_min = (1 - side) * 144 + (144 - sector * 12) - 12 - 6; - int subsurfkey_max = subsurfkey_min + 12; - for (int subsurfkey = subsurfkey_min; subsurfkey < subsurfkey_max; subsurfkey++) - { - int sskey = subsurfkey; - if (sskey < 0) - { - sskey += 288; - } - - surf = surfMaps.getTpcSurface(this_hitsetkey, (unsigned int) sskey); - - Eigen::Vector3d localFrameTranslation(0, 0, 0); - use_module_tilt = false; - if (test_layer < 4 || use_module_tilt_always) - { - // get the local frame translation that puts the local surface center at the tilted position after the local rotations are applied - unsigned int this_region = (this_layer - 7) / 16; // 0-2 - Eigen::Vector3d this_center = surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1; // mm to cm - double this_radius = std::sqrt(this_center[0] * this_center[0] + this_center[1] * this_center[1]); - float moduleRadius = TpcModuleRadii[side][sector][this_region]; // radius of the center of the module in cm - localFrameTranslation = getTpcLocalFrameTranslation(moduleRadius, this_radius, sensorAngles) * 10; // cm to mm - - // set this flag for later use - use_module_tilt = true; - } - - Acts::Transform3 transform; - transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, localFrameTranslation, sensorAnglesGlobal, trkrId, false); - Acts::GeometryIdentifier id = surf->geometryId(); - - if (localVerbosity) - { - unsigned int layer = this_layer; - std::cout << " Add transform for TPC with surface GeometryIdentifier " << id << std::endl - << " trkrid " << trkrId << " hitsetkey " << this_hitsetkey << " layer " << layer << " sector " << sector << " side " << side - << " subsurfkey " << subsurfkey << std::endl; - Acts::Vector3 center = surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1; // convert to cm - std::cout << "Ideal surface center: " << std::endl - << center << std::endl; - std::cout << "transform matrix: " << std::endl - << transform.matrix() << std::endl; - } - transformMap->addTransform(id, transform); - transformMapTransient->addTransform(id, transform); - } + // Each TPC hitsetkey has 12 fake surfaces associated with it + // We want to make a transform for every fake surface in this hitsetkey + // Loop over the sector phi angles for the fake surfaces and get each surface + auto layergeom = m_tpccellgeo->GetLayerCellGeom((int) this_layer); + auto sec_min_phi = layergeom->get_sector_min_phi(); + auto min_phi = sec_min_phi[side][sector]; + auto sec_max_phi = layergeom->get_sector_max_phi(); + auto max_phi = sec_max_phi[side][sector]; + double dphi = (max_phi - min_phi)/12.0; + for(int is = 0; is < 12; ++is) + { + double phis = min_phi + is*dphi + dphi/2.0; + double radius = layergeom->get_radius(); + double zcenter = 51.0; + if(side == 0) + { + zcenter *= -1; + } + Acts::Vector3 env_pos(radius*std::cos(phis), radius * std::sin(phis), zcenter); + Acts::Vector3 world_pos = m_tGeometry->transformTpcEnvelopeToWorld(env_pos); + unsigned short sskey = 999; + Surface this_surf = m_tGeometry->get_tpc_surface_from_coords(this_hitsetkey, world_pos, sskey); + if(sskey == 999) + { + std::cout << PHWHERE << "Failed to get surface for layer " << this_layer << " side " << side << " sector " << sector << " quit!" << std::endl; + } + /* + std::cout << " layer " << this_layer << " radius " << radius << " phis " << phis << " min_phi " << min_phi << " max_phi " << max_phi + << " side " << side << " sector " << sector << " world " << world_pos.x() << " " << world_pos.y() << " " << world_pos.z() + <<" world_radius " << sqrt(world_pos.x() * world_pos.x() + world_pos.y() * world_pos.y()) + << " sskey " << sskey << std::endl; + */ + + Eigen::Vector3d localFrameTranslation(0, 0, 0); + use_module_tilt = false; + if (test_layer < 4 || use_module_tilt_always) + { + // get the local frame translation that puts the local surface center at the tilted position after the local rotations are applied + unsigned int this_region = (this_layer - 7) / 16; // 0-2 + Eigen::Vector3d this_center = this_surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1; // mm to cm + //this_center includes the TPC tilt used in PHG4TpcDetector construction, transform to tpc envelope coords + Acts::Vector3 this_center_envelope = m_tGeometry->transformTpcWorldToEnvelope(this_center); + double this_radius = std::sqrt(this_center_envelope[0] * this_center_envelope[0] + this_center_envelope[1] * this_center_envelope[1]); + float moduleRadius = TpcModuleRadii[side][sector][this_region]; // radius of the center of the module in cm + localFrameTranslation = getTpcLocalFrameTranslation(moduleRadius, this_radius, sensorAngles) * 10; // cm to mm + + // set this flag for later use + use_module_tilt = true; + } + + Acts::Transform3 transform; + transform = newMakeTransform(this_surf, millepedeTranslation, sensorAngles, localFrameTranslation, sensorAnglesGlobal, trkrId, false); + Acts::GeometryIdentifier id = this_surf->geometryId(); + + if (localVerbosity) + { + std::cout << " Add transform for TPC with surface GeometryIdentifier " << id + << " trkrid " << trkrId << " hitsetkey " << this_hitsetkey << " layer " << this_layer << " sector " << sector + << " side " << side << std::endl; + if(localVerbosity > 1) + { + Acts::Vector3 center = this_surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1; // convert to cm + std::cout << "Ideal surface center: " << std::endl + << center << std::endl; + std::cout << "transform matrix: " << std::endl + << transform.matrix() << std::endl; + } + } + transformMap->addTransform(id, transform); + transformMapTransient->addTransform(id, transform); + } } - + break; } - + case TrkrDefs::micromegasId: { if (perturbMM) @@ -332,7 +359,7 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) transform = newMakeTransform(surf, millepedeTranslation, sensorAngles, localFrameTranslation, sensorAnglesGlobal, trkrId, false); Acts::GeometryIdentifier id = surf->geometryId(); - if (localVerbosity) + if (localVerbosity > 1) { std::cout << " Add transform for Micromegas with surface GeometryIdentifier " << id << " trkrid " << trkrId << std::endl; } @@ -461,7 +488,7 @@ Acts::Transform3 AlignmentTransformation::newMakeTransform(const Surface& surf, } } - if (localVerbosity) + if (localVerbosity > 1) { Acts::Transform3 actstransform = actsTranslationAffine * actsRotationAffine; @@ -526,7 +553,7 @@ Eigen::Vector3d AlignmentTransformation::getTpcLocalFrameTranslation(float modul dy += -Rdiff * (1 - std::cos(gamma)); dz += 0.0; - if (localVerbosity) + if (localVerbosity > 1) { std::cout << " alpha, beta, gamma " << alpha << " " << beta << " " << gamma << " radius " << moduleRadius << " Rdiff " << Rdiff << " dx, dy dz " << dx << " " << dy << " " << dz << std::endl; @@ -548,6 +575,13 @@ int AlignmentTransformation::getNodes(PHCompositeNode* topNode) return Fun4AllReturnCodes::ABORTEVENT; } + m_tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); + if (!m_tpccellgeo) + { + std::cout << PHWHERE << " unable to find DST node TPCGEOMCONTAINER" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + return 0; } @@ -621,7 +655,7 @@ void AlignmentTransformation::generateRandomPerturbations(Eigen::Vector3d angleD std::normal_distribution distribution(0, transformDev(2)); perturbationTranslation(2) = distribution(generator); } - if (localVerbosity) + if (localVerbosity > 1) { std::cout << "randomperturbationAngles" << perturbationAngles << " randomperturbationTrans:" << perturbationTranslation << std::endl; } @@ -629,7 +663,7 @@ void AlignmentTransformation::generateRandomPerturbations(Eigen::Vector3d angleD void AlignmentTransformation::extractModuleCenterPositions() { - if (localVerbosity) + if (localVerbosity > 1) { std::cout << "Extracting TPC module center radii:" << std::endl; } @@ -654,7 +688,7 @@ void AlignmentTransformation::extractModuleCenterPositions() double mod_radius = (surf_rad_in + surf_rad_out) / 2.0; TpcModuleRadii[iside][isector][iregion] = mod_radius; - if (localVerbosity) + if (localVerbosity > 1) { std::cout << " hitsetkey_in " << hitsetkey_in << " lin " << lin << " sector " << isector << " side " << iside << " region " << iregion << std::endl; std::cout << " hitsetkey_out " << hitsetkey_out << " lout " << lout << " sector " << isector << " side " << iside << " region " << iregion << std::endl; diff --git a/offline/packages/trackbase/AlignmentTransformation.h b/offline/packages/trackbase/AlignmentTransformation.h index e8eb3882d1..8df8aeb2ef 100644 --- a/offline/packages/trackbase/AlignmentTransformation.h +++ b/offline/packages/trackbase/AlignmentTransformation.h @@ -11,7 +11,7 @@ #include class PHCompositeNode; - +class PHG4TpcGeomContainer; class ActsGeometry; class AlignmentTransformation @@ -148,6 +148,8 @@ class AlignmentTransformation float TpcModuleRadii[2][12][3] = {}; // module radial center in local coords unsigned int innerLayer[3] = {}; double sectorPhi[2][12] = {}; + + PHG4TpcGeomContainer *m_tpccellgeo = nullptr; }; #endif diff --git a/offline/packages/trackbase/Makefile.am b/offline/packages/trackbase/Makefile.am index 439d73256e..52c4c3e505 100644 --- a/offline/packages/trackbase/Makefile.am +++ b/offline/packages/trackbase/Makefile.am @@ -295,6 +295,7 @@ libtrack_la_LIBADD = \ -lActsPluginRoot \ -lActsExamplesDetectorTGeo \ -lffamodules \ + -lg4detectors \ -lboost_program_options libtrack_io_la_LIBADD = \ diff --git a/offline/packages/trackbase_historic/TrackAnalysisUtils.cc b/offline/packages/trackbase_historic/TrackAnalysisUtils.cc index 49bbd418bb..4d9282e1f1 100644 --- a/offline/packages/trackbase_historic/TrackAnalysisUtils.cc +++ b/offline/packages/trackbase_historic/TrackAnalysisUtils.cc @@ -309,10 +309,10 @@ namespace TrackAnalysisUtils globalWrapper.loadNodes(topNode); globalWrapper.set_suppressCrossing(true); - auto* geometry = findNode::getClass(topNode, "ActsGeometry"); - TpcClusterMover mover; + auto* geometry = findNode::getClass(topNode, "ActsGeometry"); auto* tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); + TpcClusterMover mover; mover.initialize_geometry(tpccellgeo, geometry); mover.set_verbosity(0); diff --git a/offline/packages/trackreco/MakeActsGeometry.cc b/offline/packages/trackreco/MakeActsGeometry.cc index 88319fa938..1c51d07283 100644 --- a/offline/packages/trackreco/MakeActsGeometry.cc +++ b/offline/packages/trackreco/MakeActsGeometry.cc @@ -183,8 +183,8 @@ int MakeActsGeometry::InitRun(PHCompositeNode *topNode) m_CM_halfwidth = layergeom->get_CM_halfwidth(); m_maxSurfZ = m_max_driftlength - 0.0001; // add clearance from physical TPC gas volume length to avoid overlaps - // This transform will eventually be built using the tilt and placement variables that will be in layergeom - // TPC envelope to global transformation + // Make the transform from TPC envelope to global coordinates + // This transform is built using the tilt and placement variables from layergeom double rot_x = layergeom->get_rot_x(); double rot_y = layergeom->get_rot_y(); @@ -513,7 +513,7 @@ void MakeActsGeometry::editTPCGeometry(PHCompositeNode *topNode) return; } - if (Verbosity() > 3) + if (Verbosity() > 0) { std::cout << "EditTPCGeometry - gas volume: "; tpc_gas_north_vol->Print(); @@ -559,7 +559,7 @@ void MakeActsGeometry::addActsTpcSurfaces(TGeoVolume *tpc_gas_vol, tpc_gas_measurement_vol[ilayer]->SetFillColor(kYellow); tpc_gas_measurement_vol[ilayer]->SetVisibility(kTRUE); - if (Verbosity() > 3) + if (Verbosity() > 0) { std::cout << " Made box for layer " << ilayer << " with dx " << m_layerThickness[ilayer] << " dy " @@ -899,7 +899,8 @@ void MakeActsGeometry::makeTpcMapPairs(TrackingVolumePtr &tpcVolume) TrkrDefs::hitsetkey hitsetkey = getTpcHitSetKeyFromCoords(world_center); unsigned int layer = TrkrDefs::getLayer(hitsetkey); - + // unsigned int sector = TpcDefs::getSectorId(hitsetkey); + // unsigned int side = TpcDefs::getSide(hitsetkey); // If there is already an entry for this hitsetkey, add the surface // to its corresponding vector // std::map>::iterator mapIter; @@ -909,11 +910,13 @@ void MakeActsGeometry::makeTpcMapPairs(TrackingVolumePtr &tpcVolume) if (mapIter != m_clusterSurfaceMapTpcEdit.end()) { + //std::cout << " Adding surface to map with layer " << layer << " side " << side << " sector " << sector << std::endl; mapIter->second.push_back(surf); } else { // Otherwise make a new map entry + // std::cout << "Starting new surfvec for layer " << layer << " side " << side << " sector " << sector << std::endl; std::vector dumvec; dumvec.push_back(surf); std::pair> tmp = @@ -927,7 +930,7 @@ void MakeActsGeometry::makeTpcMapPairs(TrackingVolumePtr &tpcVolume) //____________________________________________________________________________________________ void MakeActsGeometry::makeMmMapPairs(TrackingVolumePtr &mmVolume) { - if (Verbosity()) + if (Verbosity()>1) { std::cout << "MakeActsGeometry::makeMmMapPairs - mmVolume: " << mmVolume->volumeName() << std::endl; } @@ -989,7 +992,7 @@ void MakeActsGeometry::makeMmMapPairs(TrackingVolumePtr &mmVolume) continue; } - if (Verbosity()) + if (Verbosity()>1) { std::cout << "MakeActsGeometry::makeMmMapPairs - layer: " << layer << " tileid: " << tileid << std::endl; } @@ -1155,7 +1158,7 @@ void MakeActsGeometry::makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume) auto vec3d = surf->center(m_geoCtxt); std::vector world_center = {(vec3d(0) - v_globaldisplacement[0]) / 10.0, (vec3d(1) - v_globaldisplacement[1]) / 10.0, (vec3d(2) - v_globaldisplacement[2]) / 10.0}; // convert from mm to cm double layer_rad = sqrt(pow(world_center[0], 2) + pow(world_center[1], 2)); - if (Verbosity() > 0) + if (Verbosity() > 1) { std::cout << "[DEBUG] MVTX surface center (before misalignment): (x,y,z)=(" << vec3d(0) / 10. << "," << vec3d(1) / 10. << "," << vec3d(2) / 10. << "), layer_rad=" << sqrt(pow(vec3d(0) / 10., 2) + pow(vec3d(1) / 10., 2)) << std::endl; std::cout << "[DEBUG] MVTX surface center: (x,y,z)=(" << world_center[0] << "," << world_center[1] << "," << world_center[2] << "), layer_rad=" << layer_rad << std::endl; @@ -1255,7 +1258,8 @@ void MakeActsGeometry::makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume) TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector &world) { - // the input position is assumed to be in tpc envelope coords - i.e. tilt removed + // This is used only in simulations + // so the input position is assumed to be in tpc envelope coords - i.e. tilt removed // Look up TPC surface index values from tpc envelope position of surface center // layer unsigned int layer = 999; @@ -1310,8 +1314,8 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector 3 && layer == 7) - { + if (Verbosity() > 3 && layer == 15) + { std::cout << " layer_rad " << layer_rad << " m_layerRadius[layer] " << m_layerRadius[layer-7] << " found layer " << layer << " side " << side << " world " << world[0] << " " << world[1] << " " << world[2] << " phi_world " << phi_world << " readout_mod " << readout_mod << std::endl; } diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index f038e736e1..2fc97257bb 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -1005,6 +1005,7 @@ SurfacePtrVec PHActsTrkFitter::getSurfaceVector(const SourceLinkVec& sourceLinks { const ActsSourceLink asl = sl.get(); const auto* const surf = m_tGeometry->geometry().tGeometry->findSurface(asl.geometryId()); + // std::cout << "sl: " << surf->geometryId() << std::endl; surfaces.push_back(surf); } diff --git a/offline/packages/trackreco/PHCASeeding.cc b/offline/packages/trackreco/PHCASeeding.cc index abfcc6018a..a589cad71c 100644 --- a/offline/packages/trackreco/PHCASeeding.cc +++ b/offline/packages/trackreco/PHCASeeding.cc @@ -31,6 +31,7 @@ #include #include #include // for getLayer, clu... +#include #include // ROOT includes for debugging @@ -57,6 +58,8 @@ #include // for pair, make_pair #include +#include // for uint8_t, uint16_t, uint32_t + //#define _DEBUG_ #if defined(_DEBUG_) @@ -211,6 +214,17 @@ int PHCASeeding::InitializeGeometry(PHCompositeNode* topNode) Acts::Vector3 PHCASeeding::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* cluster) const { + /* + unsigned int layer = TrkrDefs::getLayer(key); + unsigned int side = TpcDefs::getSide(key); + unsigned int sector = TpcDefs::getSectorId(key); + std::cout << " _pp_mode = " << _pp_mode + << " layer " << layer + << " side " << side + << " sector " << sector + << " subsurfkey " << cluster->getSubSurfKey() + << std::endl; + */ return _pp_mode ? m_tGeometry->getGlobalPosition(key, cluster) : m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, 0); } diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeom.h b/simulation/g4simulation/g4detectors/PHG4TpcGeom.h index 521cf2fa9c..93376c1b72 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeom.h +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeom.h @@ -195,32 +195,32 @@ class PHG4TpcGeom : public PHObject virtual double get_rot_x() const { PHOOL_VIRTUAL_WARN("get_rot_x()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } virtual double get_rot_y() const { PHOOL_VIRTUAL_WARN("get_rot_y()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } virtual double get_rot_z() const { PHOOL_VIRTUAL_WARN("get_rot_z()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } virtual double get_place_x() const { PHOOL_VIRTUAL_WARN("get_place_x()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } virtual double get_place_y() const { PHOOL_VIRTUAL_WARN("get_place_y()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } virtual double get_place_z() const { PHOOL_VIRTUAL_WARN("get_place_z()"); - return std::numeric_limits::quiet_NaN(); + return 0.0; } diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h b/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h index a92b24657a..9bca152e81 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h @@ -9,7 +9,7 @@ #include #include // for make_pair, pair -class PHG4TpcGeom; +#include class PHG4TpcGeomContainer : public PHObject { diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc index 050baf86a9..e500c53117 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc @@ -193,7 +193,7 @@ void PHG4TpcGeomv2::identify(std::ostream& os) const std::pair PHG4TpcGeomv2::get_zbounds(const int ibin) const { - if (ibin < 0 || ibin > nzbins) + if (ibin < 0 || ibin >= nzbins) { std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl; exit(1); @@ -207,7 +207,7 @@ PHG4TpcGeomv2::get_zbounds(const int ibin) const std::pair PHG4TpcGeomv2::get_etabounds(const int ibin) const { - if (ibin < 0 || ibin > nzbins) + if (ibin < 0 || ibin >= nzbins) { std::cout << PHWHERE << " Asking for invalid bin in z: " << ibin << std::endl; exit(1); @@ -222,7 +222,7 @@ PHG4TpcGeomv2::get_etabounds(const int ibin) const std::pair PHG4TpcGeomv2::get_phibounds(const int ibin) const { - if (ibin < 0 || ibin > nphibins) + if (ibin < 0 || ibin >= nphibins) { std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; exit(1); @@ -235,7 +235,7 @@ PHG4TpcGeomv2::get_phibounds(const int ibin) const int PHG4TpcGeomv2::get_zbin(const double z) const { - if (z < zmin || z > (zmin + nzbins * zstep)) + if (z < zmin || z >= (zmin + nzbins * zstep)) { // cout << PHWHERE << "Asking for bin for z outside of z range: " << z << endl; return -1; @@ -247,7 +247,7 @@ int PHG4TpcGeomv2::get_zbin(const double z) const int PHG4TpcGeomv2::get_etabin(const double eta) const { - if (eta < zmin || eta > (zmin + nzbins * zstep)) + if (eta < zmin || eta >= (zmin + nzbins * zstep)) { // cout << "Asking for bin for eta outside of eta range: " << eta << endl; return -1; @@ -259,7 +259,7 @@ int PHG4TpcGeomv2::get_etabin(const double eta) const int PHG4TpcGeomv2::get_phibin_new(const double phi) const { double norm_phi = phi; - if (phi < phimin || phi > (phimin + nphibins * phistep)) + if (phi < phimin || phi >= (phimin + nphibins * phistep)) { int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); norm_phi += 2 * M_PI * nwraparound; @@ -270,8 +270,14 @@ int PHG4TpcGeomv2::get_phibin_new(const double phi) const int PHG4TpcGeomv2::find_phibin(const double phi, int side) const { + if(side < 0 || side > 1) + { + std::cout << PHWHERE << " side is not valid, have to quit!" << std::endl; + exit(1); + } + double norm_phi = phi; - if (phi < phimin || phi > (phimin + nphibins * phistep)) + if (phi < phimin || phi >= (phimin + nphibins * phistep)) { int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); norm_phi += 2 * M_PI * nwraparound; @@ -315,8 +321,14 @@ int PHG4TpcGeomv2::find_phibin(const double phi, int side) const float PHG4TpcGeomv2::get_pad_float(const double phi, int side) const { + if(side < 0 || side > 1) + { + std::cout << PHWHERE << " side is not valid, have to quit!" << std::endl; + exit(1); + } + double norm_phi = phi; - if (phi < phimin || phi > (phimin + nphibins * phistep)) + if (phi < phimin || phi >= (phimin + nphibins * phistep)) { int nwraparound = -floor((phi - phimin) * 0.5 / M_PI); norm_phi += 2 * M_PI * nwraparound; @@ -360,7 +372,7 @@ float PHG4TpcGeomv2::get_pad_float(const double phi, int side) const float PHG4TpcGeomv2::get_tbin_float(const double z) const { - if (z < zmin || z > (zmin + nzbins * zstep)) + if (z < zmin || z >= (zmin + nzbins * zstep)) { // cout << PHWHERE << "Asking for bin for z outside of z range: " << z << endl; return -1; @@ -372,6 +384,12 @@ float PHG4TpcGeomv2::get_tbin_float(const double z) const int PHG4TpcGeomv2::get_phibin(const double phi, int side) const { + if(side < 0 || side > 1) + { + std::cout << PHWHERE << " side is not valid, have to quit!" << std::endl; + exit(1); + } + double new_phi = phi; if (phi > M_PI) { @@ -431,7 +449,7 @@ int PHG4TpcGeomv2::get_phibin(const double phi, int side) const double PHG4TpcGeomv2::get_zcenter(const int ibin) const { - if (ibin < 0 || ibin > nzbins) + if (ibin < 0 || ibin >= nzbins) { std::cout << PHWHERE << "Asking for invalid bin in z: " << ibin << std::endl; exit(1); @@ -443,7 +461,7 @@ PHG4TpcGeomv2::get_zcenter(const int ibin) const double PHG4TpcGeomv2::get_etacenter(const int ibin) const { - if (ibin < 0 || ibin > nzbins) + if (ibin < 0 || ibin >= nzbins) { std::cout << PHWHERE << "Asking for invalid bin in eta: " << ibin << std::endl; std::cout << "minbin: 0, maxbin " << nzbins << std::endl; @@ -456,7 +474,7 @@ PHG4TpcGeomv2::get_etacenter(const int ibin) const double PHG4TpcGeomv2::get_phicenter_new(const int ibin) const { - if (ibin < 0 || ibin > nphibins) + if (ibin < 0 || ibin >= nphibins) { std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; exit(1); @@ -470,8 +488,14 @@ PHG4TpcGeomv2::get_phicenter_new(const int ibin) const double PHG4TpcGeomv2::get_phicenter(const int ibin, const int side) const { + if(side < 0 || side > 1) + { + std::cout << PHWHERE << " side is not valid, have to quit!" << std::endl; + exit(1); + } + // double phi_center = -999; - if (ibin < 0 || ibin > nphibins) + if (ibin < 0 || ibin >= nphibins) { std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; exit(1); @@ -493,8 +517,14 @@ PHG4TpcGeomv2::get_phicenter(const int ibin, const int side) const double PHG4TpcGeomv2::get_phi(const float ibin, const int side) const { + if(side < 0 || side > 1) + { + std::cout << PHWHERE << " side is not valid, have to quit!" << std::endl; + exit(1); + } + // double phi_center = -999; - if (ibin < 0 || ibin > nphibins) + if (ibin < 0 || ibin >= nphibins) { std::cout << PHWHERE << "Asking for invalid bin in phi: " << ibin << std::endl; exit(1); diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h index 235256b67c..773c6451b2 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h @@ -1,7 +1,7 @@ // Tell emacs that this is a C++ source // -*- C++ -*-. -#ifndef G4DETECTORS_PHG4TPCGEOMV1_H -#define G4DETECTORS_PHG4TPCGEOMV1_H +#ifndef G4DETECTORS_PHG4TPCGEOMV2_H +#define G4DETECTORS_PHG4TPCGEOMV2_H #include "PHG4TpcGeom.h" diff --git a/simulation/g4simulation/g4eval/SvtxTruthEval.cc b/simulation/g4simulation/g4eval/SvtxTruthEval.cc index f2fcd220a7..25d4cbb125 100644 --- a/simulation/g4simulation/g4eval/SvtxTruthEval.cc +++ b/simulation/g4simulation/g4eval/SvtxTruthEval.cc @@ -324,6 +324,9 @@ std::map> SvtxTruthEval::all_tru std::vector contributing_hits_energy; std::vector> contributing_hits_entry; std::vector> contributing_hits_exit; + // contributing_hits are the original g4hits in world coords + // contributing_hits_entry, contributing_hits_exit are in envelope coords, for use in G4ClusterSize() + // gx, gy, gz are the cluster position in this layer in world coords to compare with data LayerClusterG4Hits(g4hits, contributing_hits, contributing_hits_energy, contributing_hits_entry, contributing_hits_exit, layer, gx, gy, gz, gt, gedep); if (!(gedep > 0)) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc index 6f736ccc39..e75a34961c 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc @@ -131,12 +131,12 @@ void PHG4TpcDetector::ConstructMe(G4LogicalVolume *logicWorld) logicWorld, false, false, OverlapCheck()); - /* + G4ThreeVector test_env(10.0, 40.0, 80.0); std::cout << " test_env " << test_env.x() << " " << test_env.y() << " " << test_env.z() << std::endl; G4ThreeVector test_glob = test_env.transform(rot); std::cout << " test_glob " << test_glob.x() << " " << test_glob.y() << " " << test_glob.z() << std::endl; - */ + // geometry node add_geometry_node(); diff --git a/simulation/g4simulation/g4tpc/PHG4TpcEndCapDetector.cc b/simulation/g4simulation/g4tpc/PHG4TpcEndCapDetector.cc index 583934798d..3029233299 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcEndCapDetector.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcEndCapDetector.cc @@ -94,9 +94,9 @@ void PHG4TpcEndCapDetector::ConstructMe(G4LogicalVolume *logicWorld) m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm); G4RotationMatrix rotm_center; - rotm_center.rotateX(m_Params->get_double_param("rot_x") * deg); - rotm_center.rotateY(m_Params->get_double_param("rot_y") * deg); - rotm_center.rotateZ(m_Params->get_double_param("rot_z") * deg); + rotm_center.rotateX(m_Params->get_double_param("rot_x") * rad); + rotm_center.rotateY(m_Params->get_double_param("rot_y") * rad); + rotm_center.rotateZ(m_Params->get_double_param("rot_z") * rad); G4Transform3D transform_center(rotm_center, g4vec_center); int i = 0; @@ -107,6 +107,14 @@ void PHG4TpcEndCapDetector::ConstructMe(G4LogicalVolume *logicWorld) G4Transform3D transform_side2 = transform_center * rotm_otherside * g4vec_front_z; m_EndCapAssembly->MakeImprint(logicWorld, transform_side2, i++, OverlapCheck()); + + G4ThreeVector test_env(10.0, 40.0, 80.0); + std::cout << "Endcap: rot_x " << m_Params->get_double_param("rot_x")*rad << " test_env " << test_env.x() << " " << test_env.y() << " " << test_env.z() << std::endl; + G4ThreeVector test_glob1 = test_env.transform(rotm_center); + std::cout << " test_glob1 " << test_glob1.x() << " " << test_glob1.y() << " " << test_glob1.z() << std::endl; + // G4ThreeVector test_glob2 = test_env(transform_side2); + // std::cout << " test_glob2 " << test_glob2.x() << " " << test_glob2.y() << " " << test_glob2.z() << std::endl; + return; } diff --git a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc index ec00172916..e7aae73319 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc @@ -122,7 +122,8 @@ void PHG4TpcEndCapSubsystem::SetDefaultParameters() set_default_double_param("place_x", 0.); set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); - set_default_double_param("rot_x", 0.); + // set_default_double_param("rot_x", 0.); + set_default_double_param("rot_x", -0.004); // TEST! //rad set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); diff --git a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc index 8d802c7886..e1770da257 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc @@ -144,7 +144,7 @@ void PHG4TpcSubsystem::SetDefaultParameters() set_default_double_param("place_x", 0.); set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); - set_default_double_param("rot_x", 0.); + set_default_double_param("rot_x", -0.004); // TEST! // should default to zero set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); set_default_double_param("tpc_length", 205.21); // 2 * (maxdrift 102.325 + CM halfwidth 0.28) cm diff --git a/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc b/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc index 5fb3ee433b..e62f31dce7 100644 --- a/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc +++ b/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc @@ -183,6 +183,12 @@ void TpcClusterBuilder::cluster_hits(TrkrTruthTrack* track) adc_sum += adc; } + + if(adc_sum == 0) + { + continue; + } + if (mClusHitsVerbose) { if (verbosity > 10) From d97ac63e576ec6cbfe6ab124bfe08adb54fb7b56 Mon Sep 17 00:00:00 2001 From: Anthony Denis Frawley Date: Thu, 25 Jun 2026 14:39:54 -0400 Subject: [PATCH 4/7] Fix screwd up file. Implement coderabbit suggestions. Leave -4 mrad tilt hardcoded for this test only (so no macro changes are needed); --- offline/packages/trackbase/AlignmentTransformation.cc | 5 +++-- simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h | 2 +- simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc | 5 +++-- simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc | 3 ++- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/offline/packages/trackbase/AlignmentTransformation.cc b/offline/packages/trackbase/AlignmentTransformation.cc index 8cfaa2bf88..8f25dab812 100644 --- a/offline/packages/trackbase/AlignmentTransformation.cc +++ b/offline/packages/trackbase/AlignmentTransformation.cc @@ -288,9 +288,10 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) Acts::Vector3 world_pos = m_tGeometry->transformTpcEnvelopeToWorld(env_pos); unsigned short sskey = 999; Surface this_surf = m_tGeometry->get_tpc_surface_from_coords(this_hitsetkey, world_pos, sskey); - if(sskey == 999) + if(sskey == 999 || !this_surf) { std::cout << PHWHERE << "Failed to get surface for layer " << this_layer << " side " << side << " sector " << sector << " quit!" << std::endl; + exit(1); } /* std::cout << " layer " << this_layer << " radius " << radius << " phis " << phis << " min_phi " << min_phi << " max_phi " << max_phi @@ -579,7 +580,7 @@ int AlignmentTransformation::getNodes(PHCompositeNode* topNode) if (!m_tpccellgeo) { std::cout << PHWHERE << " unable to find DST node TPCGEOMCONTAINER" << std::endl; - return Fun4AllReturnCodes::ABORTRUN; + exit(1); } return 0; diff --git a/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h b/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h index 9bca152e81..a92b24657a 100644 --- a/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomContainer.h @@ -9,7 +9,7 @@ #include #include // for make_pair, pair -#include +class PHG4TpcGeom; class PHG4TpcGeomContainer : public PHObject { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc index e7aae73319..b340c82202 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc @@ -116,14 +116,15 @@ void PHG4TpcEndCapSubsystem::SetDefaultParameters() { set_default_int_param("construction_verbosity", 0); // sizes are in cm - // angles are in deg + // angles are in rad // units should be converted to G4 units when used // implement your own here// set_default_double_param("place_x", 0.); set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); // set_default_double_param("rot_x", 0.); - set_default_double_param("rot_x", -0.004); // TEST! //rad + // angles are in rad + set_default_double_param("rot_x", -0.004); // TEMPORARY TEST! return to 0.0 set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); diff --git a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc index e1770da257..e9f491b951 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc @@ -144,7 +144,8 @@ void PHG4TpcSubsystem::SetDefaultParameters() set_default_double_param("place_x", 0.); set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); - set_default_double_param("rot_x", -0.004); // TEST! // should default to zero + // angles are in radians + set_default_double_param("rot_x", -0.004); // TEMPORARY TEST! return to 0.0 set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); set_default_double_param("tpc_length", 205.21); // 2 * (maxdrift 102.325 + CM halfwidth 0.28) cm From 45967b3f9334134cbc4e0458b22f6da5ead1e1e8 Mon Sep 17 00:00:00 2001 From: Anthony Denis Frawley Date: Thu, 25 Jun 2026 15:48:26 -0400 Subject: [PATCH 5/7] Modify check of surface vectors to use radii instead of layer from volume id. --- offline/packages/trackreco/PHActsTrkFitter.cc | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index 2fc97257bb..f355b17aa3 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -1014,25 +1014,32 @@ SurfacePtrVec PHActsTrkFitter::getSurfaceVector(const SourceLinkVec& sourceLinks void PHActsTrkFitter::checkSurfaceVec(SurfacePtrVec& surfaces) const { + // Do not assume volume id has the correct layer, check the surface radius + for (unsigned int i = 0; i < surfaces.size() - 1; i++) { const auto& surface = surfaces.at(i); const auto thisVolume = surface->geometryId().volume(); - const auto thisLayer = surface->geometryId().layer(); + const Acts::Vector3 this_center = surface->center(m_tGeometry->geometry().getGeoContext()); + double thisRadius = sqrt(this_center.x()*this_center.x()+this_center.y()*this_center.y()); + const auto nextSurface = surfaces.at(i + 1); const auto nextVolume = nextSurface->geometryId().volume(); - const auto nextLayer = nextSurface->geometryId().layer(); + const Acts::Vector3 next_center = surface->center(m_tGeometry->geometry().getGeoContext()); + double nextRadius = sqrt(next_center.x()*next_center.x()+next_center.y()*next_center.y()); + /// Implement a check to ensure surfaces are sorted if (nextVolume == thisVolume) { - if (nextLayer < thisLayer) + // if (nextLayer < thisLayer) + if (nextRadius < thisRadius) { std::cout << "PHActsTrkFitter::checkSurfaceVec - " << "Surface not in order... removing surface" - << surface->geometryId() << std::endl; + << surface->geometryId() << " with radius " << thisRadius << std::endl; surfaces.erase(surfaces.begin() + i); From 275ed89cf99c0a8e7e72e6f67c78f2c4ece864ac Mon Sep 17 00:00:00 2001 From: Anthony Denis Frawley Date: Thu, 25 Jun 2026 23:46:17 -0400 Subject: [PATCH 6/7] Set default TPC tilt angles all to zero after testing. --- simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc | 3 +-- simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc index b340c82202..4cd13647de 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc @@ -122,9 +122,8 @@ void PHG4TpcEndCapSubsystem::SetDefaultParameters() set_default_double_param("place_x", 0.); set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); - // set_default_double_param("rot_x", 0.); // angles are in rad - set_default_double_param("rot_x", -0.004); // TEMPORARY TEST! return to 0.0 + set_default_double_param("rot_x", 0.); set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); diff --git a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc index e9f491b951..80a04450f4 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc @@ -145,7 +145,7 @@ void PHG4TpcSubsystem::SetDefaultParameters() set_default_double_param("place_y", 0.); set_default_double_param("place_z", 0.); // angles are in radians - set_default_double_param("rot_x", -0.004); // TEMPORARY TEST! return to 0.0 + set_default_double_param("rot_x", 0.); set_default_double_param("rot_y", 0.); set_default_double_param("rot_z", 0.); set_default_double_param("tpc_length", 205.21); // 2 * (maxdrift 102.325 + CM halfwidth 0.28) cm From 4b72c343e8c72fdc4083ee2ce48aafa04e82db3d Mon Sep 17 00:00:00 2001 From: Joe Osborn Date: Fri, 26 Jun 2026 14:33:05 -0400 Subject: [PATCH 7/7] clang-tidy fixes --- offline/packages/trackbase/ActsGeometry.cc | 4 ++-- offline/packages/trackbase/ActsGeometry.h | 4 ++-- offline/packages/trackbase/AlignmentTransformation.cc | 2 +- offline/packages/trackreco/PHActsTrkFitter.cc | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/offline/packages/trackbase/ActsGeometry.cc b/offline/packages/trackbase/ActsGeometry.cc index b57a8f1b79..2e0f3f75c7 100644 --- a/offline/packages/trackbase/ActsGeometry.cc +++ b/offline/packages/trackbase/ActsGeometry.cc @@ -345,14 +345,14 @@ Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* c return local; } - Acts::Vector3 ActsGeometry::transformTpcWorldToEnvelope(Acts::Vector3 world) const + Acts::Vector3 ActsGeometry::transformTpcWorldToEnvelope(const Acts::Vector3& world) const { Acts::Vector3 envelope = m_tpc_world_envelope_transform * world; return envelope; } - Acts::Vector3 ActsGeometry::transformTpcEnvelopeToWorld(Acts::Vector3 envelope) const + Acts::Vector3 ActsGeometry::transformTpcEnvelopeToWorld(const Acts::Vector3& envelope) const { Acts::Vector3 world = m_tpc_world_envelope_transform.inverse() * envelope; diff --git a/offline/packages/trackbase/ActsGeometry.h b/offline/packages/trackbase/ActsGeometry.h index d92bf4af38..be9b4e01cd 100644 --- a/offline/packages/trackbase/ActsGeometry.h +++ b/offline/packages/trackbase/ActsGeometry.h @@ -78,8 +78,8 @@ class ActsGeometry Acts::Transform3 makeAffineTransform(Acts::Vector3 rotation, Acts::Vector3 translation) const; - Acts::Vector3 transformTpcWorldToEnvelope(Acts::Vector3 vin) const ; - Acts::Vector3 transformTpcEnvelopeToWorld(Acts::Vector3 vin) const ; + Acts::Vector3 transformTpcWorldToEnvelope(const Acts::Vector3& world) const ; + Acts::Vector3 transformTpcEnvelopeToWorld(const Acts::Vector3& envelope) const ; Acts::Vector2 getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster) const; Acts::Vector2 getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* cluster, short int crossing) const; diff --git a/offline/packages/trackbase/AlignmentTransformation.cc b/offline/packages/trackbase/AlignmentTransformation.cc index 8f25dab812..0e9556866e 100644 --- a/offline/packages/trackbase/AlignmentTransformation.cc +++ b/offline/packages/trackbase/AlignmentTransformation.cc @@ -269,7 +269,7 @@ void AlignmentTransformation::createMap(PHCompositeNode* topNode) // Each TPC hitsetkey has 12 fake surfaces associated with it // We want to make a transform for every fake surface in this hitsetkey // Loop over the sector phi angles for the fake surfaces and get each surface - auto layergeom = m_tpccellgeo->GetLayerCellGeom((int) this_layer); + auto* layergeom = m_tpccellgeo->GetLayerCellGeom((int) this_layer); auto sec_min_phi = layergeom->get_sector_min_phi(); auto min_phi = sec_min_phi[side][sector]; auto sec_max_phi = layergeom->get_sector_max_phi(); diff --git a/offline/packages/trackreco/PHActsTrkFitter.cc b/offline/packages/trackreco/PHActsTrkFitter.cc index f355b17aa3..4738725af7 100644 --- a/offline/packages/trackreco/PHActsTrkFitter.cc +++ b/offline/packages/trackreco/PHActsTrkFitter.cc @@ -1024,7 +1024,7 @@ void PHActsTrkFitter::checkSurfaceVec(SurfacePtrVec& surfaces) const const Acts::Vector3 this_center = surface->center(m_tGeometry->geometry().getGeoContext()); double thisRadius = sqrt(this_center.x()*this_center.x()+this_center.y()*this_center.y()); - const auto nextSurface = surfaces.at(i + 1); + const auto* nextSurface = surfaces.at(i + 1); const auto nextVolume = nextSurface->geometryId().volume(); const Acts::Vector3 next_center = surface->center(m_tGeometry->geometry().getGeoContext());