diff --git a/event/include/Track.h b/event/include/Track.h index af7ab1e0..16e51695 100644 --- a/event/include/Track.h +++ b/event/include/Track.h @@ -291,7 +291,17 @@ class Track : public TObject { double getP() const {return sqrt(px_*px_ + py_*py_ + pz_*pz_);}; double getPt() const {return sqrt(px_*px_ + py_*py_);} + + /** + * Get the track residual of the given layer. + * + * @param layer The SVT layer of interest. + * @return The track residual value of the given layer. + */ + void setTrackResid(const int layer, const double track_resid) { track_residuals_[layer]=track_resid; } + double getTrackResid(const int layer) const { return track_residuals_[layer]; } + /** * Set the lambda kink of the given layer. * @@ -386,7 +396,7 @@ class Track : public TObject { /** hit layers */ std::vector hit_layers_; - /** truth mcp hits */ + /** truth mcp hits */ std::vector> mcp_hits_; /** The distance of closest approach to the reference point. */ @@ -446,6 +456,9 @@ class Track : public TObject { /** The z position track. */ double z_{-999.}; + /** Array used to store track residuals. */ + double track_residuals_[14]; + /** Array used to store the lambda kinks for each of the sensor layers. */ double lambda_kinks_[14]; diff --git a/processors/include/TrackingProcessor.h b/processors/include/TrackingProcessor.h index 9cd43e27..ab1c2b93 100644 --- a/processors/include/TrackingProcessor.h +++ b/processors/include/TrackingProcessor.h @@ -36,7 +36,7 @@ #include "Event.h" #include "RawSvtHit.h" #include "TrackHistos.h" - +#include "MCParticle.h" // Forward declarations class TTree; @@ -117,14 +117,18 @@ class TrackingProcessor : public Processor { int debug_{false}; //!< Debug Level int doResiduals_{0}; //!< do Residuals - std::string trackResDataLcio_{""}; //!< description - TrackHistos* trkResHistos_{nullptr}; //!< description + int doHistograms_{0}; + std::string trackResDataLcio_{""}; //!< description + std::string trackXKinkDataLcio_{""}; + std::string trackZKinkDataLcio_{""}; + + TrackHistos* trkResHistos_{nullptr}; //!< description std::string resCfgFilename_{""}; //!< description std::string resoutname_{""}; //!< description double bfield_{-1.}; //!< magnetic field - int useTrackerHits_{1}; //!< Load hit collections, otherwise get from getSubdetectorHitNumbers + int useTrackerHits_{1};//!< Load hit collections, otherwise get from getSubdetectorHitNumbsrs }; // Tracking Processor diff --git a/processors/src/TrackingProcessor.cxx b/processors/src/TrackingProcessor.cxx index 8543daca..46ffdb9b 100644 --- a/processors/src/TrackingProcessor.cxx +++ b/processors/src/TrackingProcessor.cxx @@ -25,13 +25,15 @@ void TrackingProcessor::configure(const ParameterSet& parameters) { truthTracksCollLcio_ = parameters.getString("truthTrackCollLcio",truthTracksCollLcio_); truthTracksCollRoot_ = parameters.getString("truthTrackCollRoot",truthTracksCollRoot_); trackStateLocation_ = parameters.getString("trackStateLocation",trackStateLocation_); - useTrackerHits_ = parameters.getInteger("useTrackerHits",useTrackerHits_); - bfield_ = parameters.getDouble("bfield",bfield_); - + useTrackerHits_ = parameters.getInteger("useTrackerHits",useTrackerHits_); + bfield_ = parameters.getDouble("bfield",bfield_); //Residual plotting is done in this processor for the moment. doResiduals_ = parameters.getInteger("doResiduals",doResiduals_); - trackResDataLcio_ = parameters.getString("trackResDataLcio",trackResDataLcio_); - resCfgFilename_ = parameters.getString("resPlots",resCfgFilename_); + doHistograms_ = parameters.getInteger("doHistograms",doHistograms_); + trackResDataLcio_ = parameters.getString("trackResDataLcio",trackResDataLcio_); + trackXKinkDataLcio_ = parameters.getString("trackXKinkDataLcio",trackXKinkDataLcio_); + trackZKinkDataLcio_ = parameters.getString("trackZKinkDataLcio",trackZKinkDataLcio_); + resCfgFilename_ = parameters.getString("resPlots",resCfgFilename_); resoutname_ = parameters.getString("resoutname",resoutname_); } @@ -56,7 +58,7 @@ void TrackingProcessor::initialize(TTree* tree) { //Residual plotting - if (doResiduals_) { + if (doResiduals_ and doHistograms_) { trkResHistos_ = new TrackHistos(trkCollLcio_); trkResHistos_->debugMode(debug_); trkResHistos_->loadHistoConfig(resCfgFilename_); @@ -136,8 +138,7 @@ bool TrackingProcessor::process(IEvent* ievent) { std::map > SharedHits; //TODO: can we do better? (innermost) std::map SharedHitsLy0; - std::map SharedHitsLy1; - + std::map SharedHitsLy1; for (int itrack = 0; itrack < tracks->getNumberOfElements();++itrack) { SharedHits[itrack] = {}; SharedHitsLy0[itrack] = false; @@ -178,8 +179,9 @@ bool TrackingProcessor::process(IEvent* ievent) { if (bfield_>0) track->setMomentum(bfield_); - - int nHits = 0; + // Get the collection of hits associated with a LCIO Track + + int nHits = 0; if(!useTrackerHits_){ auto hitPattern = lc_track->getSubdetectorHitNumbers(); @@ -196,67 +198,64 @@ bool TrackingProcessor::process(IEvent* ievent) { track->setTrackerHitCount(nHits); } else{ + // Iterate through the collection of 3D hits (TrackerHit objects) + // associated with a track, find the corresponding hits in the HPS + // event and add references to the track + EVENT::TrackerHitVec lc_tracker_hits = lc_track->getTrackerHits(); + bool rotateHits = true; + int hitType = 0; + if (track->isKalmanTrack()) + hitType=1; //SiClusters + + for (auto lc_tracker_hit : lc_tracker_hits) { + + TrackerHit* tracker_hit = utils::buildTrackerHit(static_cast(lc_tracker_hit),rotateHits,hitType); + + std::vector rawSvthitsOn3d; + utils::addRawInfoTo3dHit(tracker_hit,static_cast(lc_tracker_hit), + raw_svt_hit_fits,&rawSvthitsOn3d,hitType); + + for (auto rhit : rawSvthitsOn3d) + rawhits_.push_back(rhit); - // Get the collection of hits associated with a LCIO Track - EVENT::TrackerHitVec lc_tracker_hits = lc_track->getTrackerHits(); - - // Iterate through the collection of 3D hits (TrackerHit objects) - // associated with a track, find the corresponding hits in the HPS - // event and add references to the track - bool rotateHits = true; - int hitType = 0; - if (track->isKalmanTrack()) - hitType=1; //SiClusters - - for (auto lc_tracker_hit : lc_tracker_hits) { - - TrackerHit* tracker_hit = utils::buildTrackerHit(static_cast(lc_tracker_hit),rotateHits,hitType); - - std::vector rawSvthitsOn3d; - utils::addRawInfoTo3dHit(tracker_hit,static_cast(lc_tracker_hit), - raw_svt_hit_fits,&rawSvthitsOn3d,hitType); - - for (auto rhit : rawSvthitsOn3d) - rawhits_.push_back(rhit); + rawSvthitsOn3d.clear(); - rawSvthitsOn3d.clear(); + if (debug_) + std::cout<getRawHits().GetEntries()<addHit(tracker_hit); + track->addHitLayer(tracker_hit->getLayer()); + hits_.push_back(tracker_hit); + + //Get shared Hits information + for (int jtrack = itrack+1; jtrack < tracks->getNumberOfElements(); ++jtrack) { + + EVENT::Track* j_lc_track = static_cast(tracks->getElementAt(jtrack)); + if (utils::isUsedByTrack(tracker_hit,j_lc_track)) { + //The hit is not already in the shared list + if (std::find(SharedHits[itrack].begin(), SharedHits[itrack].end(),tracker_hit->getID()) == SharedHits[itrack].end()) { + SharedHits[itrack].push_back(tracker_hit->getID()); + if (tracker_hit->getLayer() == 0 ) + SharedHitsLy0[itrack] = true; + if (tracker_hit->getLayer() == 1 ) + SharedHitsLy1[itrack] = true; + } + if (std::find(SharedHits[jtrack].begin(), SharedHits[jtrack].end(),tracker_hit->getID()) == SharedHits[jtrack].end()) { + SharedHits[jtrack].push_back(tracker_hit->getID()); + if (tracker_hit->getLayer() == 0 ) + SharedHitsLy0[jtrack] = true; + if (tracker_hit->getLayer() == 1 ) + SharedHitsLy1[jtrack] = true; + } + } // found shared hit + } // loop on j>i tracks + }//tracker hits + } + track->setNShared(SharedHits[itrack].size()); + track->setSharedLy0(SharedHitsLy0[itrack]); + track->setSharedLy1(SharedHitsLy1[itrack]); + - if (debug_) - std::cout<getRawHits().GetEntries()<addHit(tracker_hit); - track->addHitLayer(tracker_hit->getLayer()); - hits_.push_back(tracker_hit); - - //Get shared Hits information - for (int jtrack = itrack+1; jtrack < tracks->getNumberOfElements(); ++jtrack) { - - EVENT::Track* j_lc_track = static_cast(tracks->getElementAt(jtrack)); - if (utils::isUsedByTrack(tracker_hit,j_lc_track)) { - //The hit is not already in the shared list - if (std::find(SharedHits[itrack].begin(), SharedHits[itrack].end(),tracker_hit->getID()) == SharedHits[itrack].end()) { - SharedHits[itrack].push_back(tracker_hit->getID()); - if (tracker_hit->getLayer() == 0 ) - SharedHitsLy0[itrack] = true; - if (tracker_hit->getLayer() == 1 ) - SharedHitsLy1[itrack] = true; - } - if (std::find(SharedHits[jtrack].begin(), SharedHits[jtrack].end(),tracker_hit->getID()) == SharedHits[jtrack].end()) { - SharedHits[jtrack].push_back(tracker_hit->getID()); - if (tracker_hit->getLayer() == 0 ) - SharedHitsLy0[jtrack] = true; - if (tracker_hit->getLayer() == 1 ) - SharedHitsLy1[jtrack] = true; - } - } // found shared hit - } // loop on j>i tracks - }//tracker hits - - track->setNShared(SharedHits[itrack].size()); - track->setSharedLy0(SharedHitsLy0[itrack]); - track->setSharedLy1(SharedHitsLy1[itrack]); - } - //Get the truth tracks relations: // Get the collection of LCRelations between GBL kink data and track data variables @@ -298,27 +297,37 @@ bool TrackingProcessor::process(IEvent* ievent) { } } - tracks_.push_back(track); //Do the residual plots -- should be in another function if (doResiduals_) { EVENT::LCCollection* trackRes_data_rel{nullptr}; + EVENT::LCCollection* trackXKink_data_rel{nullptr}; + EVENT::LCCollection* trackZKink_data_rel{nullptr}; try { if (!trackResDataLcio_.empty()) trackRes_data_rel = static_cast(event->getLCCollection(trackResDataLcio_.c_str())); - } + if (!trackXKinkDataLcio_.empty()) + trackXKink_data_rel = static_cast(event->getLCCollection(trackXKinkDataLcio_.c_str())); + if (!trackZKinkDataLcio_.empty()) + trackZKink_data_rel = static_cast(event->getLCCollection(trackZKinkDataLcio_.c_str())); + + } catch (EVENT::DataNotAvailableException e) { std::cout< trackRes_data_nav = std::make_shared(trackRes_data_rel); + std::shared_ptr trackXKink_data_nav = std::make_shared(trackXKink_data_rel); + std::shared_ptr trackZKink_data_nav = std::make_shared(trackZKink_data_rel); EVENT::LCObjectVec trackRes_data_vec = trackRes_data_nav->getRelatedFromObjects(lc_track); - IMPL::LCGenericObjectImpl* trackRes_data = static_cast(trackRes_data_vec.at(0)); - + EVENT::LCObjectVec trackXKink_data_vec = trackXKink_data_nav->getRelatedFromObjects(lc_track); + EVENT::LCObjectVec trackZKink_data_vec = trackZKink_data_nav->getRelatedFromObjects(lc_track); + IMPL::LCGenericObjectImpl* trackRes_data = static_cast(trackRes_data_vec.at(0)); + IMPL::LCGenericObjectImpl* trackXKink_data = static_cast(trackXKink_data_vec.at(0)); + IMPL::LCGenericObjectImpl* trackZKink_data = static_cast(trackZKink_data_vec.at(0)); /* //Some of the residuals do not get saved because sigma is negative. Will be fixed. if (track->getTrackerHitCount() != trackRes_data->getNDouble()) { @@ -326,18 +335,30 @@ bool TrackingProcessor::process(IEvent* ievent) { std::cout<<"Hits::"<getTrackerHitCount()<<" Residuals:"<getNDouble()<getNInt()-1;i_res++) { //std::cout<<"Residual ly " << trackRes_data->getIntVal(i_res)<<" res="<< trackRes_data->getDoubleVal(i_res)<<" sigma="<getFloatVal(i_res)<getIntVal(i_res); double res = trackRes_data->getDoubleVal(i_res); double sigma = trackRes_data->getFloatVal(i_res); - trkResHistos_->FillResidualHistograms(track,ly,res,sigma); - } + track->setTrackResid(ly,res); + if(doHistograms_){trkResHistos_->FillResidualHistograms(track,ly,res,sigma);} + } + for (int i_res = 0; i_res < trackXKink_data->getNInt()-1;i_res++) { + int ly = trackXKink_data->getIntVal(i_res); + double Xkink = trackXKink_data->getDoubleVal(i_res); + double sigma = trackXKink_data->getFloatVal(i_res); + track->setLambdaKink(ly,Xkink); + } + for (int i_res = 0; i_res < trackZKink_data->getNInt()-1;i_res++) { + int ly = trackZKink_data->getIntVal(i_res); + double Zkink = trackZKink_data->getDoubleVal(i_res); + double sigma = trackZKink_data->getFloatVal(i_res); + track->setPhiKink(ly,Zkink); + } }//trackResData exists }//doResiduals - + tracks_.push_back(track); }// tracks @@ -353,8 +374,7 @@ bool TrackingProcessor::process(IEvent* ievent) { } void TrackingProcessor::finalize() { - - if (doResiduals_) { + if (doResiduals_ and doHistograms_) { TFile* outfile = new TFile(resoutname_.c_str(),"RECREATE"); trkResHistos_->saveHistos(outfile,trkCollLcio_); if (trkResHistos_) delete trkResHistos_;