diff --git a/src/core/Corryvreckan.cpp b/src/core/Corryvreckan.cpp index 50705e63da8113e990441e617aa2c98fbeefb7ba..61817f97bb7ce180b73b18e8932ea861483c886e 100644 --- a/src/core/Corryvreckan.cpp +++ b/src/core/Corryvreckan.cpp @@ -238,10 +238,11 @@ void Corryvreckan::add_units() { // NOTE: fixed by above Units::add("V", 1e-6); Units::add("kV", 1e-3); + Units::add("MV", 1); // MAGNETIC FIELD - Units::add("T", 1e-3); - Units::add("mT", 1e-6); + Units::add("T", 1); + Units::add("mT", 1e-3); // ANGLES // NOTE: these are fake units diff --git a/src/core/detector/PixelDetector.cpp b/src/core/detector/PixelDetector.cpp index 90a28d6ad943e0e7107e1585f553494eb7bd9e1c..4f4f34af9aaa958e81eb58c2db8bc85f0c7e238f 100644 --- a/src/core/detector/PixelDetector.cpp +++ b/src/core/detector/PixelDetector.cpp @@ -208,7 +208,7 @@ void PixelDetector::configure_pos_and_orientation(Configuration& config) const { PositionVector3D> PixelDetector::getIntercept(const Track* track) const { // FIXME: this is else statement can only be temporary - if(track->getType() == "GblTrack") { + if(track->getType() == "GblTrack" || track->getType() == "HelixTrack") { return track->getState(getName()); } else { // Get the distance from the plane to the track initial state diff --git a/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp b/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp index 5e752386838003b21f270fbd91bec563a12e5b64..21579c778e7f2d46f1dfb8bcb159a38d085b6687 100644 --- a/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp +++ b/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp @@ -133,7 +133,7 @@ void AlignmentTrackChi2::MinimiseTrackChi2(Int_t&, Double_t*, Double_t& result, // Refit the track track->registerPlane(AlignmentTrackChi2::globalDetector->getName(), - AlignmentTrackChi2::globalDetector->displacement().z(), + AlignmentTrackChi2::globalDetector->origin(), AlignmentTrackChi2::globalDetector->materialBudget(), AlignmentTrackChi2::globalDetector->toLocal()); LOG(DEBUG) << "Updated transformations for detector " << AlignmentTrackChi2::globalDetector->getName(); diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp index 67b50d38ac3fb8b1301b6ecfa57cecc046f7b9bb..e388ce874d0f32d21802bc971b6da03fea182ccc 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp @@ -38,27 +38,35 @@ void AnalysisTelescope::initialize() { if(detector->isDUT()) { std::string title = detector->getName() + " Telescope resolution X;x_{track}-x_{MC} [mm];events"; - telescopeResolutionX[detector->getName()] = new TH1F("telescopeResolutionX", title.c_str(), 600, -0.2, 0.2); + telescopeResolutionX[detector->getName()] = new TH1F("telescopeResolutionX", title.c_str(), 600, -3*detector->getPitch().X(), 3*detector->getPitch().X()); title = detector->getName() + " Telescope resolution Y;y_{track}-y_{MC} [mm];events"; - telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -0.2, 0.2); + telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); } else { std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{cluster} [mm];events"; - telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -0.2, 0.2); + telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); title = detector->getName() + " Biased residual Y (local);y_{track}-y_{cluster} [mm];events"; - telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -0.2, 0.2); + telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); title = detector->getName() + " Biased residual X (global);x_{track}-x_{cluster} [mm];events"; telescopeResidualsX[detector->getName()] = new TH1F("residualX_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " Biased residual Y (global);y_{track}-y_{cluster} [mm];events"; telescopeResidualsY[detector->getName()] = new TH1F("residualY_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual X (local);x_{track}-x_{MC} [mm];events"; - telescopeMCresidualsLocalX[detector->getName()] = new TH1F("residualX_MC_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual Y (local);y_{track}-y_{MC} [mm];events"; - telescopeMCresidualsLocalY[detector->getName()] = new TH1F("residualY_MC_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual X (global);x_{track}-x_{MC} [mm];events"; - telescopeMCresidualsX[detector->getName()] = new TH1F("residualX_MC_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual Y (global);y_{track}-y_{MC} [mm];events"; - telescopeMCresidualsY[detector->getName()] = new TH1F("residualY_MC_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC cluster residual X (local);x_{cluster}-x_{MC} [mm];events"; + telescopeMCresidualsClusterLocalX[detector->getName()] = new TH1F("residualX_MC_cluster_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); + title = detector->getName() + " Biased MC cluster residual Y (local);y_{cluster}-y_{MC} [mm];events"; + telescopeMCresidualsClusterLocalY[detector->getName()] = new TH1F("residualY_MC_cluster_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); + title = detector->getName() + " Biased MC track residual X (local);x_{track}-x_{MC} [mm];events"; + telescopeMCresidualsTrackLocalX[detector->getName()] = new TH1F("residualX_MC_track_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); + title = detector->getName() + " Biased MC track residual Y (local);y_{track}-y_{MC} [mm];events"; + telescopeMCresidualsTrackLocalY[detector->getName()] = new TH1F("residualY_MC_track_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); + title = detector->getName() + " Biased MC cluster residual X (global);x_{cluster}-x_{MC} [mm];events"; + telescopeMCresidualsClusterX[detector->getName()] = new TH1F("residualX_MC_cluster_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC cluster residual Y (global);y_{cluster}-y_{MC} [mm];events"; + telescopeMCresidualsClusterY[detector->getName()] = new TH1F("residualY_MC_cluster_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC track residual X (global);x_{track}-x_{MC} [mm];events"; + telescopeMCresidualsTrackX[detector->getName()] = new TH1F("residualX_MC_track_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC track residual Y (global);y_{track}-y_{MC} [mm];events"; + telescopeMCresidualsTrackY[detector->getName()] = new TH1F("residualY_MC_track_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " pixel - seed timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events"; @@ -145,7 +153,7 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { if(track->getType() != "GblTrack") { intercept = track->getIntercept(cluster->global().z()); } - auto interceptLocal = detector->globalToLocal(intercept); + auto interceptLocal = detector->getLocalIntercept(track.get()); telescopeResidualsLocalX[name]->Fill(cluster->local().x() - interceptLocal.X()); telescopeResidualsLocalY[name]->Fill(cluster->local().y() - interceptLocal.Y()); telescopeResidualsX[name]->Fill(cluster->global().x() - intercept.X()); @@ -179,13 +187,24 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { continue; } - ROOT::Math::XYZPoint particlePosition = closestApproach(cluster->local(), mcParticles); - telescopeMCresidualsLocalX[name]->Fill(cluster->local().x() + detector->getSize().X() / 2 - - particlePosition.X()); - telescopeMCresidualsLocalY[name]->Fill(cluster->local().y() + detector->getSize().Y() / 2 - - particlePosition.Y()); - telescopeMCresidualsX[name]->Fill(interceptLocal.X() + detector->getSize().X() / 2 - particlePosition.X()); - telescopeMCresidualsY[name]->Fill(interceptLocal.Y() + detector->getSize().Y() / 2 - particlePosition.Y()); + ROOT::Math::XYZPoint clusterPosLocalConv( + cluster->local().X()+detector->getSize().X()/2-detector->getPitch().X()/2, + cluster->local().Y()+detector->getSize().Y()/2-detector->getPitch().Y()/2, + cluster->local().Z()); + ROOT::Math::XYZPoint particlePositionConv = closestApproach(clusterPosLocalConv, mcParticles); + ROOT::Math::XYZPoint particlePosition( + particlePositionConv.X() - detector->getSize().X()/2 + detector->getPitch().X()/2, + particlePositionConv.Y() - detector->getSize().Y()/2 + detector->getPitch().Y()/2, + particlePositionConv.Z()); + ROOT::Math::XYZPoint particlePositionGlobal = detector->localToGlobal(particlePosition); + telescopeMCresidualsClusterLocalX[name]->Fill(cluster->local().X() - particlePosition.X()); + telescopeMCresidualsClusterLocalY[name]->Fill(cluster->local().Y() - particlePosition.Y()); + telescopeMCresidualsTrackLocalX[name]->Fill(interceptLocal.X()- particlePosition.X()); + telescopeMCresidualsTrackLocalY[name]->Fill(interceptLocal.Y()- particlePosition.Y()); + telescopeMCresidualsClusterX[name]->Fill(cluster->global().X() - particlePositionGlobal.X()); + telescopeMCresidualsClusterY[name]->Fill(cluster->global().X() - particlePositionGlobal.X()); + telescopeMCresidualsTrackX[name]->Fill(intercept.X() - particlePositionGlobal.X()); + telescopeMCresidualsTrackY[name]->Fill(intercept.X() - particlePositionGlobal.X()); } // Calculate telescope resolution at DUTs @@ -196,14 +215,19 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { continue; } - auto intercept = detector->getIntercept(track.get()); - auto interceptLocal = detector->globalToLocal(intercept); - auto particlePosition = closestApproach(interceptLocal, mcParticles); - - telescopeResolutionX[detector->getName()]->Fill(interceptLocal.X() + detector->getSize().X() / 2 - - particlePosition.X()); - telescopeResolutionY[detector->getName()]->Fill(interceptLocal.Y() + detector->getSize().Y() / 2 - - particlePosition.Y()); + auto interceptLocal = detector->getLocalIntercept(track.get()); + ROOT::Math::XYZPoint interceptLocalConv( + interceptLocal.X()+detector->getSize().X()/2-detector->getPitch().X()/2, + interceptLocal.Y()+detector->getSize().Y()/2-detector->getPitch().Y()/2, + interceptLocal.Z()); + auto particlePositionConv = closestApproach(interceptLocalConv, mcParticles); + ROOT::Math::XYZPoint particlePosition( + particlePositionConv.X() - detector->getSize().X()/2 + detector->getPitch().X()/2, + particlePositionConv.Y() - detector->getSize().Y()/2 + detector->getPitch().Y()/2, + particlePositionConv.Z()); + + telescopeResolutionX[detector->getName()]->Fill(interceptLocal.X() - particlePosition.X()); + telescopeResolutionY[detector->getName()]->Fill(interceptLocal.Y() - particlePosition.Y()); } } diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.h b/src/modules/AnalysisTelescope/AnalysisTelescope.h index 30e52fa456edb3980dbe500137509cac5205708b..671ac86ab19d5da2cc9d263fe0f10689bf5c9416 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.h +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.h @@ -35,10 +35,14 @@ namespace corryvreckan { ROOT::Math::XYZPoint closestApproach(ROOT::Math::XYZPoint position, const MCParticleVector& particles); // Histograms for each of the devices - std::map telescopeMCresidualsLocalX; - std::map telescopeMCresidualsLocalY; - std::map telescopeMCresidualsX; - std::map telescopeMCresidualsY; + std::map telescopeMCresidualsClusterLocalX; + std::map telescopeMCresidualsClusterLocalY; + std::map telescopeMCresidualsTrackLocalX; + std::map telescopeMCresidualsTrackLocalY; + std::map telescopeMCresidualsClusterX; + std::map telescopeMCresidualsClusterY; + std::map telescopeMCresidualsTrackX; + std::map telescopeMCresidualsTrackY; std::map telescopeResidualsLocalX; std::map telescopeResidualsLocalY; diff --git a/src/modules/Prealignment/Prealignment.cpp b/src/modules/Prealignment/Prealignment.cpp index 262ba4b96ec7d085e92fd48c11f315df8c1e4dda..de8cdc645f823b1de80a08e90915d511ec4afd1e 100644 --- a/src/modules/Prealignment/Prealignment.cpp +++ b/src/modules/Prealignment/Prealignment.cpp @@ -9,6 +9,7 @@ */ #include "Prealignment.h" +#include #include "tools/cuts.h" using namespace corryvreckan; @@ -25,6 +26,10 @@ Prealignment::Prealignment(Configuration& config, std::shared_ptr dete config_.setDefault("method", PrealignMethod::MEAN); config_.setDefault("fit_range_rel", 500); config_.setDefault("range_abs", Units::get(10, "mm")); + config_.setDefault("magnetic_field", ROOT::Math::XYZVector()); + config_.setDefault("momentum", Units::get(5, "GeV")); + config_.setDefault("charge", Units::get(1, "e")); + config_.setDefault("initial_angle", Units::get(0, "deg")); if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { config_.setDefault("time_cut_rel", 3.0); @@ -38,6 +43,10 @@ Prealignment::Prealignment(Configuration& config, std::shared_ptr dete range_abs = config_.get("range_abs"); method = config_.get("method"); fit_range_rel = config_.get("fit_range_rel"); + magnetic_field = config_.get("magnetic_field"); + momentum = config_.get("momentum"); + charge = config_.get("charge"); + initial_angle = config_.get("initial_angle"); LOG(DEBUG) << "Setting max_correlation_rms to : " << max_correlation_rms; LOG(DEBUG) << "Setting damping_factor to : " << damping_factor; @@ -167,6 +176,44 @@ void Prealignment::finalize(const std::shared_ptr&) { shift_Y = correlationY->GetXaxis()->GetBinCenter(binMaxY); } + if(magnetic_field.R() != 0) { + ROOT::Math::XYZVector globalZ(0, 0, 1); + auto helixZ = magnetic_field.unit(); + auto helixY = helixZ.Theta() != 0 ? helixZ.Cross(globalZ).unit() : ROOT::Math::XYZVector(0, 1, 0); + auto helixX = helixY.Cross(helixZ).unit(); + + ROOT::Math::Rotation3D toGlobal(helixX, helixY, helixZ); + auto toHelix = toGlobal.Inverse(); + + // calculate bending radius (in mm) from magnetic field, momentum and charge + double R = momentum / charge / magnetic_field.R() * globalZ.Cross(magnetic_field.unit()).R() * 10. / 3.; + + ROOT::Math::XYZVector planePosGlobal(0., 0., m_detector->displacement().Z()); + ROOT::Math::XYZVector refPosGlobal(0., 0., get_reference()->displacement().Z()); + auto planePosHelix = toHelix * planePosGlobal; + auto refPosHelix = toHelix * refPosGlobal; + planePosHelix.SetX(planePosHelix.x() + R * sin(initial_angle)); + refPosHelix.SetX(refPosHelix.x() + R * sin(initial_angle)); + + double initialShift = -(charge / abs(charge)) * R * (1 - cos(initial_angle)); + double planeShiftHelix = + -(charge / abs(charge)) * (R - sqrt(R * R - planePosHelix.x() * planePosHelix.x())) - initialShift; + double refShiftHelix = + -(charge / abs(charge)) * (R - sqrt(R * R - refPosHelix.x() * refPosHelix.x())) - initialShift; + auto planeShiftGlobal = toGlobal * ROOT::Math::XYZVector(0., planeShiftHelix, 0.); + auto refShiftGlobal = toGlobal * ROOT::Math::XYZVector(0., refShiftHelix, 0.); + + auto magneticShift_X = planeShiftGlobal.x() - refShiftGlobal.x(); + auto magneticShift_Y = planeShiftGlobal.y() - refShiftGlobal.y(); + + shift_X += magneticShift_X; + shift_Y += magneticShift_Y; + + LOG(INFO) << "Shift (from magnetic field)" << m_detector->getName() + << ": x = " << Units::display(magneticShift_X, {"mm", "um"}) + << " , y = " << Units::display(magneticShift_Y, {"mm", "um"}); + } + LOG(DEBUG) << "Shift (without damping factor)" << m_detector->getName() << ": x = " << Units::display(shift_X, {"mm", "um"}) << " , y = " << Units::display(shift_Y, {"mm", "um"}); diff --git a/src/modules/Prealignment/Prealignment.h b/src/modules/Prealignment/Prealignment.h index 4a584f97d0c623ea798cb38666dba16a90a063b8..9ebc935a22f0b80ea195b52e71429e16999c3cd4 100644 --- a/src/modules/Prealignment/Prealignment.h +++ b/src/modules/Prealignment/Prealignment.h @@ -11,6 +11,7 @@ #ifndef PREALIGNMENT_H #define PREALIGNMENT_H 1 +#include #include #include #include @@ -59,6 +60,10 @@ namespace corryvreckan { double range_abs; PrealignMethod method; int fit_range_rel; + ROOT::Math::XYZVector magnetic_field; + double momentum; + double charge; + double initial_angle; }; } // namespace corryvreckan #endif // PREALIGNMENT_H diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index 71ac995d9f7beb75e961e4e6c219feaf61520bc0..bd99a2cf7f8a54b4fd207afb832847aa488f54a1 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -24,16 +24,20 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("min_hits_on_track", 6); config_.setDefault("exclude_dut", true); - config_.setDefault("track_model", "straightline"); + config_.setDefault("track_model", "helix"); config_.setDefault("momentum", Units::get(5, "GeV")); config_.setDefault("max_plot_chi2", 50.0); config_.setDefault("volume_radiation_length", Units::get(304.2, "m")); config_.setDefault("volume_scattering", false); config_.setDefault("reject_by_roi", false); config_.setDefault("unique_cluster_usage", false); + config_.setDefault("magnetic_field", ROOT::Math::XYZVector()); + config_.setDefault("radius_tol_min", 1.0); + config_.setDefault("radius_tol_rel", 0.1); if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { config_.setDefault("time_cut_rel", 3.0); @@ -67,6 +71,28 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("volume_scattering"); reject_by_ROI_ = config_.get("reject_by_roi"); unique_cluster_usage_ = config_.get("unique_cluster_usage"); + magnetic_field_ = config_.get("magnetic_field"); + radius_tol_min_ = config_.get("radius_tol_min"); + radius_tol_rel_ = config_.get("radius_tol_rel"); + + expected_radius_ = 0.; + if(magnetic_field_ != ROOT::Math::XYZVector()) { + expected_radius_ = momentum_ * + sqrt(magnetic_field_.x() * magnetic_field_.x() + magnetic_field_.y() * magnetic_field_.y()) / + (magnetic_field_.R() * magnetic_field_.R() * 0.3); // R[mm] = pT[MeV]/(B[T]*c)*1000 + } + radius_tol_ = std::max(radius_tol_min_, expected_radius_ * radius_tol_rel_); + + // the folowing combinations of track model and magnetic field are not supported: + // straightline with magnetic field + // helix without magnetic field + if(track_model_ == "straightline" && magnetic_field_ != ROOT::Math::XYZVector()) { + throw InvalidCombinationError( + config_, {"track_model", "magnetic_field"}, track_model_ + " does not support magnetic fields."); + } else if(track_model_ == "helix" && magnetic_field_ == ROOT::Math::XYZVector()) { + throw InvalidCombinationError( + config_, {"track_model", "magnetic_field"}, track_model_ + " can only be used with magnetic fields."); + } // print a warning if volumeScatterer are used as this causes fit failures // that are still not understood @@ -96,6 +122,12 @@ void Tracking4D::initialize() { trackAngleX = new TH1F("trackAngleX", title.c_str(), 2000, -0.01, 0.01); title = "Track angle Y;angle_{y} [rad];events"; trackAngleY = new TH1F("trackAngleY", title.c_str(), 2000, -0.01, 0.01); + title = "Track radius;radius [m];events"; + trackRadius = new TH1F("trackRadius", + title.c_str(), + 500, + std::max(0., (expected_radius_ - radius_tol_) / 1000), + (expected_radius_ + radius_tol_) / 1000); title = "Track time within event;track time - event start;events"; trackTime = new TH1F("trackTime", title.c_str(), 1000, 0, 460.8); title = "Track time with respect to first trigger;track time - trigger;events"; @@ -262,7 +294,16 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { // Container for all clusters, and detectors in tracking map, KDTree> trees; - std::shared_ptr reference_first, reference_last; + std::vector> seed_detectors; + // Number of detectors needed for seeding (2 without magnetic field, 3 with magnetic field) + size_t num_seed_detectors; + if(magnetic_field_ == ROOT::Math::XYZVector()) { + LOG(DEBUG) << "No magnetic field. Seeding tracks from two detectors."; + num_seed_detectors = 2; + } else { + LOG(DEBUG) << "Mangetic field. Seeding track from three detectors."; + num_seed_detectors = 3; + } for(auto& detector : get_regular_detectors(!exclude_DUT_)) { // Get the clusters auto tempClusters = clipboard->getData(detector->getName()); @@ -277,10 +318,11 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { // Get first and last detectors with clusters on them: if(std::find(exclude_from_seed_.begin(), exclude_from_seed_.end(), detector->getName()) == exclude_from_seed_.end()) { - if(!reference_first) { - reference_first = detector; + if(seed_detectors.size() < num_seed_detectors) { + seed_detectors.push_back(detector); + } else { + seed_detectors.back() = detector; } - reference_last = detector; } else { LOG(DEBUG) << "Not using " << detector->getName() << " as seed as chosen by config file."; } @@ -288,7 +330,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { } // If there are no detectors then stop trying to track - if(trees.size() < 2) { + if(seed_detectors.size() < num_seed_detectors) { // Fill histogram tracksPerEvent->Fill(0); @@ -300,201 +342,286 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { TrackVector tracks; // Time cut for combinations of reference clusters and for reference track with additional detector - auto time_cut_ref = std::max(time_cuts_[reference_first], time_cuts_[reference_last]); - auto time_cut_ref_track = std::min(time_cuts_[reference_first], time_cuts_[reference_last]); - for(auto& clusterFirst : trees[reference_first].getAllElements()) { - for(auto& clusterLast : trees[reference_last].getAllElements()) { - LOG(DEBUG) << "Looking at next reference cluster pair"; - - if(std::fabs(clusterFirst->timestamp() - clusterLast->timestamp()) > time_cut_ref) { - LOG(DEBUG) << "Reference clusters not within time cuts."; - continue; + double time_cut_ref = std::numeric_limits::min(); + double time_cut_ref_track = std::numeric_limits::max(); + std::map, std::vector>> seed_clusters; + std::map, std::vector>::iterator> seed_clusters_it; + for(auto seedDetector : seed_detectors) { + time_cut_ref = std::max(time_cut_ref, time_cuts_[seedDetector]); + time_cut_ref_track = std::min(time_cut_ref_track, time_cuts_[seedDetector]); + seed_clusters[seedDetector] = trees[seedDetector].getAllElements(); + seed_clusters_it[seedDetector] = seed_clusters[seedDetector].begin(); + } + + // lambda to iterate through all combinations of elements of variable number of vectors + auto iterate = [&seed_detectors, &seed_clusters, &seed_clusters_it]() { + ++seed_clusters_it[seed_detectors.front()]; + for(auto seedDetectorsIt = seed_detectors.begin(); seedDetectorsIt != seed_detectors.end() - 1; seedDetectorsIt++) { + if(seed_clusters_it[*seedDetectorsIt] != seed_clusters[*seedDetectorsIt].end()) { + break; } + seed_clusters_it[*seedDetectorsIt] = seed_clusters[*seedDetectorsIt].begin(); + ++seed_clusters_it[*(seedDetectorsIt + 1)]; + } + }; - // The track finding is based on a straight line. Therefore a refTrack to extrapolate to the next plane is used - StraightLineTrack refTrack; - refTrack.addCluster(clusterFirst.get()); - refTrack.addCluster(clusterLast.get()); - auto averageTimestamp = calculate_average_timestamp(&refTrack); - refTrack.setTimestamp(averageTimestamp); + for(; seed_clusters_it[seed_detectors.back()] != seed_clusters[seed_detectors.back()].end(); iterate()) { + LOG(DEBUG) << "Looking at next reference cluster combination"; - // Make a new track - auto track = Track::Factory(track_model_); - track->addCluster(clusterFirst.get()); - track->addCluster(clusterLast.get()); + std::map, std::shared_ptr> seedClusters; + for(auto& seedDetector : seed_detectors) { + seedClusters[seedDetector] = *seed_clusters_it[seedDetector]; + } - track->setTimestamp(averageTimestamp); - if(use_volume_scatterer_) { - track->setVolumeScatter(volume_radiation_length_); + auto checkSeedTimeCuts = [&seed_detectors, + time_cut_ref](std::map, std::shared_ptr>& clusters) { + auto seedDetector1 = seed_detectors.begin(); + while(seedDetector1 != seed_detectors.end()) { + auto seedDetector2 = seedDetector1 + 1; + while(seedDetector2 != seed_detectors.end()) { + if(std::fabs(clusters[*seedDetector1]->timestamp() - clusters[*seedDetector2]->timestamp()) > + time_cut_ref) { + LOG(DEBUG) << "Reference clusters on " << (*seedDetector1)->getName() << " and " + << (*seedDetector2)->getName() << " not within time cuts."; + return false; + } + seedDetector2++; + } + seedDetector1++; } - track->setParticleMomentum(momentum_); + return true; + }; - // Loop over each subsequent plane and look for a cluster within the timing cuts - size_t detector_nr = 2; - // Get all detectors here to also include passive layers which might contribute to scattering - for(auto& detector : get_detectors()) { - if(detector->isAuxiliary()) { - continue; - } - auto detectorID = detector->getName(); - LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); - track->registerPlane( - detectorID, detector->displacement().z(), detector->materialBudget(), detector->toLocal()); + if(!checkSeedTimeCuts(seedClusters)) { + continue; + } - if(detector == reference_first || detector == reference_last) { - continue; - } + // The track finding is based on a straight line or a helix. Therefore a seedTrack to extrapolate to the next plane + // is used + std::string seedTrackModel; + if(magnetic_field_ == ROOT::Math::XYZVector()) { + LOG(DEBUG) << "Seeding tracks based on straight line model."; + seedTrackModel = "straightline"; + } else { + LOG(DEBUG) << "Seeding tracks based on helix model."; + seedTrackModel = "helix"; + } - if(exclude_DUT_ && detector->isDUT()) { - LOG(DEBUG) << "Skipping DUT plane."; - continue; - } + // Make new seed track and new track + auto seedTrack = Track::Factory(seedTrackModel); + auto track = Track::Factory(track_model_); - if(detector->isPassive()) { - LOG(DEBUG) << "Skipping passive plane."; - continue; - } + for(auto& seedDetector : seed_detectors) { + seedTrack->addCluster(seedClusters[seedDetector].get()); + seedTrack->registerPlane( + seedDetector->getName(), seedDetector->origin(), seedDetector->materialBudget(), seedDetector->toLocal()); - // Determine whether a track can still be assembled given the number of current hits and the number of - // detectors to come. Reduces computing time. - detector_nr++; - if(refTrack.getNClusters() + (trees.size() - detector_nr + 1) < min_hits_on_track_) { - LOG(DEBUG) << "No chance to find a track - too few detectors left: " << refTrack.getNClusters() << " + " - << trees.size() << " - " << detector_nr << " < " << min_hits_on_track_; - continue; - } + track->addCluster(seedClusters[seedDetector].get()); + } - if(trees.count(detector) == 0) { - LOG(TRACE) << "Skipping detector " << detector->getName() << " as it has 0 clusters."; - continue; - } + auto averageTimestamp = calculate_average_timestamp(seedTrack.get()); + seedTrack->setTimestamp(averageTimestamp); - // Get all neighbors within the timing cut - LOG(DEBUG) << "Searching for neighboring cluster on device " << detector->getName(); - LOG(DEBUG) << "- reference time is " << Units::display(refTrack.timestamp(), {"ns", "us", "s"}); - Cluster* closestCluster = nullptr; + track->setTimestamp(averageTimestamp); + if(use_volume_scatterer_) { + track->setVolumeScatter(volume_radiation_length_); + } + seedTrack->setParticleMomentum(momentum_); + seedTrack->setMagneticField(magnetic_field_); + track->setParticleMomentum(momentum_); + track->setMagneticField(magnetic_field_); + + // Loop over each subsequent plane and look for a cluster within the timing cuts + size_t detector_nr = seed_detectors.size(); + // Get all detectors here to also include passive layers which might contribute to scattering + for(auto& detector : get_detectors()) { + if(detector->isAuxiliary()) { + continue; + } + auto detectorID = detector->getName(); + LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); + track->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); - // Use spatial cut only as initial value (check if cluster is ellipse defined by cuts is done below): - double closestClusterDistance = sqrt(spatial_cuts_[detector].x() * spatial_cuts_[detector].x() + - spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + if(std::find(seed_detectors.begin(), seed_detectors.end(), detector) != seed_detectors.end()) { + continue; + } - double timeCut = std::max(time_cut_ref_track, time_cuts_[detector]); - LOG(DEBUG) << "Using timing cut of " << Units::display(timeCut, {"ns", "us", "s"}); + // register plane for seed track here to avoid doubly registering seed detectors + seedTrack->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); - auto neighbors = trees[detector].getAllElementsInTimeWindow(refTrack.timestamp(), timeCut); + if(exclude_DUT_ && detector->isDUT()) { + LOG(DEBUG) << "Skipping DUT plane."; + continue; + } - LOG(DEBUG) << "- found " << neighbors.size() << " neighbors within the correct time window"; + if(detector->isPassive()) { + LOG(DEBUG) << "Skipping passive plane."; + continue; + } - // Now look for the spatially closest cluster on the next plane - refTrack.fit(); + // Determine whether a track can still be assembled given the number of current hits and the number of + // detectors to come. Reduces computing time. If this check fails, one can directly continue with the next seed + detector_nr++; + if(seedTrack->getNClusters() + (trees.size() - detector_nr + 1) < min_hits_on_track_) { + LOG(DEBUG) << "No chance to find a track - too few detectors left: " << seedTrack->getNClusters() << " + " + << trees.size() << " - " << detector_nr << " < " << min_hits_on_track_; + break; + } - PositionVector3D> interceptPoint = detector->getLocalIntercept(&refTrack); - double interceptX = interceptPoint.X(); - double interceptY = interceptPoint.Y(); + if(trees.count(detector) == 0) { + LOG(TRACE) << "Skipping detector " << detector->getName() << " as it has 0 clusters."; + continue; + } - for(size_t ne = 0; ne < neighbors.size(); ne++) { - auto newCluster = neighbors[ne].get(); + // Get all neighbors within the timing cut + LOG(DEBUG) << "Searching for neighboring cluster on device " << detector->getName(); + LOG(DEBUG) << "- reference time is " << Units::display(seedTrack->timestamp(), {"ns", "us", "s"}); + Cluster* closestCluster = nullptr; - // Calculate the distance to the previous plane's cluster/intercept - double distanceX = interceptX - newCluster->local().x(); - double distanceY = interceptY - newCluster->local().y(); - double distance = sqrt(distanceX * distanceX + distanceY * distanceY); + // Use spatial cut only as initial value (check if cluster is ellipse defined by cuts is done below): + double closestClusterDistance = sqrt(spatial_cuts_[detector].x() * spatial_cuts_[detector].x() + + spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); - // Check if newCluster lies within ellipse defined by spatial cuts around intercept, - // following this example: - // https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/ - // - // ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse, - // > 1: outside, - // < 1: inside - // Continue if outside of ellipse: + double timeCut = std::max(time_cut_ref_track, time_cuts_[detector]); + LOG(DEBUG) << "Using timing cut of " << Units::display(timeCut, {"ns", "us", "s"}); - double norm = (distanceX * distanceX) / (spatial_cuts_[detector].x() * spatial_cuts_[detector].x()) + - (distanceY * distanceY) / (spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + auto neighbors = trees[detector].getAllElementsInTimeWindow(seedTrack->timestamp(), timeCut); - if(norm > 1) { - LOG(DEBUG) << "Cluster outside the cuts. Normalized distance: " << norm; - continue; - } + LOG(DEBUG) << "- found " << neighbors.size() << " neighbors within the correct time window"; - // If this is the closest keep it for now - if(distance < closestClusterDistance) { - closestClusterDistance = distance; - closestCluster = newCluster; - } - } + // Now look for the spatially closest cluster on the next plane + seedTrack->fit(); - if(closestCluster == nullptr) { - LOG(DEBUG) << "No cluster within spatial cut"; - continue; + // check whether the seed track has a radius within the limits of the expectation from momentum and magnetic + // field. Reduces computing time + if(magnetic_field_ != ROOT::Math::XYZVector()) { + auto fittedRadius = seedTrack->getRadius(); + if(std::fabs(fittedRadius - expected_radius_) > radius_tol_) { + LOG(DEBUG) << "Radius of seed track outside of tolerance. Looking at next seed. Expected radius: " + << Units::display(expected_radius_, {"mm", "m", "km", "um"}) + << ", tolerance: " << Units::display(radius_tol_, {"mm", "m", "km", "um"}) + << ", fitted radius: " << Units::display(fittedRadius, {"mm", "m", "km", "um"}); + break; } + } - // Add the cluster to the track - refTrack.addCluster(closestCluster); - track->addCluster(closestCluster); - averageTimestamp = calculate_average_timestamp(&refTrack); - refTrack.setTimestamp(averageTimestamp); - track->setTimestamp(averageTimestamp); + PositionVector3D> interceptPoint = detector->getLocalIntercept(seedTrack.get()); + double interceptX = interceptPoint.X(); + double interceptY = interceptPoint.Y(); - LOG(DEBUG) << "- added cluster to track"; - } + for(size_t ne = 0; ne < neighbors.size(); ne++) { + auto newCluster = neighbors[ne].get(); - // check if track has required detector(s): - auto foundRequiredDetector = [this](Track* t) { - for(auto& requireDet : require_detectors_) { - if(!requireDet.empty() && !t->hasDetector(requireDet)) { - LOG(DEBUG) << "No cluster from required detector " << requireDet << " on the track."; - return false; - } + // Calculate the distance to the previous plane's cluster/intercept + double distanceX = interceptX - newCluster->local().x(); + double distanceY = interceptY - newCluster->local().y(); + double distance = sqrt(distanceX * distanceX + distanceY * distanceY); + + // Check if newCluster lies within ellipse defined by spatial cuts around intercept, + // following this example: + // https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/ + // + // ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse, + // > 1: outside, + // < 1: inside + // Continue if outside of ellipse: + + double norm = (distanceX * distanceX) / (spatial_cuts_[detector].x() * spatial_cuts_[detector].x()) + + (distanceY * distanceY) / (spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + + if(norm > 1) { + LOG(DEBUG) << "Cluster outside the cuts. Normalized distance: " << norm; + continue; + } + + // If this is the closest keep it for now + if(distance < closestClusterDistance) { + closestClusterDistance = distance; + closestCluster = newCluster; } - return true; - }; - if(!foundRequiredDetector(track.get())) { - continue; } - // Now should have a track with one cluster from each plane - if(track->getNClusters() < min_hits_on_track_) { - LOG(DEBUG) << "Not enough clusters on the track, found " << track->getNClusters() << " but " - << min_hits_on_track_ << " required."; + if(closestCluster == nullptr) { + LOG(DEBUG) << "No cluster within spatial cut"; continue; } - // Fit the track - track->fit(); + // Add the cluster to the track + seedTrack->addCluster(closestCluster); + track->addCluster(closestCluster); + averageTimestamp = calculate_average_timestamp(seedTrack.get()); + seedTrack->setTimestamp(averageTimestamp); + track->setTimestamp(averageTimestamp); + + LOG(DEBUG) << "- added cluster to track"; + } - if(reject_by_ROI_ && track->isFitted()) { - // check if the track is within ROI for all detectors - auto ds = get_regular_detectors(!exclude_DUT_); - auto out_of_roi = - std::find_if(ds.begin(), ds.end(), [track](const auto& d) { return !d->isWithinROI(track.get()); }); - if(out_of_roi != ds.end()) { - LOG(DEBUG) << "Rejecting track outside of ROI of detector " << out_of_roi->get()->getName(); - continue; + // check if track has required detector(s): + auto foundRequiredDetector = [this](Track* t) { + for(auto& requireDet : require_detectors_) { + if(!requireDet.empty() && !t->hasDetector(requireDet)) { + LOG(DEBUG) << "No cluster from required detector " << requireDet << " on the track."; + return false; } } - // save the track - if(track->isFitted()) { - tracks.push_back(track); - } else { - LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting"; + return true; + }; + if(!foundRequiredDetector(track.get())) { + continue; + } + + // Now should have a track with one cluster from each plane + if(track->getNClusters() < min_hits_on_track_) { + LOG(DEBUG) << "Not enough clusters on the track, found " << track->getNClusters() << " but " + << min_hits_on_track_ << " required."; + continue; + } + + // Fit the track + track->fit(); + + if(reject_by_ROI_ && track->isFitted()) { + // check if the track is within ROI for all detectors + auto ds = get_regular_detectors(!exclude_DUT_); + auto out_of_roi = + std::find_if(ds.begin(), ds.end(), [track](const auto& d) { return !d->isWithinROI(track.get()); }); + if(out_of_roi != ds.end()) { + LOG(DEBUG) << "Rejecting track outside of ROI of detector " << out_of_roi->get()->getName(); continue; } - - if(timestamp_from_.empty()) { - // Improve the track timestamp by taking the average of all planes - auto timestamp = calculate_average_timestamp(track.get()); - track->setTimestamp(timestamp); - LOG(DEBUG) << "Using average cluster timestamp of " << Units::display(timestamp, "us") - << " as track timestamp."; - } else { - // use timestamp of required detector: - double track_timestamp = track->getClusterFromDetector(timestamp_from_)->timestamp(); - LOG(DEBUG) << "Using timestamp of detector " << timestamp_from_ - << " as track timestamp: " << Units::display(track_timestamp, "us"); - track->setTimestamp(track_timestamp); + } + // check whether the track has a radius within the limits of the expectation from momentum and magnetic field. + // Reduces computing time + if(magnetic_field_ != ROOT::Math::XYZVector() && track->isFitted() && track_model_ == "helix") { + auto fittedRadius = track->getRadius(); + if(std::fabs(fittedRadius - expected_radius_) > radius_tol_) { + LOG(DEBUG) << "Rejecting track with radius outside of tolerance. Expected radius: " + << Units::display(expected_radius_, {"mm", "m", "km", "um"}) + << ", tolerance: " << Units::display(radius_tol_, {"mm", "m", "km", "um"}) + << ", fitted radius: " << Units::display(fittedRadius, {"mm", "m", "km", "um"}); + continue; } } + // save the track + if(track->isFitted()) { + tracks.push_back(track); + } else { + LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting"; + continue; + } + + if(timestamp_from_.empty()) { + // Improve the track timestamp by taking the average of all planes + auto timestamp = calculate_average_timestamp(track.get()); + track->setTimestamp(timestamp); + LOG(DEBUG) << "Using average cluster timestamp of " << Units::display(timestamp, "us") << " as track timestamp."; + } else { + // use timestamp of required detector: + double track_timestamp = track->getClusterFromDetector(timestamp_from_)->timestamp(); + LOG(DEBUG) << "Using timestamp of detector " << timestamp_from_ + << " as track timestamp: " << Units::display(track_timestamp, "us"); + track->setTimestamp(track_timestamp); + } } auto duplicated_hit = [this](const Track* a, const Track* b) { @@ -555,6 +682,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { trackAngleX->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).X())); trackAngleY->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).Y())); } + trackRadius->Fill(track->getRadius() / 1000); // Make residuals auto trackClusters = track->getClusters(); for(auto& trackCluster : trackClusters) { diff --git a/src/modules/Tracking4D/Tracking4D.h b/src/modules/Tracking4D/Tracking4D.h index 68994b158754673b91d951fe9702dd489af28c0c..b2417ed62d26b3dd1b3542a32fb0b3d4dcb32fef 100644 --- a/src/modules/Tracking4D/Tracking4D.h +++ b/src/modules/Tracking4D/Tracking4D.h @@ -11,6 +11,7 @@ #ifndef TRACKING4D_H #define TRACKING4D_H 1 +#include #include #include #include @@ -45,6 +46,7 @@ namespace corryvreckan { TH1F* tracksPerEvent; TH1F* trackAngleX; TH1F* trackAngleY; + TH1F* trackRadius; TH1F* tracksVsTime; std::map residualsX_local; std::map residualsXwidth1_local; @@ -93,6 +95,11 @@ namespace corryvreckan { std::map, XYVector> spatial_cuts_; std::string timestamp_from_; std::string track_model_; + ROOT::Math::XYZVector magnetic_field_; + double expected_radius_; + double radius_tol_; + double radius_tol_min_; + double radius_tol_rel_; // Function to calculate the weighted average timestamp from the clusters of a track double calculate_average_timestamp(const Track* track); diff --git a/src/modules/TrackingMultiplet/TrackingMultiplet.cpp b/src/modules/TrackingMultiplet/TrackingMultiplet.cpp index 5fda24f476ff76e1397c0cb87ae040c629304ee5..1fcbe31a441cebee8cb89fe57c07a16a8bafa142 100644 --- a/src/modules/TrackingMultiplet/TrackingMultiplet.cpp +++ b/src/modules/TrackingMultiplet/TrackingMultiplet.cpp @@ -294,7 +294,7 @@ TrackVector TrackingMultiplet::find_multiplet_tracklets(const streams& stream, for(auto& det : cluster_trees) { auto detector = det.first.get(); trackletCandidate->registerPlane( - detector->getName(), detector->displacement().x(), detector->materialBudget(), detector->toLocal()); + detector->getName(), detector->origin(), detector->materialBudget(), detector->toLocal()); } trackletCandidate->addCluster(clusterFirst.get()); trackletCandidate->addCluster(clusterLast.get()); diff --git a/src/objects/CMakeLists.txt b/src/objects/CMakeLists.txt index 7951e664cb8c1b457e11dbcb4e73835f20fc53a3..7ae87ecdea516b21aa0ec28246f666c8a696a4a1 100644 --- a/src/objects/CMakeLists.txt +++ b/src/objects/CMakeLists.txt @@ -27,6 +27,7 @@ ROOT_GENERATE_DICTIONARY(CorryvreckanObjectsDictionary ${CMAKE_CURRENT_SOURCE_DIR}/MCParticle.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Event.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Multiplet.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/HelixTrack.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Waveform.hpp LINKDEF ${CMAKE_CURRENT_SOURCE_DIR}/Linkdef.h @@ -53,6 +54,7 @@ ADD_LIBRARY(CorryvreckanObjects SHARED GblTrack.cpp Event.cpp Multiplet.cpp + HelixTrack.cpp ${CMAKE_CURRENT_BINARY_DIR}/CorryvreckanObjectsDictionary.cxx ) diff --git a/src/objects/GblTrack.cpp b/src/objects/GblTrack.cpp index db0076f4c86277c68fdba0cdab2f5cf1f187d09d..f7212c8ff649b0a26ca169f36c5d1e96ea4161dc 100644 --- a/src/objects/GblTrack.cpp +++ b/src/objects/GblTrack.cpp @@ -56,6 +56,10 @@ void GblTrack::add_plane(std::vector::iterator& plane, double total_material) { // lambda to add plane (not the first one) and air scatterers auto globalTangent = Vector4d(0, 0, 1, 0); + auto cosLamb = std::sqrt(globalTangent.x() * globalTangent.x() + globalTangent.y() * globalTangent.y()) / + std::sqrt(globalTangent.x() * globalTangent.x() + globalTangent.y() * globalTangent.y() + + globalTangent.z() * globalTangent.z()); + auto lamb = std::acos(cosLamb); // extract the rotation from an ROOT::Math::Transfrom3D, store it in 4x4 matrix to match proteus format auto getRotation = [](Transform3D in) { @@ -69,13 +73,13 @@ void GblTrack::add_plane(std::vector::iterator& plane, auto tmp_local = plane->getToLocal() * globalTrackPos; Vector4d localPosTrack = Vector4d{tmp_local.x(), tmp_local.y(), tmp_local.z(), 1}; Vector4d localTangent = getRotation(plane->getToLocal()) * globalTangent; - double dist = localPosTrack[2]; + double dist = -localPosTrack[2]; // Minus sign to get a positive value LOG(TRACE) << "Rotation: " << getRotation(plane->getToLocal()); LOG(TRACE) << "Local tan before normalization: " << localTangent; LOG(TRACE) << "Distance: " << dist; localTangent /= localTangent.z(); - localPosTrack -= dist * localTangent; + localPosTrack += dist * localTangent; // add the local track pos for future reference - e.g. dut position: local_track_points_[plane->getName()] = ROOT::Math::XYPoint(localPosTrack(0), localPosTrack(1)); @@ -83,11 +87,15 @@ void GblTrack::add_plane(std::vector::iterator& plane, Matrix4d toTarget = getRotation(plane->getToLocal()) * getRotation(prevToGlobal); // lambda Jacobian from one scatter to the next - auto jac = [](const Vector4d& tangent, const Matrix4d& target, double distance) { + auto jac = [this](const Vector4d& tangent, const Matrix4d& target, double distance, double lambda) { + auto magneticFieldFactor = magneticField_ * 0.3; // Magnetic field factor: B*c/MeV*mm + auto BfacX = (-std::cos(lambda) * magneticFieldFactor.z() - std::sin(lambda) * magneticFieldFactor.y()) * distance; + auto BfacY = magneticFieldFactor.x() * distance; + LOG(TRACE) << "Magnetic field factor: " << magneticFieldFactor; Matrix R; R.col(0) = target.col(0); R.col(1) = target.col(1); - R.col(2) = target.col(3); + R.col(2) = target.col(0) * BfacX + target.col(1) * BfacY; Vector4d S = target * tangent * (1 / tangent[2]); Matrix F = Matrix::Zero(); @@ -99,11 +107,13 @@ void GblTrack::add_plane(std::vector::iterator& plane, F(2, 2) = -S[3] / S[2]; Matrix jaco; - jaco << F * R, (-distance / S[2]) * F * R, Matrix3d::Zero(), (1 / S[2]) * F * R; - jaco(5, 5) = 1; // a future time component + jaco << F * R, (distance / S[2]) * F * R, Matrix3d::Zero(), (1 / S[2]) * F * R; + jaco(0, 5) /= 2; + jaco(1, 5) /= 2; + jaco(5, 5) = 1; return jaco; }; - auto myjac = jac(prevTan, toTarget, dist); + auto myjac = jac(prevTan, toTarget, dist, lamb); // Layout if volume scattering active // | | | | @@ -142,16 +152,16 @@ void GblTrack::add_plane(std::vector::iterator& plane, myjac(0, 0) = 1; // Adding volume scattering if requested } else if(use_volume_scatter_) { - myjac = jac(prevTan, toTarget, frac1 * dist); + myjac = jac(prevTan, toTarget, frac1 * dist, lamb); GblPoint pVolume(toGbl * myjac * toProt); addScattertoGblPoint(pVolume, fabs(dist) / 2. / scattering_length_volume_); gblpoints_.push_back(pVolume); // We have already rotated to the next local coordinate system - myjac = jac(prevTan, Matrix4d::Identity(), frac2 * dist); + myjac = jac(prevTan, Matrix4d::Identity(), frac2 * dist, lamb); GblPoint pVolume2(toGbl * myjac * toProt); addScattertoGblPoint(pVolume2, fabs(dist) / 2. / scattering_length_volume_); gblpoints_.push_back(pVolume2); - myjac = jac(prevTan, Matrix4d::Identity(), frac1 * dist); + myjac = jac(prevTan, Matrix4d::Identity(), frac1 * dist, lamb); } auto transformedJac = toGbl * myjac * toProt; GblPoint point(transformedJac); @@ -211,7 +221,7 @@ void GblTrack::prepare_gblpoints() { for(auto& l : planes_) { total_material += l.getMaterialBudget(); if(clusters.count(l.getName()) == 1) { - l.setPosition(clusters.at(l.getName())->global().z()); + l.setPosition(clusters.at(l.getName())->global()); l.setCluster(clusters.at(l.getName())); } } @@ -221,7 +231,7 @@ void GblTrack::prepare_gblpoints() { // add volume scattering length - for now simply the distance between first and last plane if(use_volume_scatter_) { - total_material += (planes_.back().getPosition() - planes_.front().getPosition()) / scattering_length_volume_; + total_material += (planes_.back().getPosition().z() - planes_.front().getPosition().z()) / scattering_length_volume_; } // Prepare transformations for first plane: @@ -267,7 +277,7 @@ void GblTrack::fit() { prepare_gblpoints(); // perform fit - GblTrajectory traj(gblpoints_, false); // false = no magnetic field + GblTrajectory traj(gblpoints_, magneticField_.R() != 0); // true = magnetic field, false = no magnetic field double lostWeight = 0; int ndf = 0; @@ -337,7 +347,7 @@ void GblTrack::fit() { LOG(TRACE) << "Results for detector " << name << std::endl << "Fitted residual local:\t" << residual_local_.at(name) << std::endl << "Seed residual:\t" << initital_residual_.at(name) << std::endl - << "Ditted residual global:\t" << ROOT::Math::XYPoint(clusterPos - corPos); + << "Fitted residual global:\t" << ROOT::Math::XYPoint(clusterPos - corPos); } LOG(DEBUG) << "Plane: " << name << ": residual " << residual_local_[name] << ", kink: " << kink_[name]; } @@ -355,7 +365,7 @@ ROOT::Math::XYZPoint GblTrack::getIntercept(double z) const { } for(const auto& l : planes_) { - if(l.getPosition() >= z) { + if(l.getPosition().z() >= z) { found = true; break; } @@ -408,7 +418,7 @@ XYZPoint GblTrack::get_position_outside_telescope(double z) const { auto last_plane = planes_.end(); last_plane--; // No direct iterator for last_plane element // check if z is up or downstream - bool upstream = (z < first_plane->getPosition()); + bool upstream = (z < first_plane->getPosition().z()); auto outerPlane = (upstream ? first_plane->getName() : last_plane->getName()); // inner neighbour of plane - simply adjust the iterators @@ -447,7 +457,8 @@ ROOT::Math::XYZVector GblTrack::getDirection(const std::string& detectorID) cons } XYZVector GblTrack::getDirection(const double& z) const { - auto planeUpstream = std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition() > z; }); + auto planeUpstream = + std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition().z() > z; }); if(planeUpstream != planes_.end()) { return getDirection(planeUpstream->getName()); } else { diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp new file mode 100644 index 0000000000000000000000000000000000000000..69ebc8b592a124b97880f23e01dabbd36b007c1a --- /dev/null +++ b/src/objects/HelixTrack.cpp @@ -0,0 +1,317 @@ +/** + * @file + * @brief Implementation of StraightLine track object + * + * @copyright Copyright (c) 2017-2022 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#include "Eigen/Dense" +#include "StraightLineTrack.hpp" +#include "Track.hpp" +#include "core/utils/log.h" +#include "core/utils/unit.h" +#include "exceptions.h" + +using namespace corryvreckan; +using namespace Eigen; +using namespace ROOT::Math; + +ROOT::Math::XYPoint HelixTrack::getKinkAt(const std::string&) const { + return ROOT::Math::XYPoint(0, 0); +} + +double HelixTrack::getRadius() const { + if(!isFitted_) { + throw TrackError(typeid(HelixTrack), " has no defined radius before fitting"); + } + return fabs(1 / curvature_); +} + +ROOT::Math::XYZPoint HelixTrack::getState(const std::string& detectorID) const { + const Plane* plane = get_plane(detectorID); + XYZVector planeU, planeV, planeN; + plane->getToGlobal().Rotation().GetComponents(planeU, planeV, planeN); + return getIntercept(plane->getPosition(), planeN); +} + +ROOT::Math::XYZVector HelixTrack::getDirection(const std::string& detectorID) const { + const Plane* plane = get_plane(detectorID); + XYZVector planeU, planeV, planeN; + plane->getToGlobal().Rotation().GetComponents(planeU, planeV, planeN); + return getDirection(plane->getPosition(), planeN); +} + +ROOT::Math::XYZVector HelixTrack::getDirection(const double& z) const { + XYZPoint origin = XYZPoint(0, 0, z); + XYZVector normal = XYZVector(0, 0, 1); + return getDirection(origin, normal); +} + +void HelixTrack::calculate_transformations() { + // compute transformation to helix coordinate system with the direction of the magnetic field as the z-direction + XYZVector globalZ = XYZVector(0, 0, 1); + XYZVector helixZ = magneticField_.unit(); + XYZVector helixY = helixZ.Theta() != 0 ? helixZ.Cross(globalZ).unit() : XYZVector(0, 1, 0); + XYZVector helixX = helixY.Cross(helixZ).unit(); + + toGlobal_ = Rotation3D(helixX, helixY, helixZ); + toHelix_ = toGlobal_.Inverse(); +} + +void HelixTrack::calculate_chi2() { + + // Each hit provides two measurements and there are five fitting paramters. + ndof_ = 2 * track_clusters_.size() - 5; + + chi2_ = 0.; + chi2ndof_ = 0.; + + // Loop over all clusters + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + // Get the distance and the error + auto intercept = getState(cluster->detectorID()); + auto interceptLocal = get_plane(cluster->detectorID())->getToLocal() * intercept; + auto residual = cluster->local() - interceptLocal; + LOG(TRACE) << "Residual on detector " << cluster->detectorID() << ": " << residual * 1000 << " um"; + + double ex2 = cluster->errorX() * cluster->errorX(); + double ey2 = cluster->errorY() * cluster->errorY(); + chi2_ += ((residual.x() * residual.x() / ex2) + (residual.y() * residual.y() / ey2)); + } + + // Store also the chi2/degrees of freedom + chi2ndof_ = (ndof_ <= 0) ? -1 : (chi2_ / static_cast(ndof_)); +} + +void HelixTrack::calculate_residuals() { + for(const auto& c : track_clusters_) { + auto* cluster = c.get(); + + auto intercept = getState(cluster->detectorID()); + auto interceptLocal = get_plane(cluster->detectorID())->getToLocal() * intercept; + residual_global_[cluster->detectorID()] = cluster->global() - intercept; + residual_local_[cluster->detectorID()] = cluster->local() - interceptLocal; + } +} + +double HelixTrack::operator()(const double* parameters) { + + // Update the HelixTrack parameters + this->curvature_ = parameters[0]; + this->phi_ = parameters[1]; + this->dca_ = parameters[2]; + this->tanLambda_ = parameters[3]; + this->z0_ = parameters[4]; + + // Calculate the chi2 + this->calculate_chi2(); + + // Return this to minuit + return chi2_; +} + +void HelixTrack::fit() { + + if(magneticField_ == XYZVector()) { + throw TrackError(typeid(HelixTrack), " can only be used with magnetic fields"); + } + + isFitted_ = false; + calculate_transformations(); + + LOG(DEBUG) << "Starting helix fit."; + // Fit points to helix with fast helix fit described in: + // https://www.physi.uni-heidelberg.de/Publications/DiplomaKiehn.pdf + + // 1. Fit circle in transverse plane + // + // mean of positons and radius in transverse plane (x, y, r^2) + Matrix mean = Matrix::Zero(); + // mean of squared positions and radius in transverse plane + // (x^2, x*y, x*r^2 ) + // ( 0 , y^2, y*r^2 ) + // ( 0 , 0 , r^4 ) + Matrix mean2 = Matrix::Zero(); + auto n = track_clusters_.size(); + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + auto position = toHelix_ * cluster->global(); + long double x = position.x(); + long double y = position.y(); + long double r2 = x * x + y * y; + + mean(0) += x; + mean(1) += y; + mean(2) += r2; + mean2(0, 0) += x * x; + mean2(0, 1) += x * y; + mean2(0, 2) += x * r2; + mean2(1, 1) += y * y; + mean2(1, 2) += y * r2; + mean2(2, 2) += r2 * r2; + } + mean /= static_cast(n); + mean2 /= static_cast(n); + LOG(TRACE) << "mean: " << mean; + + Matrix C = Matrix::Zero(); + C(0, 0) = mean2(0, 0) - mean(0) * mean(0); + C(0, 1) = mean2(0, 1) - mean(0) * mean(1); + C(0, 2) = mean2(0, 2) - mean(0) * mean(2); + C(1, 1) = mean2(1, 1) - mean(1) * mean(1); + C(1, 2) = mean2(1, 2) - mean(1) * mean(2); + C(2, 2) = mean2(2, 2) - mean(2) * mean(2); + + LOG(TRACE) << "C: " << C; + + long double q1 = (C(2, 2) * C(0, 1)) - (C(0, 2) * C(1, 2)); + long double q2 = (C(2, 2) * (C(0, 0) - C(1, 1))) - (C(0, 2) * C(0, 2)) + (C(1, 2) * C(1, 2)); + LOG(TRACE) << "q1: " << q1 << ", q2: " << q2; + + phi_ = static_cast(0.5 * atan(2. * q1 / q2)); + long double kappa = (sin(phi_) * C(0, 2) - cos(phi_) * C(1, 2)) / C(2, 2); + long double delta = -kappa * mean(2) + sin(phi_) * mean(0) - cos(phi_) * mean(1); + + curvature_ = static_cast(2. * kappa / sqrt(1. - 4. * delta * kappa)); + dca_ = static_cast(2. * delta / (1. + sqrt(1. - 4. * delta * kappa))); + + LOG(DEBUG) << "Circle fit finished. Radius: " << Units::display(1 / curvature_, {"m", "km"}) + << ", DCA: " << Units::display(dca_, {"mm", "um"}) << ", Phi: " << Units::display(phi_, "deg"); + + // 2. Calculate arc length in transverse plane + // 3. Fit in longitudinal (a,z)-plane + Matrix2d mat(Matrix2d::Zero()); + Vector2d vec(Vector2d::Zero()); + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + auto position = toHelix_ * cluster->global(); + double a = getArcLength(position.x(), position.y()); + + double z = position.z(); + + vec += Vector2d(z, a * z); + Matrix2d err; + err << 1., a, a, a * a; + mat += err; + } + + // Check for singularities. + if(fabs(mat.determinant()) < std::numeric_limits::epsilon()) { + throw TrackFitError(typeid(this), "Martix inversion in straight line fit failed"); + } + + // Get the longitudinal parameters + Eigen::Vector2d res = mat.inverse() * vec; + + z0_ = res[0]; + tanLambda_ = res[1]; + LOG(DEBUG) << "Longitudinal fit finished. Z0: " << Units::display(z0_, {"mm", "cm"}) << ", tan(lambda): " << tanLambda_; + + // Calculate the chi2 + this->calculate_chi2(); + this->calculate_residuals(); + isFitted_ = true; +} + +ROOT::Math::XYZPoint HelixTrack::getIntercept(double z) const { + XYZPoint origin = XYZPoint(0., 0., z); + XYZVector normal = XYZVector(0., 0., 1); + return getIntercept(origin, normal); +} + +double HelixTrack::getArcLength(double x, double y) const { + // Transform to cooridnate system with origin in center of track circle + Vector2d posOrigin; + posOrigin << x, y; + Matrix2d toCenter; + toCenter << sin(phi_), -cos(phi_), cos(phi_), sin(phi_); + Vector2d posCenter = toCenter * posOrigin; + posCenter[0] -= (dca_ + 1. / curvature_); + + double a = abs(1. / curvature_) * + (curvature_ <= 0. ? std::atan2(posCenter[1], posCenter[0]) : std::atan2(posCenter[1], -posCenter[0])); + return a; +} + +ROOT::Math::XYZPoint HelixTrack::getIntercept(XYZPoint origin, XYZVector normal) const { + // Rotate plane origin and normal to helix coordinates + auto normalHelix = toHelix_ * normal; + auto originHelix = toHelix_ * origin; + + double cosLambda = 1. / sqrt(1. + tanLambda_ * tanLambda_); + double sinLambda = sqrt(1. - cosLambda * cosLambda); + XYZVector distance, trackDirection; + XYZPoint position; + // initial estimate of arc length + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(dPhi - phi_); + const double sinPhi = sin(dPhi - phi_); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) + cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; + } + } + + return toGlobal_ * position; +} + +ROOT::Math::XYZVector HelixTrack::getDirection(XYZPoint origin, XYZVector normal) const { + // Rotate plane normal vector to helix coordinates + auto normalHelix = toHelix_ * normal; + auto originHelix = toHelix_ * origin; + + double cosLambda = 1. / sqrt(1. + tanLambda_ * tanLambda_); + double sinLambda = sqrt(1. - cosLambda * cosLambda); + XYZVector distance, trackDirection; + XYZPoint position; + // initial estimate of arc length + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(phi_ + dPhi); + const double sinPhi = sin(phi_ + dPhi); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) - cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; + } + } + + return toGlobal_ * trackDirection; +} + +void HelixTrack::print(std::ostream& out) const { + out << "HelixTrack " << this->curvature_ << ", " << this->phi_ << ", " << this->dca_ << ", " << this->tanLambda_ << ", " + << this->z0_ << ", " << this->chi2_ << ", " << this->ndof_ << ", " << this->chi2ndof_ << ", " << this->timestamp(); +} diff --git a/src/objects/HelixTrack.hpp b/src/objects/HelixTrack.hpp new file mode 100644 index 0000000000000000000000000000000000000000..06c40bba46235a2a6ef1ed78b6c35559c6dbdb34 --- /dev/null +++ b/src/objects/HelixTrack.hpp @@ -0,0 +1,130 @@ +/** + * @file + * @brief Definition of Helix track object + * + * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#ifndef CORRYVRECKAN_HELIXTRACK_H +#define CORRYVRECKAN_HELIXTRACK_H 1 + +#include "Track.hpp" + +namespace corryvreckan { + /** + * @ingroup Objects + * @brief HelixTrack object + * + * This class is a simple track class which knows how to fit itself. It is derived from Track + */ + + class HelixTrack : public Track { + + public: + // Minimisation operator used by Minuit. Minuit passes the current iteration of the parameters and checks if the chi2 + // is better or worse + double operator()(const double* parameters); + + void print(std::ostream& out) const override; + + /** + * @brief The fiting routine + */ + void fit() override; + + /** + * @brief Get the track position for a certain z position + * @param z Global z position + * @return ROOT::Math::XYZPoint at z position + */ + ROOT::Math::XYZPoint getIntercept(double z) const override; + + /** + * @brief Get the track position on a plane specified by a point and the normal vector in global coordinates + * @param origin Point on the plane (in global coordinates) + * @param normal Normal vector to the plane (in global coordinates) + * @return ROOT::Math::XYZPoint at intersection with plane + */ + ROOT::Math::XYZPoint getIntercept(ROOT::Math::XYZPoint origin, ROOT::Math::XYZVector normal) const; + + /** + * @brief Get the track state at a detector + * @param detectorID Name of detector + * @return ROOT::Math::XYZPoint state at detetcor layer + */ + ROOT::Math::XYZPoint getState(const std::string& detectorID) const override; + + /** + * @brief Get the track direction at a detector + * @param detectorID Name of detector + * @return ROOT::Math::XYZPoint direction at detetcor layer + */ + ROOT::Math::XYZVector getDirection(const std::string& detectorID) const override; + + /** + * @brief Get the track direction at any given z position + * @param z Global z position + * @return ROOT::Math::XYZPoint direction at z + */ + ROOT::Math::XYZVector getDirection(const double& z) const override; + + /** + * @brief Get the track direction on a plane specified by a point and the normal vector in global coordinates + * @param origin Point on the plane (in global coordinates) + * @param normal Normal vector to the plane (in global coordinates) + * @return ROOT::Math::XYZVector at intersection with plane + */ + ROOT::Math::XYZVector getDirection(ROOT::Math::XYZPoint origin, ROOT::Math::XYZVector normal) const; + + /** + * @brief This track model does not support kinks, it therefore is always zero + * @param detectorID Detector ID at which the kink should be evaluated + * @return Kink at given detector + */ + ROOT::Math::XYPoint getKinkAt(const std::string& detectorID) const override; + + /** + * @brief Get the radius in the transverse plane at a given detector layer. + * @param detectorID Detector ID at which the kink should be evaluated + * @return radius in transverse plane at a given detector + */ + double getRadius() const override; + + /** + * @brief Get the arc length in the transverse plane of the helix + * @param x x-coordinate in transverse plane + * @param y y-coordinate in transverse plane + * @return double arc length + */ + double getArcLength(double x, double y) const; + + void setVolumeScatter(double) override{}; + + private: + /** + * @brief calculate the chi2 of the linear regression + */ + void calculate_chi2(); + + void calculate_residuals(); + + void calculate_transformations(); + + // Member variables (Karimaeki parametrization) + double curvature_; + double phi_; + double dca_; + double tanLambda_; + double z0_; + ROOT::Math::Rotation3D toHelix_; + ROOT::Math::Rotation3D toGlobal_; + + // ROOT I/O class definition - update version number when you change this class! + ClassDefOverride(HelixTrack, 1) + }; +} // namespace corryvreckan + +#endif // CORRYVRECKAN_HELIXTRACK_H diff --git a/src/objects/Linkdef.h b/src/objects/Linkdef.h index 6aa860322180de75ccc7b822027d50502b9d602e..020d13967c16e494cb3b1f7789a18b95fabe4ec5 100644 --- a/src/objects/Linkdef.h +++ b/src/objects/Linkdef.h @@ -23,6 +23,7 @@ #pragma link C++ class corryvreckan::StraightLineTrack + ; #pragma link C++ class corryvreckan::GblTrack + ; #pragma link C++ class corryvreckan::Multiplet + ; +#pragma link C++ class corryvreckan::HelixTrack + ; #pragma link C++ class corryvreckan::MCParticle + ; #pragma link C++ class corryvreckan::Event + ; #pragma link C++ class corryvreckan::Track::Plane + ; diff --git a/src/objects/Track.cpp b/src/objects/Track.cpp index 7aeb68a045a6c9d89bd96cd5915bb5e8bb67dcec..2e618c04b6abe39fd970ea008af343c341cc7111 100644 --- a/src/objects/Track.cpp +++ b/src/objects/Track.cpp @@ -15,11 +15,11 @@ using namespace corryvreckan; -Track::Plane::Plane(std::string name, double z, double x_x0, Transform3D to_local) - : z_(z), x_x0_(x_x0), name_(std::move(name)), to_local_(to_local) {} +Track::Plane::Plane(std::string name, ROOT::Math::XYZPoint origin, double x_x0, Transform3D to_local) + : origin_(origin), x_x0_(x_x0), name_(std::move(name)), to_local_(to_local) {} -double Track::Plane::getPosition() const { - return z_; +ROOT::Math::XYZPoint Track::Plane::getPosition() const { + return origin_; } double Track::Plane::getMaterialBudget() const { @@ -51,7 +51,7 @@ Transform3D Track::Plane::getToGlobal() const { } bool Track::Plane::operator<(const Plane& pl) const { - return z_ < pl.z_; + return origin_.z() < pl.origin_.z(); } void Track::Plane::setCluster(const Cluster* cluster) { @@ -59,7 +59,7 @@ void Track::Plane::setCluster(const Cluster* cluster) { } void Track::Plane::print(std::ostream& os) const { - os << "Plane at " << z_ << " with rad. length " << x_x0_ << ", name " << name_ << " and"; + os << "Plane at " << origin_ << " with rad. length " << x_x0_ << ", name " << name_ << " and"; if(hasCluster()) { os << "cluster with global pos: " << getCluster()->global(); } else { @@ -126,6 +126,10 @@ void Track::setParticleMomentum(double p) { momentum_ = p; } +void Track::setMagneticField(ROOT::Math::XYZVector B) { + magneticField_ = B; +} + double Track::getChi2() const { if(!isFitted_) { throw RequestParameterBeforeFitError(this, "chi2"); @@ -227,8 +231,12 @@ double Track::getMaterialBudget(const std::string& detectorID) const { return budget; } -void Track::registerPlane(const std::string& name, double z, double x0, Transform3D g2l) { - Plane p(name, z, x0, g2l); +double Track::getRadius() const { + return 0.; +} + +void Track::registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, double x0, Transform3D g2l) { + Plane p(name, origin, x0, g2l); auto pl = std::find_if(planes_.begin(), planes_.end(), [&p](const Plane& plane) { return plane.getName() == p.getName(); }); if(pl == planes_.end()) { @@ -238,7 +246,7 @@ void Track::registerPlane(const std::string& name, double z, double x0, Transfor } } -Track::Plane* Track::get_plane(std::string detetorID) { +const Track::Plane* Track::get_plane(const std::string detetorID) const { auto plane = std::find_if(planes_.begin(), planes_.end(), [&detetorID](Plane const& p) { return p.getName() == detetorID; }); if(plane == planes_.end()) { @@ -252,6 +260,8 @@ std::shared_ptr corryvreckan::Track::Factory(const std::string& trackMode return std::make_shared(); } else if(trackModel == "gbl") { return std::make_shared(); + } else if(trackModel == "helix") { + return std::make_shared(); } else { throw UnknownTrackModel(typeid(Track), trackModel); } diff --git a/src/objects/Track.hpp b/src/objects/Track.hpp index ae669c8225e5158b9a3846e09207cfca047a2ac0..8e738c2a1b76e7d7a571af3d24f7cbe5318c493b 100644 --- a/src/objects/Track.hpp +++ b/src/objects/Track.hpp @@ -104,6 +104,12 @@ namespace corryvreckan { */ void setParticleMomentum(double momentum); + /** + * @brief Set the magnetic field + * @param magneticField + */ + void setMagneticField(ROOT::Math::XYZVector magneticField); + /** * @brief Get the chi2 of the track fit * @return chi2 @@ -223,6 +229,13 @@ namespace corryvreckan { */ virtual ROOT::Math::XYPoint getKinkAt(const std::string& detectorID) const = 0; + /** + * @brief Get the radius in the transverse plane at a given detector layer. + * @param detectorID Detector ID at which the kink should be evaluated + * @return radius in transverse plane at a given detector + */ + virtual double getRadius() const; + /** * @brief Get the materialBudget of a detector layer * @param detectorID @@ -232,12 +245,12 @@ namespace corryvreckan { virtual void setVolumeScatter(double length) = 0; - void registerPlane(const std::string& name, double z, double x0, Transform3D g2l); + void registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, double x0, Transform3D g2l); class Plane { public: Plane() = default; - Plane(std::string name, double z, double x_x0, Transform3D to_local); + Plane(std::string name, ROOT::Math::XYZPoint origin, double x_x0, Transform3D to_local); /** * @brief Required virtual destructor */ @@ -259,7 +272,7 @@ namespace corryvreckan { Plane& operator=(Plane&& rhs) = default; /// @} - double getPosition() const; + ROOT::Math::XYZPoint getPosition() const; double getMaterialBudget() const; bool hasCluster() const; @@ -271,7 +284,7 @@ namespace corryvreckan { // sorting overload bool operator<(const Plane& pl) const; // set elements that might be unknown at construction - void setPosition(double z) { z_ = z; } + void setPosition(ROOT::Math::XYZPoint origin) { origin_ = origin; } void setCluster(const Cluster* cluster); void print(std::ostream& os) const; @@ -282,7 +295,8 @@ namespace corryvreckan { void petrifyHistory(); private: - double z_, x_x0_; + ROOT::Math::XYZPoint origin_; + double x_x0_; std::string name_; PointerWrapper cluster_; Transform3D to_local_; @@ -291,7 +305,7 @@ namespace corryvreckan { void loadHistory() override; void petrifyHistory() override; - Plane* get_plane(std::string detetorID); + const Plane* get_plane(const std::string detetorID) const; std::vector> track_clusters_; std::vector> associated_clusters_; std::map residual_local_; @@ -305,6 +319,7 @@ namespace corryvreckan { double chi2ndof_; bool isFitted_{}; double momentum_{-1}; + ROOT::Math::XYZVector magneticField_; // ROOT I/O class definition - update version number when you change this class! ClassDefOverride(Track, 11) @@ -315,6 +330,7 @@ namespace corryvreckan { // include all tracking methods here to have one header to be include everywhere #include "GblTrack.hpp" +#include "HelixTrack.hpp" #include "Multiplet.hpp" #include "StraightLineTrack.hpp" #endif // CORRYVRECKAN_TRACK_H diff --git a/src/objects/objects.h b/src/objects/objects.h index 47b054f2470be16e585919e6ce6e1c2dcee65fda..d51aa1517931ac9ed045972674c37f4c76543071 100644 --- a/src/objects/objects.h +++ b/src/objects/objects.h @@ -19,5 +19,6 @@ namespace corryvreckan { /** * @brief Tuple containing all objects */ - using OBJECTS = std::tuple; + using OBJECTS = + std::tuple; } // namespace corryvreckan