Skip to content
23 changes: 13 additions & 10 deletions DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc
Original file line number Diff line number Diff line change
Expand Up @@ -812,6 +812,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))
Expand All @@ -826,12 +831,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)
Expand All @@ -849,8 +855,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
Expand Down Expand Up @@ -1028,7 +1034,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);

Expand All @@ -1042,9 +1047,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));
Expand All @@ -1053,10 +1060,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)
Expand Down
76 changes: 69 additions & 7 deletions RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -37,22 +39,27 @@ class CTPPSDiamondLocalTrackFitter : public edm::stream::EDProducer<> {
edm::EDGetTokenT<edm::DetSetVector<CTPPSDiamondRecHit> > recHitsToken_;
const edm::ParameterSet trk_algo_params_;
std::unordered_map<CTPPSDetId, std::unique_ptr<CTPPSDiamondTrackRecognition> > trk_algo_;
edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> ctppsGeometryEventToken_;
};

CTPPSDiamondLocalTrackFitter::CTPPSDiamondLocalTrackFitter(const edm::ParameterSet& iConfig)
: recHitsToken_(
consumes<edm::DetSetVector<CTPPSDiamondRecHit> >(iConfig.getParameter<edm::InputTag>("recHitsTag"))),
trk_algo_params_(iConfig.getParameter<edm::ParameterSet>("trackingAlgorithmParams")) {
trk_algo_params_(iConfig.getParameter<edm::ParameterSet>("trackingAlgorithmParams")),
ctppsGeometryEventToken_(esConsumes<CTPPSGeometry, VeryForwardRealGeometryRecord>()){
produces<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
}

void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
// prepare the output
auto pOut = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
auto pOutLocal = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
auto pOutGlobal = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();

edm::Handle<edm::DetSetVector<CTPPSDiamondRecHit> > 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();
Expand All @@ -65,17 +72,72 @@ void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventS
trk_algo_[detid] = std::make_unique<CTPPSDiamondTrackRecognition>(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) {
Expand Down