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(); diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.cc b/offline/packages/TrackingDiagnostics/TrackResiduals.cc index 4b613c5dd5..bd12b6b877 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.cc +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.cc @@ -118,7 +118,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..a7663349d5 100644 --- a/offline/packages/tpc/TpcClusterMover.cc +++ b/offline/packages/tpc/TpcClusterMover.cc @@ -45,13 +45,21 @@ 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; + 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; + int layer = 0; PHG4TpcGeomContainer::ConstRange layerrange = cellgeo->get_begin_end(); for (PHG4TpcGeomContainer::ConstIterator layeriter = layerrange.first; @@ -67,8 +75,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 +89,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 +150,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 32479b6e5f..0752b557d0 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -30,7 +30,8 @@ #include #include // for SubsysReco -#include +//#include +#include #include #include @@ -806,9 +807,9 @@ namespace // This is local position double clusiphi = iphi_sum / adc_sum; double clusit = it_sum / adc_sum; - - // This is the global position - double clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side); + + // this is the phi position in the TPC envelope + double env_clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side); double clust = t_sum / adc_sum; // ADC of centroid bin @@ -835,7 +836,7 @@ namespace if (my_data.layergeom->get_phistep() > 0) { - padphase = (clusphi - maxphi) / my_data.layergeom->get_phistep(); + padphase = (env_clusphi - maxphi) / my_data.layergeom->get_phistep(); } if (my_data.layergeom->get_zstep() > 0) @@ -843,26 +844,25 @@ namespace tbinphase = (clust - maxt) / my_data.layergeom->get_zstep(); } - 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); // 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, @@ -875,6 +875,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 @@ -961,6 +963,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 @@ -978,9 +981,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)); @@ -1524,16 +1529,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..2e0f3f75c7 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,43 @@ 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; + + 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_envelope.z() < 0 && side != 0) { continue; } + if(surf_center_envelope.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 +215,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 +259,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 +276,9 @@ Surface ActsGeometry::get_tpc_surface_from_coords( } return nullptr; } + */ + - return surf_vec[surf_index]; } //________________________________________________________________________________________________ @@ -291,3 +344,17 @@ Acts::Vector2 ActsGeometry::getLocalCoords(TrkrDefs::cluskey key, TrkrCluster* c return local; } + + Acts::Vector3 ActsGeometry::transformTpcWorldToEnvelope(const Acts::Vector3& world) const + { + Acts::Vector3 envelope = m_tpc_world_envelope_transform * world; + + return envelope; + } + + Acts::Vector3 ActsGeometry::transformTpcEnvelopeToWorld(const 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..be9b4e01cd 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(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; 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..0e9556866e 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; } @@ -254,67 +257,92 @@ 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; - - // 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); - } + // std::cout << " *** module hitsetkey " << hitsetkey << " this_hitsetkey " << this_hitsetkey << " this layer " << this_layer << " side " << side << " sector " << sector << std::endl; + + // 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 || !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 + << " 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 +360,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 +489,7 @@ Acts::Transform3 AlignmentTransformation::newMakeTransform(const Surface& surf, } } - if (localVerbosity) + if (localVerbosity > 1) { Acts::Transform3 actstransform = actsTranslationAffine * actsRotationAffine; @@ -526,7 +554,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 +576,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; + exit(1); + } + return 0; } @@ -621,7 +656,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 +664,7 @@ void AlignmentTransformation::generateRandomPerturbations(Eigen::Vector3d angleD void AlignmentTransformation::extractModuleCenterPositions() { - if (localVerbosity) + if (localVerbosity > 1) { std::cout << "Extracting TPC module center radii:" << std::endl; } @@ -643,26 +678,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 > 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; + std::cout << " module radius " << mod_radius << std::endl; + } + } } } @@ -670,7 +703,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 +713,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 +727,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/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 f7dbbb2b80..4d9282e1f1 100644 --- a/offline/packages/trackbase_historic/TrackAnalysisUtils.cc +++ b/offline/packages/trackbase_historic/TrackAnalysisUtils.cc @@ -308,12 +308,13 @@ namespace TrackAnalysisUtils TpcGlobalPositionWrapper globalWrapper; globalWrapper.loadNodes(topNode); globalWrapper.set_suppressCrossing(true); - TpcClusterMover mover; - auto* tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); - mover.initialize_geometry(tpccellgeo); - mover.set_verbosity(0); + auto* geometry = findNode::getClass(topNode, "ActsGeometry"); + auto* tpccellgeo = findNode::getClass(topNode, "TPCGEOMCONTAINER"); + TpcClusterMover mover; + mover.initialize_geometry(tpccellgeo, geometry); + mover.set_verbosity(0); std::vector> global_raw; diff --git a/offline/packages/trackreco/MakeActsGeometry.cc b/offline/packages/trackreco/MakeActsGeometry.cc index 7f2deec166..1c51d07283 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 - + + // 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(); + 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); @@ -304,6 +341,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) @@ -475,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(); @@ -521,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 " @@ -852,15 +890,17 @@ 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); - + // 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; @@ -870,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 = @@ -888,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; } @@ -950,7 +992,7 @@ void MakeActsGeometry::makeMmMapPairs(TrackingVolumePtr &mmVolume) continue; } - if (Verbosity()) + if (Verbosity()>1) { std::cout << "MakeActsGeometry::makeMmMapPairs - layer: " << layer << " tileid: " << tileid << std::endl; } @@ -1116,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; @@ -1216,7 +1258,9 @@ void MakeActsGeometry::makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume) TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector &world) { - // Look up TPC surface index values from world position of surface center + // 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; double layer_rad = sqrt(pow(world[0], 2) + pow(world[1], 2)); @@ -1226,6 +1270,7 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector= tpc_ref_radius_low && layer_rad < tpc_ref_radius_high) { layer = ilayer; @@ -1254,7 +1299,7 @@ TrkrDefs::hitsetkey MakeActsGeometry::getTpcHitSetKeyFromCoords(std::vector 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; + } + if (readout_mod >= m_nTpcModulesPerLayer) { std::cout << PHWHERE @@ -1277,16 +1328,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 caf7183281..4738725af7 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); @@ -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); } @@ -1013,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 auto nextSurface = surfaces.at(i + 1); + 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); 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/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..93376c1b72 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 0.0; + } + virtual double get_rot_y() const + { + PHOOL_VIRTUAL_WARN("get_rot_y()"); + return 0.0; + } + virtual double get_rot_z() const + { + PHOOL_VIRTUAL_WARN("get_rot_z()"); + return 0.0; + } + virtual double get_place_x() const + { + PHOOL_VIRTUAL_WARN("get_place_x()"); + return 0.0; + } + virtual double get_place_y() const + { + PHOOL_VIRTUAL_WARN("get_place_y()"); + return 0.0; + } + virtual double get_place_z() const + { + PHOOL_VIRTUAL_WARN("get_place_z()"); + return 0.0; + } + + 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..e500c53117 --- /dev/null +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.cc @@ -0,0 +1,623 @@ +#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 +{ + 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)) + { + 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 +{ + 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)) + { + 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 +{ + 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) + { + 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 +{ + 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) + { + 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 +{ + 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) + { + 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..773c6451b2 --- /dev/null +++ b/simulation/g4simulation/g4detectors/PHG4TpcGeomv2.h @@ -0,0 +1,158 @@ +// Tell emacs that this is a C++ source +// -*- C++ -*-. +#ifndef G4DETECTORS_PHG4TPCGEOMV2_H +#define G4DETECTORS_PHG4TPCGEOMV2_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..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)) { @@ -457,8 +460,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 +479,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 +670,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..e75a34961c 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/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..4cd13647de 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcEndCapSubsystem.cc @@ -116,13 +116,14 @@ 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.); + // angles are in rad + 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 8d802c7886..80a04450f4 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcSubsystem.cc @@ -144,6 +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.); + // angles are in radians 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/TpcClusterBuilder.cc b/simulation/g4simulation/g4tpc/TpcClusterBuilder.cc index eb39ed8033..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) @@ -330,9 +336,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);