Skip to content

This is a hpstr pull request meant to correspond to the Java Issue 1079 #206

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rodwyer100 wants to merge 13 commits into
base: master
Choose a base branch
from
Open

This is a hpstr pull request meant to correspond to the Java Issue 1079 #206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
814b6d7
Add files via upload
rodwyer100 Feb 5, 2025
4716502
Add files via upload
rodwyer100 Feb 5, 2025
18f957c
Add files via upload
rodwyer100 Feb 5, 2025
1aa5412
Add files via upload
rodwyer100 Feb 5, 2025
d086556
Update utilities.cxx
rodwyer100 Feb 5, 2025
a326975
Update TrackingProcessor.h
rodwyer100 Feb 5, 2025
cda5882
Update TrackingProcessor.h
rodwyer100 Feb 5, 2025
a12fe24
Update TrackingProcessor.h
rodwyer100 Feb 5, 2025
d45adab
Update TrackingProcessor.cxx
rodwyer100 Feb 5, 2025
6359cee
Update TrackingProcessor.cxx
rodwyer100 Feb 6, 2025
1ae07bb
Add files via upload
rodwyer100 Feb 6, 2025
b548f7b
Add files via upload
rodwyer100 Feb 6, 2025
269552e
Update TrackingProcessor.cxx
rodwyer100 Feb 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 14 additions & 1 deletion event/include/Track.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*
Expand Down Expand Up @@ -386,7 +396,7 @@ class Track : public TObject {
/** hit layers */
std::vector<int> hit_layers_;

/** truth mcp hits */
/** truth mcp hits */
std::vector<std::pair<int,int>> mcp_hits_;

/** The distance of closest approach to the reference point. */
Expand Down Expand Up @@ -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];

Expand Down
12 changes: 8 additions & 4 deletions processors/include/TrackingProcessor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
#include "Event.h"
#include "RawSvtHit.h"
#include "TrackHistos.h"

#include "MCParticle.h"

// Forward declarations
class TTree;
Expand Down Expand Up @@ -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

Expand Down
180 changes: 100 additions & 80 deletions processors/src/TrackingProcessor.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_);

}
Expand All @@ -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_);
Expand Down Expand Up @@ -136,8 +138,7 @@ bool TrackingProcessor::process(IEvent* ievent) {
std::map <int, std::vector<int> > SharedHits;
//TODO: can we do better? (innermost)
std::map <int, bool> SharedHitsLy0;
std::map <int, bool> SharedHitsLy1;

std::map <int, bool> SharedHitsLy1;
for (int itrack = 0; itrack < tracks->getNumberOfElements();++itrack) {
SharedHits[itrack] = {};
SharedHitsLy0[itrack] = false;
Expand Down Expand Up @@ -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();
Expand All @@ -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<IMPL::TrackerHitImpl*>(lc_tracker_hit),rotateHits,hitType);

std::vector<RawSvtHit*> rawSvthitsOn3d;
utils::addRawInfoTo3dHit(tracker_hit,static_cast<IMPL::TrackerHitImpl*>(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<IMPL::TrackerHitImpl*>(lc_tracker_hit),rotateHits,hitType);

std::vector<RawSvtHit*> rawSvthitsOn3d;
utils::addRawInfoTo3dHit(tracker_hit,static_cast<IMPL::TrackerHitImpl*>(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<<tracker_hit->getRawHits().GetEntries()<<std::endl;
// Add a reference to the hit
track->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<EVENT::Track*>(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<<tracker_hit->getRawHits().GetEntries()<<std::endl;
// Add a reference to the hit
track->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<EVENT::Track*>(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
Expand Down Expand Up @@ -298,46 +297,68 @@ 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::LCCollection*>(event->getLCCollection(trackResDataLcio_.c_str()));
}
if (!trackXKinkDataLcio_.empty())
trackXKink_data_rel = static_cast<EVENT::LCCollection*>(event->getLCCollection(trackXKinkDataLcio_.c_str()));
if (!trackZKinkDataLcio_.empty())
trackZKink_data_rel = static_cast<EVENT::LCCollection*>(event->getLCCollection(trackZKinkDataLcio_.c_str()));

}
catch (EVENT::DataNotAvailableException e)
{
std::cout<<e.what()<<std::endl;
}

if (trackRes_data_rel) {
if ((trackRes_data_rel)and(trackXKink_data_rel)and(trackZKink_data_rel)) {
std::shared_ptr<UTIL::LCRelationNavigator> trackRes_data_nav = std::make_shared<UTIL::LCRelationNavigator>(trackRes_data_rel);
std::shared_ptr<UTIL::LCRelationNavigator> trackXKink_data_nav = std::make_shared<UTIL::LCRelationNavigator>(trackXKink_data_rel);
std::shared_ptr<UTIL::LCRelationNavigator> trackZKink_data_nav = std::make_shared<UTIL::LCRelationNavigator>(trackZKink_data_rel);
EVENT::LCObjectVec trackRes_data_vec = trackRes_data_nav->getRelatedFromObjects(lc_track);
IMPL::LCGenericObjectImpl* trackRes_data = static_cast<IMPL::LCGenericObjectImpl*>(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<IMPL::LCGenericObjectImpl*>(trackRes_data_vec.at(0));
IMPL::LCGenericObjectImpl* trackXKink_data = static_cast<IMPL::LCGenericObjectImpl*>(trackXKink_data_vec.at(0));
IMPL::LCGenericObjectImpl* trackZKink_data = static_cast<IMPL::LCGenericObjectImpl*>(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()) {
std::cout<<"WARNING-Different number of hit on track residuals wrt measurements"<<std::endl;
std::cout<<"Hits::"<<track->getTrackerHitCount()<<" Residuals:"<<trackRes_data->getNDouble()<<std::endl;
}
*/

//Last int is the volume
for (int i_res = 0; i_res < trackRes_data->getNInt()-1;i_res++) {
//std::cout<<"Residual ly " << trackRes_data->getIntVal(i_res)<<" res="<< trackRes_data->getDoubleVal(i_res)<<" sigma="<<trackRes_data->getFloatVal(i_res)<<std::endl;
int ly = trackRes_data->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

Expand All @@ -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_;
Expand Down