diff --git a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
index d3bcd2ef85b2371a70c758f1117753761526280e..9881d83a24e0933f025c9563f86619718d064902 100644
--- a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
+++ b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
@@ -28,22 +28,29 @@ inline edm4hep::TrackState getTrackStateAt(edm4hep::ConstTrack track, int locati
}
inline std::array<float,15> getCovMatrix(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).covMatrix;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).covMatrix;
+ std::array<float,15> dummy{};
+ return dummy;
}
inline float getTanLambda(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).tanLambda;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).tanLambda;
+ return 0;
}
inline float getOmega(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).omega;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).omega;
+ return 0;
}
inline float getD0(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).D0;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).D0;
+ return 0;
}
inline float getZ0(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).Z0;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).Z0;
+ return 0;
}
inline float getPhi(const edm4hep::ConstTrack &track) {
- return track.getTrackStates(0).phi;
+ if(track.trackStates_size()>0) return track.getTrackStates(0).phi;
+ return 0;
}
diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
index ed76016f0f6cb516e6e47515edfcae6a07f35863..aed3b55e5f148ff98cc77b64b61bee3fdba42bd3 100644
--- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp
@@ -233,7 +233,7 @@ StatusCode ClupatraAlg::initialize() {
StatusCode ClupatraAlg::execute() {
- debug() << "Clupatra Algorithm started" << endmsg;
+ info() << "Clupatra Algorithm started" << endmsg;
// clock_t start = clock() ;
Timer timer ;
@@ -537,7 +537,7 @@ StatusCode ClupatraAlg::execute() {
debug() << " call fitter for seed cluster with " << (*icv)->size() << " hits " << endmsg;
int counter = 0;
for( Clusterer::cluster_type::iterator ci=(*icv)->begin(), end= (*icv)->end() ; ci!=end; ++ci ) {
- debug() << counter++ << " " << (*ci)->first->edm4hepHit << " \nlayer " << (*ci)->first->layer << endmsg;
+ debug() << counter++ << " " << (*ci)->first->edm4hepHit.id() << " layer " << (*ci)->first->layer << endmsg;
}
@@ -562,9 +562,8 @@ StatusCode ClupatraAlg::execute() {
ConstTrack edm4hepTrk( converter( *icv ) ) ;
// debug() << "Goes goes here" << endmsg;
- // FIXME Mingrui
- debug() << "============= poor seed cluster - no hits added - started from row " << outerRow << "\n"
- << edm4hepTrk << endmsg;
+ debug() << "============= poor seed cluster - no hits added - started from row " << outerRow << " track id="
+ << edm4hepTrk.id() << endmsg;
for( Clusterer::cluster_type::iterator ci=(*icv)->begin(), end= (*icv)->end() ; ci!=end; ++ci ) {
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
index 400774f26b5aa2c19da438f47fbce05d0fbf2c20..3cffd0f0719998ec93f88631c6e50b9f794a81cd 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
@@ -288,7 +288,7 @@ namespace clupatra_new{
encoder[UTIL::ILDCellID0::subdet] = UTIL::ILDDetID::TPC ;
edm4hep::ConstTrackerHit firstHit; // = 0 ;
- firstHit.unlink();
+ //firstHit.unlink();
IMarlinTrack* bwTrk = 0 ;
@@ -1354,11 +1354,18 @@ start:
//IMPL::TrackerHitImpl* thi = dynamic_cast<IMPL::TrackerHitImpl*> ( hitsInFit[i].first ) ;
//thi->setQualityBit( UTIL::ILDTrkHitQualityBit::USED_IN_FIT , 1 ) ;
}
-
- edm4hep::TrackState tsIP;
- edm4hep::TrackState tsFH;
- edm4hep::TrackState tsLH;
- edm4hep::TrackState tsCA;
+ edm4hep::TrackState tsBase;
+ tsBase.D0 = 0;
+ tsBase.phi = 0;
+ tsBase.omega = 0;
+ tsBase.Z0 = 0;
+ tsBase.tanLambda = 0;
+ tsBase.referencePoint = edm4hep::Vector3f(0,0,0);
+ tsBase.covMatrix = std::array<float, 15>{};
+ edm4hep::TrackState tsIP(tsBase);
+ edm4hep::TrackState tsFH(tsBase);
+ edm4hep::TrackState tsLH(tsBase);
+ edm4hep::TrackState tsCA(tsBase);
tsIP.location = lcio::TrackState::AtIP;
tsFH.location = lcio::TrackState::AtFirstHit;
@@ -1385,10 +1392,9 @@ start:
code = mtrk->getTrackState( fHit, tsFH, chi2, ndf ) ;
if( code != MarlinTrk::IMarlinTrack::success ){
-
- std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at first Hit !!?? "
- << " error code : " << MarlinTrk::errorCode( code )
- << std::endl ;
+ //std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at first Hit !!?? "
+ // << " error code : " << MarlinTrk::errorCode( code )
+ // << std::endl ;
}
// ======= get TrackState at last hit ========================
@@ -1398,17 +1404,17 @@ start:
#if use_fit_at_last_hit
code = mtrk->getTrackState( lHit, tsLH, chi2, ndf ) ;
#else // get the track state at the last hit by propagating from the last(first) constrained fit position (a la MarlinTrkUtils)
- edm4hep::ConstTrackerHit last_constrained_hit;
+ edm4hep::ConstTrackerHit last_constrained_hit(0);
code = mtrk->getTrackerHitAtPositiveNDF( last_constrained_hit );
- mtrk->smooth() ;
+ //code = mtrk->smooth() ;
if( code != MarlinTrk::IMarlinTrack::success ){
- std::cout << "last_constrained_hit is unavaibale" << std::endl;
+ //std::cout << "last_constrained_hit is avaibale? " << last_constrained_hit.isAvailable() << std::endl;
}
- else{
+ //else{
// std::cout << "the hit" << std::endl;
// std::cout << lHit.getCellID0() << std::endl;
// std::cout << last_constrained_hit << std::endl;
- //code = mtrk->smooth() ;
+ code = mtrk->smooth() ;
// std::cout << "Smooth success" << std::endl;
// gear::Vector3D last_hit_pos( lHit->getPosition()[0], lHit->getPosition()[1], lHit->getPosition()[2] );
// std::cout << "lHit = " << lHit << std::endl;
@@ -1420,14 +1426,13 @@ start:
edm4hep::Vector3d last_hit_pos( lHit.getPosition() );
// std::cout << "lHit = " << lHit << std::endl;
code = mtrk->propagate( last_hit_pos, last_constrained_hit, tsLH, chi2, ndf);
- }
+ //}
// std::cout << "Propagate success" << std::endl;
#endif
// std::cout << "We are here" << std::endl;
if( code != MarlinTrk::IMarlinTrack::success ){
-
- std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at last Hit !!?? " << std::endl ;
+ //std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at last Hit !!?? " << std::endl ;
}
// ======= get TrackState at calo face ========================
@@ -1459,8 +1464,7 @@ start:
}
if ( code !=MarlinTrk::IMarlinTrack::success ) {
-
- // streamlog_out( DEBUG6 ) << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at calo face !!?? " << std::endl ;
+ //std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not get TrackState at calo face !!?? " << std::endl ;
}
// ======= get TrackState at IP ========================
@@ -1475,10 +1479,9 @@ start:
code = ( UsePropagate ? mtrk->propagate( ipv, fHit, tsIP, chi2, ndf ) : mtrk->extrapolate( ipv, tsIP, chi2, ndf ) ) ;
if( code != MarlinTrk::IMarlinTrack::success ){
-
- std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not extrapolate TrackState to IP !!?? " << std::endl ;
+ //std::cout << " >>>>>>>>>>> PLCIOTrackConverter : could not extrapolate TrackState to IP !!?? " << std::endl ;
}
-
+ //std::cout << "D0 = " << tsLH.D0 << " phi = " << tsLH.phi << " omega = " << tsLH.omega << " Z0 = " << tsLH.Z0 << " tanLambda = " << tsLH.tanLambda << std::endl;
trk.addToTrackStates( tsIP ) ;
trk.addToTrackStates( tsFH ) ;
trk.addToTrackStates( tsLH ) ;
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
index e8aecc0c78ef9209c52fdffe5d1d2b9d9cd4bfbd..793cc8d1eb4772611963d00a90c324ae5454051e 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
@@ -156,8 +156,8 @@ namespace clupatra_new{
// FIXME debug Mingrui
//streamlog_out( DEBUG ) << " add hit << " << *first << " to layer " << (*first)->first->layer << std::endl ;
-
- hLV[ (*first)->first->layer ].push_back( *first ) ;
+ //std::cout << "fucd debug: " << " add hit << " << *first << " to layer " << (*first)->first->layer << std::endl ;
+ hLV[ (*first)->first->layer ].push_back( *first ) ;
++first ;
}
}
@@ -437,6 +437,7 @@ namespace clupatra_new{
// track state at last hit migyt be rubish....
// const lcio::TrackState* ts = ( outward ? trk->getTrackState( lcio::TrackState::AtLastHit ) : trk->getTrackState( lcio::TrackState::AtFirstHit ) ) ;
+ if(!hasTrackStateAt(trk, lcio::TrackState::AtFirstHit)) return false; //equivalent to pointer ts ==0 in lcio, fucd
const edm4hep::TrackState ts = getTrackStateAt(trk, lcio::TrackState::AtFirstHit ) ;
@@ -455,7 +456,7 @@ namespace clupatra_new{
int nHit = trk.trackerHits_size() ;
- if( nHit == 0 /*|| ts ==0*/ ) // FIXME Mingrui Since it is no pointer
+ if( nHit == 0 /*|| ts ==0*/ ) // FIXME Mingrui Since it is no pointer, fixed previous by hasTrackStateAt, fucd
return false ;
// float initial_chi2 = trk->getChi2() ;
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
index 223bc9379a8f928395babd7b2968b7fabbf0cbb7..38c54b400e37111151c69a5833e8bd5708abed06 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -92,8 +92,8 @@ std::string toString( int iTrk, edm4hep::ConstTrack tpcTrack, float bField=3.5 )
}
DECLARE_COMPONENT( FullLDCTrackingAlg )
-FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
- // _description = "Performs full tracking in ILD detector" ;
+FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+ // _description = "Performs full tracking in ILD detector" ;
_encoder = new UTIL::BitField64(lcio::ILDCellID0::encoder_string);
@@ -111,7 +111,7 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc
//declareProperty("VTXRawHits", _inVXDRawColHdl, "VXD SimTrackerHit collection");
// Input track collections
- declareProperty("TPCTracks", _TPCTrackColHdl, "TPC Track Collection");
+ declareProperty("TPCTracks", _TPCTrackColHdl, "TPC Track Collection");
declareProperty("SiTracks", _SiTrackColHdl, "Si Track Collection");
// Input relation collections
@@ -135,10 +135,10 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc
-StatusCode FullLDCTrackingAlg::initialize() {
+StatusCode FullLDCTrackingAlg::initialize() {
// usually a good idea to
- // printParameters();
+ // printParameters();
_nRun = -1 ;
_nEvt = 0 ;
PI = acos(-1.);
@@ -168,7 +168,7 @@ StatusCode FullLDCTrackingAlg::initialize() {
}
/*
-void FullLDCTrackingAlg::processRunHeader( LCRunHeader* run) {
+void FullLDCTrackingAlg::processRunHeader( LCRunHeader* run) {
_nRun++ ;
_nEvt = 0;
@@ -178,18 +178,18 @@ void FullLDCTrackingAlg::processRunHeader( LCRunHeader* run) {
}
*/
-StatusCode FullLDCTrackingAlg::execute() {
+StatusCode FullLDCTrackingAlg::execute() {
// debug() << endmsg;
- debug() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << " event = " << _nEvt << endmsg;
+ info() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << " event = " << _nEvt << endmsg;
// debug() << endmsg;
auto outCol = _OutputTrackColHdl.createAndPut();
prepareVectors();
debug() << "************************************PrepareVectors done..." << endmsg;
-
- debug() << "************************************Merge TPC/Si ..." << endmsg;
+ debug() << "************************************Merge TPC/Si ..." << endmsg;
+
MergeTPCandSiTracks();
debug() << "************************************Merging done ..." << endmsg;
@@ -206,20 +206,24 @@ StatusCode FullLDCTrackingAlg::execute() {
debug() << "************************************Trying non combined tracks ..." << endmsg;
AddNotCombinedTracks( );
debug() << "************************************Non combined tracks added ..." << endmsg;
- CheckTracks( );
+ //CheckTracks( );
debug() << "************************************Add Non assigned hits ..." << endmsg;
AddNotAssignedHits();
debug() << "************************************Non assigned hits added ..." << endmsg;
- AddTrackColToEvt(_trkImplVec, outCol);
- debug() << "Collections added to event ..." << endmsg;
+ try{
+ AddTrackColToEvt(_trkImplVec, outCol);
+ debug() << "Collections added to event ..." << endmsg;
+ }
+ catch(std::runtime_error& e){
+ error() << e.what() << std::endl;
+ }
CleanUp();
debug() << "Cleanup is done." << endmsg;
_nEvt++;
// getchar();
- // streamlog_out(DEBUG5) << endmsg;
- // streamlog_out(DEBUG5) << endmsg;
+ debug() << "FullLDCTrackingAlg execute() finished" << endmsg;
return StatusCode::SUCCESS;
}
@@ -230,7 +234,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
// if we want to point back to the hits we need to set the flag
//LCFlagImpl trkFlag(0) ;
//trkFlag.setBit( LCIO::TRBIT_HITS ) ;
- //colTRK->setFlag( trkFlag.getFlag() ) ;
+ //colTRK->setFlag( trkFlag.getFlag() ) ;
// LCCollectionVec * colRel = NULL;
@@ -244,7 +248,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
//SJA:FIXME: So here we are going to do one final refit. This can certainly be optimised, but rather than worry about the mememory management right now lets make it work, and optimise it later ...
-
+ debug() << "Total " << nTrkCand << " trkCand to deal" << endmsg;
for (int iTRK=0;iTRK<nTrkCand;++iTRK) {
TrackExtended * trkCand = trkVec[iTRK];
@@ -252,11 +256,11 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
std::vector<edm4hep::ConstTrackerHit> trkHits;
- debug() << " Trying to add track " << trkCand << " to final lcio collection " << endmsg;
+ //debug() << " Trying to add track " << trkCand << "(" << iTRK << ")" << " to final lcio collection " << endmsg;
int nHits = int(hitVec.size());
- debug() << " Trying to add track " << trkCand << " to final lcio collection nHits = " << nHits << endmsg;
+ debug() << " Trying to add track " << trkCand << "(" << iTRK << ")" <<" to final lcio collection nHits = " << nHits << endmsg;
for (int ihit=0;ihit<nHits;++ihit) {
@@ -268,10 +272,10 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
edm4hep::ConstTrackerHit trkHit = hitVec[ihit]->getTrackerHit();
if(trkHit.isAvailable()) {
- trkHits.push_back(trkHit);
+ trkHits.push_back(trkHit);
}
else{
- throw EVENT::Exception( std::string("FullLDCTrackingAlg::AddTrackColToEvt: TrackerHit pointer == NULL ") ) ;
+ throw std::runtime_error("FullLDCTrackingAlg::AddTrackColToEvt: TrackerHit pointer == NULL ");
}
}
@@ -279,7 +283,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
if( trkHits.size() < 3 ) {
debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Cannot fit less than 3 hits. Number of hits = " << trkHits.size() << endmsg;
- continue ;
+ continue ;
}
edm4hep::Track track;// = new edm4hep::Track;
@@ -322,12 +326,11 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
}
if(hasTrackStateAt(te->getTrack(), 3/*lcio::TrackState::AtLastHit*/)){
-
- debug() << "Initialise Fit with trackstate from last hit" << group << endmsg;
-
+ debug() << "Initialise Fit with trackstate from last hit " << group << " " << te << " " << te->getTrack().id() << endmsg;
+ if(te->getTrack().trackStates_size()>4) warning() << "TrackState size > 4" << endmsg;
ts_initial = getTrackStateAt(te->getTrack(), 3/*lcio::TrackState::AtLastHit*/);
-
prefit_set = true;
+ debug() << "ts_initial " << ts_initial << endmsg;
}
}
@@ -345,6 +348,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
ts_initial.referencePoint = ref;
ts_initial.location = 1/*lcio::TrackState::AtIP*/;
+ debug() << "ts_initial " << ts_initial << endmsg;
}
ts_initial.covMatrix = covMatrix;
@@ -367,41 +371,47 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
for (std::vector< std::pair<float, edm4hep::ConstTrackerHit> >::iterator it=r2_values.begin(); it!=r2_values.end(); ++it) {
trkHits.push_back(it->second);
}
-
+
+ //for(int iHit=0;iHit<trkHits.size();iHit++){
+ //const edm4hep::Vector3d& pos = trkHits[iHit].getPosition();
+ //debug() << "Hit " << iHit << ": r= " << sqrt(pos.x*pos.x+pos.y*pos.y) << " id = " << trkHits[iHit].id() << endmsg;
+ //}
+
+ debug() << "trkHits ready" << endmsg;
bool fit_backwards = IMarlinTrack::backward;
MarlinTrk::IMarlinTrack* marlinTrk = _trksystem->createTrack();
+ debug() << "marlinTrk created " << endmsg;
-
- int error = 0;
+ int error_code = 0;
try {
- error = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit);
+ error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit);
} catch (...) {
// delete track;
// delete marlinTrk;
-
- throw ;
+ error() << "exception happened while createFinalisedLCIOTrack!" << endmsg;
+ throw std::runtime_error("Need more check");
}
-
+ debug() << "createFinalisedLCIOTrack finished" << endmsg;
#ifdef MARLINTRK_DIAGNOSTICS_ON
- if ( error != IMarlinTrack::success && _runMarlinTrkDiagnostics ) {
+ if ( error_code != IMarlinTrack::success && _runMarlinTrkDiagnostics ) {
void * dcv = _trksystem->getDiagnositicsPointer();
DiagnosticsController* dc = static_cast<DiagnosticsController*>(dcv);
dc->skip_current_track();
- }
+ }
#endif
- std::vector<std::pair<edm4hep::ConstTrackerHit , double> > hits_in_fit ;
+ std::vector<std::pair<edm4hep::ConstTrackerHit , double> > hits_in_fit ;
std::vector<std::pair<edm4hep::ConstTrackerHit , double> > outliers ;
- std::vector<edm4hep::ConstTrackerHit> all_hits;
- all_hits.reserve(300);
+ std::vector<edm4hep::ConstTrackerHit> all_hits;
+ all_hits.reserve(hits_in_fit.size());
marlinTrk->getHitsInFit(hits_in_fit);
@@ -409,7 +419,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
all_hits.push_back(hits_in_fit[ihit].first);
}
- UTIL::BitField64 cellID_encoder( lcio::ILDCellID0::encoder_string ) ;
+ UTIL::BitField64 cellID_encoder( lcio::ILDCellID0::encoder_string ) ;
MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder);
@@ -424,15 +434,14 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
delete marlinTrk;
- if( error != IMarlinTrack::success ) {
- debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit failed with error code " << error << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
-
- //delete Track;
+ if( error_code != IMarlinTrack::success ) {
+ debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit failed with error code " << error_code << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
+ //delete Track;
continue ;
}
- if( track.getNdf() < 0) {
- debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
+ if( track.getNdf() < 0) {
+ debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg;
//delete Track;
continue ;
}
@@ -457,17 +466,20 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
// check if it is a tpc looper ...
if( UTIL::BitSet32( subTrack->getTrack().getType() )[ lcio::ILDDetID::TPC ] ) {
-
+ int nseg = subTrack->getTrack().tracks_size();
+ for(unsigned iSeg=0;iSeg<nseg;iSeg++){
+ track.addToTracks(subTrack->getTrack().getTracks(iSeg));
+ }
//const TrackVec segments = subTrack->getTrack().getTracks();
- std::vector<edm4hep::ConstTrack> segments;
- std::copy(subTrack->getTrack().tracks_begin(), subTrack->getTrack().tracks_end(), std::back_inserter(segments));
- if ( segments.empty() == false ) {
+ //std::vector<edm4hep::ConstTrack> segments;
+ //std::copy(subTrack->getTrack().tracks_begin(), subTrack->getTrack().tracks_end(), std::back_inserter(segments));
+ //if ( segments.empty() == false ) {
- for (unsigned iSeg=0;iSeg<segments.size();++iSeg) {
- track.addToTracks(segments[iSeg]);
- }
+ //for (unsigned iSeg=0;iSeg<segments.size();++iSeg) {
+ // track.addToTracks(segments[iSeg]);
+ //}
- }
+ //}
}
}
@@ -536,7 +548,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
eTot += trkP;
pxTot += trkPx;
pyTot += trkPy;
- pzTot += trkPz;
+ pzTot += trkPz;
nTotTracks++;
colTRK->push_back(track);
@@ -569,9 +581,8 @@ void FullLDCTrackingAlg::prepareVectors() {
_trkImplVec.clear();
_candidateCombinedTracks.clear();
-
std::map <edm4hep::ConstTrackerHit,TrackerHitExtended*> mapTrackerHits;
-
+
// Reading TPC hits
const edm4hep::TrackerHitCollection* hitTPCCol = nullptr;
try {
@@ -585,7 +596,7 @@ void FullLDCTrackingAlg::prepareVectors() {
debug() << "Number of TPC hits = " << nelem << endmsg;
for (edm4hep::ConstTrackerHit hit : *hitTPCCol) {
TrackerHitExtended * hitExt = new TrackerHitExtended(hit);
-
+ //info() << "TPC hit " << hit.id() << " " << hitExt << endmsg;
// Covariance Matrix in LCIO is defined in XYZ convert to R-Phi-Z
// For no error in r
@@ -609,7 +620,8 @@ void FullLDCTrackingAlg::prepareVectors() {
unsigned int layer = static_cast<unsigned int>(getLayerID(hit));
- debug() << " TPC Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << tpcRPhiRes << " dz " << tpcZRes << " layer = " << layer << endmsg;
+ debug() << " TPC Hit added : @ " << pos[0] << " " << pos[1] << " (r=" << sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << ") " << pos[2]
+ << " drphi " << tpcRPhiRes << " dz " << tpcZRes << " layer = " << layer << endmsg;
}
}
@@ -641,7 +653,7 @@ void FullLDCTrackingAlg::prepareVectors() {
hitExt->setResolutionZ(0.1);
// type and det are no longer used, set to INT_MAX to try and catch any missuse
- hitExt->setType(int(INT_MAX));
+ hitExt->setType(int(INT_MAX));
hitExt->setDet(int(INT_MAX));
_allFTDHits.push_back( hitExt );
@@ -673,7 +685,7 @@ void FullLDCTrackingAlg::prepareVectors() {
for (int i=0; i<3; ++i) {
pos[i] = hit.getPosition()[i];
- }
+ }
debug() << " FTD Pixel Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " layer = " << layer << endmsg;
@@ -749,7 +761,7 @@ void FullLDCTrackingAlg::prepareVectors() {
for (int i=0; i<3; ++i) {
pos[i] = hit.getPosition()[i];
- }
+ }
debug() << " FTD SpacePoint Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " layer = " << layer << endmsg;
@@ -810,15 +822,14 @@ void FullLDCTrackingAlg::prepareVectors() {
if ( UTIL::BitSet32( trkhit.getType() )[ UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ] ) {
fatal() << "FullLDCTrackingAlg: SIT Hit cannot be of type UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL COMPOSITE SPACEPOINTS must be use instead. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg;
exit(1);
- }
+ }
// most likely case: COMPOSITE_SPACEPOINT hits formed from stereo strip hits
else if ( UTIL::BitSet32( trkhit.getType() )[ UTIL::ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT ] ) {
// SJA:FIXME: fudge for now by a factor of two and ignore covariance
- drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]);
- dz = sqrt(trkhit.getCovMatrix()[5]);
-
- }
+ drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]);
+ dz = sqrt(trkhit.getCovMatrix()[5]);
+ }
// or a PIXEL based SIT, using 2D TrackerHitPlane like the VXD above
else if ( UTIL::BitSet32( trkhit.getType() )[3]) {
// FIXME Should make it correct
@@ -842,10 +853,10 @@ void FullLDCTrackingAlg::prepareVectors() {
// FIXME should make it correct
// drphi = trkhit_P->getdU();
- // dz = trkhit_P->getdV();
+ // dz = trkhit_P->getdV();
drphi = trkhit.getCovMatrix()[2];
dz = trkhit.getCovMatrix()[5];
- }
+ }
// or a simple cylindrical design, as used in the LOI
/* FIXME, fucd
else if ( true ) {
@@ -854,16 +865,15 @@ void FullLDCTrackingAlg::prepareVectors() {
// drphi = trkhit_C->getdRPhi();
// dz = trkhit_C->getdZ();
drphi = 1.0;
- dz = 1.0;
- }
+ dz = 1.0;
+ }
*/
// this would be very unlikely, but who knows ... just an ordinary TrackerHit, which is not a COMPOSITE_SPACEPOINT
else {
// SJA:FIXME: fudge for now by a factor of two and ignore covariance
- drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]);
- dz = sqrt(trkhit.getCovMatrix()[5]);
-
+ drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]);
+ dz = sqrt(trkhit.getCovMatrix()[5]);
}
// now that the hit type has been established carry on and create a
@@ -1052,7 +1062,7 @@ void FullLDCTrackingAlg::prepareVectors() {
hitExt->setResolutionZ(trkhit.getCovMatrix()[5]);
// type and det are no longer used, set to INT_MAX to try and catch any missuse
- hitExt->setType(int(INT_MAX));
+ hitExt->setType(int(INT_MAX));
hitExt->setDet(int(INT_MAX));
_allVTXHits.push_back( hitExt );
mapTrackerHits[trkhit] = hitExt;
@@ -1084,26 +1094,27 @@ void FullLDCTrackingAlg::prepareVectors() {
int iTrk = -1;
for(edm4hep::ConstTrack tpcTrack : *tpcTrackCol){
iTrk++;
- ConstTrackerHitVec hitVec(tpcTrack.trackerHits_begin(), tpcTrack.trackerHits_end());
- int nHits = int(hitVec.size());
+ //ConstTrackerHitVec hitVec(tpcTrack.trackerHits_begin(), tpcTrack.trackerHits_end());
+ int nHits = tpcTrack.trackerHits_size();//int(hitVec.size());
debug() << toString( iTrk, tpcTrack , _bField ) << endmsg;
-
+ //edm4hep::TrackState ts = getTrackStateAt(tpcTrack, 3);
TrackExtended * trackExt = new TrackExtended( tpcTrack );
-
+ //debug() << "TPC TrackExtended " << iTrk << " = " << trackExt << " " << tpcTrack.id() << " D0=" << getD0(tpcTrack) << " phi=" << getPhi(tpcTrack) << " " << ts.phi << endmsg;
trackExt->setOmega(getOmega(tpcTrack));
trackExt->setTanLambda(getTanLambda(tpcTrack));
trackExt->setPhi(getPhi(tpcTrack));
trackExt->setD0(getD0(tpcTrack));
trackExt->setZ0(getZ0(tpcTrack));
+
float cov[15];
- float param[5];
+ //float param[5];
// float reso[4];
- param[0] = getOmega(tpcTrack);
- param[1] = getTanLambda(tpcTrack);
- param[2] = getPhi(tpcTrack);
- param[3] = getD0(tpcTrack);
- param[4] = getZ0(tpcTrack);
+ //param[0] = getOmega(tpcTrack);
+ //param[1] = getTanLambda(tpcTrack);
+ //param[2] = getPhi(tpcTrack);
+ //param[3] = getD0(tpcTrack);
+ //param[4] = getZ0(tpcTrack);
std::array<float, 15> Cov = getCovMatrix(tpcTrack);
int NC = int(Cov.size());
@@ -1113,17 +1124,23 @@ void FullLDCTrackingAlg::prepareVectors() {
trackExt->setCovMatrix(cov);
trackExt->setNDF(tpcTrack.getNdf());
- trackExt->setChi2(tpcTrack.getChi2());
+ trackExt->setChi2(tpcTrack.getChi2());
for (int iHit=0;iHit<nHits;++iHit) {
- edm4hep::ConstTrackerHit hit = hitVec[iHit];
- TrackerHitExtended * hitExt = mapTrackerHits[hit];
+ edm4hep::ConstTrackerHit hit = tpcTrack.getTrackerHits(iHit);//hitVec[iHit];
+ if(!hit.isAvailable()) error() << "Tracker hit not available" << endmsg;
+ //info() << "hit " << hit.id() << " " << hit.getCellID() << " " << hit.getPosition()[0] << " " << hit.getPosition()[1] << " " << hit.getPosition()[2] << endmsg;
+ auto it = mapTrackerHits.find(hit);
+ if(it==mapTrackerHits.end()) error() << "Cannot find hit " << hit.id() << " in map" << endmsg;
+ else continue;
+ TrackerHitExtended * hitExt = it->second;
+ //info() << hit.id() << " " << hitExt << endmsg;
hitExt->setTrackExtended( trackExt );
- trackExt->addTrackerHitExtended( hitExt );
- }
- _allTPCTracks.push_back( trackExt );
- }
+ trackExt->addTrackerHitExtended( hitExt );
+ }
+ _allTPCTracks.push_back( trackExt );
+ }
}
-
+ //checkTrackState(3);
// Reading Si Tracks
const edm4hep::TrackCollection* siTrackCol = nullptr;
try {
@@ -1141,43 +1158,44 @@ void FullLDCTrackingAlg::prepareVectors() {
iTrk++;
double prob = ( siTrack.getNdf() > 0 ? gsl_cdf_chisq_Q( siTrack.getChi2() , (double) siTrack.getNdf() ) : 0. ) ;
if( prob < _minChi2ProbForSiliconTracks ) {
- debug() << "Si Tracks " << siTrack << " id : " << siTrack.id() << " rejected with prob " << prob << " < " << _minChi2ProbForSiliconTracks << endmsg;
+ debug() << "Si Tracks " << iTrk << " id : " << siTrack.id() << " rejected with prob " << prob << " < " << _minChi2ProbForSiliconTracks << endmsg;
continue;
}
TrackExtended * trackExt = new TrackExtended( siTrack );
- ConstTrackerHitVec hitVec(siTrack.trackerHits_begin(), siTrack.trackerHits_end());
- int nHits = int(hitVec.size());
+ //info() << "Si TrackExtended " << iTrk << " = " << trackExt << endmsg;
+ //ConstTrackerHitVec hitVec(siTrack.trackerHits_begin(), siTrack.trackerHits_end());
+ int nHits = siTrack.trackerHits_size();//int(hitVec.size());
trackExt->setOmega(getOmega(siTrack));
trackExt->setTanLambda(getTanLambda(siTrack));
trackExt->setPhi(getPhi(siTrack));
trackExt->setD0(getD0(siTrack));
trackExt->setZ0(getZ0(siTrack));
float cov[15];
- float param[5];
+ //float param[5];
- param[0] = getOmega(siTrack);
- param[1] = getTanLambda(siTrack);
- param[2] = getPhi(siTrack);
- param[3] = getD0(siTrack);
- param[4] = getZ0(siTrack);
+ //param[0] = getOmega(siTrack);
+ //param[1] = getTanLambda(siTrack);
+ //param[2] = getPhi(siTrack);
+ //param[3] = getD0(siTrack);
+ //param[4] = getZ0(siTrack);
std::array<float, 15> Cov = getCovMatrix(siTrack);
int NC = int(Cov.size());
for (int ic=0;ic<NC;ic++) {
cov[ic] = Cov[ic];
- }
+ }
trackExt->setCovMatrix(cov);
trackExt->setNDF(siTrack.getNdf());
- trackExt->setChi2(siTrack.getChi2());
+ trackExt->setChi2(siTrack.getChi2());
char strg[200];
HelixClass helixSi;
for (int iHit=0;iHit<nHits;++iHit) {
- edm4hep::ConstTrackerHit hit = hitVec[iHit];
+ edm4hep::ConstTrackerHit hit = siTrack.getTrackerHits(iHit);//hitVec[iHit];
TrackerHitExtended * hitExt = mapTrackerHits[hit];
hitExt->setTrackExtended( trackExt );
- trackExt->addTrackerHitExtended( hitExt );
+ trackExt->addTrackerHitExtended( hitExt );
}
float d0Si = trackExt->getD0();
@@ -1207,7 +1225,7 @@ void FullLDCTrackingAlg::prepareVectors() {
void FullLDCTrackingAlg::CleanUp(){
- int nNonCombTpc = int(_allNonCombinedTPCTracks.size());
+ int nNonCombTpc = int(_allNonCombinedTPCTracks.size());
for (int i=0;i<nNonCombTpc;++i) {
TrackExtended * trkExt = _allNonCombinedTPCTracks[i];
GroupTracks * group = trkExt->getGroupTracks();
@@ -1215,7 +1233,7 @@ void FullLDCTrackingAlg::CleanUp(){
}
_allNonCombinedTPCTracks.clear();
- int nNonCombSi = int(_allNonCombinedSiTracks.size());
+ int nNonCombSi = int(_allNonCombinedSiTracks.size());
for (int i=0;i<nNonCombSi;++i) {
TrackExtended * trkExt = _allNonCombinedSiTracks[i];
GroupTracks * group = trkExt->getGroupTracks();
@@ -1284,15 +1302,16 @@ void FullLDCTrackingAlg::CleanUp(){
TrackExtended * trkExt = _allCombinedTracks[i];
GroupTracks * group = trkExt->getGroupTracks();
delete group;
- delete trkExt;
+ delete trkExt;
}
_allCombinedTracks.clear();
- // int nImplTrk = int(_trkImplVec.size());
- // for (int i=0;i<nImplTrk;++i) {
- // TrackExtended * trkImpl = _trkImplVec[i];
- // delete trkImpl;
- // }
+ // fucd: why comment in Marlin? if open, crash while save collection
+ //int nImplTrk = int(_trkImplVec.size());
+ //for (int i=0;i<nImplTrk;++i) {
+ // TrackExtended * trkImpl = _trkImplVec[i];
+ // delete trkImpl;
+ //}
_trkImplVec.clear();
//AS: Dont delete the individual entries, some of them are cleared elsewhere, I think
@@ -1483,28 +1502,28 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
for (int ih=0;ih<nSiHits;++ih) {
edm4hep::ConstTrackerHit trkHit = siHitVec[ih]->getTrackerHit();
- if(trkHit.isAvailable()) {
- trkHits.push_back(trkHit);
+ if(trkHit.isAvailable()) {
+ trkHits.push_back(trkHit);
}
else{
throw EVENT::Exception( std::string("FullLDCTrackingAlg::CombineTracks: TrackerHit pointer == NULL ") ) ;
}
- }
+ }
for (int ih=0;ih<nTPCHits;++ih) {
edm4hep::ConstTrackerHit trkHit = tpcHitVec[ih]->getTrackerHit();
- if(trkHit.isAvailable()) {
- trkHits.push_back(trkHit);
+ if(trkHit.isAvailable()) {
+ trkHits.push_back(trkHit);
}
else{
throw EVENT::Exception( std::string("FullLDCTrackingAlg::CombineTracks: TrackerHit pointer == NULL ") ) ;
}
- }
+ }
double chi2_D;
int ndf;
- if( trkHits.size() < 3 ) {
+ if( trkHits.size() < 3 ) {
return 0 ;
@@ -1543,7 +1562,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
pre_fit = getTrackStateAt(tpcTrack->getTrack(), 3/*lcio::TrackState::AtLastHit*/);
}
catch(std::runtime_error& e){
- error() << e.what() << " should be checked (TPC track)" << endmsg;
+ error() << e.what() << " should be checked (TPC track) " << tpcTrack << endmsg;
try{
pre_fit = getTrackStateAt(siTrack->getTrack(), 3/*lcio::TrackState::AtLastHit*/);
}
@@ -1573,14 +1592,14 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
debug() << "FullLDCTrackingAlg::CombineTracks: creation of fit fails with error " << errorCode << endmsg;
return 0;
}
-
+ debug() << "createFit finished" << endmsg;
edm4hep::Vector3d point(0.,0.,0.); // nominal IP
edm4hep::TrackState trkState ;
errorCode = marlin_trk.propagate(point, trkState, chi2_D, ndf ) ;
if ( errorCode != IMarlinTrack::success ) {
- debug() << "FullLDCTrackingAlg::CombineTracks: propagate to IP fails with error " << errorCode << endmsg;
+ debug() << "FullLDCTrackingAlg::CombineTracks: propagate to IP fails with error " << errorCode << endmsg;
return 0;
}
@@ -1636,7 +1655,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
edm4hep::ConstTrackerHit rawHit = Navigation::Instance()->GetTrackerHit(hit.getRawHits(ihit));
hits.push_back(rawHit);
}
- else debug() << "not space point, id=" << hit.id() << endmsg;
+ else debug() << "not space point, id=" << hit.id() << endmsg;
}
catch(std::runtime_error& e){
debug() << e.what() << " raw hit id = " << hit.getRawHits(ihit) << " for TrackerHit id = " << hit.id() << endmsg;
@@ -1719,7 +1738,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
OutputTrack->setD0(d0);
OutputTrack->setChi2(chi2_D);
OutputTrack->setNDF(ndf);
-
+ debug() << "create TrackExtended " << OutputTrack << endmsg;
float cov[15];
for (int i = 0 ; i<15 ; ++i) {
@@ -1781,7 +1800,6 @@ void FullLDCTrackingAlg::SortingTrackHitPairs(TrackHitPairVec & trackHitPairVec)
}
}
-
}
/*
@@ -1806,7 +1824,7 @@ void FullLDCTrackingAlg::Sorting(TrackExtendedVec & trackVec) {
trackVec[j] = trackVec[j+1];
trackVec[j+1] = Temp;
}
- }
+ }
}
/*
@@ -1849,7 +1867,7 @@ void FullLDCTrackingAlg::SelectCombinedTracks() {
// associate the current group to the two sub tracks
firstTrack->setGroupTracks(group);
- secondTrack->setGroupTracks(group);
+ secondTrack->setGroupTracks(group);
// get the tracker hits ...
TrackerHitExtendedVec firstVec = firstTrack->getTrackerHitExtendedVec();
@@ -1872,7 +1890,7 @@ void FullLDCTrackingAlg::SelectCombinedTracks() {
if (zpos>edges[1])
edges[1] = zpos;
if (zpos<edges[0])
- edges[0] = zpos;
+ edges[0] = zpos;
}
// get min and max z for the second sub track
@@ -1896,7 +1914,7 @@ void FullLDCTrackingAlg::SelectCombinedTracks() {
int iopt = 1;
// here it is assumed that the tpc tracks is the secondTrack ...
// PrintOutMerging(secondTrack,firstTrack,iopt);
- }
+ }
}
}
else { // if(nTracks>2)
@@ -1905,7 +1923,7 @@ void FullLDCTrackingAlg::SelectCombinedTracks() {
}
}
-void FullLDCTrackingAlg::AddNotCombinedTracks() {
+void FullLDCTrackingAlg::AddNotCombinedTracks() {
int nTPCTrk = int(_allTPCTracks.size());
int nSiTrk = int(_allSiTracks.size());
@@ -1915,7 +1933,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
allMergedTracks.clear();
// forcing merging of Si and TPC track segments
- if (_forceMerging==1) {
+ if (_forceMerging==1) {
// loop over all TPC tracks
for (int i=0;i<nTPCTrk;++i) {
@@ -1957,14 +1975,15 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
float chi2O = dOmega/_dOmegaForForcedMerging;
float chi2A = angle/_angleForForcedMerging;
- float deltaP = chi2O*chi2O + chi2A*chi2A;
+ float deltaP = chi2O*chi2O + chi2A*chi2A;
// if this is the best match (diffMin) set the possible merger
if (deltaP<diffMin) {
diffMin = deltaP;
siTrkToAttach = trkExtSi;
}
- } else {
+ }
+ else {
if (_debug >= 3) {
int iopt = 7;
debug() << significance << " " << angleSignificance << endmsg;
@@ -1985,20 +2004,20 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
OutputTrack->setGroupTracks(group);
// trkExtSi->setGroupTracks(group);
- // trkExtTPC->setGroupTracks(group);
+ // trkExtTPC->setGroupTracks(group);
OutputTrack->setOmega(trkExtTPC->getOmega());
OutputTrack->setTanLambda(trkExtSi->getTanLambda());
OutputTrack->setPhi(trkExtSi->getPhi());
OutputTrack->setZ0(trkExtSi->getZ0());
OutputTrack->setD0(trkExtSi->getD0());
-
+ debug() << "create TrackExtended " << OutputTrack << endmsg;
float covMatTPC[15];
float covMatSi[15];
float covMat[15];
for (int iCov=0;iCov<15;++iCov) {
covMatTPC[iCov] = trkExtTPC->getCovMatrix()[iCov];
- covMatSi[iCov] = trkExtSi->getCovMatrix()[iCov];
+ covMatSi[iCov] = trkExtSi->getCovMatrix()[iCov];
covMat[iCov] = covMatSi[iCov];
}
@@ -2011,10 +2030,10 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
OutputTrack->setCovMatrix(covMat);
TrackerHitExtendedVec tpcHitVec = trkExtTPC->getTrackerHitExtendedVec();
- TrackerHitExtendedVec siHitVec = trkExtSi->getTrackerHitExtendedVec();
+ TrackerHitExtendedVec siHitVec = trkExtSi->getTrackerHitExtendedVec();
int nTPCHits = int( tpcHitVec.size());
- int nSiHits = int( siHitVec.size());
+ int nSiHits = int( siHitVec.size());
float edges[2];
edges[0] = 1.0e+20;
@@ -2031,11 +2050,11 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
edges[0] = zpos;
if (zpos>edges[1])
edges[1] = zpos;
- }
+ }
for (int iH=0;iH<nTPCHits;++iH) {
TrackerHitExtended * hitExt = tpcHitVec[iH];
OutputTrack->addTrackerHitExtended(hitExt);
- hitExt->setUsedInFit(true);
+ hitExt->setUsedInFit(true);
edm4hep::ConstTrackerHit hit = hitExt->getTrackerHit();
float zpos = float(hit.getPosition()[2]);
if (zpos<edges[0])
@@ -2243,8 +2262,8 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
for (int iTrk=0;iTrk<nTrk;++iTrk) {
TrackExtended * trkInGroup = segVec[iTrk];
int iComp = 1;
- float dPt = CompareTrk(trkExt,trkInGroup,_d0CutToMergeTPC,_z0CutToMergeTPC,iComp);
- if (dPt >= dPtMin) {
+ float dPt = CompareTrk(trkExt,trkInGroup,_d0CutToMergeTPC,_z0CutToMergeTPC,iComp);
+ if (dPt >= dPtMin) {
// PrintOutMerging(trkExt,trkInGroup,iopt);
}
}
@@ -2384,7 +2403,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
}
if (CombTrkToAttach != NULL) { // attach TPC segment to existing Comb Track
- GroupTracks * groupToAttach = CombTrkToAttach->getGroupTracks();
+ GroupTracks * groupToAttach = CombTrkToAttach->getGroupTracks();
TrackExtended * SiCombTrk = groupToAttach->getTrackExtendedVec()[0];
TrackExtended * TpcCombTrk = groupToAttach->getTrackExtendedVec()[1];
@@ -2392,7 +2411,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
int iopt = 4;
// PrintOutMerging(keyTrack,SiCombTrk,iopt);
iopt = 5;
- // PrintOutMerging(keyTrack,TpcCombTrk,iopt);
+ // PrintOutMerging(keyTrack,TpcCombTrk,iopt);
}
for (int iTrk=0;iTrk<nTrk;iTrk++) {
@@ -2441,16 +2460,16 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
// get the lcio track which is behind this segemnt
edm4hep::ConstTrack track = segment->getTrack();
- ConstTrackerHitVec hitVec(track.trackerHits_begin(), track.trackerHits_end());
+ //ConstTrackerHitVec hitVec(track.trackerHits_begin(), track.trackerHits_end());
debug() << "Group of orphaned TPC tracks: trying track " << track.id() << endmsg;
- int nHits = int(hitVec.size());
+ int nHits = track.trackerHits_size();//int(hitVec.size());
// loop over it's hits
for (int iH=0;iH<nHits;++iH) {
- float zPosi = fabs(hitVec[iH].getPosition()[2]);
+ float zPosi = fabs(track.getTrackerHits(iH).getPosition()[2]);//fabs(hitVec[iH].getPosition()[2]);
// if this segment has the hit closest to the IP so far
if (zPosi<zMin) {
@@ -2478,7 +2497,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
TrackExtended * segment = segVec[iTrk];
TrackerHitExtendedVec hitVecS = segment->getTrackerHitExtendedVec();
- int nHitS = int(hitVecS.size());
+ int nHitS = int(hitVecS.size());
// loop over the hits for the current segment
for (int iH=0;iH<nHitS;++iH) {
@@ -2489,7 +2508,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
segment->setGroupTracks( newGroup );
newGroup->addTrackExtended( segment );
trkHitExt->setUsedInFit( false );
- chosenTrack->addTrackerHitExtended( trkHitExt );
+ chosenTrack->addTrackerHitExtended( trkHitExt );
}
else {
trkHitExt->setUsedInFit( true );
@@ -2531,7 +2550,7 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
trkExt->setGroupTracks( newGrp );
}
- }
+ }
}
for (int i=0;i<nSiTrk;++i) { // adding left-over Si segments to the list of tracks
@@ -2554,14 +2573,14 @@ void FullLDCTrackingAlg::AddNotCombinedTracks() {
GroupTracks * newGrp = new GroupTracks();
newGrp->addTrackExtended( trkExt );
- trkExt->setGroupTracks( newGrp );
+ trkExt->setGroupTracks( newGrp );
}
}
}
-void FullLDCTrackingAlg::CheckTracks() {
+void FullLDCTrackingAlg::CheckTracks() {
for(unsigned int i = 0; i< _trkImplVec.size();i++){
TrackExtended *first = _trkImplVec[i];
@@ -2688,7 +2707,7 @@ float FullLDCTrackingAlg::CompareTrkII(TrackExtended * first, TrackExtended * se
float result = 1.0e+20;
- Angle = 1.0e+20;
+ Angle = 1.0e+20;
float omegaFirst = first->getOmega();
float omegaSecond = second->getOmega();
float deltaOmega = fabs((omegaFirst-omegaSecond)/omegaSecond);
@@ -2955,7 +2974,7 @@ float FullLDCTrackingAlg::CompareTrk(TrackExtended * first, TrackExtended * seco
if (momSecond<momFirst)
den = momSecond;
- result = nom/den;
+ result = nom/den;
}
else {
@@ -3005,7 +3024,7 @@ float FullLDCTrackingAlg::CompareTrk(TrackExtended * first, TrackExtended * seco
TrackerHitExtendedVec hitVec = trkGrp->getTrackerHitExtendedVec();
for(unsigned int i =0; i<hitVec.size(); ++i){
- hitvecFirst.push_back(hitVec[i]->getTrackerHit());
+ hitvecFirst.push_back(hitVec[i]->getTrackerHit());
}
}
}
@@ -3126,7 +3145,7 @@ float FullLDCTrackingAlg::CompareTrk(TrackExtended * first, TrackExtended * seco
}
delete combinedTrack->getGroupTracks();
delete combinedTrack;
- }
+ }
else {
//std::cout << "Could not combine track " << endmsg;
// if close in pt and the fraction of hits matching the helix extraolations is greater than _maxFractionOfOutliersCutHighPtMerge
@@ -3146,7 +3165,7 @@ float FullLDCTrackingAlg::CompareTrk(TrackExtended * first, TrackExtended * seco
(fcloseSecond > _maxFractionOfOutliersCutHighPtMerge || fcloseFirst > _maxFractionOfOutliersCutHighPtMerge)
){
- split = true;
+ split = true;
}
@@ -3176,7 +3195,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
// if (_assignSETHits>0) { // Assignment of SET Hits
//
- // const gear::GearParameters& pSETDet = Global::GEAR->getGearParameters("SET");
+ // const gear::GearParameters& pSETDet = Global::GEAR->getGearParameters("SET");
// int nLayersSET = int(pSETDet.getDoubleVals("SETLayerRadius").size());
//
// int nSETHits = _allSETHits.size();
@@ -3188,7 +3207,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
// edm4hep::ConstTrackerHit hit = trkHit->getTrackerHit();
// int layer = getLayerID(trkHit);
// if (layer>=0&&layer<nLayersSET)
- // SETHits[layer].push_back(trkHit);
+ // SETHits[layer].push_back(trkHit);
// }
// for (int iL=0; iL<nLayersSET; ++iL) { // loop over SET layers
// TrackerHitExtendedVec hitVec = SETHits[iL];
@@ -3199,7 +3218,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
// if (_assignETDHits>0) { // Assignment of ETD Hits
//
- // const gear::GearParameters& pETDDet = Global::GEAR->getGearParameters("ETD");
+ // const gear::GearParameters& pETDDet = Global::GEAR->getGearParameters("ETD");
// int nLayersETD = int(pETDDet.getDoubleVals("ETDLayerZ").size());
//
// int nETDHits = _allETDHits.size();
@@ -3261,10 +3280,10 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
if (_assignSITHits>0) { // Treatment of left-over SIT hits
debug() << "Assign SIT hits *********************************" << endmsg;
- std::vector<TrackerHitExtendedVec> nonAssignedSITHits;
+ std::vector<TrackerHitExtendedVec> nonAssignedSITHits;
nonAssignedSITHits.resize(_nLayersSIT);
- int nSITHits = int(_allSITHits.size());
+ int nSITHits = int(_allSITHits.size());
// loop over all SIT hits ...
for (int iH=0;iH<nSITHits;++iH) {
@@ -3282,7 +3301,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
nonAssignedSITHits[layer].push_back(trkHitExt);
}
}
- }
+ }
for (int iL=_nLayersSIT-1;iL>=0;--iL) { // reverse loop over layers in Si
@@ -3341,8 +3360,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
TrackerHitExtendedVec hitVec = nonAssignedFTDHits[iL];
debug() << "AddNotAssignedHits : Try to assign hits from layer " << iL << " : Number of hits = " << hitVec.size() << endmsg;
- AssignSiHitsToTracks(hitVec, _distCutForFTDHits);
-
+ AssignSiHitsToTracks(hitVec, _distCutForFTDHits);
}
}
}
@@ -3377,7 +3395,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
if ( nonAssignedVTXHits[iL].empty() == false ) {
TrackerHitExtendedVec hitVec = nonAssignedVTXHits[iL];
debug() << "AddNotAssignedHits : Try to assign hits from layer " << iL << " : Number of hits = " << hitVec.size() << endmsg;
- AssignSiHitsToTracks(hitVec, _distCutForVTXHits);
+ AssignSiHitsToTracks(hitVec, _distCutForVTXHits);
}
}
}
@@ -3422,7 +3440,7 @@ void FullLDCTrackingAlg::CleanUpExtrapolations() {
TrackExtended * trk = _trkImplVec[iTrk];
HelixClass * helix = _trackExtrapolatedHelix[trk];
delete helix;
- }
+ }
_trackExtrapolatedHelix.clear();
}
@@ -3462,7 +3480,7 @@ void FullLDCTrackingAlg::AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, f
for (int iT=0;iT<nTrk;++iT) {
TrackExtended * trkExt = _trkImplVec[iT];
- float tanLambda = trkExt->getTanLambda();
+ float tanLambda = trkExt->getTanLambda();
float product = pos[2]*tanLambda;
// check that the hit and track are in the same z-half, which won't work for the rare cases of something going backwards ...
@@ -3609,7 +3627,7 @@ void FullLDCTrackingAlg::AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, f
pre_fit.referencePoint = ref;
- pre_fit.location = 1/*lcio::TrackState::AtIP*/;
+ pre_fit.location = 1/*lcio::TrackState::AtIP*/;
// setup initial dummy covariance matrix
std::array<float,15> covMatrix;
@@ -3799,11 +3817,11 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
TrackerHitExtended * trkHitExt = hitVecGrp[iH];
float pos[3] = {0.,0.,0.};
for (int iC=0;iC<3;++iC)
- pos[iC] = float(trkHitExt->getTrackerHit().getPosition()[iC]);
+ pos[iC] = float(trkHitExt->getTrackerHit().getPosition()[iC]);
if (pos[2]>zMax) {
zMax = pos[2];
for (int iC=0;iC<3;++iC)
- endPoint[iC] = pos[iC];
+ endPoint[iC] = pos[iC];
}
if (pos[2]<zMin) {
zMin = pos[2];
@@ -3832,7 +3850,7 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
HitPositions[iH].push_back(float(temppos[1]));
HitPositions[iH].push_back(float(temppos[2]));
HitSign[iH]=std::signbit(temppos[2]);
- }
+ }
debug() << "AssignTPCHitsToTracks: Starting loop " << nTrk << " tracks and " << nHits << " hits" << endmsg;
@@ -3984,7 +4002,7 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec,
TrackExtended * trkExt = _allNonCombinedTPCTracks[iT];
- float tanLambda = trkExt->getTanLambda();
+ float tanLambda = trkExt->getTanLambda();
float product = pos[2]*tanLambda;
debug() << "AssignSiHitsToTracks : product = " << product << " z hit = " << pos[2] << endmsg;
@@ -4030,8 +4048,7 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec,
trkHitPair->getTrackerHitExtended();
- if (flagTrack[trkExt] && flagHit[trkHitExt]) {
-
+ if (flagTrack[trkExt] && flagHit[trkHitExt]) {
// get the hits already assigned to the track
TrackerHitExtendedVec hitsInTrack = trkExt->getTrackerHitExtendedVec();
@@ -4044,7 +4061,7 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec,
// count the number of hits used in the fit
if (hitInTrack->getUsedInFit()) {
- nHitsInFit++;
+ nHitsInFit++;
}
}
@@ -4060,8 +4077,8 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec,
if (hitInTrack->getUsedInFit()) {
edm4hep::ConstTrackerHit hit = hitInTrack->getTrackerHit();
iHitInFit++;
- if(hit.isAvailable()) {
- trkHits.push_back(hit);
+ if(hit.isAvailable()) {
+ trkHits.push_back(hit);
}
else{
throw EVENT::Exception( std::string("FullLDCTrackingAlg::AssignSiHitsToTracks: TrackerHit pointer == NULL ") ) ;
@@ -4262,7 +4279,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
}
else if (iopt==8) {
secondColName = _SiTrackMCPCollName;
- }
+ }
else {
secondColName = _TPCTrackMCPCollName;
}
@@ -4320,19 +4337,19 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
float z0Second = secondTrackExt->getZ0();
float omegaSecond = secondTrackExt->getOmega();
float tanLSecond = secondTrackExt->getTanLambda();
- float phi0Second = secondTrackExt->getPhi();
+ float phi0Second = secondTrackExt->getPhi();
HelixClass firstHelix;
firstHelix.Initialize_Canonical(phi0First,d0First,z0First,omegaFirst,tanLFirst,_bField);
float pxFirst = firstHelix.getMomentum()[0];
float pyFirst = firstHelix.getMomentum()[1];
- float pzFirst = firstHelix.getMomentum()[2];
+ float pzFirst = firstHelix.getMomentum()[2];
HelixClass secondHelix;
secondHelix.Initialize_Canonical(phi0Second,d0Second,z0Second,omegaSecond,tanLSecond,_bField);
float pxSecond = secondHelix.getMomentum()[0];
float pySecond = secondHelix.getMomentum()[1];
- float pzSecond = secondHelix.getMomentum()[2];
+ float pzSecond = secondHelix.getMomentum()[2];
// get the momentum differences
@@ -4369,7 +4386,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
den = pSecond;
dPplus = dPplus/den;
- dPminus = dPminus/den;
+ dPminus = dPminus/den;
// now check if this was a Erroneous merger ...
if (firstMCP!=secondMCP && iopt < 6) {
@@ -4481,7 +4498,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
streamlog_out(DEBUG4) << endmsg;
- }
+ }
// ... or if it was an incorrect TPC to Si rejection ...
else if (firstMCP==secondMCP && ( (iopt == 6) || (iopt == 7) ) ) {
@@ -4528,7 +4545,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
sprintf(strg," %5.3f\n",secondWght);
streamlog_out(DEBUG4) << strg;
- streamlog_out(DEBUG4) << "Difference in +P = " << dPplus << " -P = " << dPminus << endmsg;
+ streamlog_out(DEBUG4) << "Difference in +P = " << dPplus << " -P = " << dPminus << endmsg;
TrackExtended * combinedTrack = CombineTracks(firstTrackExt,secondTrackExt, _maxAllowedPercentageOfOutliersForTrackCombination, true);
@@ -4540,7 +4557,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
streamlog_out(DEBUG4) << "Could not combine track " << endmsg;
}
streamlog_out(DEBUG4) << " Overlap = " << SegmentRadialOverlap(firstTrackExt, secondTrackExt) << " veto = " << VetoMerge(firstTrackExt, secondTrackExt) << endmsg;
- streamlog_out(DEBUG4) << endmsg;
+ streamlog_out(DEBUG4) << endmsg;
}
// ... or if it was an correct merger ...
@@ -4598,7 +4615,7 @@ void FullLDCTrackingAlg::PrintOutMerging(TrackExtended * firstTrackExt, TrackExt
catch( ... ){};
-}
+}
*/
/*
@@ -4612,7 +4629,7 @@ void FullLDCTrackingAlg::GeneralSorting(int * index, float * val, int direct, in
}
for (int i = 0 ; i < nVal-1; i++) {
- for (int j = 0; j < nVal-i-1; j++) {
+ for (int j = 0; j < nVal-i-1; j++) {
valOne = val[j];
valTwo = val[j+1];
bool order = valOne > valTwo;
@@ -4627,7 +4644,7 @@ void FullLDCTrackingAlg::GeneralSorting(int * index, float * val, int direct, in
index[j] = index[j+1];
index[j+1] = indTemp;
}
- }
+ }
}
@@ -4655,7 +4672,7 @@ int FullLDCTrackingAlg::SegmentRadialOverlap(TrackExtended* first, TrackExtended
TrackerHitExtendedVec hitVec = trkGrp->getTrackerHitExtendedVec();
for(unsigned int i =0; i<hitVec.size(); ++i){
- hitvecFirst.push_back(hitVec[i]->getTrackerHit());
+ hitvecFirst.push_back(hitVec[i]->getTrackerHit());
}
}
}
@@ -4690,7 +4707,7 @@ int FullLDCTrackingAlg::SegmentRadialOverlap(TrackExtended* first, TrackExtended
float rj = sqrt(xj*xj+yj*yj);
if(fabs(ri-rj)<_tpc_pad_height/2.0)count++;
}
- }
+ }
return count;
}
@@ -4717,19 +4734,19 @@ bool FullLDCTrackingAlg::VetoMerge(TrackExtended* firstTrackExt, TrackExtended*
const float z0Second = secondTrackExt->getZ0();
const float omegaSecond = secondTrackExt->getOmega();
const float tanLSecond = secondTrackExt->getTanLambda();
- const float phi0Second = secondTrackExt->getPhi();
+ const float phi0Second = secondTrackExt->getPhi();
HelixClass firstHelix;
firstHelix.Initialize_Canonical(phi0First,d0First,z0First,omegaFirst,tanLFirst,_bField);
const float pxFirst = firstHelix.getMomentum()[0];
const float pyFirst = firstHelix.getMomentum()[1];
- const float pzFirst = firstHelix.getMomentum()[2];
+ const float pzFirst = firstHelix.getMomentum()[2];
HelixClass secondHelix;
secondHelix.Initialize_Canonical(phi0Second,d0Second,z0Second,omegaSecond,tanLSecond,_bField);
const float pxSecond = secondHelix.getMomentum()[0];
const float pySecond = secondHelix.getMomentum()[1];
- const float pzSecond = secondHelix.getMomentum()[2];
+ const float pzSecond = secondHelix.getMomentum()[2];
const float pFirst = sqrt(pxFirst*pxFirst+pyFirst*pyFirst+pzFirst*pzFirst);
const float pSecond = sqrt(pxSecond*pxSecond+pySecond*pySecond+pzSecond*pzSecond);
@@ -4767,7 +4784,7 @@ bool FullLDCTrackingAlg::VetoMerge(TrackExtended* firstTrackExt, TrackExtended*
return veto;
}
-StatusCode FullLDCTrackingAlg::finalize() {
+StatusCode FullLDCTrackingAlg::finalize() {
delete _encoder ;
return StatusCode::SUCCESS;
}
@@ -4822,7 +4839,7 @@ void FullLDCTrackingAlg::setupGearGeom(){
const EVENT::DoubleVec& SIT_r = pSIT.getDoubleVals("SITLayerRadius" ) ;
const EVENT::DoubleVec& SIT_hl = pSIT.getDoubleVals("SITSupportLayerHalfLength" ) ;
- _nLayersSIT = SIT_r.size() ;
+ _nLayersSIT = SIT_r.size();
if (_nLayersSIT != SIT_r.size() || _nLayersSIT != SIT_hl.size()) {
fatal() << "FullLDCTrackingAlg Miss-match between DoubleVec and nlayers exit(1) called from file " << __FILE__ << " line " << __LINE__ << endmsg;
@@ -4831,7 +4848,7 @@ void FullLDCTrackingAlg::setupGearGeom(){
}
catch( ... ){
debug() << " ### gear::SIT Parameters from as defined in SSit03 Not Present in GEAR FILE" << endmsg;
- }
+ }
}
//-- SET Parameters--
@@ -4874,7 +4891,7 @@ void FullLDCTrackingAlg::setupGearGeom(){
}
//-- FTD Parameters--
- _petalBasedFTDWithOverlaps = false;
+ _petalBasedFTDWithOverlaps = false;
_nLayersFTD = 0;
try{
@@ -4897,12 +4914,12 @@ void FullLDCTrackingAlg::setupGearGeom(){
}
// SJA: Here we increase the size of _nlayersFTD as we are treating the
- _nLayersFTD =_zLayerFTD.size() ;
+ _nLayersFTD =_zLayerFTD.size() ;
}
catch( ... ){
debug() << " ### gear::FTDParameters Not Present in GEAR FILE" << endmsg;
- }
+ }
if( _nLayersFTD == 0 ){
// FTD
@@ -4923,6 +4940,18 @@ void FullLDCTrackingAlg::setupGearGeom(){
}
catch( ... ){
debug() << " ### gear::FTD Parameters as defined in SFtd05 Not Present in GEAR FILE" << endmsg;
- }
+ }
+ }
+}
+
+void FullLDCTrackingAlg::checkTrackState(int type){
+ debug() << "check TrackState " << endmsg;
+ for(int i=0;i<_allTPCTracks.size();i++){
+ TrackExtended* tpcTrack=_allTPCTracks[i];
+ if(hasTrackStateAt(tpcTrack->getTrack(), type/*lcio::TrackState::AtLastHit*/)){
+ //info() << tpcTrack->getTrack().id() << endmsg;
+ edm4hep::TrackState ts = getTrackStateAt(tpcTrack->getTrack(), type/*lcio::TrackState::AtLastHit*/);
+ debug() << i << " " << tpcTrack->getTrack().id() << ": ts " << ts << endmsg;
+ }
}
}
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
index cf3e174a81988133c63b645eff5b6c5b6d2ccffd..aa73c98c75e1dcce87be47f6a2364ae571250b3c 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
@@ -307,7 +307,8 @@ protected:
int SegmentRadialOverlap(TrackExtended* pTracki, TrackExtended* pTrackj);
bool VetoMerge(TrackExtended* firstTrackExt, TrackExtended* secondTrackExt);
-
+
+ void checkTrackState(int type=0);
int _nRun ;
int _nEvt ;
diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h b/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
index 901b5a276def17154b5c4a51afcb5c6c836a199b..150dd24bf2f87f2e30504e172729e9d2f713c377 100644
--- a/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
+++ b/Service/TrackSystemSvc/include/TrackSystemSvc/IMarlinTrack.h
@@ -55,7 +55,7 @@ namespace MarlinTrk{
/** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
* this order will define the direction of the energy loss used in the fit
*/
- virtual int addHit(edm4hep::ConstTrackerHit hit) = 0 ;
+ virtual int addHit(edm4hep::ConstTrackerHit& hit) = 0 ;
/** initialise the fit using the hits added up to this point -
* the fit direction has to be specified using IMarlinTrack::backward or IMarlinTrack::forward.
@@ -82,12 +82,12 @@ namespace MarlinTrk{
/** update the current fit using the supplied hit, return code via int. Provides the Chi2 increment to the fit from adding the hit via reference.
* the given hit will not be added if chi2increment > maxChi2Increment.
*/
- virtual int addAndFit( edm4hep::ConstTrackerHit hit, double& chi2increment, double maxChi2Increment=DBL_MAX ) = 0 ;
+ virtual int addAndFit( edm4hep::ConstTrackerHit& hit, double& chi2increment, double maxChi2Increment=DBL_MAX ) = 0 ;
/** obtain the chi2 increment which would result in adding the hit to the fit. This method will not alter the current fit, and the hit will not be stored in the list of hits or outliers
*/
- virtual int testChi2Increment( edm4hep::ConstTrackerHit hit, double& chi2increment ) = 0 ;
+ virtual int testChi2Increment( edm4hep::ConstTrackerHit& hit, double& chi2increment ) = 0 ;
/** smooth all track states
@@ -97,7 +97,7 @@ namespace MarlinTrk{
/** smooth track states from the last filtered hit back to the measurement site associated with the given hit
*/
- virtual int smooth( edm4hep::ConstTrackerHit hit ) = 0 ;
+ virtual int smooth( edm4hep::ConstTrackerHit& hit ) = 0 ;
@@ -110,7 +110,7 @@ namespace MarlinTrk{
/** get track state at measurement associated with the given hit, returning TrackState, chi2 and ndf via reference
*/
- virtual int getTrackState( edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
+ virtual int getTrackState( edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
/** get the list of hits included in the fit, together with the chi2 contributions of the hits.
* Pointers to the hits together with their chi2 contribution will be filled into a vector of
@@ -144,7 +144,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
* returning TrackState, chi2 and ndf via reference
*/
- virtual int propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
+ virtual int propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
/** propagate fit to numbered sensitive layer, returning TrackState, chi2, ndf and integer ID of the intersected sensitive detector element via reference
@@ -154,7 +154,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of the intersected sensitive detector element via reference
*/
- virtual int propagateToLayer(int layerID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) = 0;
+ virtual int propagateToLayer(int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) = 0;
/** propagate the fit to sensitive detector element, returning TrackState, chi2 and ndf via reference
*/
@@ -163,7 +163,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning TrackState, chi2 and ndf via reference
*/
- virtual int propagateToDetElement( int detEementID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) = 0 ;
+ virtual int propagateToDetElement( int detEementID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) = 0 ;
@@ -176,7 +176,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
* returning TrackState, chi2 and ndf via reference
*/
- virtual int extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
+ virtual int extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) = 0 ;
/** extrapolate the fit to numbered sensitive layer, returning TrackState, chi2, ndf and integer ID of the intersected sensitive detector element via reference
*/
@@ -185,7 +185,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of the intersected sensitive detector element via reference
*/
- virtual int extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) = 0 ;
+ virtual int extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) = 0 ;
/** extrapolate the fit to sensitive detector element, returning TrackState, chi2 and ndf via reference
*/
@@ -194,7 +194,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning TrackState, chi2 and ndf via reference
*/
- virtual int extrapolateToDetElement( int detEementID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) = 0 ;
+ virtual int extrapolateToDetElement( int detEementID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) = 0 ;
// INTERSECTORS
@@ -207,7 +207,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning intersection point in global coordinates and integer ID of the intersected sensitive detector element via reference
*/
- virtual int intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit hit, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) = 0 ;
+ virtual int intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) = 0 ;
/** extrapolate the fit to numbered sensitive detector element, returning intersection point in global coordinates via reference
@@ -217,7 +217,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning intersection point in global coordinates via reference
*/
- virtual int intersectionWithDetElement( int detEementID, edm4hep::ConstTrackerHit hit, edm4hep::Vector3d& point, int mode=modeClosest ) = 0 ;
+ virtual int intersectionWithDetElement( int detEementID, edm4hep::ConstTrackerHit& hit, edm4hep::Vector3d& point, int mode=modeClosest ) = 0 ;
protected:
diff --git a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
index f41413eacaa43c691780c23da5154ed8bfce0a22..0929b8f8ce7f4adb2eab26726d3d7d3df0deb004 100644
--- a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
+++ b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc
@@ -76,7 +76,7 @@ namespace MarlinTrk {
MarlinKalTestTrack::MarlinKalTestTrack(MarlinKalTest* ktest)
- : _ktest(ktest) {
+ : _ktest(ktest), _trackHitAtPositiveNDF(edm4hep::ConstTrackerHit(0)) {
_kaltrack = new TKalTrack() ;
_kaltrack->SetOwner() ;
@@ -88,7 +88,7 @@ namespace MarlinTrk {
_fitDirection = false ;
_smoothed = false ;
- _trackHitAtPositiveNDF = 0;
+ //_trackHitAtPositiveNDF = 0;
_hitIndexAtPositiveNDF = 0;
#ifdef MARLINTRK_DIAGNOSTICS_ON
@@ -109,27 +109,26 @@ namespace MarlinTrk {
- int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit trkhit) {
+ int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit& trkhit) {
return this->addHit( trkhit, _ktest->findMeasLayer( trkhit )) ;
}
- int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit trkhit, const ILDVMeasLayer* ml) {
+ int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit& trkhit, const ILDVMeasLayer* ml) {
//std::cout << "MarlinKalTestTrack::addHit: trkhit = " << trkhit.id() << " addr: " << trkhit << " ml = " << ml << std::endl ;
if( trkhit.isAvailable() && ml ) {
//if(ml){
return this->addHit( trkhit, ml->ConvertLCIOTrkHit(trkhit), ml) ;
}
else {
- //std::cout << "MarlinKalTestTrack::addHit: trkhit = " << trkhit.id() << " addr: " << trkhit << " ml = " << ml << std::endl ;
- //streamlog_out( ERROR ) << " MarlinKalTestTrack::addHit - bad inputs " << trkhit << " ml : " << ml << std::endl ;
+ //std::cout << "MarlinKalTestTrack::addHit: - bad inputs trkhit = " << trkhit.id() << " addr: " << trkhit << " ml = " << ml << std::endl ;
return bad_intputs ;
}
return bad_intputs ;
}
- int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit trkhit, ILDVTrackHit* kalhit, const ILDVMeasLayer* ml) {
+ int MarlinKalTestTrack::addHit( edm4hep::ConstTrackerHit& trkhit, ILDVTrackHit* kalhit, const ILDVMeasLayer* ml) {
//std::cout << "MarlinKalTestTrack::addHit: trkhit = " << trkhit.id() << " ILDVTrackHit: " << kalhit << " ml = " << ml << std::endl ;
if( kalhit && ml ) {
//if(ml){
@@ -230,7 +229,7 @@ namespace MarlinTrk {
}
/*
- streamlog_out(DEBUG2) << "MarlinKalTestTrack::initialise Create initial helix from hits: \n "
+ std::cout << "DEBUG<<<<MarlinKalTestTrack::initialise Create initial helix from hits: \n "
<< "P1 x = " << x1.x() << " y = " << x1.y() << " z = " << x1.z() << " r = " << x1.Perp() << "\n "
<< "P2 x = " << x2.x() << " y = " << x2.y() << " z = " << x2.z() << " r = " << x2.Perp() << "\n "
<< "P3 x = " << x3.x() << " y = " << x3.y() << " z = " << x3.z() << " r = " << x3.Perp() << "\n "
@@ -289,13 +288,12 @@ namespace MarlinTrk {
_initialised = true ;
/*
- streamlog_out( DEBUG2 ) << " track parameters used for init : " << std::scientific << std::setprecision(6)
+ streamlog_out( DEBUG2 ) << " track parameters used for init : " << std::scientific << std::setprecision(6)
<< "\t D0 " << 0.0
<< "\t Phi :" << toBaseRange( helstart.GetPhi0() + M_PI/2. )
- << "\t Omega " << 1. /helstart.GetRho()
+ << "\t Omega " << 1. /helstart.GetRho()
<< "\t Z0 " << 0.0
<< "\t tan(Lambda) " << helstart.GetTanLambda()
-
<< "\t pivot : [" << helstart.GetPivot().X() << ", " << helstart.GetPivot().Y() << ", " << helstart.GetPivot().Z()
<< " - r: " << std::sqrt( helstart.GetPivot().X()*helstart.GetPivot().X()+helstart.GetPivot().Y()*helstart.GetPivot().Y() ) << "]"
<< std::endl ;
@@ -335,7 +333,7 @@ namespace MarlinTrk {
if (_kalhits->GetEntries() == 0) {
- //streamlog_out( ERROR) << "<<<<<< MarlinKalTestTrack::Initialise: Number of Hits is Zero. Cannot Initialise >>>>>>>" << std::endl;
+ std::cout << "ERROR<<<<<< MarlinKalTestTrack::Initialise: Number of Hits is Zero. Cannot Initialise >>>>>>>" << std::endl;
return error ;
}
@@ -348,27 +346,25 @@ namespace MarlinTrk {
}
/*
- streamlog_out( DEBUG2 ) << "MarlinKalTestTrack::initialise using TrackState: track parameters used for init : "
- << "\t D0 " << ts.getD0()
- << "\t Phi :" << ts.getPhi()
- << "\t Omega " << ts.getOmega()
- << "\t Z0 " << ts.getZ0()
- << "\t tan(Lambda) " << ts.getTanLambda()
-
- << "\t pivot : [" << ts.getReferencePoint()[0] << ", " << ts.getReferencePoint()[1] << ", " << ts.getReferencePoint()[2]
- << " - r: " << std::sqrt( ts.getReferencePoint()[0]*ts.getReferencePoint()[0]+ts.getReferencePoint()[1]*ts.getReferencePoint()[1] ) << "]"
- << std::endl ;
+ std::cout << "MarlinKalTestTrack::initialise using TrackState: track parameters used for init : "
+ << " D0 " << ts.D0
+ << " Phi :" << ts.phi
+ << " Omega " << ts.omega
+ << " Z0 " << ts.Z0
+ << " tan(Lambda) " << ts.tanLambda
+ << " pivot : [" << ts.referencePoint[0] << ", " << ts.referencePoint[1] << ", " << ts.referencePoint[2]
+ << " - r: " << std::sqrt( ts.referencePoint[0]*ts.referencePoint[0]+ts.referencePoint[1]*ts.referencePoint[1] ) << "]"
+ << std::endl ;
*/
-
_fitDirection = fitDirection ;
-
+
// for GeV, Tesla, R in mm
double alpha = bfield_z * 2.99792458E-4 ;
double kappa;
if ( bfield_z == 0.0 )
kappa = DBL_MAX;
else kappa = ts.omega / alpha ;
-
+
THelicalTrack helix( -ts.D0,
toBaseRange( ts.phi - M_PI/2. ) ,
kappa,
@@ -380,7 +376,6 @@ namespace MarlinTrk {
bfield_z );
TMatrixD cov(5,5) ;
-
std::array<float, 15> covLCIO = ts.covMatrix;
cov( 0 , 0 ) = covLCIO[ 0] ; // d0, d0
@@ -416,14 +411,14 @@ namespace MarlinTrk {
// cov.Print();
// move the helix to either the position of the last hit or the first depending on initalise_at_end
-
+
// default case initalise_at_end
int index = _kalhits->GetEntries() - 1 ;
// or initialise at start
if( _fitDirection == IMarlinTrack::forward ){
index = 0 ;
}
-
+
TVTrackHit* kalhit = dynamic_cast<TVTrackHit *>(_kalhits->At(index));
double dphi;
@@ -439,7 +434,6 @@ namespace MarlinTrk {
initial_pivot = TVector3(0.0,0.0,0.0);
}
-
// ---------------------------
// Create an initial start site for the track using the hit
// ---------------------------
@@ -468,17 +462,16 @@ namespace MarlinTrk {
surf->CalcXingPointWith(helix, initial_pivot, phi);
} else {
- //streamlog_out( ERROR) << "<<<<<<<<< MarlinKalTestTrack::initialise: dynamic_cast failed for TVSurface >>>>>>>" << std::endl;
+ //std::cout << "ERROR<<<<<<<<< MarlinKalTestTrack::initialise: dynamic_cast failed for TVSurface >>>>>>>" << std::endl;
return error ;
}
}
else {
- //streamlog_out( ERROR) << "<<<<<<<<< MarlinKalTestTrack::initialise: dynamic_cast failed for hit type >>>>>>>" << std::endl;
+ //std::cout << "ERROR<<<<<<<<< MarlinKalTestTrack::initialise: dynamic_cast failed for hit type >>>>>>>" << std::endl;
return error ;
}
-
TVTrackHit& dummyHit = *pDummyHit;
//SJA:FIXME: this constants should go in a header file
@@ -489,14 +482,12 @@ namespace MarlinTrk {
// use dummy hit to create initial site
TKalTrackSite& initialSite = *new TKalTrackSite(dummyHit);
-
initialSite.SetHitOwner();// site owns hit
initialSite.SetOwner(); // site owns states
-
+
// ---------------------------
// Set up initial track state
// ---------------------------
-
helix.MoveTo( initial_pivot, dphi, 0, &cov );
static TKalMatrix initialState(kSdim,1) ;
@@ -523,17 +514,15 @@ namespace MarlinTrk {
if (kSdim == 6) covK(5,5) = 1.e6; // t0
// covK.Print();
-
+
// Add initial states to the site
initialSite.Add(new TKalTrackState(initialState,covK,initialSite,TVKalSite::kPredicted));
initialSite.Add(new TKalTrackState(initialState,covK,initialSite,TVKalSite::kFiltered));
-
+
// add the initial site to the track: that is, give the track initial parameters and covariance
// matrix at the starting measurement layer
_kaltrack->Add(&initialSite);
-
_initialised = true ;
-
#ifdef MARLINTRK_DIAGNOSTICS_ON
@@ -603,7 +592,7 @@ namespace MarlinTrk {
// although calling smooth will natrually update delta chi2.
- if (!_kaltrack->AddAndFilter(*temp_site)) {
+ if (!_kaltrack->AddAndFilter(*temp_site)) {
chi2increment = temp_site->GetDeltaChi2() ;
// get the measurement layer of the current hit
@@ -670,16 +659,16 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::addAndFit( edm4hep::ConstTrackerHit trkhit, double& chi2increment, double maxChi2Increment) {
+ int MarlinKalTestTrack::addAndFit( edm4hep::ConstTrackerHit& trkhit, double& chi2increment, double maxChi2Increment) {
- if( ! trkhit.isAvailable() ) {
+ if( ! trkhit.isAvailable() ) {
std::cout << "Error: MarlinKalTestTrack::addAndFit(edm4hep::TrackerHit trkhit, double& chi2increment, double maxChi2Increment): trkhit == 0" << std::endl;
return bad_intputs ;
}
const ILDVMeasLayer* ml = _ktest->findMeasLayer( trkhit ) ;
- if( ml == 0 ){
+ if( ml == 0 ){
// fg: not sure if ml should ever be 0 - but it seems to happen,
// if point is not on surface and more than one surface exists ...
@@ -735,6 +724,10 @@ namespace MarlinTrk {
<< " NDF = " << _kaltrack->GetNDF()
<< std::endl;
*/
+ //std::cout << "ERROR<<<addAndFit(" << _trackHitAtPositiveNDF.isAvailable() << " " << _kaltrack->GetNDF() << std::endl;
+ }
+ else{
+ //std::cout << "ERROR<<<addAndFit(" << _trackHitAtPositiveNDF.isAvailable() << " " << _kaltrack->GetNDF() << std::endl;
}
return success ;
@@ -743,16 +736,16 @@ namespace MarlinTrk {
- int MarlinKalTestTrack::testChi2Increment( edm4hep::ConstTrackerHit trkhit, double& chi2increment ) {
+ int MarlinKalTestTrack::testChi2Increment( edm4hep::ConstTrackerHit& trkhit, double& chi2increment ) {
- //if( ! trkhit ) {
+ //if( ! trkhit ) {
// streamlog_out( ERROR) << "MarlinKalTestTrack::addAndFit(edm4hep::TrackerHit trkhit, double& chi2increment, double maxChi2Increment): trkhit == 0" << std::endl;
// return IMarlinTrack::bad_intputs ;
//}
const ILDVMeasLayer* ml = _ktest->findMeasLayer( trkhit ) ;
- if( ml == 0 ){
+ if( ml == 0 ){
// fg: not sure if ml should ever be 0 - but it seems to happen,
// if point is not on surface and more than one surface exists ...
@@ -829,7 +822,7 @@ namespace MarlinTrk {
_trackHitAtPositiveNDF = trkhit;
_hitIndexAtPositiveNDF = _kaltrack->IndexOf( site );
/*
- streamlog_out( DEBUG2 ) << ">>>>>>>>>>> Fit is now constrained at : "
+ streamlog_out( DEBUG2 ) << ">>>>>>>>>>> Fit is now constrained at : "
<< decodeILD( trkhit.getCellID() )
<< " pos " << trkhit.getPosition()
<< " trkhit = " << _trackHitAtPositiveNDF
@@ -837,7 +830,11 @@ namespace MarlinTrk {
<< " NDF = " << _kaltrack->GetNDF()
<< std::endl;
*/
+ //std::cout << "ERROR<<<fit( double maxChi2Increment )" << _trackHitAtPositiveNDF.isAvailable() << " " << _kaltrack->GetNDF() << std::endl;
}
+ else{
+ //std::cout << "ERROR<<<fit( double maxChi2Increment )" << _trackHitAtPositiveNDF.isAvailable() << " " << _kaltrack->GetNDF() << std::endl;
+ }
}
else { // hit rejected by the filter, so store in the list of rejected hits
@@ -903,7 +900,7 @@ namespace MarlinTrk {
/** smooth track states from the last filtered hit back to the measurement site associated with the given hit
*/
- int MarlinKalTestTrack::smooth( edm4hep::ConstTrackerHit trkhit ) {
+ int MarlinKalTestTrack::smooth( edm4hep::ConstTrackerHit& trkhit ) {
//streamlog_out( DEBUG2 ) << "MarlinKalTestTrack::smooth( edm4hep::TrackerHit " << trkhit << " ) " << std::endl ;
@@ -944,7 +941,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::getTrackState( edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
+ int MarlinKalTestTrack::getTrackState( edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
//streamlog_out( DEBUG2 ) << "MarlinKalTestTrack::getTrackState(edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts ) using hit: " << trkhit << " with cellID0 = " << trkhit.getCellID() << std::endl ;
@@ -973,7 +970,7 @@ namespace MarlinTrk {
// as they will be added to _hit_chi2_values in the order of fitting
// not in the order of time
//
-// if( _fitDirection == IMarlinTrack::backward ){
+// if( _fitDirection == IMarlinTrack::backward ){
// std::reverse_copy( _hit_chi2_values.begin() , _hit_chi2_values.end() , std::back_inserter( hits ) ) ;
// } else {
// std::copy( _hit_chi2_values.begin() , _hit_chi2_values.end() , std::back_inserter( hits ) ) ;
@@ -992,7 +989,7 @@ namespace MarlinTrk {
// // as they will be added to _hit_chi2_values in the order of fitting
// // not in the order of time
//
-// if( _fitDirection == IMarlinTrack::backward ){
+// if( _fitDirection == IMarlinTrack::backward ){
// std::reverse_copy( _outlier_chi2_values.begin() , _outlier_chi2_values.end() , std::back_inserter( hits ) ) ;
// } else {
// std::copy( _outlier_chi2_values.begin() , _outlier_chi2_values.end() , std::back_inserter( hits ) ) ;
@@ -1003,7 +1000,7 @@ namespace MarlinTrk {
int MarlinKalTestTrack::getNDF( int& ndf ){
- if( _initialised == false ) {
+ if( _initialised == false ) {
return error;
} else {
@@ -1020,12 +1017,13 @@ namespace MarlinTrk {
return success;
}
else{
+ //std::cout << "WARNING<<<<<MarlinKalTestTrack::getTrackerHitAtPositiveNDF>>>>_trackHitAtPositiveNDF not available" << std::endl;
return error;
}
}
- int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ){
+ int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ){
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1033,7 +1031,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
+ int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ) {
TKalTrackSite* site = 0 ;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1044,7 +1042,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, const TKalTrackSite& site ,edm4hep::TrackState& ts, double& chi2, int& ndf ){
+ int MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, const TKalTrackSite& site ,edm4hep::TrackState& ts, double& chi2, int& ndf ){
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::extrapolate( const edm4hep::Vector3d& point, edm4hep::TrackState& ts, double& chi2, int& ndf ) called " << std::endl ;
@@ -1076,7 +1074,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1085,7 +1083,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1097,7 +1095,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::extrapolateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::extrapolateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID ) called " << std::endl ;
@@ -1113,7 +1111,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1122,7 +1120,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1134,7 +1132,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::extrapolateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::extrapolateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) called " << std::endl ;
@@ -1163,7 +1161,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ){
+ int MarlinKalTestTrack::propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf ){
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1264,7 +1262,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::propagateToLayer( int layerID, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1273,7 +1271,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToLayer( int layerID, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::propagateToLayer( int layerID, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1285,7 +1283,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode ) {
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::propagateToLayer( int layerID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID ) called " << std::endl;
@@ -1302,7 +1300,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1311,7 +1309,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToDetElement( int detElementID, edm4hep::ConstTrackerHit trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::propagateToDetElement( int detElementID, edm4hep::ConstTrackerHit& trkhit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1323,7 +1321,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::propagateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
+ int MarlinKalTestTrack::propagateToDetElement( int detElementID, const TKalTrackSite& site, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) {
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::propagateToDetElement( int detElementID, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode ) called " << std::endl ;
@@ -1339,7 +1337,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::Vector3d& point, int mode ) {
+ int MarlinKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::Vector3d& point, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
@@ -1349,7 +1347,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::ConstTrackerHit trkhit, edm4hep::Vector3d& point, int mode ) {
+ int MarlinKalTestTrack::intersectionWithDetElement( int detElementID, edm4hep::ConstTrackerHit& trkhit, edm4hep::Vector3d& point, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1407,7 +1405,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::intersectionWithLayer( int layerID, edm4hep::Vector3d& point, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::intersectionWithLayer( int layerID, edm4hep::Vector3d& point, int& detElementID, int mode ) {
const TKalTrackSite& site = *(dynamic_cast<const TKalTrackSite*>(_kaltrack->Last())) ;
const ILDVMeasLayer* ml = 0;
@@ -1416,7 +1414,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit trkhit, edm4hep::Vector3d& point, int& detElementID, int mode ) {
+ int MarlinKalTestTrack::intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit& trkhit, edm4hep::Vector3d& point, int& detElementID, int mode ) {
TKalTrackSite* site = 0;
int error_code = getSiteFromLCIOHit(trkhit, site);
@@ -1429,7 +1427,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const ILDVMeasLayer*& ml, int mode ) {
+ int MarlinKalTestTrack::intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, const ILDVMeasLayer*& ml, int mode ) {
//streamlog_out(DEBUG2) << "MarlinKalTestTrack::intersectionWithLayer( int layerID, const TKalTrackSite& site, edm4hep::Vector3d& point, int& detElementID, int mode) called layerID = " << layerID << std::endl;
@@ -1508,7 +1506,7 @@ namespace MarlinTrk {
//streamlog_out(DEBUG1) << std::endl ;
- if( crossing_exist == 0 ) {
+ if( crossing_exist == 0 ) {
return no_intersection ;
}
else {
@@ -1545,8 +1543,7 @@ namespace MarlinTrk {
if( error_code == success ) {
// make sure we get the next crossing
- if( fabs(dphi) < dphi_min ) {
-
+ if( fabs(dphi) < dphi_min ) {
dphi_min = fabs(dphi) ;
surf_found = true ;
ml = meas_modules[i];
@@ -1631,24 +1628,22 @@ namespace MarlinTrk {
pivot[2] = helix.GetPivot().Z() ;
ts.referencePoint = pivot;
/*
- streamlog_out( DEBUG2 ) << " kaltest track parameters: "
- << " chi2/ndf " << chi2 / ndf
- << " chi2 " << chi2
- << " ndf " << ndf
- << " prob " << TMath::Prob(chi2, ndf)
- << std::endl
-
- << "\t D0 " << d0 << "[+/-" << sqrt( covLCIO[0] ) << "] "
- << "\t Phi :" << phi << "[+/-" << sqrt( covLCIO[2] ) << "] "
- << "\t Omega " << omega << "[+/-" << sqrt( covLCIO[5] ) << "] "
- << "\t Z0 " << z0 << "[+/-" << sqrt( covLCIO[9] ) << "] "
- << "\t tan(Lambda) " << tanLambda << "[+/-" << sqrt( covLCIO[14]) << "] "
-
- << "\t pivot : [" << pivot[0] << ", " << pivot[1] << ", " << pivot[2]
- << " - r: " << std::sqrt( pivot[0]*pivot[0]+pivot[1]*pivot[1] ) << "]"
- << std::endl ;
+ std::cout << " kaltest track parameters: "
+ << " chi2/ndf " << chi2 / ndf
+ << " chi2 " << chi2
+ << " ndf " << ndf
+ << " prob " << TMath::Prob(chi2, ndf)
+ << std::endl
+ << "\t D0 " << d0 << "[+/-" << sqrt( covLCIO[0] ) << "] "
+ << "\t Phi :" << phi << "[+/-" << sqrt( covLCIO[2] ) << "] "
+ << "\t Omega " << omega << "[+/-" << sqrt( covLCIO[5] ) << "] "
+ << "\t Z0 " << z0 << "[+/-" << sqrt( covLCIO[9] ) << "] "
+ << "\t tan(Lambda) " << tanLambda << "[+/-" << sqrt( covLCIO[14]) << "] "
+
+ << "\t pivot : [" << pivot[0] << ", " << pivot[1] << ", " << pivot[2]
+ << " - r: " << std::sqrt( pivot[0]*pivot[0]+pivot[1]*pivot[1] ) << "]"
+ << std::endl ;
*/
-
}
@@ -1667,7 +1662,7 @@ namespace MarlinTrk {
}
- int MarlinKalTestTrack::getSiteFromLCIOHit( edm4hep::ConstTrackerHit trkhit, TKalTrackSite*& site ) const {
+ int MarlinKalTestTrack::getSiteFromLCIOHit( edm4hep::ConstTrackerHit& trkhit, TKalTrackSite*& site ) const {
std::map<edm4hep::ConstTrackerHit,TKalTrackSite*>::const_iterator it;
diff --git a/Service/TrackSystemSvc/src/MarlinKalTestTrack.h b/Service/TrackSystemSvc/src/MarlinKalTestTrack.h
index 1df08e4972dc66d53050a9f1fbad5a94ff6af2a2..f54afc9b509b8d96d54e313d521d12bb9a36a2bd 100644
--- a/Service/TrackSystemSvc/src/MarlinKalTestTrack.h
+++ b/Service/TrackSystemSvc/src/MarlinKalTestTrack.h
@@ -49,17 +49,17 @@ namespace MarlinTrk{
/** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
* this order will define the direction of the energy loss used in the fit
*/
- int addHit(edm4hep::ConstTrackerHit hit) ;
+ int addHit(edm4hep::ConstTrackerHit& hit) ;
/** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
* this order will define the direction of the energy loss used in the fit
*/
- int addHit(edm4hep::ConstTrackerHit trkhit, const ILDVMeasLayer* ml) ;
+ int addHit(edm4hep::ConstTrackerHit& trkhit, const ILDVMeasLayer* ml) ;
/** add hit to track - the hits have to be added ordered in time ( i.e. typically outgoing )
* this order will define the direction of the energy loss used in the fit
*/
- int addHit( edm4hep::ConstTrackerHit trkhit, ILDVTrackHit* kalhit, const ILDVMeasLayer* ml) ;
+ int addHit( edm4hep::ConstTrackerHit& trkhit, ILDVTrackHit* kalhit, const ILDVMeasLayer* ml) ;
/** initialise the fit using the hits added up to this point -
* the fit direction has to be specified using IMarlinTrack::backward or IMarlinTrack::forward.
@@ -90,13 +90,13 @@ namespace MarlinTrk{
/** smooth track states from the last filtered hit back to the measurement site associated with the given hit
*/
- int smooth( edm4hep::ConstTrackerHit hit ) ;
+ int smooth( edm4hep::ConstTrackerHit& hit ) ;
/** update the current fit using the supplied hit, return code via int. Provides the Chi2 increment to the fit from adding the hit via reference.
* the given hit will not be added if chi2increment > maxChi2Increment.
*/
- int addAndFit( edm4hep::ConstTrackerHit hit, double& chi2increment, double maxChi2Increment=DBL_MAX ) ;
+ int addAndFit( edm4hep::ConstTrackerHit& hit, double& chi2increment, double maxChi2Increment=DBL_MAX ) ;
/** update the current fit using the supplied hit, return code via int. Provides the Chi2 increment to the fit from adding the hit via reference.
* the given hit will not be added if chi2increment > maxChi2Increment.
@@ -106,7 +106,7 @@ namespace MarlinTrk{
/** obtain the chi2 increment which would result in adding the hit to the fit. This method will not alter the current fit, and the hit will not be stored in the list of hits or outliers
*/
- int testChi2Increment( edm4hep::ConstTrackerHit hit, double& chi2increment ) ;
+ int testChi2Increment( edm4hep::ConstTrackerHit& hit, double& chi2increment ) ;
// Track State Accessesors
@@ -118,7 +118,7 @@ namespace MarlinTrk{
/** get track state at measurement associated with the given hit, returning TrackState, chi2 and ndf via reference
*/
- int getTrackState( edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+ int getTrackState( edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
/** get the list of hits included in the fit, together with the chi2 contributions of the hits.
@@ -154,7 +154,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
* returning TrackState, chi2 and ndf via reference
*/
- int propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+ int propagate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
/** propagate the fit at the provided measurement site, to the point of closest approach to the given point,
@@ -170,7 +170,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
*/
- int propagateToLayer( int layerID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+ int propagateToLayer( int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
/** propagate the fit at the measurement site, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
@@ -184,7 +184,7 @@ namespace MarlinTrk{
/** propagate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning TrackState, chi2 and ndf via reference
*/
- int propagateToDetElement( int detEementID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+ int propagateToDetElement( int detEementID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
/** propagate the fit at the measurement site, to sensitive detector element,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
@@ -202,7 +202,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to the point of closest approach to the given point,
* returning TrackState, chi2 and ndf via reference
*/
- int extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
+ int extrapolate( const edm4hep::Vector3d& point, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf ) ;
/** extrapolate the fit at the measurement site, to the point of closest approach to the given point,
* returning TrackState, chi2 and ndf via reference
@@ -216,7 +216,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
*/
- int extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
+ int extrapolateToLayer( int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int& detElementID, int mode=modeClosest ) ;
/** extrapolate the fit at the measurement site, to numbered sensitive layer,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
@@ -230,7 +230,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning TrackState, chi2 and ndf via reference
*/
- int extrapolateToDetElement( int detEementID, edm4hep::ConstTrackerHit hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
+ int extrapolateToDetElement( int detEementID, edm4hep::ConstTrackerHit& hit, edm4hep::TrackState& ts, double& chi2, int& ndf, int mode=modeClosest ) ;
/** extrapolate the fit at the measurement site, to sensitive detector element,
* returning TrackState, chi2, ndf and integer ID of sensitive detector element via reference
@@ -250,7 +250,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to numbered sensitive layer,
* returning intersection point in global coordinates and integer ID of the intersected sensitive detector element via reference
*/
- int intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit hit, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) ;
+ int intersectionWithLayer( int layerID, edm4hep::ConstTrackerHit& hit, edm4hep::Vector3d& point, int& detElementID, int mode=modeClosest ) ;
/** extrapolate the fit at the measurement site, to numbered sensitive layer,
* returning intersection point in global coordinates and integer ID of the intersected sensitive detector element via reference
@@ -265,7 +265,7 @@ namespace MarlinTrk{
/** extrapolate the fit at the measurement site associated with the given hit, to sensitive detector element,
* returning intersection point in global coordinates via reference
*/
- int intersectionWithDetElement( int detElementID, edm4hep::ConstTrackerHit hit, edm4hep::Vector3d& point, int mode=modeClosest ) ;
+ int intersectionWithDetElement( int detElementID, edm4hep::ConstTrackerHit& hit, edm4hep::Vector3d& point, int mode=modeClosest ) ;
/** extrapolate the fit at the measurement site, to sensitive detector element,
* returning intersection point in global coordinates via reference
@@ -297,12 +297,13 @@ namespace MarlinTrk{
/** get the measurement site associated with the given lcio TrackerHit trkhit
*/
- int getSiteFromLCIOHit( edm4hep::ConstTrackerHit trkhit, TKalTrackSite*& site ) const ;
+ int getSiteFromLCIOHit( edm4hep::ConstTrackerHit& trkhit, TKalTrackSite*& site ) const ;
/** helper function to restrict the range of the azimuthal angle to ]-pi,pi]*/
inline double toBaseRange( double phi) const {
+ phi = fmod(phi, 2.*M_PI);
while( phi <= -M_PI ){ phi += 2. * M_PI ; }
while( phi > M_PI ){ phi -= 2. * M_PI ; }
return phi ;
diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
index a3654ec7f8dfec1c2b472b5b507bc9c67a6f8ab8..41723ee622010dcaa175551b62d4e5a903f3a6af 100644
--- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
+++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
@@ -37,7 +37,7 @@ namespace MarlinTrk {
//
// int nrows = 5;
//
-// TMatrixD cov(nrows,nrows) ;
+// TMatrixD cov(nrows,nrows) ;
//
// bool matrix_is_positive_definite = true;
//
@@ -148,16 +148,14 @@ namespace MarlinTrk {
int fit_status = createFit(hit_list, marlinTrk, pre_fit, bfield_z, fit_backwards, maxChi2Increment);
-
- if( fit_status != IMarlinTrack::success ){
-
- //std::cout << "MarlinTrk::createFinalisedLCIOTrack fit failed: fit_status = " << fit_status << std::endl;
-
+ //std::cout << "createFit finished. status = " << fit_status << std::endl;
+ if( fit_status != IMarlinTrack::success ){
+ //std::cout << "MarlinTrk::createFinalisedLCIOTrack fit failed: fit_status = " << fit_status << std::endl;
return fit_status;
}
int error = finaliseLCIOTrack(marlinTrk, track, hit_list);
-
+ //std::cout << "finaliseLCIOTrack. status = " << error << std::endl;
return error;
}
@@ -189,7 +187,7 @@ namespace MarlinTrk {
//
// if (( fit_backwards == IMarlinTrack::backward && pre_fit->getLocation() != lcio::TrackState::AtLastHit )
// ||
-// ( fit_backwards == IMarlinTrack::forward && pre_fit->getLocation() != lcio::TrackState::AtFirstHit )) {
+// ( fit_backwards == IMarlinTrack::forward && pre_fit->getLocation() != lcio::TrackState::AtFirstHit )) {
// std::stringstream ss ;
//
// ss << "MarlinTrk::createFinalisedLCIOTrack track state must be set at either first or last hit. Location = ";
@@ -214,7 +212,7 @@ namespace MarlinTrk {
for( it = hit_list.begin() ; it != hit_list.end() ; ++it ) {
edm4hep::ConstTrackerHit trkHit = *it;
- bool isSuccessful = false;
+ bool isSuccessful = false;
//std::cout << "debug: TrackerHit pointer " << trkHit << std::endl;
if( UTIL::BitSet32( trkHit.getType() )[ UTIL::ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT ] ){ //it is a composite spacepoint
//Split it up and add both hits to the MarlinTrk
@@ -227,7 +225,7 @@ namespace MarlinTrk {
if( marlinTrk->addHit( rawHit ) == IMarlinTrack::success ){
isSuccessful = true; //if at least one hit from the spacepoint gets added
++ndof_added;
-// streamlog_out(DEBUG4) << "MarlinTrk::createFit ndof_added = " << ndof_added << std::endl;
+ //std::cout << "DEBUG<<<<<MarlinTrk::createFit ndof_added = " << ndof_added << std::endl;
}
}
}
@@ -235,7 +233,7 @@ namespace MarlinTrk {
if (marlinTrk->addHit( trkHit ) == IMarlinTrack::success ) {
isSuccessful = true;
ndof_added += 2;
-// streamlog_out(DEBUG4) << "MarlinTrk::createFit ndof_added = " << ndof_added << std::endl;
+ //std::cout << "DEBUG<<<<<MarlinTrk::createFit ndof_added = " << ndof_added << std::endl;
}
}
@@ -243,13 +241,13 @@ namespace MarlinTrk {
added_hits.push_back(trkHit);
}
else{
- //std::cout << "MarlinTrkUtils::createFit Hit " << it - hit_list.begin() << " Dropped " << std::endl;
+ //std::cout << "DEBUG<<<<<MarlinTrkUtils::createFit Hit " << it - hit_list.begin() << " Dropped " << std::endl;
}
}
if( ndof_added < MIN_NDF ) {
- //streamlog_out(DEBUG2) << "MarlinTrk::createFit : Cannot fit less with less than " << MIN_NDF << " degrees of freedom. Number of hits = " << added_hits.size() << " ndof = " << ndof_added << std::endl;
+ //std::cout << "ERROR<<<<<MarlinTrk::createFit : Cannot fit less with less than " << MIN_NDF << " degrees of freedom. Number of hits = " << added_hits.size() << " ndof = " << ndof_added << std::endl;
return IMarlinTrack::bad_intputs;
}
@@ -262,12 +260,12 @@ namespace MarlinTrk {
return_error = marlinTrk->initialise( *pre_fit, bfield_z, IMarlinTrack::backward ) ;
if (return_error != IMarlinTrack::success) {
- //streamlog_out(DEBUG5) << "MarlinTrk::createFit Initialisation of track fit failed with error : " << return_error << std::endl;
+ std::cout << "MarlinTrk::createFit Initialisation of track fit failed with error : " << return_error << std::endl;
return return_error;
}
-
+ //std::cout << "debug:===================initialise " << return_error << std::endl;
///////////////////////////////////////////////////////
// try fit and return error
///////////////////////////////////////////////////////
@@ -329,10 +327,10 @@ namespace MarlinTrk {
if ( fit_backwards == IMarlinTrack::backward ) {
pre_fit->location = MarlinTrk::Location::AtLastHit;
- helixTrack.moveRefPoint(hit_list.back().getPosition()[0], hit_list.back().getPosition()[1], hit_list.back().getPosition()[2]);
+ helixTrack.moveRefPoint(hit_list.back().getPosition()[0], hit_list.back().getPosition()[1], hit_list.back().getPosition()[2]);
} else {
pre_fit->location = MarlinTrk::Location::AtFirstHit;
- helixTrack.moveRefPoint(hit_list.front().getPosition()[0], hit_list.front().getPosition()[1], hit_list.front().getPosition()[2]);
+ helixTrack.moveRefPoint(hit_list.front().getPosition()[0], hit_list.front().getPosition()[1], hit_list.front().getPosition()[2]);
}
@@ -435,8 +433,8 @@ namespace MarlinTrk {
bool is_outlier = false;
// here we loop over outliers as this will be faster than looping over the used hits
for ( unsigned ohit = 0; ohit < outliers.size(); ++ohit) {
- if ( rawHit.id() == outliers[ohit].first.id() ) {
- is_outlier = true;
+ if ( rawHit.id() == outliers[ohit].first.id() ) {
+ is_outlier = true;
break; // break out of loop over outliers
}
}
@@ -450,8 +448,8 @@ namespace MarlinTrk {
bool is_outlier = false;
// here we loop over outliers as this will be faster than looping over the used hits
for ( unsigned ohit = 0; ohit < outliers.size(); ++ohit) {
- if ( trkHit == outliers[ohit].first ) {
- is_outlier = true;
+ if ( trkHit == outliers[ohit].first ) {
+ is_outlier = true;
break; // break out of loop over outliers
}
}
@@ -484,12 +482,12 @@ namespace MarlinTrk {
edm4hep::TrackState* trkStateAtLastHit = new edm4hep::TrackState() ;
edm4hep::ConstTrackerHit lastHit = hits_in_fit.front().first;
- edm4hep::ConstTrackerHit last_constrained_hit = 0 ;
+ edm4hep::ConstTrackerHit last_constrained_hit(0);// = 0 ;
marlintrk->getTrackerHitAtPositiveNDF(last_constrained_hit);
return_error = marlintrk->smooth(lastHit);
- if ( return_error != IMarlinTrack::success ) {
+ if ( return_error != IMarlinTrack::success ) {
//streamlog_out(DEBUG4) << "MarlinTrk::finaliseLCIOTrack: return_code for smoothing to " << lastHit << " = " << return_error << " NDF = " << ndf << std::endl;
delete trkStateAtFirstHit;
delete trkStateAtLastHit;
@@ -511,7 +509,7 @@ namespace MarlinTrk {
return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
- if ( return_error != IMarlinTrack::success ) {
+ if ( return_error != IMarlinTrack::success ) {
//streamlog_out(DEBUG4) << "MarlinTrk::finaliseLCIOTrack: return_code for propagation = " << return_error << " NDF = " << ndf << std::endl;
delete trkStateIP;
delete trkStateAtFirstHit;
@@ -569,7 +567,7 @@ namespace MarlinTrk {
trkStateAtLastHit->location = MarlinTrk::Location::AtLastHit;
track->addToTrackStates(*trkStateAtLastHit);
} else {
- //streamlog_out( DEBUG5 ) << " >>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Last Hit " << last_constrained_hit << std::endl ;
+ std::cout << "ERROR>>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Last Hit " << last_constrained_hit.id() << std::endl ;
delete trkStateAtLastHit;
}
@@ -630,7 +628,7 @@ namespace MarlinTrk {
double chi2 = -DBL_MAX;
int ndf = 0;
- UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ;
+ UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string );
encoder.reset() ; // reset to 0
encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::ECAL ;
@@ -669,7 +667,7 @@ namespace MarlinTrk {
if( track == 0 ){
throw std::runtime_error( std::string("MarlinTrk::addHitsToTrack: TrackImpl == NULL ") ) ;
}
- std::map<int, int> hitNumbers;
+ std::map<int, int> hitNumbers;
hitNumbers[UTIL::ILDDetID::VXD] = 0;
hitNumbers[UTIL::ILDDetID::SIT] = 0;
@@ -682,7 +680,7 @@ namespace MarlinTrk {
cellID_encoder.setValue(hit_list.at(j).getCellID()) ;
int detID = cellID_encoder[UTIL::ILDCellID0::subdet];
++hitNumbers[detID];
- //std::cout << "debug: " << "Hit from Detector " << detID << std::endl;
+ //std::cout << "debug: " << "Hit from Detector " << detID << std::endl;
}
int offset = 2 ;
@@ -713,7 +711,7 @@ namespace MarlinTrk {
throw std::runtime_error( std::string("MarlinTrk::addHitsToTrack: TrackImpl == NULL ") ) ;
}
- std::map<int, int> hitNumbers;
+ std::map<int, int> hitNumbers;
hitNumbers[lcio::ILDDetID::VXD] = 0;
hitNumbers[lcio::ILDDetID::SIT] = 0;
@@ -726,7 +724,7 @@ namespace MarlinTrk {
cellID_encoder.setValue(hit_list.at(j).first.getCellID()) ;
int detID = cellID_encoder[UTIL::ILDCellID0::subdet];
++hitNumbers[detID];
- //std::cout << "debug: Hit from Detector " << detID << std::endl;
+ //std::cout << "debug: Hit from Detector " << detID << std::endl;
}
int offset = 2 ;