diff --git a/Tracking/TrkEventCnv/TrkJiveXML/CMakeLists.txt b/Tracking/TrkEventCnv/TrkJiveXML/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..fdf9349f49d5d74612b6b159d6c6df15e03e5642 --- /dev/null +++ b/Tracking/TrkEventCnv/TrkJiveXML/CMakeLists.txt @@ -0,0 +1,42 @@ +################################################################################ +# Package: TrkJiveXML +################################################################################ + +# Declare the package name: +atlas_subdir( TrkJiveXML ) + +# Declare the package's dependencies: +atlas_depends_on_subdirs( PUBLIC + Control/AthenaBaseComps + GaudiKernel + Tracking/TrkEvent/TrkParameters + graphics/JiveXML + PRIVATE + Control/DataModel + Control/StoreGate + DetectorDescription/AtlasDetDescr + Reconstruction/Particle + Tracking/TrkEvent/TrkCompetingRIOsOnTrack + Tracking/TrkEvent/TrkEventPrimitives + Tracking/TrkEvent/TrkEventUtils + Tracking/TrkEvent/TrkMeasurementBase + Tracking/TrkEvent/TrkParticleBase + Tracking/TrkEvent/TrkRIO_OnTrack + Tracking/TrkEvent/TrkSegment + Tracking/TrkEvent/TrkTrack + Tracking/TrkEvent/TrkTrackLink + Tracking/TrkEvent/TrkTrackSummary + Tracking/TrkEvent/TrkTruthData + Tracking/TrkEvent/VxVertex + Tracking/TrkTools/TrkToolInterfaces ) + +# Component(s) in the package: +atlas_add_component( TrkJiveXML + src/*.cxx + src/components/*.cxx + LINK_LIBRARIES AthenaBaseComps GaudiKernel TrkParameters JiveXMLLib DataModel StoreGateLib SGtests AtlasDetDescr Particle TrkCompetingRIOsOnTrack TrkEventPrimitives TrkEventUtils TrkMeasurementBase TrkParticleBase TrkRIO_OnTrack TrkSegment TrkTrack TrkTrackSummary TrkTruthData VxVertex TrkToolInterfaces ) + +# Install files from the package: +atlas_install_headers( TrkJiveXML ) +atlas_install_joboptions( share/*.py ) + diff --git a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/SegmentRetriever.h b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/SegmentRetriever.h index 9830a0cf7ff4e737fc6b2f9dbb273c97c908a78a..30651bdf2abb416beaaf50ff7d63793df880f347 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/SegmentRetriever.h +++ b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/SegmentRetriever.h @@ -33,7 +33,7 @@ namespace JiveXML { SegmentRetriever(const std::string& type,const std::string& name,const IInterface* parent); /// Retrieve all the data - virtual StatusCode retrieve(ToolHandle<IFormatTool> FormatTool); + virtual StatusCode retrieve(ToolHandle<IFormatTool> &FormatTool); /// Return the name of the data type virtual std::string dataTypeName() const { return typeName; }; diff --git a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/TrackRetriever.h b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/TrackRetriever.h index c3f4ac1e4db8d5794778f595f5598a4f10c9df4c..eac5a6085d384cb9e632a3a75b3f32d06f7a5120 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/TrackRetriever.h +++ b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/TrackRetriever.h @@ -67,7 +67,7 @@ namespace JiveXML{ TrackRetriever(const std::string& type,const std::string& name,const IInterface* parent); /// Retrieve all the data - virtual StatusCode retrieve(ToolHandle<IFormatTool> FormatTool); + virtual StatusCode retrieve(ToolHandle<IFormatTool> &FormatTool); /// Return the name of the data type virtual std::string dataTypeName() const { return typeName; }; diff --git a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/VertexRetriever.h b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/VertexRetriever.h index c1262399b93b19e15c6879ebbecd98f55ffcf3c8..2b22ce364d7eb0ac0cfcba5bfa2add0404f3925f 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/VertexRetriever.h +++ b/Tracking/TrkEventCnv/TrkJiveXML/TrkJiveXML/VertexRetriever.h @@ -39,7 +39,7 @@ namespace JiveXML{ VertexRetriever(const std::string& t,const std::string& n,const IInterface* p); /// Retrieve all the data - virtual StatusCode retrieve(ToolHandle<IFormatTool> FormatTool); + virtual StatusCode retrieve(ToolHandle<IFormatTool> &FormatTool); /// Return the name of the data type virtual std::string dataTypeName() const { return typeName; }; diff --git a/Tracking/TrkEventCnv/TrkJiveXML/share/TrkJiveXML_DataTypes.py b/Tracking/TrkEventCnv/TrkJiveXML/share/TrkJiveXML_DataTypes.py index f6eb139055ab27d1313d24d3ac054268783aa2e8..c6e8ce46426d3ad2c3234be19724cfa3a221f76d 100755 --- a/Tracking/TrkEventCnv/TrkJiveXML/share/TrkJiveXML_DataTypes.py +++ b/Tracking/TrkEventCnv/TrkJiveXML/share/TrkJiveXML_DataTypes.py @@ -19,7 +19,7 @@ theTrackRetriever = JiveXML__TrackRetriever (name = "TrackRetriever") ### Switch to the general tracks: theTrackRetriever.PriorityTrackCollection = "Tracks" ### Select only assorted muon tracks: -theTrackRetriever.OtherTrackCollections = ["CombinedFitMuonTracks","MuonSpectrometerTracks","ConvertedStacoTracks","ConvertedMuIdCBTracks","CombinedInDetTracks","GSFTracks"] +theTrackRetriever.OtherTrackCollections = ["CombinedMuonTracks","MuonSpectrometerTracks","ConvertedStacoTracks","ConvertedMuIdCBTracks","CombinedInDetTracks","GSFTracks"] ### ### The Event Filter track collections are not written to XML by default. ### To write them out, you must uncomment the following line: diff --git a/Tracking/TrkEventCnv/TrkJiveXML/src/SegmentRetriever.cxx b/Tracking/TrkEventCnv/TrkJiveXML/src/SegmentRetriever.cxx index 1f8e2c49c456f2c6fe5172f6509d29db1f822118..b0b9fa3723c66bc4a322a1760d2dbb19efce7cb1 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/src/SegmentRetriever.cxx +++ b/Tracking/TrkEventCnv/TrkJiveXML/src/SegmentRetriever.cxx @@ -35,7 +35,7 @@ namespace JiveXML { * - for each segment try to obtain hits and add as multiple collection * @param FormatTool the tool that will create formated output from the DataMap */ - StatusCode SegmentRetriever::retrieve(ToolHandle<IFormatTool> FormatTool) { + StatusCode SegmentRetriever::retrieve(ToolHandle<IFormatTool> &FormatTool) { //be verbose if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieving " << dataTypeName() <<endreq; diff --git a/Tracking/TrkEventCnv/TrkJiveXML/src/TrackRetriever.cxx b/Tracking/TrkEventCnv/TrkJiveXML/src/TrackRetriever.cxx index 711e8328e7d3fcb0257d05b58dad40613bd3d53c..cbc7f80e9d0a67221245620d277d1f4706782bf1 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/src/TrackRetriever.cxx +++ b/Tracking/TrkEventCnv/TrkJiveXML/src/TrackRetriever.cxx @@ -13,13 +13,10 @@ #include "TrkTrackSummary/TrackSummary.h" #include "TrkToolInterfaces/ITrackSummaryTool.h" -//#include "Particle/TrackParticleContainer.h" -//#include "TrkParameters/Perigee.h" #include "TrkParameters/TrackParameters.h" #include "TrkEventPrimitives/FitQuality.h" #include "TrkEventPrimitives/TrackStateDefs.h" -//ova#include "TrkEventPrimitives/CovarianceMatrix.h" #include "TrkEventPrimitives/JacobianThetaPToCotThetaPt.h" #include "TrkTruthData/TrackTruthCollection.h" #include "TrkTruthData/TrackTruthKey.h" @@ -43,731 +40,668 @@ namespace JiveXML { ///Namespace for all helper functions namespace TrackRetrieverHelpers { - /** - * Obtain the perigee paramets for a given track, if available, and fill - * them in the corresponding data vectors. Perigee parameters are written - * out in the old format using \cot\theta and q/p_T - * @return the perigee parameter object for further use - */ - const Trk::Perigee* getPerigeeParameters( const Trk::Track* track, - DataVect& pt, DataVect& d0, DataVect& z0, DataVect& phi0, DataVect& cotTheta, DataVect& covMatrix){ - - /** - * Get perigee parameters in old format (d_0, z_0, \phi, \cot\theta, q/p_T), - * whereas tracking uses (d_0, z_0, \phi, \theta, q/p), - * therefore a transformation of the covariance matrix is needed - */ - const Trk::Perigee *perigee = track->perigeeParameters(); - - //return immediately if there is no perigee information - if (! perigee) return NULL ; - - //write out p_T - if ((perigee->parameters())[Trk::qOverP]==0) pt.push_back(DataType(9999.)); - else pt.push_back( (perigee->charge() > 0) ? DataType(perigee->pT()/CLHEP::GeV) : DataType((-perigee->pT())/CLHEP::GeV)); - - d0.push_back(DataType((perigee->parameters())[Trk::d0]/CLHEP::cm)); - z0.push_back(DataType(perigee->parameters()[Trk::z0]/CLHEP::cm)); - phi0.push_back(DataType(perigee->parameters()[Trk::phi0])); - - if (perigee->parameters()[Trk::theta] == 0.) cotTheta.push_back(DataType(9999.)); - else cotTheta.push_back(DataType(1./tan(perigee->parameters()[Trk::theta]))); - - // CLHEP->Eigen migration. jpt Dec'13 - // https://twiki.cern.ch/twiki/bin/viewauth/Atlas/MigrationCLHEPtoEigen - // https://twiki.cern.ch/twiki/bin/viewauth/Atlas/MigrationToUpdatedEDM#Changes_to_TrkParameters - - /// get transformed covariance matrix - AmgSymMatrix(5) covVert; - //covVert.clear(); - //const Trk::Perigee* per = - // dynamic_cast<const Trk::Perigee*>(track->perigee()); - const AmgSymMatrix(5)* covariance = perigee ? perigee->covariance() : NULL; - if (perigee && covariance) { - // do trafo to old format - double measuredTheta = perigee->parameters()[Trk::theta]; - double measuredQoverp = perigee->parameters()[Trk::qOverP]; - const Trk::JacobianThetaPToCotThetaPt theJac( measuredTheta, measuredQoverp ); - covVert = covariance->similarity(theJac); - }else{ - for ( int ii=0; ii<20; ii++){ // placeholder. Do this nicer. - covVert(ii) = 0.; - } - } - //Scale covariance matrix values to get good resolution with fixed - //precision in JiveXML data - - const long scale = 10000; -// Migration: Now only has diagonal elements from covariance matrix ? - covMatrix.push_back(DataType(covVert(0)*scale/100.)); // 5 elements - covMatrix.push_back(DataType(covVert(1)*scale/100.)); - covMatrix.push_back(DataType(covVert(2)*scale/100.)); - covMatrix.push_back(DataType(covVert(3)*scale/100.)); - covMatrix.push_back(DataType(covVert(4)*scale/100.)); - -// Used to be 15 elements before migration, so need to put 10 placeholders - for ( int i=0; i<10; i++){ - covMatrix.push_back(DataType( 0. )); - } - - // get transformed covariance matrix -//-----------------old code, before migration --------------- -/* - Trk::CovarianceMatrix covVert(5); - const Trk::Perigee* perigee = dynamic_cast<const Trk::Perigee*>(perigee); - - if (perigee) { - // do trafo to old format - Trk::JacobianThetaPToCotThetaPt theJac(perigee->parameters()[Trk::theta], perigee->parameters()[Trk::qOverP]); - covVert = perigee->localErrorMatrix().covariance().similarity(theJac); - } -*/ - //Scale covariance matrix values to get good resolution with fixed - //precision in JiveXML data -/* - long scale = 10000; - covMatrix.push_back(DataType(covVert[0][0]*scale/100.)); - covMatrix.push_back(DataType(covVert[1][0]*scale/100.)); - covMatrix.push_back(DataType(covVert[1][1]*scale/100.)); - covMatrix.push_back(DataType(covVert[2][0]*scale/10.)); - covMatrix.push_back(DataType(covVert[2][1]*scale/10.)); - covMatrix.push_back(DataType(covVert[2][2]*scale/1.)); - covMatrix.push_back(DataType(covVert[3][0]*scale/10.)); - covMatrix.push_back(DataType(covVert[3][1]*scale/10.)); - covMatrix.push_back(DataType(covVert[3][2]*scale/1.)); - covMatrix.push_back(DataType(covVert[3][3]*scale/1.)); - covMatrix.push_back(DataType(covVert[4][0]*scale/0.01)); - covMatrix.push_back(DataType(covVert[4][1]*scale/0.01)); - covMatrix.push_back(DataType(covVert[4][2]*scale/0.001)); - covMatrix.push_back(DataType(covVert[4][3]*scale/0.001)); - covMatrix.push_back(DataType(covVert[4][4]*scale/0.000001)); -*/ -//-----------------end of old code, before migration --------------- - //All for perigee, return object for use by other functions - return perigee ; - } - - /** - * Get a list of track-State on Surfaces for measurement and - * outlier hits, sorted using the perigee comparison functions - * @return a std::vector of Trk::TrackStateOnSurface* - */ - std::vector<const Trk::TrackStateOnSurface*> getTrackStateOnSurfaces( const Trk::Track* track, const Trk::Perigee* perigee, bool doHitsSorting){ - - // vector for the return object - std::vector<const Trk::TrackStateOnSurface*> TSoSVec; - - // loop over TrackStateOnSurfaces to extract interesting ones - DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos = track->trackStateOnSurfaces()->begin(); - for (; tsos!=track->trackStateOnSurfaces()->end(); ++tsos) { - // include measurements AND outliers: - if ((*tsos)->type(Trk::TrackStateOnSurface::Measurement) || (*tsos)->type(Trk::TrackStateOnSurface::Outlier) ) { - // add to temp vector - TSoSVec.push_back(*tsos); - } // end if TSoS is measurement or outlier - } // end loop over TSoS - - // sort the selected TSoS, if not already sorted using a comparison functor - - if (perigee) { - //Get hold of the comparison functor - Trk::TrackStateOnSurfaceComparisonFunction *compFunc - = new Trk::TrackStateOnSurfaceComparisonFunction(perigee->position(), perigee->momentum(\ -)); - if (compFunc){ - if (doHitsSorting) { - //Sort track state on surface if needed - if (TSoSVec.size() > 2 && !is_sorted(TSoSVec.begin(), TSoSVec.end(), *compFunc)) - sort(TSoSVec.begin(), TSoSVec.end(), *compFunc); - } - } - delete compFunc; - } // end if compFunc - - //Now return that vector - return TSoSVec; - } - - /** - * Get polyline hits if available. Polyline tracks that have less than 2 points are not useful - skip - */ - void getPolylineFromHits( const std::vector<const Trk::TrackStateOnSurface*>& TSoSVec, - DataVect& polylineX, DataVect& polylineY, DataVect& polylineZ, DataVect& numPolyline){ - - int numPoly = 0 ; - if (TSoSVec.size() > 1) { - //Loop over track state on surfaces - std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter; - for (tsosIter=TSoSVec.begin(); tsosIter!=TSoSVec.end(); ++tsosIter) { - // get global track position - if (!(*tsosIter)->trackParameters()) continue ; - const Amg::Vector3D& pos = (*tsosIter)->trackParameters()->position(); -// const Trk::GlobalPosition& pos = (*tsosIter)->trackParameters()->position(); - polylineX.push_back(DataType(pos.x()/10.)); - polylineY.push_back(DataType(pos.y()/10.)); - polylineZ.push_back(DataType(pos.z()/10.)); - ++numPoly; - } - } - //Store counter as well - numPolyline.push_back(DataType(numPoly)); - } - - - /** - * Retrieve all the basic hit information from the Trk::TrackStateOnSurface - * @return the reconstruction input object (RIO) for further use - */ - const Trk::RIO_OnTrack* getBaseInfoFromHit( const Trk::TrackStateOnSurface* tsos, const AtlasDetectorID* idHelper, - DataVect& isOutlier, DataVect& hits, DataVect& driftSign, DataVect& tsosDetType){ - - //Get corresponding measurement - const Trk::MeasurementBase *measurement = tsos->measurementOnTrack(); - - //If measurement is invalid, return imediately - if ( ! measurement ) return NULL ; - - // get RIO_OnTrack - const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(measurement); - - // check for competing RIO_OnTracks - if (!rot) { - // try to get identifier by CompetingROT: - const Trk::CompetingRIOsOnTrack* comprot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(measurement); - //Get the input object with highest probability - if (comprot) rot = &(comprot->rioOnTrack(comprot->indexOfMaxAssignProb())); - } - - //If there is still no RIO_onTrack, return Null - if (!rot) return NULL ; - - // Now start writing out values: - // Check if this is an outlier hit - isOutlier.push_back(DataType(tsos->type(Trk::TrackStateOnSurface::Outlier))); - - //Now try to get the identifier, create an empty invalid one if no rot - Identifier hitId (rot->identify()); - //Check if it is valid, othwerise store 0 - hits.push_back( hitId.is_valid()?DataType( hitId.get_compact() ):DataType(0) ); - - // get sign of drift radius for TRT measurements - int theDriftSign = 0; - if (idHelper->is_trt(hitId)) { - // get local parameters - theDriftSign = measurement->localParameters()[Trk::driftRadius] > 0. ? 1 : -1; - } - driftSign.push_back(DataType(theDriftSign)); - - //Now get the detector type of the hit - if ( !hitId.is_valid() ){ - tsosDetType.push_back(DataType("unident")); - } else if (idHelper->is_pixel(hitId) ) { - tsosDetType.push_back(DataType("PIX")); // is PIX in Atlantis - } else if (idHelper->is_sct(hitId)) { - tsosDetType.push_back(DataType("SIL")); // is SIL in Atlantis - } else if (idHelper->is_trt(hitId)) { - tsosDetType.push_back(DataType("TRT")); - } else if (idHelper->is_mdt(hitId)) { - tsosDetType.push_back(DataType("MDT")); - } else if (idHelper->is_csc(hitId)) { - tsosDetType.push_back(DataType("CSC")); - } else if (idHelper->is_rpc(hitId)) { - tsosDetType.push_back(DataType("RPC")); - } else if (idHelper->is_tgc(hitId)) { - tsosDetType.push_back(DataType("TGC")); - } else { - tsosDetType.push_back(DataType("unident")); - } - - //Return the reco input object - return rot; - } - - - /** - * Get the residual pull information from the Trk::TrackStateOnSurface hit - */ - void getResidualPullFromHit( const Trk::TrackStateOnSurface* tsos, const Trk::RIO_OnTrack* rot, - ToolHandle<Trk::IResidualPullCalculator> residualPullCalculator, - DataVect& tsosResLoc1, DataVect& tsosResLoc2, DataVect& tsosPullLoc1, DataVect& tsosPullLoc2 ){ - - //Define default return values for invalid states - double ResLoc1 = -99.; - double ResLoc2 = -99.; - double PullLoc1 = -99.; - double PullLoc2 = -99.; - - // get TrackParameters on the surface to calculate residual - const Trk::TrackParameters* tsosParameters = tsos->trackParameters(); - - //Check we got the parameters - if (tsosParameters){ - - /** - * Using track residual tool: ResidualPullCalculator - * -excerpt from Tracking/TrkValidation/TrkValTools/src/BasicValidationNtupleTool.cxx - */ - - //Get the residualPull object - const Trk::ResidualPull* residualPull = residualPullCalculator->residualPull(rot, tsosParameters,Trk::ResidualPull::Biased); - - if (residualPull) { - //Get the first residual - ResLoc1 = residualPull->residual()[Trk::loc1]; - //Get the second residual for more than 1 dimension - if (residualPull->dimension() >= 2) ResLoc2 = residualPull->residual()[Trk::loc2]; - - if ((residualPull->isPullValid()) ) { - //Get the first residual - PullLoc1 = residualPull->pull()[Trk::loc1]; - //Get the second residual for more than 1 dimension - if (residualPull->dimension() >= 2) PullLoc2 = residualPull->pull()[Trk::loc2]; - } - } // end if resdiualPull - } // end if tsosParameters - - //Finally store the values - tsosResLoc1.push_back( ResLoc1 ); - tsosResLoc2.push_back( ResLoc2 ); - tsosPullLoc1.push_back( PullLoc1 ); - tsosPullLoc2.push_back( PullLoc2 ); - } - - /** - * Get the barcode of the associated truth track - */ - void getTruthFromTrack( const Trk::Track* track, const TrackCollection* trackCollection, - const TrackTruthCollection* truthCollection, DataVect& barcode){ - - if (!truthCollection) { - //Fill with zero if none found - barcode.push_back(DataType(0)); - return ; - } - - //Get the element link to this track collection - ElementLink<TrackCollection> tracklink; - tracklink.setElement(track); - tracklink.setStorableObject(*trackCollection); - - //try to find it in the truth collection - std::map<Trk::TrackTruthKey,TrackTruth>::const_iterator tempTrackTruthItr = truthCollection->find(tracklink); - - //Fill with zero if none found - if (tempTrackTruthItr == truthCollection->end()){ - //Fill with zero if none found - barcode.push_back(DataType(0)); - return ; - } - - //if found, Store the barcode - barcode.push_back(DataType((*tempTrackTruthItr).second.particleLink().barcode())); - } - - } //namespace TrackRetrieverHelpers - - /** - * This is the standard AthAlgTool constructor - * @param type AlgTool type name - * @param name AlgTool instance name - * @param parent AlgTools parent owning this tool - **/ - TrackRetriever::TrackRetriever(const std::string& type,const std::string& name,const IInterface* parent): - AthAlgTool(type,name,parent), - typeName("Track"), - m_residualPullCalculator("Trk::ResidualPullCalculator/ResidualPullCalculator") { - - //Declare the interface - declareInterface<IDataRetriever>(this); - - m_trackSumTool = ToolHandle<Trk::ITrackSummaryTool>("Trk::TrackSummaryTool/InDetTrackSummaryTool"); - - //Properties - declareProperty("ResidualPullCalculator" , m_residualPullCalculator, "ToolHandle to ResidualPullCaclulator" ); - declareProperty("PriorityTrackCollection", m_PriorityTrackCollection = "CombinedInDetTracks", "Track collections to retrieve first, shown as default in Atlantis"); - declareProperty("OtherTrackCollections" , m_OtherTrackCollections , "Track collections to retrieve, all if empty"); - declareProperty("TrackTruthColName" , m_TrackTruthCollection = "TrackTruthCollection", - "Track collection from which to retrieve the truth associations for the priority track collection"); - declareProperty("DoWriteHLT" , m_doWriteHLT = false, "Whether to write HLTAutoKey objects"); - declareProperty("DoWriteResiduals" , m_doWriteResiduals = true, "Whether to write TrackResiduals"); - declareProperty("DoHitsSorting" , m_doHitsSorting = false, "Whether to sort hits (TrackStateOnSurfaces)"); - declareProperty("DoHitsDetails" , m_doHitsDetails = true, "Whether to write hits details (TrackStateOnSurfaces)"); - } - - /** - * Initialize before event loop - * - retrieve the residual-pull tool - * - setup the ID helper - */ - StatusCode TrackRetriever::initialize() { - - - //Set up ATLAS ID helper to be able to identify the RIO's det-subsystem. - if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) { - msg(MSG::FATAL) << "Could not get AtlasDetectorID helper" << endreq; - return StatusCode::FAILURE; - } - - // try to retrieve residual-pull calculation only if requested - if (m_doWriteResiduals) - if ( m_residualPullCalculator.retrieve().isFailure()) { - //Give a warning - if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Cannot retrieve ResidualPullCalculator tool " << m_residualPullCalculator << endreq; - if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << " -> will not write out residuals and pulls! " << endreq; - //switch of residual writing - m_doWriteResiduals = false ; - } - //All fine - - // get TrackSummaryTool - if ( m_trackSumTool.retrieve().isFailure() ) { - if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Failed to retrieve tool " << m_trackSumTool << endreq; - } else { - if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieved tool " << m_trackSumTool << endreq; - } - - - return StatusCode::SUCCESS; - } - - - - /** - * For each track collection retrieve all data - * - loop over tracks in all collections - * - for each track get basic parameters - * @param FormatTool the tool that will create formated output from the DataMap - */ - StatusCode TrackRetriever::retrieve(ToolHandle<IFormatTool> FormatTool) { - - //be verbose - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieving " << dataTypeName() << endreq; - - //Generate a list of requested track collections - typedef std::pair< TrackCollection , std::string > tracksNamePair; - std::vector< tracksNamePair > requestedTrackColls; - - //First try to get hold of the priorty track collection - const TrackCollection* tracks = NULL ; - if (evtStore()->retrieve(tracks, m_PriorityTrackCollection).isFailure()){ - if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Unable to retrieve requested priority track collection " - << m_PriorityTrackCollection << endreq; - } else { - //Add this to the list of requested collections - requestedTrackColls.push_back(tracksNamePair(*tracks,m_PriorityTrackCollection)); - } - - //If we have been given an explicit list, try to retrieve these in the order - //they were given - std::vector<std::string>::iterator CollNameItr = m_OtherTrackCollections.begin(); - for ( ; CollNameItr != m_OtherTrackCollections.end(); ++CollNameItr){ - const TrackCollection* tracks = NULL ; - if (evtStore()->retrieve(tracks, (*CollNameItr)).isFailure()){ - if (msgLvl(MSG::WARNING)) msg(MSG::WARNING) << "Unable to retrieve requested track collection " << (*CollNameItr) << endreq; - continue ; - } - - //skip if it's priority collection again - if ((*CollNameItr) == m_PriorityTrackCollection){ continue; } - - //Add them to the list of requested collections - requestedTrackColls.push_back(tracksNamePair(*tracks,(*CollNameItr))); - - } // end OtherCollections loop - - //If no collections had been requested explicitly, loop over all of them - if (m_OtherTrackCollections.empty()) { - - //Get an iterator over all other track collections - const DataHandle<TrackCollection> trackCollIter, trackCollEnd; - if ((evtStore()->retrieve(trackCollIter, trackCollEnd)).isFailure()){ - if (msgLvl(MSG::ERROR)) msg(MSG::ERROR) << "Unable to retrieve track collection iterator" << endreq; - return StatusCode::RECOVERABLE; - } - - //Next loop over all collections - for (; trackCollIter!=trackCollEnd; trackCollIter++) { - - //Check if this is an HLT-AutoKey collection - if ((trackCollIter.key().find("HLT",0) != std::string::npos) && (!m_doWriteHLT)){ - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Ignoring HLT-AutoKey collection " << trackCollIter.key() << endreq; - continue ; - } - - //Ignore TRTSeed and MuonSlimmedTrackCollections - if ( (trackCollIter.key()=="TRTSeeds") || (trackCollIter.key() =="MuonSlimmedTrackCollection")) { - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Always ignoring collection " << trackCollIter.key() << endreq; - continue ; - } - - //Next try to retrieve the actual collection - if (evtStore()->retrieve(tracks,trackCollIter.key()).isFailure()){ - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Unable to retrieve collection " << trackCollIter.key() << endreq; - continue ; - } - - //Check if this collection is not already in the list - if (std::find(requestedTrackColls.begin(), requestedTrackColls.end(),tracksNamePair(*tracks,trackCollIter.key())) != requestedTrackColls.end()) - continue ; - - //Add this to the list of requested collections - requestedTrackColls.push_back(tracksNamePair(*tracks,trackCollIter.key())); - } //loop over all track collections - } - - /** - * Second step: - * Now loop over all collections a retrieve the actual values - **/ - - //Loop over the collection list we have assembled above - std::vector<tracksNamePair>::iterator tracksNamePairItr = requestedTrackColls.begin(); - for ( ; tracksNamePairItr != requestedTrackColls.end(); ++tracksNamePairItr){ - - //save some typing by getting a handle on the collection pointer - const TrackCollection* trackCollection =&((*tracksNamePairItr).first); - std::string collectionName = (*tracksNamePairItr).second; - - //Some sanity checks - if ( trackCollection->size() == 0){ - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Ignoring empty track collection " << collectionName << endreq; - } else { - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieving data for track collection " << collectionName << endreq; - } - - /** - * Try to find the appropiate truth collection - * StoreGate key can be either 'TruthCollection' or just 'Truth' suffix. Check both. - * Using 'SG::contains' check first, avoids spitting out SG warnings - */ - const TrackTruthCollection *truthCollection = NULL ; - //Check for given truth collection - if ( collectionName == m_PriorityTrackCollection ){ - if ( evtStore()->contains<TrackTruthCollection>(m_TrackTruthCollection) ){ - evtStore()->retrieve(truthCollection, m_TrackTruthCollection).ignore(); - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Found TrackTruthCollection with key " << m_TrackTruthCollection << endreq; - } - //Check for name+TruthCollection - } else if ( evtStore()->contains<TrackTruthCollection>(collectionName+"TruthCollection") ){ - evtStore()->retrieve(truthCollection, collectionName+"TruthCollection").ignore(); - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Found TrackTruthCollection with key " << collectionName << "TruthCollection" << endreq; - //Check for name+Truth - } else if ( evtStore()->contains<TrackTruthCollection>(collectionName+"Truth") ){ - evtStore()->retrieve(truthCollection, collectionName+"Truth").ignore(); - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Found TrackTruthCollection with key " << collectionName << "Truth" << endreq; - // No matching track truth collection found at all - } else { - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Could not find matching TrackTruthCollection for " << collectionName << endreq; - truthCollection = NULL ; - } - - // Make a list of track-wise entries and reserve enough space - DataVect id; id.reserve(trackCollection->size()); - DataVect chi2; chi2.reserve(trackCollection->size()); - DataVect numDoF; numDoF.reserve(trackCollection->size()); - DataVect trackAuthor; trackAuthor.reserve(trackCollection->size()); - DataVect barcode; barcode.reserve(trackCollection->size()); - DataVect numHits; numHits.reserve(trackCollection->size()); - DataVect numPolyline; numPolyline.reserve(trackCollection->size()); - DataVect nBLayerHits; nBLayerHits.reserve(trackCollection->size()); - DataVect nPixHits; nPixHits.reserve(trackCollection->size()); - DataVect nSCTHits; nSCTHits.reserve(trackCollection->size()); - DataVect nTRTHits; nTRTHits.reserve(trackCollection->size()); - - // vectors with measurement- or TrackStateOnSurface-wise entries - // reserve space later on - DataVect polylineX; - DataVect polylineY; - DataVect polylineZ; - DataVect tsosResLoc1; - DataVect tsosResLoc2; - DataVect tsosPullLoc1; - DataVect tsosPullLoc2; - DataVect tsosDetType; - DataVect isOutlier; - DataVect driftSign; - DataVect hits; - - //Store wether this collection has perigee parameters - DataVect pt; pt.reserve(trackCollection->size()); - DataVect d0; d0.reserve(trackCollection->size()); - DataVect z0; z0.reserve(trackCollection->size()); - DataVect phi0; phi0.reserve(trackCollection->size()); - DataVect cotTheta; cotTheta.reserve(trackCollection->size()); - //Covariance matrix has 15 entries per track - DataVect covMatrix; covMatrix.reserve(trackCollection->size() * 15 ); - - // Now loop over all tracks in the collection - TrackCollection::const_iterator track; - for (track=trackCollection->begin(); track!=trackCollection->end(); ++track) { - - /** - * General track fit info - */ - id.push_back(DataType(id.size())); //<! simple counter starting from 0 - chi2.push_back(DataType((*track)->fitQuality()->chiSquared())); - numDoF.push_back(DataType((*track)->fitQuality()->numberDoF())); - trackAuthor.push_back(DataType((*track)->info().trackFitter())); - - /** - * Get truth Information - */ - TrackRetrieverHelpers::getTruthFromTrack(*track,trackCollection,truthCollection,barcode); - - /** - * Get Perigee parameters (if available) - */ - const Trk::Perigee* perigee = TrackRetrieverHelpers::getPerigeeParameters(*track, pt, d0, z0, phi0, cotTheta, covMatrix); - - /** - * Get number of Pix/SCT/TRT hits - */ - const Trk::TrackSummary* summary = NULL; - summary = m_trackSumTool->createSummary(**track); - - if(summary==0){ - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Track summary is NULL " << endreq; - nBLayerHits.push_back(DataType(0)); - nPixHits.push_back(DataType(0)); - nSCTHits.push_back(DataType(0)); - nTRTHits.push_back(DataType(0)); - } - else{ - nBLayerHits.push_back(DataType(summary->get(Trk::numberOfBLayerHits))); - nPixHits.push_back(DataType(summary->get(Trk::numberOfPixelHits))); - nSCTHits.push_back(DataType(summary->get(Trk::numberOfSCTHits))); - nTRTHits.push_back(DataType(summary->get(Trk::numberOfTRTHits))); - } - - /** - * Get sorted list of track state on surfaces - */ - // Vector of interesting TrackStateOnSurfaces - std::vector<const Trk::TrackStateOnSurface*> TSoSVec = TrackRetrieverHelpers::getTrackStateOnSurfaces(*track,perigee,m_doHitsSorting); - - /** - * Get polyline information - */ - // Reserving some space for polyline hits - polylineX.reserve(polylineX.size()+TSoSVec.size()); - polylineY.reserve(polylineY.size()+TSoSVec.size()); - polylineZ.reserve(polylineZ.size()+TSoSVec.size()); - - //And fill them - TrackRetrieverHelpers::getPolylineFromHits(TSoSVec,polylineX,polylineY,polylineZ,numPolyline); - - /** - * RIO association and outlier id - */ - //Reserve some space for resPull and other hit info - isOutlier.reserve(isOutlier.size()+TSoSVec.size()); - hits.reserve(hits.size()+TSoSVec.size()); - driftSign.reserve(driftSign.size()+TSoSVec.size()); - tsosResLoc1.reserve(tsosResLoc1.size()+TSoSVec.size()); - tsosResLoc2.reserve(tsosResLoc2.size()+TSoSVec.size()); - tsosPullLoc1.reserve(tsosPullLoc1.size()+TSoSVec.size()); - tsosPullLoc2.reserve(tsosPullLoc2.size()+TSoSVec.size()); - tsosDetType.reserve(tsosDetType.size()+TSoSVec.size()); - - //Now loop over tracks and fill them - std::vector< const Trk::TrackStateOnSurface* >::const_iterator TSoSItr = TSoSVec.begin(); - //Count number of hits stored in this loop - long nHits = 0; - if (m_doHitsDetails){ // disable only for HeavyIons ! - for (; TSoSItr != TSoSVec.end(); ++TSoSItr){ - - // This produces the full long debug dump for TSoS: - if (msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << (**TSoSItr) << endreq; - - //Get the basic hit information - const Trk::RIO_OnTrack* rot = TrackRetrieverHelpers::getBaseInfoFromHit(*TSoSItr, m_idHelper, isOutlier, hits, driftSign, tsosDetType ); - - //tell if this didn't work out - if (!rot){ - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Could not obtain RIO for TSoS of type " << (*TSoSItr)->dumpType() << endreq; - continue ; - } - - //count this as a hit - ++nHits; - - //if we shell retrieve residuals, also get those - if (m_doWriteResiduals) - TrackRetrieverHelpers::getResidualPullFromHit( *TSoSItr, rot, m_residualPullCalculator, tsosResLoc1, tsosResLoc2, tsosPullLoc1, tsosPullLoc2); - - } - } // end hits details - - //Store number of retrieved hits for which we have retrieved information - numHits.push_back(nHits); - - } // end loop over tracks in collection - - //Now fill everything in a datamap - DataMap m_DataMap; - // Start with mandatory entries - m_DataMap["id"] = id; - m_DataMap["chi2"] = chi2; - m_DataMap["numDoF"] = numDoF; - m_DataMap["trackAuthor"] = trackAuthor; - m_DataMap["barcode"] = barcode; - m_DataMap["numHits"] = numHits; - m_DataMap["nBLayerHits"] = nBLayerHits; - m_DataMap["nPixHits"] = nPixHits; - m_DataMap["nSCTHits"] = nSCTHits; - m_DataMap["nTRTHits"] = nTRTHits; - m_DataMap["numPolyline"] = numPolyline; - - // if perigee parameters are not available, leave the corresponding subtags empty. - // This way atlantis knows that such tracks can only be displayed as polylines. - if (pt.size() > 0){ - m_DataMap["pt"] = pt; - m_DataMap["d0"] = d0; - m_DataMap["z0"] = z0; - m_DataMap["phi0"] = phi0; - m_DataMap["cotTheta"] = cotTheta; - m_DataMap["covMatrix multiple=\"15\""] = covMatrix; - } - - // vectors with measurement- or TrackStateOnSurface-wise entries - if ( polylineX.size() > 0){ - std::string numPolyPerTrack = DataType(polylineX.size()/((double)id.size())).toString(); - m_DataMap["polylineX multiple=\"" + numPolyPerTrack + "\""] = polylineX; - m_DataMap["polylineY multiple=\"" + numPolyPerTrack + "\""] = polylineY; - m_DataMap["polylineZ multiple=\"" + numPolyPerTrack + "\""] = polylineZ; - } - - if ( hits.size() > 0){ - std::string numHitsPerTrack = DataType(hits.size()/((double)id.size())).toString(); - m_DataMap["hits multiple=\"" + numHitsPerTrack + "\""] = hits; - m_DataMap["isOutlier multiple=\""+numHitsPerTrack+"\""] = isOutlier; - m_DataMap["driftSign multiple=\""+numHitsPerTrack+"\""] = driftSign; - - if (m_doWriteResiduals){ - // hits counter in principle not needed anymore: - m_DataMap["numTsos"] = numHits; - m_DataMap["tsosResLoc1 multiple=\""+numHitsPerTrack+"\""] = tsosResLoc1; - m_DataMap["tsosResLoc2 multiple=\""+numHitsPerTrack+"\""] = tsosResLoc2; - m_DataMap["tsosPullLoc1 multiple=\""+numHitsPerTrack+"\""] = tsosPullLoc1; - m_DataMap["tsosPullLoc2 multiple=\""+numHitsPerTrack+"\""] = tsosPullLoc2; - m_DataMap["tsosDetType multiple=\""+numHitsPerTrack+"\""] = tsosDetType; - } - } - - //forward data to formating tool - if ( FormatTool->AddToEvent(dataTypeName(), collectionName, &m_DataMap).isFailure()) - return StatusCode::RECOVERABLE; - - //Be verbose - if (msgLvl(MSG::DEBUG)) { - msg(MSG::DEBUG) << dataTypeName() << " collection " << collectionName; - msg(MSG::DEBUG) << " retrieved with " << id.size() << " entries"<< endreq; - } - - } //loop over track collections - - //All collections retrieved okay - return StatusCode::SUCCESS; - - } //retrieve + /** + * Obtain the perigee paramets for a given track, if available, and fill + * them in the corresponding data vectors. Perigee parameters are written + * out in the old format using \cot\theta and q/p_T + * @return the perigee parameter object for further use + */ + const Trk::Perigee* getPerigeeParameters( const Trk::Track* track, + DataVect& pt, DataVect& d0, DataVect& z0, DataVect& phi0, + DataVect& cotTheta, DataVect& covMatrix){ + + /** + * Get perigee parameters in old format (d_0, z_0, \phi, \cot\theta, q/p_T), + * whereas tracking uses (d_0, z_0, \phi, \theta, q/p), + * therefore a transformation of the covariance matrix is needed + */ + const Trk::Perigee *perigee = track->perigeeParameters(); + + //return immediately if there is no perigee information + if (! perigee) return nullptr ; + + //write out p_T + if ((perigee->parameters())[Trk::qOverP]==0) pt.push_back(DataType(9999.)); + else pt.push_back( (perigee->charge() > 0) ? DataType(perigee->pT()/CLHEP::GeV) : DataType((-perigee->pT())/CLHEP::GeV)); + + d0.push_back(DataType((perigee->parameters())[Trk::d0]/CLHEP::cm)); + z0.push_back(DataType(perigee->parameters()[Trk::z0]/CLHEP::cm)); + phi0.push_back(DataType(perigee->parameters()[Trk::phi0])); + + if (perigee->parameters()[Trk::theta] == 0.) cotTheta.push_back(DataType(9999.)); + else cotTheta.push_back(DataType(1./tan(perigee->parameters()[Trk::theta]))); + + // CLHEP->Eigen migration. jpt Dec'13 + // https://twiki.cern.ch/twiki/bin/viewauth/Atlas/MigrationCLHEPtoEigen + // https://twiki.cern.ch/twiki/bin/viewauth/Atlas/MigrationToUpdatedEDM#Changes_to_TrkParameters + + /// get transformed covariance matrix + AmgSymMatrix(5) covVert; + + const AmgSymMatrix(5)* covariance = perigee->covariance(); //perigee cannot be null here + if (perigee && covariance) { + // do trafo to old format + double measuredTheta = perigee->parameters()[Trk::theta]; + double measuredQoverp = perigee->parameters()[Trk::qOverP]; + const Trk::JacobianThetaPToCotThetaPt theJac( measuredTheta, measuredQoverp ); + covVert = covariance->similarity(theJac); + }else{ + for ( int ii=0; ii<20; ii++){ // placeholder. Do this nicer. + covVert(ii) = 0.; + } + } + //Scale covariance matrix values to get good resolution with fixed + //precision in JiveXML data + + const long scale = 10000; + const double thisScale(scale/100.); + // Migration: Now only has diagonal elements from covariance matrix ? + covMatrix.push_back(DataType(covVert(0)*thisScale)); // 5 elements + covMatrix.push_back(DataType(covVert(1)*thisScale)); + covMatrix.push_back(DataType(covVert(2)*thisScale)); + covMatrix.push_back(DataType(covVert(3)*thisScale)); + covMatrix.push_back(DataType(covVert(4)*thisScale)); + + // Used to be 15 elements before migration, so need to put 10 placeholders + for ( int i=0; i<10; i++){ + covMatrix.push_back(DataType( 0. )); + } + + //All for perigee, return object for use by other functions + return perigee ; + } + + /** + * Get a list of track-State on Surfaces for measurement and + * outlier hits, sorted using the perigee comparison functions + * @return a std::vector of Trk::TrackStateOnSurface* + */ + std::vector<const Trk::TrackStateOnSurface*> getTrackStateOnSurfaces( const Trk::Track* track, const Trk::Perigee* perigee, bool doHitsSorting){ + // vector for the return object + std::vector<const Trk::TrackStateOnSurface*> TSoSVec; + + // loop over TrackStateOnSurfaces to extract interesting ones + DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos = track->trackStateOnSurfaces()->begin(); + for (; tsos!=track->trackStateOnSurfaces()->end(); ++tsos) { + // include measurements AND outliers: + if ((*tsos)->type(Trk::TrackStateOnSurface::Measurement) || (*tsos)->type(Trk::TrackStateOnSurface::Outlier) ) { + // add to temp vector + TSoSVec.push_back(*tsos); + } // end if TSoS is measurement or outlier + } // end loop over TSoS + + // sort the selected TSoS, if not already sorted using a comparison functor + + if (perigee) { + //Get hold of the comparison functor + Trk::TrackStateOnSurfaceComparisonFunction *compFunc + = new Trk::TrackStateOnSurfaceComparisonFunction(perigee->position(), perigee->momentum(\ + )); + if (compFunc){ + if (doHitsSorting) { + //Sort track state on surface if needed + if (TSoSVec.size() > 2 && !is_sorted(TSoSVec.begin(), TSoSVec.end(), *compFunc)) + sort(TSoSVec.begin(), TSoSVec.end(), *compFunc); + } + } + delete compFunc; + } // end if compFunc + + //Now return that vector + return TSoSVec; + } + + /** + * Get polyline hits if available. Polyline tracks that have less than 2 points are not useful - skip + */ + void getPolylineFromHits( const std::vector<const Trk::TrackStateOnSurface*>& TSoSVec, + DataVect& polylineX, DataVect& polylineY, DataVect& polylineZ, DataVect& numPolyline){ + int numPoly = 0 ; + if (TSoSVec.size() > 1) { + //Loop over track state on surfaces + std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter; + const double onetenth(0.1); + for (tsosIter=TSoSVec.begin(); tsosIter!=TSoSVec.end(); ++tsosIter) { + // get global track position + if (!(*tsosIter)->trackParameters()) continue ; + const Amg::Vector3D& pos = (*tsosIter)->trackParameters()->position(); + polylineX.push_back(DataType(pos.x()*onetenth)); + polylineY.push_back(DataType(pos.y()*onetenth)); + polylineZ.push_back(DataType(pos.z()*onetenth)); + ++numPoly; + } + } + //Store counter as well + numPolyline.push_back(DataType(numPoly)); + } + + + /** + * Retrieve all the basic hit information from the Trk::TrackStateOnSurface + * @return the reconstruction input object (RIO) for further use + */ + const Trk::RIO_OnTrack* getBaseInfoFromHit( const Trk::TrackStateOnSurface* tsos, const AtlasDetectorID* idHelper, + DataVect& isOutlier, DataVect& hits, DataVect& driftSign, DataVect& tsosDetType){ + + //Get corresponding measurement + const Trk::MeasurementBase *measurement = tsos->measurementOnTrack(); + + //If measurement is invalid, return imediately + if ( ! measurement ) return nullptr ; + + // get RIO_OnTrack + const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(measurement); + + // check for competing RIO_OnTracks + if (!rot) { + // try to get identifier by CompetingROT: + const Trk::CompetingRIOsOnTrack* comprot = dynamic_cast<const Trk::CompetingRIOsOnTrack*>(measurement); + //Get the input object with highest probability + if (comprot) rot = &(comprot->rioOnTrack(comprot->indexOfMaxAssignProb())); + } + + //If there is still no RIO_onTrack, return Null + if (!rot) return nullptr ; + + // Now start writing out values: + // Check if this is an outlier hit + isOutlier.push_back(DataType(tsos->type(Trk::TrackStateOnSurface::Outlier))); + + //Now try to get the identifier, create an empty invalid one if no rot + Identifier hitId (rot->identify()); + //Check if it is valid, othwerise store 0 + hits.push_back( hitId.is_valid()?DataType( hitId.get_compact() ):DataType(0) ); + + // get sign of drift radius for TRT measurements + int theDriftSign = 0; + if (idHelper->is_trt(hitId)) { + // get local parameters + theDriftSign = measurement->localParameters()[Trk::driftRadius] > 0. ? 1 : -1; + } + driftSign.push_back(DataType(theDriftSign)); + + //Now get the detector type of the hit + if ( !hitId.is_valid() ){ + tsosDetType.push_back(DataType("unident")); + } else if (idHelper->is_pixel(hitId) ) { + tsosDetType.push_back(DataType("PIX")); // is PIX in Atlantis + } else if (idHelper->is_sct(hitId)) { + tsosDetType.push_back(DataType("SIL")); // is SIL in Atlantis + } else if (idHelper->is_trt(hitId)) { + tsosDetType.push_back(DataType("TRT")); + } else if (idHelper->is_mdt(hitId)) { + tsosDetType.push_back(DataType("MDT")); + } else if (idHelper->is_csc(hitId)) { + tsosDetType.push_back(DataType("CSC")); + } else if (idHelper->is_rpc(hitId)) { + tsosDetType.push_back(DataType("RPC")); + } else if (idHelper->is_tgc(hitId)) { + tsosDetType.push_back(DataType("TGC")); + } else { + tsosDetType.push_back(DataType("unident")); + } + + //Return the reco input object + return rot; + } + + + /** + * Get the residual pull information from the Trk::TrackStateOnSurface hit + */ + void getResidualPullFromHit( const Trk::TrackStateOnSurface* tsos, const Trk::RIO_OnTrack* rot, + const ToolHandle<Trk::IResidualPullCalculator> & residualPullCalculator, + DataVect& tsosResLoc1, DataVect& tsosResLoc2, DataVect& tsosPullLoc1, DataVect& tsosPullLoc2 ){ + + //Define default return values for invalid states + double ResLoc1 = -99.; + double ResLoc2 = -99.; + double PullLoc1 = -99.; + double PullLoc2 = -99.; + + // get TrackParameters on the surface to calculate residual + const Trk::TrackParameters* tsosParameters = tsos->trackParameters(); + + //Check we got the parameters + if (tsosParameters){ + + /** + * Using track residual tool: ResidualPullCalculator + * -excerpt from Tracking/TrkValidation/TrkValTools/src/BasicValidationNtupleTool.cxx + */ + + //Get the residualPull object + const Trk::ResidualPull* residualPull = residualPullCalculator->residualPull(rot, tsosParameters,Trk::ResidualPull::Biased); + + if (residualPull) { + //Get the first residual + ResLoc1 = residualPull->residual()[Trk::loc1]; + //Get the second residual for more than 1 dimension + if (residualPull->dimension() >= 2) ResLoc2 = residualPull->residual()[Trk::loc2]; + + if ((residualPull->isPullValid()) ) { + //Get the first residual + PullLoc1 = residualPull->pull()[Trk::loc1]; + //Get the second residual for more than 1 dimension + if (residualPull->dimension() >= 2) PullLoc2 = residualPull->pull()[Trk::loc2]; + } + } // end if residualPull + } // end if tsosParameters + + //Finally store the values + tsosResLoc1.push_back( ResLoc1 ); + tsosResLoc2.push_back( ResLoc2 ); + tsosPullLoc1.push_back( PullLoc1 ); + tsosPullLoc2.push_back( PullLoc2 ); + } + + /** + * Get the barcode of the associated truth track + */ + void getTruthFromTrack( const Trk::Track* track, const TrackCollection* trackCollection, + const TrackTruthCollection* truthCollection, DataVect& barcode){ + + if (!truthCollection) { + //Fill with zero if none found + barcode.push_back(DataType(0)); + return ; + } + + //Get the element link to this track collection + ElementLink<TrackCollection> tracklink; + tracklink.setElement(track); + tracklink.setStorableObject(*trackCollection); + + //try to find it in the truth collection + std::map<Trk::TrackTruthKey,TrackTruth>::const_iterator tempTrackTruthItr = truthCollection->find(tracklink); + + //Fill with zero if none found + if (tempTrackTruthItr == truthCollection->end()){ + //Fill with zero if none found + barcode.push_back(DataType(0)); + return ; + } + + //if found, Store the barcode + barcode.push_back(DataType((*tempTrackTruthItr).second.particleLink().barcode())); + } + + } //namespace TrackRetrieverHelpers + + /** + * This is the standard AthAlgTool constructor + * @param type AlgTool type name + * @param name AlgTool instance name + * @param parent AlgTools parent owning this tool + **/ + TrackRetriever::TrackRetriever(const std::string& type,const std::string& name,const IInterface* parent): + AthAlgTool(type,name,parent), + typeName("Track"), + m_residualPullCalculator("Trk::ResidualPullCalculator/ResidualPullCalculator") { + //Declare the interface + declareInterface<IDataRetriever>(this); + m_trackSumTool = ToolHandle<Trk::ITrackSummaryTool>("Trk::TrackSummaryTool/InDetTrackSummaryTool"); + //Properties + declareProperty("ResidualPullCalculator" , m_residualPullCalculator, "ToolHandle to ResidualPullCaclulator" ); + declareProperty("PriorityTrackCollection", m_PriorityTrackCollection = "CombinedInDetTracks", "Track collections to retrieve first, shown as default in Atlantis"); + declareProperty("OtherTrackCollections" , m_OtherTrackCollections , "Track collections to retrieve, all if empty"); + declareProperty("TrackTruthColName" , m_TrackTruthCollection = "TrackTruthCollection", + "Track collection from which to retrieve the truth associations for the priority track collection"); + declareProperty("DoWriteHLT" , m_doWriteHLT = false, "Whether to write HLTAutoKey objects"); + declareProperty("DoWriteResiduals" , m_doWriteResiduals = true, "Whether to write TrackResiduals"); + declareProperty("DoHitsSorting" , m_doHitsSorting = false, "Whether to sort hits (TrackStateOnSurfaces)"); + declareProperty("DoHitsDetails" , m_doHitsDetails = true, "Whether to write hits details (TrackStateOnSurfaces)"); + } + + /** + * Initialize before event loop + * - retrieve the residual-pull tool + * - setup the ID helper + */ + StatusCode TrackRetriever::initialize() { + //Set up ATLAS ID helper to be able to identify the RIO's det-subsystem. + if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) { + msg(MSG::FATAL) << "Could not get AtlasDetectorID helper" << endreq; + return StatusCode::FAILURE; + } + // try to retrieve residual-pull calculation only if requested + if (m_doWriteResiduals){ + if ( m_residualPullCalculator.retrieve().isFailure()) { + //Give a warning + ATH_MSG_WARNING( "Cannot retrieve ResidualPullCalculator tool " << m_residualPullCalculator); + ATH_MSG_WARNING( " -> will not write out residuals and pulls! " ); + //switch of residual writing + m_doWriteResiduals = false ; + } + } + //All fine + + // get TrackSummaryTool + if ( m_trackSumTool.retrieve().isFailure() ) { + ATH_MSG_WARNING( "Failed to retrieve tool " << m_trackSumTool ); + } else { + if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieved tool " << m_trackSumTool << endreq; + } + return StatusCode::SUCCESS; + } + + + + /** + * For each track collection retrieve all data + * - loop over tracks in all collections + * - for each track get basic parameters + * @param FormatTool the tool that will create formated output from the DataMap + */ + StatusCode TrackRetriever::retrieve(ToolHandle<IFormatTool> &FormatTool) { + //be verbose + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Retrieving " << dataTypeName() << endreq; + //Generate a list of requested track collections + typedef std::pair< TrackCollection , std::string > tracksNamePair; + std::vector< tracksNamePair > requestedTrackColls; + //First try to get hold of the priorty track collection + const TrackCollection* tracks = NULL ; + if (evtStore()->retrieve(tracks, m_PriorityTrackCollection).isFailure()){ + ATH_MSG_WARNING( "Unable to retrieve requested priority track collection "<< m_PriorityTrackCollection); + } else { + //Add this to the list of requested collections + requestedTrackColls.push_back(tracksNamePair(*tracks,m_PriorityTrackCollection)); + } + + //If we have been given an explicit list, try to retrieve these in the order + //they were given + std::vector<std::string>::iterator CollNameItr = m_OtherTrackCollections.begin(); + for ( ; CollNameItr != m_OtherTrackCollections.end(); ++CollNameItr){ + const TrackCollection* tracks = NULL ; + if (evtStore()->retrieve(tracks, (*CollNameItr)).isFailure()){ + ATH_MSG_WARNING( "Unable to retrieve requested track collection " << (*CollNameItr) ); + continue ; + } + //skip if it's priority collection again + if ((*CollNameItr) == m_PriorityTrackCollection){ continue; } + //Add them to the list of requested collections + requestedTrackColls.push_back(tracksNamePair(*tracks,(*CollNameItr))); + } // end OtherCollections loop + //If no collections had been requested explicitly, loop over all of them + if (m_OtherTrackCollections.empty()) { + //Get an iterator over all other track collections + const DataHandle<TrackCollection> trackCollIter, trackCollEnd; + if ((evtStore()->retrieve(trackCollIter, trackCollEnd)).isFailure()){ + ATH_MSG_ERROR( "Unable to retrieve track collection iterator" ); + return StatusCode::RECOVERABLE; + } + + //Next loop over all collections + for (; trackCollIter!=trackCollEnd; trackCollIter++) { + //Check if this is an HLT-AutoKey collection + if ((trackCollIter.key().find("HLT",0) != std::string::npos) && (!m_doWriteHLT)){ + ATH_MSG_DEBUG( "Ignoring HLT-AutoKey collection " << trackCollIter.key()); + continue ; + } + //Ignore TRTSeed and MuonSlimmedTrackCollections + if ( (trackCollIter.key()=="TRTSeeds") || (trackCollIter.key() =="MuonSlimmedTrackCollection")) { + ATH_MSG_DEBUG( "Always ignoring collection " << trackCollIter.key() ); + continue ; + } + //Next try to retrieve the actual collection + if (evtStore()->retrieve(tracks,trackCollIter.key()).isFailure()){ + ATH_MSG_DEBUG( "Unable to retrieve collection " << trackCollIter.key()); + continue ; + } + //Check if this collection is not already in the list + if (std::find(requestedTrackColls.begin(), requestedTrackColls.end(),tracksNamePair(*tracks,trackCollIter.key())) != requestedTrackColls.end()) + continue ; + + //Add this to the list of requested collections + requestedTrackColls.push_back(tracksNamePair(*tracks,trackCollIter.key())); + } //loop over all track collections + } + + /** + * Second step: + * Now loop over all collections a retrieve the actual values + **/ + + //Loop over the collection list we have assembled above + std::vector<tracksNamePair>::iterator tracksNamePairItr = requestedTrackColls.begin(); + for ( ; tracksNamePairItr != requestedTrackColls.end(); ++tracksNamePairItr){ + //save some typing by getting a handle on the collection pointer + const TrackCollection* trackCollection =&((*tracksNamePairItr).first); + std::string collectionName = (*tracksNamePairItr).second; + //Some sanity checks + if ( trackCollection->size() == 0){ + ATH_MSG_DEBUG( "Ignoring empty track collection " << collectionName ); + } else { + ATH_MSG_DEBUG( "Retrieving data for track collection " << collectionName); + } + + /** + * Try to find the appropiate truth collection + * StoreGate key can be either 'TruthCollection' or just 'Truth' suffix. Check both. + * Using 'SG::contains' check first, avoids spitting out SG warnings + */ + const TrackTruthCollection *truthCollection = NULL ; + //Check for given truth collection + if ( collectionName == m_PriorityTrackCollection ){ + if ( evtStore()->contains<TrackTruthCollection>(m_TrackTruthCollection) ){ + evtStore()->retrieve(truthCollection, m_TrackTruthCollection).ignore(); + ATH_MSG_DEBUG( "Found TrackTruthCollection with key " << m_TrackTruthCollection ); + } + //Check for name+TruthCollection + } else if ( evtStore()->contains<TrackTruthCollection>(collectionName+"TruthCollection") ){ + evtStore()->retrieve(truthCollection, collectionName+"TruthCollection").ignore(); + ATH_MSG_DEBUG( "Found TrackTruthCollection with key " << collectionName << "TruthCollection" ); + //Check for name+Truth + } else if ( evtStore()->contains<TrackTruthCollection>(collectionName+"Truth") ){ + evtStore()->retrieve(truthCollection, collectionName+"Truth").ignore(); + ATH_MSG_DEBUG( "Found TrackTruthCollection with key " << collectionName << "Truth" ); + // No matching track truth collection found at all + } else { + ATH_MSG_DEBUG( "Could not find matching TrackTruthCollection for " << collectionName ); + truthCollection = nullptr ; + } + + // Make a list of track-wise entries and reserve enough space + DataVect id; id.reserve(trackCollection->size()); + DataVect chi2; chi2.reserve(trackCollection->size()); + DataVect numDoF; numDoF.reserve(trackCollection->size()); + DataVect trackAuthor; trackAuthor.reserve(trackCollection->size()); + DataVect barcode; barcode.reserve(trackCollection->size()); + DataVect numHits; numHits.reserve(trackCollection->size()); + DataVect numPolyline; numPolyline.reserve(trackCollection->size()); + DataVect nBLayerHits; nBLayerHits.reserve(trackCollection->size()); + DataVect nPixHits; nPixHits.reserve(trackCollection->size()); + DataVect nSCTHits; nSCTHits.reserve(trackCollection->size()); + DataVect nTRTHits; nTRTHits.reserve(trackCollection->size()); + + // vectors with measurement- or TrackStateOnSurface-wise entries + // reserve space later on + DataVect polylineX; + DataVect polylineY; + DataVect polylineZ; + DataVect tsosResLoc1; + DataVect tsosResLoc2; + DataVect tsosPullLoc1; + DataVect tsosPullLoc2; + DataVect tsosDetType; + DataVect isOutlier; + DataVect driftSign; + DataVect hits; + + //Store wether this collection has perigee parameters + DataVect pt; pt.reserve(trackCollection->size()); + DataVect d0; d0.reserve(trackCollection->size()); + DataVect z0; z0.reserve(trackCollection->size()); + DataVect phi0; phi0.reserve(trackCollection->size()); + DataVect cotTheta; cotTheta.reserve(trackCollection->size()); + //Covariance matrix has 15 entries per track + DataVect covMatrix; covMatrix.reserve(trackCollection->size() * 15 ); + + // Now loop over all tracks in the collection + TrackCollection::const_iterator track; + for (track=trackCollection->begin(); track!=trackCollection->end(); ++track) { + /** + * General track fit info + */ + id.push_back(DataType(id.size())); //<! simple counter starting from 0 + chi2.push_back(DataType((*track)->fitQuality()->chiSquared())); + numDoF.push_back(DataType((*track)->fitQuality()->numberDoF())); + trackAuthor.push_back(DataType((*track)->info().trackFitter())); + + /** + * Get truth Information + */ + TrackRetrieverHelpers::getTruthFromTrack(*track,trackCollection,truthCollection,barcode); + + /** + * Get Perigee parameters (if available) + */ + const Trk::Perigee* perigee = TrackRetrieverHelpers::getPerigeeParameters(*track, pt, d0, z0, phi0, cotTheta, covMatrix); + + /** + * Get number of Pix/SCT/TRT hits + */ + const Trk::TrackSummary* summary = nullptr; + summary = m_trackSumTool->createSummary(**track); + + if(not summary){ + ATH_MSG_DEBUG( "Track summary is NULL " ); + nBLayerHits.push_back(DataType(0)); + nPixHits.push_back(DataType(0)); + nSCTHits.push_back(DataType(0)); + nTRTHits.push_back(DataType(0)); + }else{ + nBLayerHits.push_back(DataType(summary->get(Trk::numberOfBLayerHits))); + nPixHits.push_back(DataType(summary->get(Trk::numberOfPixelHits))); + nSCTHits.push_back(DataType(summary->get(Trk::numberOfSCTHits))); + nTRTHits.push_back(DataType(summary->get(Trk::numberOfTRTHits))); + } + + /** + * Get sorted list of track state on surfaces + */ + // Vector of interesting TrackStateOnSurfaces + std::vector<const Trk::TrackStateOnSurface*> TSoSVec = TrackRetrieverHelpers::getTrackStateOnSurfaces(*track,perigee,m_doHitsSorting); + + /** + * Get polyline information + */ + // Reserving some space for polyline hits + polylineX.reserve(polylineX.size()+TSoSVec.size()); + polylineY.reserve(polylineY.size()+TSoSVec.size()); + polylineZ.reserve(polylineZ.size()+TSoSVec.size()); + + //And fill them + TrackRetrieverHelpers::getPolylineFromHits(TSoSVec,polylineX,polylineY,polylineZ,numPolyline); + + /** + * RIO association and outlier id + */ + //Reserve some space for resPull and other hit info + isOutlier.reserve(isOutlier.size()+TSoSVec.size()); + hits.reserve(hits.size()+TSoSVec.size()); + driftSign.reserve(driftSign.size()+TSoSVec.size()); + tsosResLoc1.reserve(tsosResLoc1.size()+TSoSVec.size()); + tsosResLoc2.reserve(tsosResLoc2.size()+TSoSVec.size()); + tsosPullLoc1.reserve(tsosPullLoc1.size()+TSoSVec.size()); + tsosPullLoc2.reserve(tsosPullLoc2.size()+TSoSVec.size()); + tsosDetType.reserve(tsosDetType.size()+TSoSVec.size()); + + //Now loop over tracks and fill them + std::vector< const Trk::TrackStateOnSurface* >::const_iterator TSoSItr = TSoSVec.begin(); + //Count number of hits stored in this loop + long nHits = 0; + if (m_doHitsDetails){ // disable only for HeavyIons ! + for (; TSoSItr != TSoSVec.end(); ++TSoSItr){ + // This produces the full long debug dump for TSoS: + ATH_MSG_VERBOSE( (**TSoSItr) ); + //Get the basic hit information + const Trk::RIO_OnTrack* rot = TrackRetrieverHelpers::getBaseInfoFromHit(*TSoSItr, m_idHelper, isOutlier, hits, driftSign, tsosDetType ); + //tell if this didn't work out + if (!rot){ + ATH_MSG_VERBOSE( "Could not obtain RIO for TSoS of type " << (*TSoSItr)->dumpType() ); + continue ; + } + //count this as a hit + ++nHits; + + //if we shell retrieve residuals, also get those + if (m_doWriteResiduals){ + TrackRetrieverHelpers::getResidualPullFromHit( *TSoSItr, rot, m_residualPullCalculator, tsosResLoc1, tsosResLoc2, tsosPullLoc1, tsosPullLoc2); + } + } + } // end hits details + + //Store number of retrieved hits for which we have retrieved information + numHits.push_back(nHits); + + } // end loop over tracks in collection + + //Now fill everything in a datamap + DataMap m_DataMap; + // Start with mandatory entries + m_DataMap["id"] = id; + m_DataMap["chi2"] = chi2; + m_DataMap["numDoF"] = numDoF; + m_DataMap["trackAuthor"] = trackAuthor; + m_DataMap["barcode"] = barcode; + m_DataMap["numHits"] = numHits; + m_DataMap["nBLayerHits"] = nBLayerHits; + m_DataMap["nPixHits"] = nPixHits; + m_DataMap["nSCTHits"] = nSCTHits; + m_DataMap["nTRTHits"] = nTRTHits; + m_DataMap["numPolyline"] = numPolyline; + + // if perigee parameters are not available, leave the corresponding subtags empty. + // This way atlantis knows that such tracks can only be displayed as polylines. + if (pt.size() > 0){ + m_DataMap["pt"] = pt; + m_DataMap["d0"] = d0; + m_DataMap["z0"] = z0; + m_DataMap["phi0"] = phi0; + m_DataMap["cotTheta"] = cotTheta; + m_DataMap["covMatrix multiple=\"15\""] = covMatrix; + } + + // vectors with measurement- or TrackStateOnSurface-wise entries + if ( polylineX.size() > 0){ + std::string numPolyPerTrack = DataType(polylineX.size()/((double)id.size())).toString(); + m_DataMap["polylineX multiple=\"" + numPolyPerTrack + "\""] = polylineX; + m_DataMap["polylineY multiple=\"" + numPolyPerTrack + "\""] = polylineY; + m_DataMap["polylineZ multiple=\"" + numPolyPerTrack + "\""] = polylineZ; + } + + if ( hits.size() > 0){ + std::string numHitsPerTrack = DataType(hits.size()/((double)id.size())).toString(); + m_DataMap["hits multiple=\"" + numHitsPerTrack + "\""] = hits; + m_DataMap["isOutlier multiple=\""+numHitsPerTrack+"\""] = isOutlier; + m_DataMap["driftSign multiple=\""+numHitsPerTrack+"\""] = driftSign; + + if (m_doWriteResiduals){ + // hits counter in principle not needed anymore: + m_DataMap["numTsos"] = numHits; + m_DataMap["tsosResLoc1 multiple=\""+numHitsPerTrack+"\""] = tsosResLoc1; + m_DataMap["tsosResLoc2 multiple=\""+numHitsPerTrack+"\""] = tsosResLoc2; + m_DataMap["tsosPullLoc1 multiple=\""+numHitsPerTrack+"\""] = tsosPullLoc1; + m_DataMap["tsosPullLoc2 multiple=\""+numHitsPerTrack+"\""] = tsosPullLoc2; + m_DataMap["tsosDetType multiple=\""+numHitsPerTrack+"\""] = tsosDetType; + } + } + + //forward data to formating tool + if ( FormatTool->AddToEvent(dataTypeName(), collectionName, &m_DataMap).isFailure()) + return StatusCode::RECOVERABLE; + + //Be verbose + if (msgLvl(MSG::DEBUG)) { + msg(MSG::DEBUG) << dataTypeName() << " collection " << collectionName; + msg(MSG::DEBUG) << " retrieved with " << id.size() << " entries"<< endreq; + } + + } //loop over track collections + + //All collections retrieved okay + return StatusCode::SUCCESS; + + } //retrieve } //namespace JiveXML diff --git a/Tracking/TrkEventCnv/TrkJiveXML/src/VertexRetriever.cxx b/Tracking/TrkEventCnv/TrkJiveXML/src/VertexRetriever.cxx index 6098607b19ca83348f20bce3002b2ca895e9bba0..5982d805ca7e996bcf93e00de5f879b49268c70f 100644 --- a/Tracking/TrkEventCnv/TrkJiveXML/src/VertexRetriever.cxx +++ b/Tracking/TrkEventCnv/TrkJiveXML/src/VertexRetriever.cxx @@ -150,7 +150,7 @@ namespace JiveXML { * - find index in track collection for each track associated with this vertex * @param FormatTool the tool that will create formated output from the DataMap */ - StatusCode VertexRetriever::retrieve(ToolHandle<IFormatTool> FormatTool){ + StatusCode VertexRetriever::retrieve(ToolHandle<IFormatTool> &FormatTool){ //Get an iterator over all vertex collections, //return if there are none