Skip to content

Commit

Permalink
Merge pull request #75 from kmcdermo/full-det-tracking
Browse files Browse the repository at this point in the history 
Move segmentMap_ outside of Event class, remove ETASEG ifdef
  • Loading branch information
@osschar
osschar authored Apr 15, 2017
2 parents 1d8db95 + ddb1d05 commit a172f45
Show file tree
Hide file tree
Showing 7 changed files with 47 additions and 101 deletions.
59 changes: 16 additions & 43 deletions Event.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ static bool tracksByPhi(const Track& t1, const Track& t2)
return t1.posPhi()<t2.posPhi();
}

#ifdef ETASEG
inline bool sortByEta(const Hit& hit1, const Hit& hit2){
return hit1.eta()<hit2.eta();
}
Expand All @@ -35,14 +34,12 @@ inline bool sortByEta(const Hit& hit1, const Hit& hit2){
inline bool sortByZ(const Hit& hit1, const Hit& hit2){
return hit1.z()<hit2.z();
}
#endif

Event::Event(const Geometry& g, Validation& v, int evtID, int threads) :
geom_(g), validation_(v),
evtID_(evtID), threads_(threads), mcHitIDCounter_(0)
{
layerHits_.resize(Config::nTotalLayers);
segmentMap_.resize(Config::nTotalLayers);

validation_.resetValidationMaps(); // need to reset maps for every event.
}
Expand All @@ -53,7 +50,6 @@ void Event::Reset(int evtID)
mcHitIDCounter_ = 0;

for (auto&& l : layerHits_) { l.clear(); }
for (auto&& l : segmentMap_) { l.clear(); }

simHitsInfo_.clear();
simTrackStates_.clear();
Expand Down Expand Up @@ -196,17 +192,18 @@ void Event::Simulate()
simTrackStates_ = tmpTSVec;
}

void Event::Segment()
void Event::Segment(BinInfoMap & segmentMap)
{
#ifdef DEBUG
bool debug=true;
#endif
segmentMap.resize(Config::nTotalLayers);

//sort in phi and dump hits per layer, fill phi partitioning
for (int ilayer=0; ilayer<layerHits_.size(); ++ilayer) {
dprint("Hits in layer=" << ilayer);

#ifdef ETASEG
segmentMap_[ilayer].resize(Config::nEtaPart);
segmentMap[ilayer].resize(Config::nEtaPart);
// eta first then phi
std::sort(layerHits_[ilayer].begin(), layerHits_[ilayer].end(), sortByZ);
std::vector<int> lay_eta_bin_count(Config::nEtaPart);
Expand Down Expand Up @@ -243,56 +240,32 @@ void Event::Segment()
int firstPhiBinIdx = lastPhiIdxFound+1;
int phiBinSize = lay_eta_phi_bin_count[phibin];
BinInfo phiBinInfo(firstPhiBinIdx,phiBinSize);
segmentMap_[ilayer][etabin].push_back(phiBinInfo);
segmentMap[ilayer][etabin].push_back(phiBinInfo);
if (phiBinSize>0){
lastPhiIdxFound+=phiBinSize;
}
#ifdef DEBUG
if ((debug) && (phiBinSize !=0)) dprintf("ilayer: %1u etabin: %1u phibin: %2u first: %2u last: %2u \n",
ilayer, etabin, phibin,
segmentMap_[ilayer][etabin][phibin].first,
segmentMap_[ilayer][etabin][phibin].second+segmentMap_[ilayer][etabin][phibin].first
segmentMap[ilayer][etabin][phibin].first,
segmentMap[ilayer][etabin][phibin].second+segmentMap[ilayer][etabin][phibin].first
);
#endif
} // end loop over storing phi index
} // end loop over storing eta index
#else
segmentMap_[ilayer].resize(1); // only one eta bin for special case, avoid ifdefs
std::sort(layerHits_[ilayer].begin(), layerHits_[ilayer].end(), sortByPhi);
std::vector<int> lay_phi_bin_count(Config::nPhiPart);//should it be 63? - yes!
for (int ihit=0;ihit<layerHits_[ilayer].size();++ihit) {
dprint("hit r/phi/eta : " << layerHits_[ilayer][ihit].r() << " "
<< layerHits_[ilayer][ihit].phi() << " " << layerHits_[ilayer][ihit].eta());

int phibin = getPhiPartition(layerHits_[ilayer][ihit].phi());
lay_phi_bin_count[phibin]++;
}

//now set index and size in partitioning map
int lastIdxFound = -1;
for (int bin=0; bin<Config::nPhiPart; ++bin) {
int binSize = lay_phi_bin_count[bin];
int firstBinIdx = lastIdxFound+1;
BinInfo binInfo(firstBinIdx, binSize);
segmentMap_[ilayer][0].push_back(binInfo); // [0] bin is just the only eta bin ... reduce ifdefs
if (binSize>0){
lastIdxFound+=binSize;
}
}
#endif
} // end loop over layers

#ifdef DEBUG
for (int ilayer = 0; ilayer < Config::nLayers; ilayer++) {
dmutex_guard;
int etahitstotal = 0;
for (int etabin = 0; etabin < Config::nEtaPart; etabin++){
int etahits = segmentMap_[ilayer][etabin][Config::nPhiPart-1].first + segmentMap_[ilayer][etabin][Config::nPhiPart-1].second - segmentMap_[ilayer][etabin][0].first;
int etahits = segmentMap[ilayer][etabin][Config::nPhiPart-1].first + segmentMap[ilayer][etabin][Config::nPhiPart-1].second - segmentMap[ilayer][etabin][0].first;
std::cout << "etabin: " << etabin << " hits in bin: " << etahits << std::endl;
etahitstotal += etahits;

for (int phibin = 0; phibin < Config::nPhiPart; phibin++){
// if (segmentMap_[ilayer][etabin][phibin].second > 3) {std::cout << " phibin: " << phibin << " hits: " << segmentMap_[ilayer][etabin][phibin].second << std::endl;}
// if (segmentMap[ilayer][etabin][phibin].second > 3) {std::cout << " phibin: " << phibin << " hits: " << segmentMap[ilayer][etabin][phibin].second << std::endl;}
}
}
std::cout << "layer: " << ilayer << " totalhits: " << etahitstotal << std::endl;
Expand All @@ -303,12 +276,12 @@ void Event::Segment()
RemapHits(simTracks_);
}

void Event::Seed()
void Event::Seed(const BinInfoMap & segmentMap)
{
#ifdef ENDTOEND
buildSeedsByRoadSearch(seedTracks_,seedTracksExtra_,layerHits_,segmentMap_,*this);
// buildSeedsByRoadTriplets(seedTracks_,seedTracksExtra_,layerHits_,segmentMap_,*this);
//buildSeedsByRZFirstRPhiSecond(seedTracks_,seedTracksExtra_,layerHits_,segmentMap_,*this);
buildSeedsByRoadSearch(seedTracks_,seedTracksExtra_,layerHits_,segmentMap,*this);
// buildSeedsByRoadTriplets(seedTracks_,seedTracksExtra_,layerHits_,segmentMap,*this);
//buildSeedsByRZFirstRPhiSecond(seedTracks_,seedTracksExtra_,layerHits_,segmentMap,*this);
#else
buildSeedsByMC(simTracks_,seedTracks_,seedTracksExtra_,*this);
simTracksExtra_ = seedTracksExtra_;
Expand All @@ -317,10 +290,10 @@ void Event::Seed()
validation_.alignTrackExtra(seedTracks_,seedTracksExtra_); // if we sort here, also have to sort seedTracksExtra and redo labels.
}

void Event::Find()
void Event::Find(const BinInfoMap & segmentMap)
{
buildTracksBySeeds(*this);
// buildTracksByLayers(*this);
buildTracksBySeeds(segmentMap,*this);
// buildTracksByLayers(segmentMap,*this);

// From CHEP-2015
// buildTestSerial(*this, Config::nlayers_per_seed, Config::maxCandsPerSeed, Config::chi2Cut, Config::nSigma, Config::minDPhi);
Expand Down
10 changes: 3 additions & 7 deletions Event.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@ class Event
void Reset(int evtID);
void RemapHits(TrackVec & tracks);
void Simulate();
void Segment();
void Seed();
void Find();
void Segment(BinInfoMap & segmentMap);
void Seed(const BinInfoMap & segmentMap);
void Find(const BinInfoMap & segmentMap);
void Fit();
void Validate();
void PrintStats(const TrackVec&, TrackExtraVec&);
Expand Down Expand Up @@ -52,10 +52,6 @@ class Event
// There should be a better way of doing that.
int seedEtaSeparators_[5];

// phi-eta partitioning map: vector of vector of vectors of std::pairs.
// vec[nLayers][nEtaBins][nPhiBins]
BinInfoMap segmentMap_;

TSVec simTrackStates_;
static std::mutex printmutex;
};
Expand Down
7 changes: 2 additions & 5 deletions Makefile.config
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,10 +109,7 @@ ENDTOEND := -DENDTOEND
# have to be adjusted.
WITH_TBB := yes

# 14. Use for eta-phi segmentation
USE_ETA_SEGMENTATION := -DETASEG

# 15. Use inward fit in Conformal fit + final KF Fit
# 14. Use inward fit in Conformal fit + final KF Fit
INWARD_FIT := -DINWARDFIT

################################################################
Expand All @@ -137,7 +134,7 @@ ifdef USE_CUDA
endif
endif

CPPFLAGS += ${USE_STATE_VALIDITY_CHECKS} ${USE_SCATTERING} ${USE_LINEAR_INTERPOLATION} ${ENDTOEND} ${USE_ETA_SEGMENTATION} ${INWARD_FIT} ${GEN_FLAT_ETA}
CPPFLAGS += ${USE_STATE_VALIDITY_CHECKS} ${USE_SCATTERING} ${USE_LINEAR_INTERPOLATION} ${ENDTOEND} ${INWARD_FIT}

ifdef USE_VTUNE_NOTIFY
ifdef VTUNE_AMPLIFIER_XE_2016_DIR
Expand Down
27 changes: 11 additions & 16 deletions buildtest.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,15 @@ static std::mutex evtlock;
#endif
typedef candvec::const_iterator canditer;

void extendCandidate(const Event& ev, const cand_t& cand, candvec& tmp_candidates, int ilay, int seedID, bool debug);
void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_t& cand, candvec& tmp_candidates, int ilay, int seedID, bool debug);

inline bool sortByHitsChi2(const cand_t& cand1, const cand_t& cand2)
{
if (cand1.nFoundHits()==cand2.nFoundHits()) return cand1.chi2()<cand2.chi2();
return cand1.nFoundHits()>cand2.nFoundHits();
}

void processCandidates(Event& ev, candvec& candidates, int ilay, int seedID, const bool debug)
void processCandidates(const BinInfoMap & segmentMap, Event& ev, candvec& candidates, int ilay, int seedID, const bool debug)
{
auto& evt_track_candidates(ev.candidateTracks_);

Expand All @@ -44,7 +44,7 @@ void processCandidates(Event& ev, candvec& candidates, int ilay, int seedID, con
{
//loop over running candidates
for (auto&& cand : candidates) {
extendCandidate(ev, cand, tmp_candidates, ilay, seedID, debug);
extendCandidate(segmentMap, ev, cand, tmp_candidates, ilay, seedID, debug);
}
if (tmp_candidates.size()>Config::maxCandsPerSeed) {
dprint("huge size=" << tmp_candidates.size() << " keeping best "<< Config::maxCandsPerSeed << " only");
Expand All @@ -69,7 +69,7 @@ void processCandidates(Event& ev, candvec& candidates, int ilay, int seedID, con
}
}

void buildTracksBySeeds(Event& ev)
void buildTracksBySeeds(const BinInfoMap & segmentMap, Event& ev)
{
auto& evt_track_candidates(ev.candidateTracks_);
const auto& evt_lay_hits(ev.layerHits_);
Expand Down Expand Up @@ -101,7 +101,7 @@ void buildTracksBySeeds(Event& ev)
auto&& candidates(track_candidates[iseed]);
for (int ilay=Config::nlayers_per_seed;ilay<evt_lay_hits.size();++ilay) {//loop over layers, starting from after the seed
dprint("going to layer #" << ilay << " with N cands=" << track_candidates.size());
processCandidates(ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
processCandidates(segmentMap, ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
}
//end of layer loop
}//end of process seeds loop
Expand All @@ -123,7 +123,7 @@ void buildTracksBySeeds(Event& ev)
}
}

void buildTracksByLayers(Event& ev)
void buildTracksByLayers(const BinInfoMap & segmentMap, Event& ev)
{
auto& evt_track_candidates(ev.candidateTracks_);
const auto& evt_lay_hits(ev.layerHits_);
Expand All @@ -147,15 +147,15 @@ void buildTracksByLayers(Event& ev)
[&](const tbb::blocked_range<size_t>& seediter) {
for (auto iseed = seediter.begin(); iseed != seediter.end(); ++iseed) {
auto&& candidates(track_candidates[iseed]);
processCandidates(ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
processCandidates(segmentMap, ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
}
}); //end of process seeds loop
#else
//process seeds
for (auto iseed = 0; iseed != evt_seeds.size(); ++iseed) {
const auto& seed(evt_seeds[iseed]);
auto&& candidates(track_candidates[iseed]);
processCandidates(ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
processCandidates(segmentMap, ev, candidates, ilay, evt_seeds_extra[iseed].seedID(), debug);
}
#endif
} //end of layer loop
Expand All @@ -174,13 +174,13 @@ void buildTracksByLayers(Event& ev)
}
}

void extendCandidate(const Event& ev, const cand_t& cand, candvec& tmp_candidates, int ilayer, int seedID, bool debug)
void extendCandidate(const BinInfoMap & segmentMap, const Event& ev, const cand_t& cand, candvec& tmp_candidates, int ilayer, int seedID, bool debug)
{
std::vector<int> branch_hit_indices; // temp variable for validation... could be used for cand hit builder engine!
const Track& tkcand = cand;
const TrackState& updatedState = cand.state();
const auto& evt_lay_hits(ev.layerHits_);
const auto& segLayMap(ev.segmentMap_[ilayer]);
const auto& segLayMap(segmentMap[ilayer]);

dprint("processing candidate with nHits=" << tkcand.nFoundHits());
#ifdef LINEARINTERP
Expand All @@ -200,17 +200,12 @@ void extendCandidate(const Event& ev, const cand_t& cand, candvec& tmp_candidate
const float predy = propState.parameters.At(1);
const float predz = propState.parameters.At(2);

#ifdef ETASEG
const float eta = getEta(std::sqrt(getRad2(predx,predy)),predz);
const float deta = std::sqrt(std::abs(getEtaErr2(predx,predy,predz,propState.errors.At(0,0),propState.errors.At(1,1),propState.errors.At(2,2),propState.errors.At(0,1),propState.errors.At(0,2),propState.errors.At(1,2))));
const float nSigmaDeta = std::min(std::max(Config::nSigma*deta,(float) Config::minDEta), (float) Config::maxDEta); // something to tune -- minDEta = 0.0
const auto etaBinMinus = getEtaPartition(eta-nSigmaDeta);
const auto etaBinPlus = getEtaPartition(eta+nSigmaDeta);
#else
const float nSigmaDeta = 0.;
const auto etaBinMinus = 0;
const auto etaBinPlus = 0;
#endif

const float phi = getPhi(predx,predy); //std::atan2(predy,predx);
const float dphi = std::sqrt(std::abs(getPhiErr2(predx,predy,propState.errors.At(0,0),propState.errors.At(1,1),propState.errors.At(0,1))));
const float nSigmaDphi = std::min(std::max(Config::nSigma*dphi,(float) Config::minDPhi), (float) Config::maxDPhi);
Expand Down
5 changes: 3 additions & 2 deletions buildtest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@
#include "Track.h"
#include "Geometry.h"
#include "Event.h"
#include "BinInfoUtils.h"

void buildTracksBySeeds(Event& ev);
void buildTracksByLayers(Event& ev);
void buildTracksBySeeds(const BinInfoMap &, Event& ev);
void buildTracksByLayers(const BinInfoMap &, Event& ev);
#endif
17 changes: 11 additions & 6 deletions main.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "Matrix.h"
#include "Event.h"
#include "Validation.h"
#include "BinInfoUtils.h"

#include "fittestEndcap.h"

Expand Down Expand Up @@ -285,12 +286,16 @@ int main(int argc, const char* argv[])
continue;
}

/*simulate time*/ ticks[0] += delta(t0);
ev.Segment(); ticks[1] += delta(t0);
ev.Seed(); ticks[2] += delta(t0);
ev.Find(); ticks[3] += delta(t0);
ev.Fit(); ticks[4] += delta(t0);
ev.Validate(); ticks[5] += delta(t0);
// phi-eta partitioning map: vector of vector of vectors of std::pairs.
// vec[nLayers][nEtaBins][nPhiBins]
BinInfoMap segmentMap;

/*simulate time*/ ticks[0] += delta(t0);
ev.Segment(segmentMap); ticks[1] += delta(t0);
ev.Seed(segmentMap); ticks[2] += delta(t0);
ev.Find(segmentMap); ticks[3] += delta(t0);
ev.Fit(); ticks[4] += delta(t0);
ev.Validate(); ticks[5] += delta(t0);

std::cout << "sim: " << ev.simTracks_.size() << " seed: " << ev.seedTracks_.size() << " found: "
<< ev.candidateTracks_.size() << " fit: " << ev.fitTracks_.size() << std::endl;
Expand Down
23 changes: 1 addition & 22 deletions seedtest.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,13 +314,8 @@ void buildHitPairs(const std::vector<HitVec>& evt_lay_hits, const BinInfoLayerMa
const float outerhitz = evt_lay_hits[1][ihit].z(); // remember, layer[0] is first layer! --> second layer = [1]
const float outerphi = evt_lay_hits[1][ihit].phi();

#ifdef ETASEG // z cut is the largest killer... up to 3 sigma is probably best
const auto etaBinMinus = getEtaPartition(getEta(Config::fRadialSpacing,(outerhitz-Config::seed_z0cut)/2.));
const auto etaBinPlus = getEtaPartition(getEta(Config::fRadialSpacing,(outerhitz+Config::seed_z0cut)/2.));
#else
const auto etaBinMinus = 0;
const auto etaBinPlus = 0;
#endif
const auto phiBinMinus = getPhiPartition(normalizedPhi(outerphi - Config::lay01angdiff));
const auto phiBinPlus = getPhiPartition(normalizedPhi(outerphi + Config::lay01angdiff));

Expand Down Expand Up @@ -382,24 +377,12 @@ void buildHitTripletsCurve(const std::vector<HitVec>& evt_lay_hits, const BinInf
intersectThirdLayer(aneg,bneg,x1,y1,negx2,negy2);
const float negPhi = getPhi(negx2,negy2);

#ifdef ETASEG
const float thirdZline = 2*hit1.z()-hit0.z(); // for dz displacements -- straight line window
const auto etaBinMinus = getEtaPartition(getEta(lay3rad,thirdZline)-Config::dEtaSeedTrip);
const auto etaBinPlus = getEtaPartition(getEta(lay3rad,thirdZline)+Config::dEtaSeedTrip);
#else
const auto etaBinMinus = 0;
const auto etaBinPlus = 0;
#endif
const auto phiBinMinus = getPhiPartition(negPhi);
const auto phiBinPlus = getPhiPartition(posPhi);

#ifdef BROKEN_DEBUG
const float lay2phi = evt_lay_hits[2][ev.simTracks_[ev.simHitsInfo_[hit0.mcHitID()].mcTrackID()].getHitIdx(2)].phi();
dprint("lay0 phi: " << hit0.phi() << " lay1 phi: " << hit1.phi() << std::endl <<
"negPhi: " << negPhi << " lay2 phi: " << lay2phi << " posPhi: " << posPhi << std::endl <<
"binM: " << phiBinMinus << " phi2Bin: " << getPhiPartition(lay2phi) << " binP: " << phiBinPlus << std::endl);
#endif

std::vector<int> cand_hit_indices = getCandHitIndices(etaBinMinus,etaBinPlus,phiBinMinus,phiBinPlus,segLayMap);
for (auto&& cand_hit_idx : cand_hit_indices){
TripletIdx hit_triplet;
Expand Down Expand Up @@ -439,14 +422,10 @@ void buildHitTripletsApproxWindow(const std::vector<HitVec>& evt_lay_hits, const
const Hit& hit0 = evt_lay_hits[0][hit_pair[0]];
const Hit& hit1 = evt_lay_hits[1][hit_pair[1]];

#ifdef ETASEG
const float thirdZline = 2*hit1.z()-hit0.z(); // for dz displacements -- straight line window
const auto etaBinMinus = getEtaPartition(getEta(thirdRad,thirdZline)-Config::dEtaSeedTrip);
const auto etaBinPlus = getEtaPartition(getEta(thirdRad,thirdZline)+Config::dEtaSeedTrip);
#else
const auto etaBinMinus = 0;
const auto etaBinPlus = 0;
#endif

const float linePhi = getPhi(hit1.position()[0] - hit0.position()[0], hit1.position()[1] - hit0.position()[1]);
float thirdPhiMinus = 0.0f;
float thirdPhiPlus = 0.0f;
Expand Down

0 comments on commit a172f45

Please sign in to comment.