Commit c0e83502 authored by Jens Kroeger's avatar Jens Kroeger
Browse files

Merge branch 'unique' into 'master'

Tracking4D: Unique cluster usage

See merge request !436
parents cb277807 81fcda81
Pipeline #2700547 passed with stages
in 20 minutes and 3 seconds
...@@ -25,6 +25,7 @@ The DUT plane can be excluded from the track finding. ...@@ -25,6 +25,7 @@ The DUT plane can be excluded from the track finding.
* `volume_scattering`: Select if volume scattering will be taken into account - defaults to false * `volume_scattering`: Select if volume scattering will be taken into account - defaults to false
* `volume_radiation_length`: Define the radiation length of the volume around the telescope. Defaults to dry air with a radiation length of`304.2 m` * `volume_radiation_length`: Define the radiation length of the volume around the telescope. Defaults to dry air with a radiation length of`304.2 m`
* `reject_by_roi`: If true, tracks intercepting any detector outside its ROI will be rejected. Defaults to `false`. * `reject_by_roi`: If true, tracks intercepting any detector outside its ROI will be rejected. Defaults to `false`.
* `uinque_cluster_usage`: Only use a cluster for one track - in the case of multiple assignments, the track with the best chi2/ndof is kept. Defaults to `false`
### Plots produced ### Plots produced
......
...@@ -32,6 +32,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector<std::shared_ptr<Detect ...@@ -32,6 +32,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector<std::shared_ptr<Detect
config_.setDefault<double>("volume_radiation_length", Units::get<double>(304.2, "m")); config_.setDefault<double>("volume_radiation_length", Units::get<double>(304.2, "m"));
config_.setDefault<bool>("volume_scattering", false); config_.setDefault<bool>("volume_scattering", false);
config_.setDefault<bool>("reject_by_roi", false); config_.setDefault<bool>("reject_by_roi", false);
config_.setDefault<bool>("unique_cluster_usage", false);
if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) {
config_.setDefault("time_cut_rel", 3.0); config_.setDefault("time_cut_rel", 3.0);
...@@ -62,6 +63,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector<std::shared_ptr<Detect ...@@ -62,6 +63,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector<std::shared_ptr<Detect
volume_radiation_length_ = config_.get<double>("volume_radiation_length"); volume_radiation_length_ = config_.get<double>("volume_radiation_length");
use_volume_scatterer_ = config_.get<bool>("volume_scattering"); use_volume_scatterer_ = config_.get<bool>("volume_scattering");
reject_by_ROI_ = config_.get<bool>("reject_by_roi"); reject_by_ROI_ = config_.get<bool>("reject_by_roi");
unique_cluster_usage_ = config_.get<bool>("unique_cluster_usage");
// print a warning if volumeScatterer are used as this causes fit failures // print a warning if volumeScatterer are used as this causes fit failures
// that are still not understood // that are still not understood
...@@ -402,70 +404,6 @@ StatusCode Tracking4D::run(const std::shared_ptr<Clipboard>& clipboard) { ...@@ -402,70 +404,6 @@ StatusCode Tracking4D::run(const std::shared_ptr<Clipboard>& clipboard) {
LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting"; LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting";
continue; continue;
} }
// Fill histograms
trackChi2->Fill(track->getChi2());
clustersPerTrack->Fill(static_cast<double>(track->getNClusters()));
trackChi2ndof->Fill(track->getChi2ndof());
if(!(track_model_ == "gbl")) {
trackAngleX->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).X()));
trackAngleY->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).Y()));
}
// Make residuals
auto trackClusters = track->getClusters();
for(auto& trackCluster : trackClusters) {
string detectorID = trackCluster->detectorID();
ROOT::Math::XYZPoint globalRes = track->getGlobalResidual(detectorID);
ROOT::Math::XYPoint localRes = track->getLocalResidual(detectorID);
residualsX_local[detectorID]->Fill(localRes.X());
residualsX_global[detectorID]->Fill(globalRes.X());
pullX_local[detectorID]->Fill(localRes.x() / track->getClusterFromDetector(detectorID)->errorX());
pullX_global[detectorID]->Fill(globalRes.x() / track->getClusterFromDetector(detectorID)->errorX());
pullY_local[detectorID]->Fill(localRes.Y() / track->getClusterFromDetector(detectorID)->errorY());
pullY_global[detectorID]->Fill(globalRes.Y() / track->getClusterFromDetector(detectorID)->errorY());
if(trackCluster->columnWidth() == 1) {
residualsXwidth1_local[detectorID]->Fill(localRes.X());
residualsXwidth1_global[detectorID]->Fill(globalRes.X());
} else if(trackCluster->columnWidth() == 2) {
residualsXwidth2_local[detectorID]->Fill(localRes.X());
residualsXwidth2_global[detectorID]->Fill(globalRes.X());
} else if(trackCluster->columnWidth() == 3) {
residualsXwidth3_local[detectorID]->Fill(localRes.X());
residualsXwidth3_global[detectorID]->Fill(globalRes.X());
}
residualsY_local[detectorID]->Fill(localRes.Y());
residualsY_global[detectorID]->Fill(globalRes.Y());
if(trackCluster->rowWidth() == 1) {
residualsYwidth1_local[detectorID]->Fill(localRes.Y());
residualsYwidth1_global[detectorID]->Fill(globalRes.Y());
} else if(trackCluster->rowWidth() == 2) {
residualsYwidth2_local[detectorID]->Fill(localRes.Y());
residualsYwidth2_global[detectorID]->Fill(globalRes.Y());
} else if(trackCluster->rowWidth() == 3) {
residualsYwidth3_local[detectorID]->Fill(localRes.Y());
residualsYwidth3_global[detectorID]->Fill(globalRes.Y());
}
residualsZ_global[detectorID]->Fill(globalRes.Z());
}
for(auto& detector : get_detectors()) {
if(detector->isAuxiliary()) {
continue;
}
auto det = detector->getName();
if(!kinkX.count(det)) {
LOG(WARNING) << "Skipping writing kinks due to missing init of histograms for " << det;
continue;
}
XYPoint kink = track->getKinkAt(det);
kinkX.at(det)->Fill(kink.x());
kinkY.at(det)->Fill(kink.y());
}
if(timestamp_from_.empty()) { if(timestamp_from_.empty()) {
// Improve the track timestamp by taking the average of all planes // Improve the track timestamp by taking the average of all planes
...@@ -483,12 +421,109 @@ StatusCode Tracking4D::run(const std::shared_ptr<Clipboard>& clipboard) { ...@@ -483,12 +421,109 @@ StatusCode Tracking4D::run(const std::shared_ptr<Clipboard>& clipboard) {
} }
} }
tracksPerEvent->Fill(static_cast<double>(tracks.size())); auto duplicated_hit = [this](const Track* a, const Track* b) {
for(auto d : get_detectors()) {
if(a->getClusterFromDetector(d->getName()) == b->getClusterFromDetector(d->getName()) &&
!(b->getClusterFromDetector(d->getName()) == nullptr)) {
LOG(DEBUG) << "Duplicated hit on " << d->getName() << ": rejecting track";
return true;
}
}
return false;
};
// Save the tracks on the clipboard // Save the tracks on the clipboard
if(tracks.size() > 0) { if(tracks.size() > 0) {
// if requested ensure unique usage of clusters
if(unique_cluster_usage_ && tracks.size() > 1) {
// sort by chi2:
LOG_ONCE(WARNING) << "Rejecting tracks with same hits";
std::sort(tracks.begin(), tracks.end(), [](const shared_ptr<Track> a, const shared_ptr<Track> b) {
return (a->getChi2() / a->getNdof()) < (b->getChi2() / b->getNdof());
});
// remove tracks with hit that is used twice
auto track1 = tracks.begin();
while(track1 != tracks.end()) {
auto track2 = track1 + 1;
while(track2 != tracks.end()) {
// if hit is used twice delete the track
if(duplicated_hit(track2->get(), track1->get())) {
track2 = tracks.erase(track2);
} else {
track2++;
}
}
track1++;
}
}
clipboard->putData(tracks); clipboard->putData(tracks);
} }
for(auto track : tracks) {
trackChi2->Fill(track->getChi2());
clustersPerTrack->Fill(static_cast<double>(track->getNClusters()));
trackChi2ndof->Fill(track->getChi2ndof());
if(!(track_model_ == "gbl")) {
trackAngleX->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).X()));
trackAngleY->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).Y()));
}
// Make residuals
auto trackClusters = track->getClusters();
for(auto& trackCluster : trackClusters) {
string detectorID = trackCluster->detectorID();
ROOT::Math::XYZPoint globalRes = track->getGlobalResidual(detectorID);
ROOT::Math::XYPoint localRes = track->getLocalResidual(detectorID);
residualsX_local[detectorID]->Fill(localRes.X());
residualsX_global[detectorID]->Fill(globalRes.X());
pullX_local[detectorID]->Fill(localRes.x() / track->getClusterFromDetector(detectorID)->errorX());
pullX_global[detectorID]->Fill(globalRes.x() / track->getClusterFromDetector(detectorID)->errorX());
pullY_local[detectorID]->Fill(localRes.Y() / track->getClusterFromDetector(detectorID)->errorY());
pullY_global[detectorID]->Fill(globalRes.Y() / track->getClusterFromDetector(detectorID)->errorY());
if(trackCluster->columnWidth() == 1) {
residualsXwidth1_local[detectorID]->Fill(localRes.X());
residualsXwidth1_global[detectorID]->Fill(globalRes.X());
} else if(trackCluster->columnWidth() == 2) {
residualsXwidth2_local[detectorID]->Fill(localRes.X());
residualsXwidth2_global[detectorID]->Fill(globalRes.X());
} else if(trackCluster->columnWidth() == 3) {
residualsXwidth3_local[detectorID]->Fill(localRes.X());
residualsXwidth3_global[detectorID]->Fill(globalRes.X());
}
residualsY_local[detectorID]->Fill(localRes.Y());
residualsY_global[detectorID]->Fill(globalRes.Y());
if(trackCluster->rowWidth() == 1) {
residualsYwidth1_local[detectorID]->Fill(localRes.Y());
residualsYwidth1_global[detectorID]->Fill(globalRes.Y());
} else if(trackCluster->rowWidth() == 2) {
residualsYwidth2_local[detectorID]->Fill(localRes.Y());
residualsYwidth2_global[detectorID]->Fill(globalRes.Y());
} else if(trackCluster->rowWidth() == 3) {
residualsYwidth3_local[detectorID]->Fill(localRes.Y());
residualsYwidth3_global[detectorID]->Fill(globalRes.Y());
}
residualsZ_global[detectorID]->Fill(globalRes.Z());
}
for(auto& detector : get_detectors()) {
if(detector->isAuxiliary()) {
continue;
}
auto det = detector->getName();
if(!kinkX.count(det)) {
LOG(WARNING) << "Skipping writing kinks due to missing init of histograms for " << det;
continue;
}
XYPoint kink = track->getKinkAt(det);
kinkX.at(det)->Fill(kink.x());
kinkY.at(det)->Fill(kink.y());
}
}
tracksPerEvent->Fill(static_cast<double>(tracks.size()));
LOG(DEBUG) << "End of event"; LOG(DEBUG) << "End of event";
return StatusCode::Success; return StatusCode::Success;
......
...@@ -74,6 +74,7 @@ namespace corryvreckan { ...@@ -74,6 +74,7 @@ namespace corryvreckan {
bool exclude_DUT_; bool exclude_DUT_;
bool use_volume_scatterer_; bool use_volume_scatterer_;
bool reject_by_ROI_; bool reject_by_ROI_;
bool unique_cluster_usage_;
std::vector<std::string> require_detectors_; std::vector<std::string> require_detectors_;
std::map<std::shared_ptr<Detector>, double> time_cuts_; std::map<std::shared_ptr<Detector>, double> time_cuts_;
std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_; std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment