diff --git a/DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc b/DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc index 5932cd09b0976..626c730c62d3a 100644 --- a/DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc +++ b/DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc @@ -42,19 +42,56 @@ //---------------------------------------------------------------------------------------------------- // Utility for efficiency computations +// NOTE: "CTPPSDiamondLocalTrack" refers to "local reconstruction", NOT coordinate system! +// - track.x0() is in GLOBAL coordinates (output from CTPPSDiamondLocalTrackFitter) +// - track.x0Sigma() is in LOCAL coordinates (track width) +// - det->translation() is in GLOBAL coordinates (channel center) +// - det->getDiamondDimensions() are in LOCAL coordinates bool channelAlignedWithTrack(const CTPPSGeometry* geom, const CTPPSDiamondDetId& detid, - const CTPPSDiamondLocalTrack& localTrack, + const CTPPSDiamondLocalTrack& track, const float tolerance = 1) { const DetGeomDesc* det = geom->sensor(detid); - const float x_pos = det->translation().x(), - x_width = 2.0 * det->getDiamondDimensions().xHalfWidth; // parameters stand for half the size - return ((x_pos + 0.5 * x_width > localTrack.x0() - localTrack.x0Sigma() - tolerance && - x_pos + 0.5 * x_width < localTrack.x0() + localTrack.x0Sigma() + tolerance) || - (x_pos - 0.5 * x_width > localTrack.x0() - localTrack.x0Sigma() - tolerance && - x_pos - 0.5 * x_width < localTrack.x0() + localTrack.x0Sigma() + tolerance) || - (x_pos - 0.5 * x_width < localTrack.x0() - localTrack.x0Sigma() - tolerance && - x_pos + 0.5 * x_width > localTrack.x0() + localTrack.x0Sigma() + tolerance)); + const auto& rot = det->rotation(); + const auto& trans = det->translation(); + + // Channel dimensions in LOCAL coordinates + const float x_half_width_local = det->getDiamondDimensions().xHalfWidth; + const float channel_x_min_local = -x_half_width_local; + const float channel_x_max_local = +x_half_width_local; + + // Transform track position from GLOBAL to LOCAL coordinates of this channel + auto trackGlobal = CTPPSGeometry::Vector(track.x0(), track.y0(), track.z0()); + auto trackLocal = rot.Inverse() * (trackGlobal - trans); + + // Now track position is in local coordinates + const float track_x_local = trackLocal.x(); + const float track_x_sigma = track.x0Sigma(); // already in local coordinates + const float track_x_min_local = track_x_local - track_x_sigma - tolerance; + const float track_x_max_local = track_x_local + track_x_sigma + tolerance; + + // Check overlap in LOCAL coordinates (new method) + bool overlaps = (channel_x_min_local <= track_x_max_local) && (track_x_min_local <= channel_x_max_local); + + // double rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz; + // rot.GetComponents(rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz); + // double angle_z = std::atan2(rot_yx, rot_xx) * 180.0 / M_PI; + + // std::cout << "=== channelAlignedWithTrack ===" << std::endl; + // std::cout << "DetId: arm=" << detid.arm() << " station=" << detid.station() << " rp=" << detid.rp() << " plane=" << detid.plane() << " channel=" << detid.channel() << std::endl; + // std::cout << "Channel position (GLOBAL): " << trans.x() << " mm" << std::endl; + // std::cout << "Channel half-width (LOCAL): " << x_half_width_local << " mm" << std::endl; + // std::cout << "Rotation around Z axis: " << angle_z << " degrees" << std::endl; + // std::cout << "Track GLOBAL: x=" << track.x0() << ", x0sigma=" << track.x0Sigma() << std::endl; + // std::cout << "Track LOCAL (after transform): x=" << track_x_local << std::endl; + // std::cout << "Channel LOCAL X range: [" << channel_x_min_local << ", " << channel_x_max_local << "]" << std::endl; + // std::cout << "Track LOCAL X range: [" << track_x_min_local << ", " << track_x_max_local << "]" << std::endl; + // std::cout << "Track X range OLD: [" << track.x0() - track.x0Sigma() - tolerance << ", " << track.x0() + track.x0Sigma() + tolerance << "]" << std::endl; + // std::cout << "Channel X range OLD: [" << trans.x() - x_half_width_local << ", " << trans.x() + x_half_width_local << "]" << std::endl; + // std::cout << "NEW method result: " << (overlaps ? "TRUE" : "FALSE") << std::endl; + // std::cout << "==========================================" << std::endl; + + return overlaps; } namespace dds { @@ -812,6 +849,11 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet // to preprare correlation plot, we need to select tracks from nr_hr pot in 220m station const CTPPSDiamondDetId detId_220nr_hr(tracks_220nr_hr.detId()); + // extract the rotation angle + auto cosRotAngle_220nr_hr = diamRotations_.at(detId_220nr_hr) * CTPPSGeometry::Vector(1, 0, 0); + + + // selecting only tracks from 220nr station, realised as skipping tracks from 220cyl station if ((detId_220nr_hr.rp() == CTPPS_DIAMOND_CYL_RP_ID) && (detId_220nr_hr.station() == CTPPS_DIAMOND_CYL_STATION_ID)) @@ -826,12 +868,13 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet continue; // get the bins from per-pot plots - int startBin_220nr_hr = trackHistoInTimeTmp->FindBin( - track_220nr_hr.x0() - diamShifts_[detId_220nr_hr.rpId()].global - track_220nr_hr.x0Sigma()); + auto localTrackX_220nr_hr = (track_220nr_hr.x0() - diamTranslations_.at(detId_220nr_hr).x() + diamHalfWidths_.at(detId_220nr_hr) - track_220nr_hr.x0Sigma()) / cosRotAngle_220nr_hr.x(); + int startBin_220nr_hr = trackHistoInTimeTmp->FindBin(localTrackX_220nr_hr); int numOfBins_220nr_hr = 2 * track_220nr_hr.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM; for (const auto& tracks_220cyl : *diamondLocalTracks) { CTPPSDiamondDetId detId_220cyl(tracks_220cyl.detId()); + auto cosRotAngle_220cyl = diamRotations_.at(detId_220cyl) * CTPPSGeometry::Vector(1, 0, 0); // select tracks in the same arm, but belonging to the cylindrical pot // that means skipping tracks from the opposite arm (skip if 220nr_hr.arm != 220cyl.arm) @@ -849,8 +892,8 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet for (const auto& track_220cyl : tracks_220cyl) { if (!track_220cyl.isValid()) continue; - int startBin_220cyl = trackHistoTmpYAxis->FindBin( - track_220cyl.x0() - diamShifts_[detId_220cyl.rpId()].global - track_220cyl.x0Sigma()); + auto localTrackX_220cyl = (track_220cyl.x0() - diamTranslations_.at(detId_220cyl).x() + diamHalfWidths_.at(detId_220cyl) - track_220cyl.x0Sigma()) / cosRotAngle_220cyl.x(); + int startBin_220cyl = trackHistoTmpYAxis->FindBin(localTrackX_220cyl); int numOfBins_220cyl = 2 * track_220cyl.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM; // fill the correlation plot @@ -1028,7 +1071,6 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet // Using CTPPSDiamondLocalTrack for (const auto& tracks : *diamondLocalTracks) { const CTPPSDiamondDetId detId(tracks.detId()), detId_pot(detId.rpId()); - const auto& x_shift = diamShifts_.at(detId_pot); // extract the rotation angle for the diamond pot auto cosRotAngle = diamRotations_.at(detId_pot) * CTPPSGeometry::Vector(1, 0, 0); @@ -1042,9 +1084,11 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet if (potPlots_.count(detId_pot) == 0) continue; + auto localTrackX = (track.x0() - diamTranslations_.at(detId_pot).x() + diamHalfWidths_.at(detId_pot) - track.x0Sigma()) / cosRotAngle.x(); + TH2F* trackHistoOOTTmp = potPlots_[detId_pot].trackDistributionOOT->getTH2F(); TAxis* trackHistoOOTTmpYAxis = trackHistoOOTTmp->GetYaxis(); - int startBin = trackHistoOOTTmpYAxis->FindBin(track.x0() - x_shift.global - track.x0Sigma()); + int startBin = trackHistoOOTTmpYAxis->FindBin((localTrackX)); int numOfBins = 2 * track.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM; for (int i = 0; i < numOfBins; ++i) trackHistoOOTTmp->Fill(track.ootIndex(), trackHistoOOTTmpYAxis->GetBinCenter(startBin + i)); @@ -1053,10 +1097,6 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet TH1F* trackHistoInTimeTmp = potPlots_[detId_pot].trackDistribution->getTH1F(); // X coordinate of the left edge of the track in the local coordinate system - auto localTrackX = - (track.x0() - diamTranslations_.at(detId_pot).x() + diamHalfWidths_.at(detId_pot) - track.x0Sigma()) / - cosRotAngle.x(); - int startBin = trackHistoInTimeTmp->FindBin((localTrackX)); int numOfBins = 2 * track.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM; for (int i = 0; i < numOfBins; ++i) diff --git a/RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc b/RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc index d9c3af54d031d..f23ff67a29f44 100644 --- a/RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc +++ b/RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc @@ -24,6 +24,8 @@ #include "DataFormats/CTPPSReco/interface/CTPPSDiamondLocalTrack.h" #include "RecoPPS/Local/interface/CTPPSDiamondTrackRecognition.h" +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" class CTPPSDiamondLocalTrackFitter : public edm::stream::EDProducer<> { public: @@ -37,22 +39,27 @@ class CTPPSDiamondLocalTrackFitter : public edm::stream::EDProducer<> { edm::EDGetTokenT > recHitsToken_; const edm::ParameterSet trk_algo_params_; std::unordered_map > trk_algo_; + edm::ESGetToken ctppsGeometryEventToken_; }; CTPPSDiamondLocalTrackFitter::CTPPSDiamondLocalTrackFitter(const edm::ParameterSet& iConfig) : recHitsToken_( consumes >(iConfig.getParameter("recHitsTag"))), - trk_algo_params_(iConfig.getParameter("trackingAlgorithmParams")) { + trk_algo_params_(iConfig.getParameter("trackingAlgorithmParams")), + ctppsGeometryEventToken_(esConsumes()){ produces >(); } void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { // prepare the output - auto pOut = std::make_unique >(); + auto pOutLocal = std::make_unique >(); + auto pOutGlobal = std::make_unique >(); edm::Handle > recHits; iEvent.getByToken(recHitsToken_, recHits); + const CTPPSGeometry* ctppsGeometry = &iSetup.getData(ctppsGeometryEventToken_); + // clear all hits possibly inherited from previous event for (auto& algo_vs_id : trk_algo_) algo_vs_id.second->clear(); @@ -65,17 +72,72 @@ void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventS trk_algo_[detid] = std::make_unique(trk_algo_params_); for (const auto& hit : vec) // skip hits without a leading edge - if (hit.ootIndex() != CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING) - trk_algo_[detid]->addHit(hit); + if (hit.ootIndex() != CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING){ + + // create a copy of the hit to transform it to the local coordinates in X and Y + CTPPSDiamondRecHit hitCopy(hit); + auto localVector = CTPPSGeometry::Vector(hit.x(), hit.y(), hit.z()); + const auto diam = ctppsGeometry->sensor(detid); + // do the global to local transformation + // Global = Rotation * Local + Translation + // Local = Rotation^-1 * (Global - Translation) + localVector -= diam->translation(); + localVector = diam->rotation().Inverse() * localVector; + hitCopy.setX(localVector.x()); + hitCopy.setY(localVector.y()); + // print X,Y,Z coordinates before and after transformation + // std::cout << "Hit before transformation: X=" << hit.x() << ", Y=" << hit.y() << ", Z=" << hit.z() << std::endl; + // std::cout << "Hit after transformation: X=" << localVector.x() << std::endl; + + trk_algo_[detid]->addHit(hitCopy); + } + } - // build the tracks for all stations + // build the local tracks for all stations for (auto& algo_vs_id : trk_algo_) { - auto& tracks = pOut->find_or_insert(algo_vs_id.first); + auto& tracks = pOutLocal->find_or_insert(algo_vs_id.first); algo_vs_id.second->produceTracks(tracks); } - iEvent.put(std::move(pOut)); + for (const auto& vec : *pOutLocal) { + const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp()); + const auto diam = ctppsGeometry->sensor(detid); + pOutGlobal->find_or_insert(detid); + for (const auto& track : vec) { + auto trackCopy = track; // make a copy of the track to transform it + auto globalVector = CTPPSGeometry::Vector(track.x0(), track.y0(), track.z0()); + + // do the global to local transformation, doing first the rotation and then the translation + globalVector = diam->rotation() * globalVector; + globalVector += diam->translation(); + auto globalPoint = math::XYZPoint(globalVector.x(), globalVector.y(), track.z0()); + + trackCopy.setPosition(globalPoint); + pOutGlobal->operator[](detid).push_back(trackCopy); + } + } + + // // print the local tracks + // std::cout << "Local tracks for event " << iEvent.id() << ":" << std::endl; + // for (const auto& vec : *pOutLocal) { + // const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp()); + // std::cout << " DetId: " << detid << std::endl; + // for (const auto& track : vec) { + // std::cout << " Track: x0=" << track.x0() << ", y0=" << track.y0() << ", z0=" << track.z0() << std::endl; + // } + // } + // // print the global tracks + // std::cout << "Global tracks for event " << iEvent.id() << ":" << std::endl; + // for (const auto& vec : *pOutGlobal) { + // const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp()); + // std::cout << " DetId: " << detid << std::endl; + // for (const auto& track : vec) { + // std::cout << " Track: x0=" << track.x0() << ", y0=" << track.y0() << ", z0=" << track.z0() << std::endl; + // } + // } + + iEvent.put(std::move(pOutGlobal)); } void CTPPSDiamondLocalTrackFitter::fillDescriptions(edm::ConfigurationDescriptions& descr) { diff --git a/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc b/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc index 43c3a0a5b2663..c3635daab958a 100644 --- a/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc +++ b/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc @@ -372,23 +372,50 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe if (timing_RP_track_multiplicity[tr_ti.rpId()] > max_n_timing_tracks_) continue; - // interpolation from tracking RPs + // Interpolation in GLOBAL coordinates from tracking RPs (pixels) const double z_ti = hGeometry->rpTranslation(tr_ti.rpId()).z(); const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j); - const double x_inter = f_i * tr_i.x() + f_j * tr_j.x(); + const double x_inter_global = f_i * tr_i.x() + f_j * tr_j.x(); + const double y_inter_global = f_i * tr_i.y() + f_j * tr_j.y(); const double x_inter_unc_sq = f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc(); - const double de_x = tr_ti.x() - x_inter; - const double de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq); + // Get timing detector geometry for coordinate transformation + const DetGeomDesc* det = hGeometry->sensor(tr_ti.rpId()); + const auto& rot = det->rotation(); + const auto& trans = det->translation(); + + // Transform interpolated point from GLOBAL to LOCAL coordinates of timing (diamond) detector + auto interGlobal = CTPPSGeometry::Vector(x_inter_global, y_inter_global, z_ti); + auto interLocal = rot.Inverse() * (interGlobal - trans); + const double x_inter_local = interLocal.x(); + + // Transform timing track from GLOBAL to LOCAL coordinates + auto trackGlobal = CTPPSGeometry::Vector(tr_ti.x(), tr_ti.y(), z_ti); + auto trackLocal = rot.Inverse() * (trackGlobal - trans); + const double x_ti_local = trackLocal.x(); + + // Compare in LOCAL coordinates (both now in same rotated frame) + const double de_x = x_ti_local - x_inter_local; + + // The track uncertainty tr_ti.xUnc() is already in LOCAL coordinates + const double x_inter_unc_local = sqrt(x_inter_unc_sq); + const double de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_local * x_inter_unc_local); const double r = (de_x_unc > 0.) ? de_x / de_x_unc : 1E100; const bool matching = (ac.getTiTrMin() < r && r < ac.getTiTrMax()); - if (verbosity_) - ssLog << "ti=" << ti << ", i=" << i << ", j=" << j << " | z_ti=" << z_ti << ", z_i=" << z_i - << ", z_j=" << z_j << " | x_ti=" << tr_ti.x() << ", x_inter=" << x_inter << ", de_x=" << de_x - << ", de_x_unc=" << de_x_unc << ", matching=" << matching << std::endl; + if (verbosity_) { + double rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz; + rot.GetComponents(rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz); + double angle_z = std::atan2(rot_yx, rot_xx) * 180.0 / M_PI; + ssLog << "ti=" << ti << ", i=" << i << ", j=" << j + << " | angle_z=" << angle_z << " deg" + << " | x_ti_global=" << tr_ti.x() << ", x_ti_local=" << x_ti_local + << " | x_inter_global=" << x_inter_global << ", x_inter_local=" << x_inter_local + << " | de_x=" << de_x << ", de_x_unc=" << de_x_unc + << ", r=" << r << ", matching=" << matching << std::endl; + } if (!matching) continue;