Commit ea7d95b0 authored by Lennart Huth's avatar Lennart Huth
Browse files

added comments from discussions with Simon

parent f04217c2
......@@ -79,7 +79,7 @@ StatusCode FileWriter::run(std::shared_ptr<Clipboard> clipboard) {
auto type_idx = block.first;
auto class_name = corryvreckan::demangle(type_idx.name());
auto class_name_full = corryvreckan::demangle(type_idx.name(), true);
LOG(TRACE) << "Received objects of type " << class_name << " in " << block.second.size() << " blocks";
LOG(TRACE) << "Received objects of type \"" << class_name << "\" in " << block.second.size() << " blocks";
// Check if these objects should be stored
if((!include_.empty() && include_.find(class_name) == include_.end()) ||
......
......@@ -18,6 +18,7 @@ Clusters from the first plane in Z (named the seed plane) are related to cluster
* `require_detectors`: Names of detectors which are required to have a cluster on the track. If a track does not have a cluster from all detectors listed here, it is rejected. If empty, no detector is required. Default is empty.
* `timestamp_from`: Defines the detector which provides the track timestamp. This detector also needs to be set as `required_detector`. If empty, the average timestamp of all clusters on the track will be used. Empty by default.
* `track_model`: Select the track model used for reconstruction. A simple line fit ignoring scattering (`straightline`) and a General-Broken-Lines (`gbl`) are currently supported. Defaults to `straightline`.
* `momentum`: Set the beam momentum. Defaults to 5 GeV
### Plots produced
......
......@@ -33,7 +33,7 @@ Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detecto
requireDetectors = m_config.getArray<std::string>("require_detectors", {""});
timestampFrom = m_config.get<std::string>("timestamp_from", {});
trackModel = m_config.get<std::string>("track_model", "straightline");
momentum = m_config.get<double>("momentum", Units::get<double>(5000, "MeV"));
momentum = m_config.get<double>("momentum", Units::get<double>(5, "GeV"));
// spatial cut, relative (x * spatial_resolution) or absolute:
if(m_config.count({"spatial_cut_rel", "spatial_cut_abs"}) > 1) {
......@@ -98,23 +98,19 @@ void Tracking4D::initialise() {
title = detectorID + " Residual X;x_{track}-x [mm];events";
residualsX[detectorID] = new TH1F("residualsX", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual X, size 1;x_{track}-x [#mu m];events";
residualsXMM[detectorID] = new TH1F("residualsXMM", title.c_str(), 500, -250, 250);
title = detectorID + " Residual X, size 1;x_{track}-x [mm];events";
title = detectorID + " Residual X, cluster column width 1;x_{track}-x [mm];events";
residualsXwidth1[detectorID] = new TH1F("residualsXwidth1", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual X, size 2;x_{track}-x [mm];events";
title = detectorID + " Residual X, cluster column width 2;x_{track}-x [mm];events";
residualsXwidth2[detectorID] = new TH1F("residualsXwidth2", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual X, size 3;x_{track}-x [mm];events";
title = detectorID + " Residual X, cluster column width 3;x_{track}-x [mm];events";
residualsXwidth3[detectorID] = new TH1F("residualsXwidth3", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual Y;y_{track}-y [mm];events";
residualsY[detectorID] = new TH1F("residualsY", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual Y;y_{track}-y [#mum];events";
residualsYMM[detectorID] = new TH1F("residualsYMM", title.c_str(), 500, -250, 250);
title = detectorID + " Residual Y, size 1;y_{track}-y [mm];events";
title = detectorID + " Residual Y, cluster row width 1;y_{track}-y [mm];events";
residualsYwidth1[detectorID] = new TH1F("residualsYwidth1", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual Y, size 2;y_{track}-y [mm];events";
title = detectorID + " Residual Y, cluster row width 2;y_{track}-y [mm];events";
residualsYwidth2[detectorID] = new TH1F("residualsYwidth2", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Residual Y, size 3;y_{track}-y [mm];events";
title = detectorID + " Residual Y, cluster row width 3;y_{track}-y [mm];events";
residualsYwidth3[detectorID] = new TH1F("residualsYwidth3", title.c_str(), 500, -0.1, 0.1);
title = detectorID + " Pull X;x_{track}-x/resolution;events";
......@@ -159,8 +155,9 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
trees[detectorID] = clusterTree;
}
// the detector always needs to be listed as we would like to add the material budget information
if(!detector->isAuxiliary())
if(!detector->isAuxiliary()) {
detectors.push_back(detectorID);
}
}
// If there are no detectors then stop trying to track
......@@ -184,6 +181,7 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
LOG(DEBUG) << "Looking at next seed cluster";
auto track = Track::Factory(trackModel);
// The track finding is based on a straight line. Therefore a refTrack to extrapolate to the next plane is used here
auto refTrack = new StraightLineTrack();
// Add the cluster to the track
track->addCluster(cluster);
......@@ -200,7 +198,7 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
// always add the material budget:
track->addMaterial(detectorID, det->materialBudget(), det->displacement().z());
LOG(TRACE) << "added MB for " << detectorID << " at z = " << det->displacement().z();
LOG(TRACE) << "added material budget for " << detectorID << " at z = " << det->displacement().z();
if(trees.count(detectorID) == 0) {
LOG(TRACE) << "Skipping detector " << det->name() << " as it has 0 clusters.";
continue;
......@@ -335,8 +333,6 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
residualsX[detectorID]->Fill(track->residual(detectorID).X());
pullX[detectorID]->Fill(track->residual(detectorID).X() / track->clusters().front()->errorX());
pullY[detectorID]->Fill(track->residual(detectorID).Y() / track->clusters().front()->errorY());
residualsXMM[detectorID]->Fill(1000 * track->residual(detectorID).X());
residualsYMM[detectorID]->Fill(1000 * track->residual(detectorID).y());
if(trackCluster->columnWidth() == 1)
residualsXwidth1[detectorID]->Fill(track->residual(detectorID).X());
if(trackCluster->columnWidth() == 2)
......@@ -353,11 +349,6 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
}
for(auto& det : detectors) {
if(!kinkX.count(det)) {
LOG(WARNING) << "Skipping writing kinks due to missing init of histograms for " << det;
continue;
}
XYPoint kink = track->kink(det);
kinkX.at(det)->Fill(kink.x());
kinkY.at(det)->Fill(kink.y());
......
......@@ -40,7 +40,6 @@ namespace corryvreckan {
std::map<std::string, TH1F*> residualsYwidth1;
std::map<std::string, TH1F*> residualsYwidth2;
std::map<std::string, TH1F*> residualsYwidth3;
std::map<std::string, TH1F*> residualsXMM, residualsYMM;
std::map<std::string, TH1F*> kinkX;
std::map<std::string, TH1F*> kinkY;
......
......@@ -11,8 +11,6 @@ using namespace gbl;
GblTrack::GblTrack() : Track() {}
GblTrack::GblTrack(const GblTrack& track) : Track(track) {
if(track.getType() != this->getType())
throw TrackModelChanged(typeid(*this), track.getType(), this->getType());
m_sorted_budgets = track.m_sorted_budgets;
}
......
......@@ -8,8 +8,6 @@ using namespace corryvreckan;
StraightLineTrack::StraightLineTrack() : Track(), m_direction(0, 0, 1.), m_state(0, 0, 0.) {}
StraightLineTrack::StraightLineTrack(const StraightLineTrack& track) : Track(track) {
if(track.getType() != this->getType())
throw TrackModelChanged(typeid(*this), track.getType(), this->getType());
m_direction = track.m_direction;
m_state = track.m_state;
}
......
......@@ -6,6 +6,10 @@ using namespace corryvreckan;
Track::Track() : m_momentum(-1) {}
Track::Track(const Track& track) : Object(track.detectorID(), track.timestamp()) {
if(track.getType() != this->getType())
throw TrackModelChanged(typeid(*this), track.getType(), this->getType());
m_isFitted = track.isFitted();
m_chi2 = track.chi2();
m_ndof = track.ndof();
......
......@@ -205,14 +205,14 @@ namespace corryvreckan {
* @param detectorID
* @return 2D residual as ROOT::Math::XYPoint
*/
ROOT::Math::XYPoint residual(std::string detectorID) { return m_residual[detectorID]; }
ROOT::Math::XYPoint residual(std::string detectorID) const { return m_residual.at(detectorID); }
/**
* @brief Get the kink at a given detector layer. This is ill defined for last and first layer
* @param detectorID
* @return 2D kink as ROOT::Math::XYPoint
*/
ROOT::Math::XYPoint kink(std::string detectorID) { return m_kink[detectorID]; }
ROOT::Math::XYPoint kink(std::string detectorID) const { return m_kink.at(detectorID); }
/**
* @brief Get the materialBudget of a detector layer
......
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