hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TrackingProcessor.cxx
Go to the documentation of this file.
1
2#include "TrackingProcessor.h"
3#include "utilities.h"
4
5TrackingProcessor::TrackingProcessor(const std::string& name, Process& process)
6 : Processor(name, process) {
7 }
8
11
13
14 std::cout << "Configuring TrackingProcessor" << std::endl;
15 try
16 {
17 debug_ = parameters.getInteger("debug", debug_);
18 trkCollLcio_ = parameters.getString("trkCollLcio", trkCollLcio_);
19 trkCollRoot_ = parameters.getString("trkCollRoot", trkCollRoot_);
20 kinkRelCollLcio_ = parameters.getString("kinkRelCollLcio", kinkRelCollLcio_);
21 trkRelCollLcio_ = parameters.getString("trkRelCollLcio", trkRelCollLcio_);
22 trkhitCollRoot_ = parameters.getString("trkhitCollRoot", trkhitCollRoot_);
23 hitFitsCollLcio_ = parameters.getString("hitFitsCollLcio", hitFitsCollLcio_);
24 rawhitCollRoot_ = parameters.getString("rawhitCollRoot", rawhitCollRoot_);
25 truthTracksCollLcio_ = parameters.getString("truthTrackCollLcio",truthTracksCollLcio_);
26 truthTracksCollRoot_ = parameters.getString("truthTrackCollRoot",truthTracksCollRoot_);
27 trackStateLocation_ = parameters.getString("trackStateLocation",trackStateLocation_);
28 useTrackerHits_ = parameters.getInteger("useTrackerHits",useTrackerHits_);
29 bfield_ = parameters.getDouble("bfield",bfield_);//this will overwrite the momentum; only use for pre-v3 slcio
30
31 //Residual plotting is done in this processor for the moment.
32 doResiduals_ = parameters.getInteger("doResiduals",doResiduals_);
33 trackResDataLcio_ = parameters.getString("trackResDataLcio",trackResDataLcio_);
34 resCfgFilename_ = parameters.getString("resPlots",resCfgFilename_);
35 resoutname_ = parameters.getString("resoutname",resoutname_);
36
37 }
38 catch (std::runtime_error& error)
39 {
40 std::cout << error.what() << std::endl;
41 }
42
43}
44
46 tree->Branch(trkCollRoot_.c_str(), &tracks_);
47
48 if (!trkhitCollRoot_.empty())
49 tree->Branch(trkhitCollRoot_.c_str(), &hits_);
50
51 if (!rawhitCollRoot_.empty())
52 tree->Branch(rawhitCollRoot_.c_str(), &rawhits_);
53
54 if (!truthTracksCollRoot_.empty())
55 tree->Branch(truthTracksCollRoot_.c_str(),&truthTracks_);
56
57
58 //Residual plotting
59 if (doResiduals_) {
65 }
66
67}
68
70
71
72 //Clean up
73 if (tracks_.size() > 0 ) {
74 for (std::vector<Track *>::iterator it = tracks_.begin(); it != tracks_.end(); ++it) {
75 delete *it;
76 }
77 tracks_.clear();
78 }
79
80 if (hits_.size() > 0) {
81 for (std::vector<TrackerHit *>::iterator it = hits_.begin(); it != hits_.end(); ++it) {
82 delete *it;
83 }
84 hits_.clear();
85 }
86
87 if (rawhits_.size() > 0) {
88 for (std::vector<RawSvtHit *>::iterator it = rawhits_.begin(); it != rawhits_.end(); ++it) {
89 delete *it;
90 }
91 rawhits_.clear();
92 }
93
94 if (truthTracks_.size() > 0) {
95 for (std::vector<Track *>::iterator it = truthTracks_.begin(); it != truthTracks_.end(); ++it) {
96 delete *it;
97 }
98 truthTracks_.clear();
99 }
100
101 Event* event = static_cast<Event*> (ievent);
102 // Get the collection of 3D hits from the LCIO event. If no such collection
103 // exist, a DataNotAvailableException is thrown
104
105 // Get decoders to read cellids
106 UTIL::BitField64 decoder("system:6,barrel:3,layer:4,module:12,sensor:1,side:32:-2,strip:12");
107
108 UTIL::LCRelationNavigator* rawTracker_hit_fits_nav = nullptr;
109 EVENT::LCCollection* raw_svt_hit_fits = nullptr;
110 //Check to see if fits are in the file
111 auto evColls = event->getLCEvent()->getCollectionNames();
112 auto it = std::find (evColls->begin(), evColls->end(), hitFitsCollLcio_.c_str());
113 bool hasFits = true;
114 if(it == evColls->end()) hasFits = false;
115 if(hasFits)
116 {
117 raw_svt_hit_fits = event->getLCCollection(hitFitsCollLcio_.c_str());
118 // Heap an LCRelation navigator which will allow faster access
119 rawTracker_hit_fits_nav = new UTIL::LCRelationNavigator(raw_svt_hit_fits);
120 }
121
122 EVENT::LCCollection* tracks{nullptr};
123 try
124 {
125 // Get all track collections from the event
126 tracks = event->getLCCollection(trkCollLcio_.c_str());
127 }
128 catch (EVENT::DataNotAvailableException e)
129 {
130 std::cout << e.what() << std::endl;
131 return false;
132 }
133
134 //Initialize map of shared hits
135 std::map <int, std::vector<int> > SharedHits;
136 //TODO: can we do better? (innermost)
137 std::map <int, bool> SharedHitsLy0;
138 std::map <int, bool> SharedHitsLy1;
139
140 for (int itrack = 0; itrack < tracks->getNumberOfElements();++itrack) {
141 SharedHits[itrack] = {};
142 SharedHitsLy0[itrack] = false;
143 SharedHitsLy1[itrack] = false;
144 }
145
146 // Loop over all the LCIO Tracks and add them to the HPS event.
147 for (int itrack = 0; itrack < tracks->getNumberOfElements(); ++itrack) {
148
149 // Get a LCIO Track from the LCIO event
150 EVENT::Track* lc_track = static_cast<EVENT::Track*>(tracks->getElementAt(itrack));
151
152 // Get the collection of LCRelations between GBL kink data and track data variables
153 // and the corresponding track.
154 EVENT::LCCollection* gbl_kink_data{nullptr};
155 EVENT::LCCollection* track_data{nullptr};
156 try
157 {
158 if (!kinkRelCollLcio_.empty())
159 gbl_kink_data = static_cast<EVENT::LCCollection*>(event->getLCCollection(kinkRelCollLcio_.c_str()));
160 if (!trkRelCollLcio_.empty())
161 track_data = static_cast<EVENT::LCCollection*>(event->getLCCollection(trkRelCollLcio_.c_str()));
162 }
163 catch (EVENT::DataNotAvailableException e)
164 {
165 std::cout << e.what() << std::endl;
166 if (!gbl_kink_data)
167 std::cout<<"TrackingProcessor::Failed retrieving " << kinkRelCollLcio_ <<std::endl;
168 if (!track_data)
169 std::cout<<"TrackingProcessor::Failed retrieving " << trkRelCollLcio_ <<std::endl;
170
171 }
172
173 // Add a track to the event
174 Track* track = utils::buildTrack(lc_track,trackStateLocation_, gbl_kink_data,track_data);
175
176 //Override the momentum of the track if the bfield_ > 0
177 if (bfield_>0)
178 track->setMomentum(bfield_); //this will overwrite the momentum; only use for pre-v3 slcio
179
180
181 int nHits = 0;
182
183 if(!useTrackerHits_){
184 auto hitPattern = lc_track->getSubdetectorHitNumbers();
185 for( int h=0; h<hitPattern.size(); h++){
186 if( hitPattern[h]>0 ){
187 track->addHitLayer(h);
188 nHits++;
189 }
190 }
191 if(nHits==0){
192 std::cout << "[TrackingProcessor::process] WARNING Track with no hits -- likely getSubdetectorHitNumbers isn't filled!" << std::endl;
193 std::cout << "[TrackingProcessor::process] WARNING If hit collections are saved, you may want to turn on useTrackerHits" << std::endl;
194 }
195 track->setTrackerHitCount(nHits);
196 }
197 else{
198
199 // Get the collection of hits associated with a LCIO Track
200 EVENT::TrackerHitVec lc_tracker_hits = lc_track->getTrackerHits();
201
202 // Iterate through the collection of 3D hits (TrackerHit objects)
203 // associated with a track, find the corresponding hits in the HPS
204 // event and add references to the track
205 bool rotateHits = true;
206 int hitType = 0;
207 if (track->isKalmanTrack())
208 hitType=1; //SiClusters
209
210 for (auto lc_tracker_hit : lc_tracker_hits) {
211
212 TrackerHit* tracker_hit = utils::buildTrackerHit(static_cast<IMPL::TrackerHitImpl*>(lc_tracker_hit),rotateHits,hitType);
213
214 std::vector<RawSvtHit*> rawSvthitsOn3d;
215 utils::addRawInfoTo3dHit(tracker_hit,static_cast<IMPL::TrackerHitImpl*>(lc_tracker_hit),
216 raw_svt_hit_fits,&rawSvthitsOn3d,hitType);
217
218 for (auto rhit : rawSvthitsOn3d)
219 rawhits_.push_back(rhit);
220
221 rawSvthitsOn3d.clear();
222
223 if (debug_)
224 std::cout<<tracker_hit->getRawHits().GetEntries()<<std::endl;
225 // Add a reference to the hit
226 track->addHit(tracker_hit);
227 track->addHitLayer(tracker_hit->getLayer());
228 hits_.push_back(tracker_hit);
229
230 //Get shared Hits information
231 for (int jtrack = itrack+1; jtrack < tracks->getNumberOfElements(); ++jtrack) {
232
233 EVENT::Track* j_lc_track = static_cast<EVENT::Track*>(tracks->getElementAt(jtrack));
234 if (utils::isUsedByTrack(tracker_hit,j_lc_track)) {
235 //The hit is not already in the shared list
236 if (std::find(SharedHits[itrack].begin(), SharedHits[itrack].end(),tracker_hit->getID()) == SharedHits[itrack].end()) {
237 SharedHits[itrack].push_back(tracker_hit->getID());
238 if (tracker_hit->getLayer() == 0 )
239 SharedHitsLy0[itrack] = true;
240 if (tracker_hit->getLayer() == 1 )
241 SharedHitsLy1[itrack] = true;
242 }
243 if (std::find(SharedHits[jtrack].begin(), SharedHits[jtrack].end(),tracker_hit->getID()) == SharedHits[jtrack].end()) {
244 SharedHits[jtrack].push_back(tracker_hit->getID());
245 if (tracker_hit->getLayer() == 0 )
246 SharedHitsLy0[jtrack] = true;
247 if (tracker_hit->getLayer() == 1 )
248 SharedHitsLy1[jtrack] = true;
249 }
250 } // found shared hit
251 } // loop on j>i tracks
252 }//tracker hits
253
254 track->setNShared(SharedHits[itrack].size());
255 track->setSharedLy0(SharedHitsLy0[itrack]);
256 track->setSharedLy1(SharedHitsLy1[itrack]);
257 }
258
259 //Get the truth tracks relations:
260
261 // Get the collection of LCRelations between GBL kink data and track data variables
262 // and the corresponding track.
263 EVENT::LCCollection* truth_tracks_rel{nullptr};
264
265 try
266 {
267 if (!truthTracksCollLcio_.empty())
268 truth_tracks_rel = static_cast<EVENT::LCCollection*>(event->getLCCollection(truthTracksCollLcio_.c_str()));
269 }
270 catch (EVENT::DataNotAvailableException e)
271 {
272 std::cout << e.what() << std::endl;
273 if (!truth_tracks_rel)
274 std::cout<<"Failed retrieving " << truthTracksCollLcio_ <<std::endl;
275 }
276
277
278 if (truth_tracks_rel) {
279
280 std::shared_ptr<UTIL::LCRelationNavigator> truth_tracks_nav = std::make_shared<UTIL::LCRelationNavigator>(truth_tracks_rel);
281 //Get the truth_track associated with the lcio_track
282 EVENT::LCObjectVec lc_truth_tracks = truth_tracks_nav->getRelatedToObjects(lc_track);
283 if (lc_truth_tracks.size() < 1) {
284 std::cout<<"Track with id "<<lc_track->id()<< " doesn't have a truth matched track "<<std::endl;
285 }
286 else {
287 EVENT::Track* lc_truth_track = static_cast<EVENT::Track*> (lc_truth_tracks.at(0));
288 Track* truth_track = utils::buildTrack(lc_truth_track,trackStateLocation_,nullptr,nullptr);
289 track->setTruthLink(truth_track);
290 if (bfield_>0)
291 truth_track->setMomentum(bfield_);//this will overwrite the momentum; only use for pre-v3 slcio
292 //truth tracks phi needs to be corrected
293 if (truth_track->getPhi() > TMath::Pi())
294 truth_track->setPhi(truth_track->getPhi() - (TMath::Pi()) * 2.);
295
296 truthTracks_.push_back(truth_track);
297 }
298
299 }
300 tracks_.push_back(track);
301
302
303
304 //Do the residual plots -- should be in another function
305 if (doResiduals_) {
306 EVENT::LCCollection* trackRes_data_rel{nullptr};
307 try {
308 if (!trackResDataLcio_.empty())
309 trackRes_data_rel = static_cast<EVENT::LCCollection*>(event->getLCCollection(trackResDataLcio_.c_str()));
310 }
311 catch (EVENT::DataNotAvailableException e)
312 {
313 std::cout<<e.what()<<std::endl;
314 }
315
316 if (trackRes_data_rel) {
317 std::shared_ptr<UTIL::LCRelationNavigator> trackRes_data_nav = std::make_shared<UTIL::LCRelationNavigator>(trackRes_data_rel);
318 EVENT::LCObjectVec trackRes_data_vec = trackRes_data_nav->getRelatedFromObjects(lc_track);
319 IMPL::LCGenericObjectImpl* trackRes_data = static_cast<IMPL::LCGenericObjectImpl*>(trackRes_data_vec.at(0));
320
321 /*
322 //Some of the residuals do not get saved because sigma is negative. Will be fixed.
323 if (track->getTrackerHitCount() != trackRes_data->getNDouble()) {
324 std::cout<<"WARNING-Different number of hit on track residuals wrt measurements"<<std::endl;
325 std::cout<<"Hits::"<<track->getTrackerHitCount()<<" Residuals:"<<trackRes_data->getNDouble()<<std::endl;
326 }
327 */
328
329 //Last int is the volume
330 for (int i_res = 0; i_res < trackRes_data->getNInt()-1;i_res++) {
331 //std::cout<<"Residual ly " << trackRes_data->getIntVal(i_res)<<" res="<< trackRes_data->getDoubleVal(i_res)<<" sigma="<<trackRes_data->getFloatVal(i_res)<<std::endl;
332 int ly = trackRes_data->getIntVal(i_res);
333 double res = trackRes_data->getDoubleVal(i_res);
334 double sigma = trackRes_data->getFloatVal(i_res);
335 trkResHistos_->FillResidualHistograms(track,ly,res,sigma);
336 }
337 }//trackResData exists
338 }//doResiduals
339
340
341 }// tracks
342
343 //delete
344 if (rawTracker_hit_fits_nav) {
345 delete rawTracker_hit_fits_nav; rawTracker_hit_fits_nav = nullptr;}
346
347 //event->addCollection("TracksInfo", &tracks_);
348 //event->addCollection("TrackerHitsInfo", &hits_);
349 //event->addCollection("TrackerHitsRawInfo", &rawhits_);
350
351 return true;
352}
353
355
356 if (doResiduals_) {
357 TFile* outfile = new TFile(resoutname_.c_str(),"RECREATE");
359 if (trkResHistos_) delete trkResHistos_;
360 trkResHistos_=nullptr;
361 }
362}
363
#define DECLARE_PROCESSOR(CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Definition Processor.h:139
Definition Event.h:35
virtual void DefineHistos()
Definition of histograms from json config.
void debugMode(bool debug)
set debug
virtual void loadHistoConfig(const std::string histoConfigFile)
load histogram config
virtual void saveHistos(TFile *outF=nullptr, std::string folder="")
save histograms
Definition IEvent.h:7
description
Base class for all event processing components.
Definition Processor.h:34
virtual bool process()
Process the histograms and generate analysis output.
Definition Processor.h:95
description
Definition TrackHistos.h:17
void FillResidualHistograms(Track *track, int ly, double res, double sigma)
Fill residual histograms.
void doTrackComparisonPlots(bool doplots)
Compare tracks.
Definition Track.h:32
double getPhi() const
Definition Track.h:88
void setMomentum(double bfield=0.52)
Definition Track.cxx:63
void setPhi(const double phi0)
Definition Track.h:89
int getLayer() const
Definition TrackerHit.h:109
int getID() const
Definition TrackerHit.h:124
TRefArray getRawHits() const
Definition TrackerHit.h:36
Insert description here. more details.
std::string truthTracksCollRoot_
description
std::string trkRelCollLcio_
collection name
std::string trackStateLocation_
Specify track state used for track collection DEFAULT AtIP.
int debug_
Debug Level.
std::string trkhitCollRoot_
description
double bfield_
magnetic field
std::string trkCollRoot_
collection name
int useTrackerHits_
Load hit collections, otherwise get from getSubdetectorHitNumbers.
virtual void configure(const ParameterSet &parameters)
Configure the Processor.
virtual void finalize()
Callback for the Processor to take any necessary action when the processing of events finishes.
std::string trackResDataLcio_
description
std::vector< Track * > tracks_
virtual void initialize(TTree *tree)
Callback for the Processor to take any necessary action when the processing of events starts.
std::string hitFitsCollLcio_
collection name
std::string resoutname_
description
std::vector< RawSvtHit * > rawhits_
std::vector< TrackerHit * > hits_
std::vector< Track * > truthTracks_
std::string truthTracksCollLcio_
description
std::string resCfgFilename_
description
std::string trkCollLcio_
collection name
int doResiduals_
do Residuals
std::string kinkRelCollLcio_
collection name
TrackHistos * trkResHistos_
description
std::string rawhitCollRoot_
collection name
TrackingProcessor(const std::string &name, Process &process)
Class constructor.
bool addRawInfoTo3dHit(TrackerHit *tracker_hit, IMPL::TrackerHitImpl *lc_tracker_hit, EVENT::LCCollection *raw_svt_fits, std::vector< RawSvtHit * > *rawHits=nullptr, int type=0, bool storeRawHit=true)
description
Track * buildTrack(EVENT::Track *lc_track, std::string trackstate_location, EVENT::LCCollection *gbl_kink_data, EVENT::LCCollection *track_data)
description
bool isUsedByTrack(IMPL::TrackerHitImpl *lc_tracker_hit, EVENT::Track *lc_track)
description
TrackerHit * buildTrackerHit(IMPL::TrackerHitImpl *lc_trackerHit, bool rotate=true, int type=0)
description