diff --git a/Pr/PrAlgorithms/doc/release.notes b/Pr/PrAlgorithms/doc/release.notes index 8f2dbf7572667853414f8549fe6f65dc6e15d2cd..9b78e544d596bd94c88791f03da089cef6bf53ee 100644 --- a/Pr/PrAlgorithms/doc/release.notes +++ b/Pr/PrAlgorithms/doc/release.notes @@ -3,6 +3,10 @@ ! Responsible : Yasmine Amhis ! Purpose : Pattern reconstruction for the upgrade !----------------------------------------------------------------------------- + +! 2016-12-14 - Jeroen van Tilburg + - Adapt to new FT geometry and numbering scheme + 2016-04-18 Renato Quagliani - Fix warnings for unused variables 2016-04-18 Renato Quagliani diff --git a/Pr/PrAlgorithms/src/PrFTHitManager.cpp b/Pr/PrAlgorithms/src/PrFTHitManager.cpp index 21d9ab1cabef943f7539cb06db2a67f0decdb9f7..c27bc56980dcf64758c28c38a0a327b9b3f372f7 100644 --- a/Pr/PrAlgorithms/src/PrFTHitManager.cpp +++ b/Pr/PrAlgorithms/src/PrFTHitManager.cpp @@ -6,8 +6,6 @@ // local #include "PrFTHitManager.h" -//#include "Kernel/FTChannelID.h" - //----------------------------------------------------------------------------- // Implementation file for class : PrFTHitManager // @@ -17,19 +15,15 @@ // Declaration of the Tool Factory DECLARE_TOOL_FACTORY( PrFTHitManager ) - //============================================================================= // Standard constructor, initializes variables //============================================================================= PrFTHitManager::PrFTHitManager( const std::string& type, - const std::string& name, - const IInterface* parent ) + const std::string& name, + const IInterface* parent ) : PrHitManager ( type, name , parent ) { declareInterface<PrFTHitManager>(this); - declareProperty( "DoSmearing" , m_doSmearing = false); - declareProperty( "XSmearing", m_xSmearing = -1. ); - declareProperty( "ZSmearing", m_zSmearing = -1. ); } //========================================================================= @@ -38,66 +32,43 @@ PrFTHitManager::PrFTHitManager( const std::string& type, void PrFTHitManager::buildGeometry ( ) { // -- only build geometry once. if(m_geometryBuilt) return; - if( m_doSmearing){ - IRndmGenSvc* randSvc = svc<IRndmGenSvc>( "RndmGenSvc", true ); - StatusCode sc = m_gauss.initialize( randSvc, Rndm::Gauss( 0., 1. ) ); - if ( sc.isFailure() ){ - error() << "Could not initialize Rndm::Gauss generator" << endmsg; - return; - } - } - m_ftDet = getDet<DeFTDetector>( DeFTDetectorLocation::Default ); + // Retrieve and initialize DeFT (no test: exception in case of failure) + m_ftDet = getDet<DeFTDetector>( DeFTDetectorLocation::Default ); if(UNLIKELY(! m_ftDet ) ){ error()<<" Cannot load the FTDetector from default location "<<endmsg; return; } + if( m_ftDet->version() < 61 ) { + error() << "This version requires FTDet v6.1 or higher" << endmsg; + return; + } - //version 20 is the detector with monolayer and fibremat structure, including the dead regions - if( msgLevel(MSG::DEBUG) ) debug() << "DETECTOR VERSION: " << m_ftDet->version() << endmsg; - //Different Detector Segments has to be used to set the geometry of zone Up and Down to be consistent with new DeFTFibreMat - DetectorSegment segUp; - DetectorSegment segDown; - for ( std::vector<DeFTLayer*>::const_iterator itL = m_ftDet->layers().begin(); //loop over layers - m_ftDet->layers().end() != itL; ++itL ) { - int id = (*itL)->params()->param<int>("layerID"); //ask layerID - - //FTChannelID (layer,module,mat,SiPMID,netCellID) - LHCb::FTChannelID ChUp( id, 2u, 0u, 0u, 0u ); //"0u" means "unsigned 0" - const DeFTFibreMat* ftMatUp = m_ftDet->findFibreMat(ChUp); - segUp = ftMatUp->createDetSegment( ChUp, 0. ); - - LHCb::FTChannelID ChDown( id, 2u, 1u, 0u, 0u ); //"0u" means "unsigned 0" - const DeFTFibreMat* ftMatDown = m_ftDet->findFibreMat(ChDown); - segDown = ftMatDown->createDetSegment( ChDown, 0. ); - - if ( msgLevel(MSG::DEBUG) ){ - debug() << "STEREO ANGLE Up: " << ftMatUp->angle() << endmsg; - debug() << "STEREO ANGLE Down: " << ftMatDown->angle() << endmsg; - } - - - //The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone. - //This is important, since these are used in the following. - //They are set once for each zone in this method. - MakeZone( 2*id, segUp , -4090., 4090., -50., 3030. ); - MakeZone( 2*id+1 , segDown, -4090., 4090., -3030., 50. ); - - //----> Debug zones - if( msgLevel(MSG::DEBUG)){ - debug() << "Layer " << id << " z " << zone(2*id).z() - << " angle " << zone(2*id).dxDy() << endmsg; + // Loop over all layers + for( auto station : m_ftDet->stations() ) { + for( auto layer : station->layers() ) { + int id = 4*(station->stationID() - 1) + layer->layerID(); + + DetectorSegment seg(0, layer->globalZ(), layer->dxdy(), layer->dzdy(), 0., 0.); + float xmax = 0.5*layer->sizeX(); + float ymax = 0.5*layer->sizeY(); + + //The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone. + //This is important, since these are used in the following. + //They are set once for each zone in this method. + MakeZone( 2*id, seg, -xmax, xmax, -25., ymax ); // Small overlap (25 mm) for stereo layers + MakeZone( 2*id+1, seg, -xmax, xmax, -ymax, 25. ); // Small overlap (25 mm) for stereo layers + + //----> Debug zones + if( msgLevel(MSG::DEBUG)){ + debug() << "Layer " << id << " z " << zone(2*id).z() + << " angle " << zone(2*id).dxDy() << endmsg; + } + //----> Debug zones } - //----> Debug zones } - - //---> Debug smearing - if( msgLevel(MSG::DEBUG)) debug() << "XSmearing " << m_xSmearing << " ZSmearing " << m_zSmearing << endmsg; - //---> Debug smearing - + m_geometryBuilt = true; - - } @@ -107,63 +78,41 @@ void PrFTHitManager::buildGeometry ( ) { void PrFTHitManager::decodeData ( ) { if(m_decodedData) return; - - if( msgLevel(MSG::DEBUG) ) debug() << "I AM IN DECODEDATA " << endmsg; - + typedef FastClusterContainer<LHCb::FTLiteCluster,int> FTLiteClusters; - FTLiteClusters* clus = get<FTLiteClusters>( LHCb::FTLiteClusterLocation::Default ); - if( msgLevel(MSG::DEBUG) ) debug() << "Retrieved " << clus->size() << " clusters" << endmsg; - const DeFTFibreMat* ftMat = nullptr; - const DeFTFibreMat* anaFtMat = nullptr; - unsigned int oldFibreMatID = 99999999; - - for ( FTLiteClusters::iterator itC = clus->begin(); clus->end() != itC; ++itC ) { - /// find fibremat to which the cluster belongs - anaFtMat = m_ftDet->findFibreMat( (*itC).channelID() ); - if(anaFtMat->FibreMatID() != oldFibreMatID) { - oldFibreMatID = anaFtMat->FibreMatID(); - ftMat = anaFtMat; - if ( nullptr == ftMat ) { - info() << "FiberMat not found for FT channelID " << (*itC).channelID() << endmsg; - oldFibreMatID = 99999999; - } + FTLiteClusters* clusters = get<FTLiteClusters>( LHCb::FTLiteClusterLocation::Default ); + if ( msgLevel( MSG::DEBUG) ) debug() << "Retrieved " << clusters->size() << " clusters" << endmsg; + + const DeFTMat* mat = nullptr; + unsigned int prevMatID = 99999999; + for ( auto clus : *clusters ) { + if( clus.channelID().uniqueMat() != prevMatID ) { + mat = m_ftDet->findMat( clus.channelID() ); + prevMatID = mat->elementID().uniqueMat(); } - double fraction = (*itC).fraction() + 0.125; // Truncated to 4 bits = 0.25. Add half of it - LHCb::FTChannelID id = (*itC).channelID(); - int lay = (*itC).channelID().layer(); - int zone = (int)((*itC).channelID().mat() > 0); //if yes, it is bottom (now bottom is 1) + double fraction = clus.fraction(); + LHCb::FTChannelID id = clus.channelID(); + + int lay = 4*(id.station()-1) + id.layer(); + int zone = int(mat->isBottom()); int code = 2*lay + zone; - //variables - const double fibreMatHalfSizeY = ftMat->getFibreMatHalfSizeY(); - const double cellSizeX = ftMat->getCellSizeX(); - const float dzDy = ftMat->getDzDy(); - const float dxDy = -ftMat->getTanAngle();// = -m_tanAngle - const double hitLocalX = ftMat->cellLocalX(id) + fraction*cellSizeX; - double xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz; - ftMat->geometry()->toGlobalMatrix().GetComponents( xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz ); - - //const float yBorder = yx* hitLocalX + yy * (2*zone-1)*fibreMatHalfSizeY + dy; - const float yMin = yx* hitLocalX - yy *fibreMatHalfSizeY + dy; - const float yMax = yx* hitLocalX + yy *fibreMatHalfSizeY + dy; - - //calculate hit x,z at y=0 in global frame (as local y=0 NOT global y=0) - const double hitGlobala_X = hitLocalX * xx + dx; - const double hitGlobalb_X = hitGlobala_X + fibreMatHalfSizeY*xy; - const double hitGlobala_Y = hitLocalX * yx+ dy; - const double hitGlobalb_Y = hitGlobala_Y + fibreMatHalfSizeY*yy; - const double hitGlobala_Z = hitLocalX * zx + dz; - const double hitGlobalb_Z = hitGlobala_Z + fibreMatHalfSizeY*zy; - const double axy = (hitGlobala_X-hitGlobalb_X)/(hitGlobala_Y-hitGlobalb_Y); - const double azy = (hitGlobala_Z-hitGlobalb_Z)/(hitGlobala_Y-hitGlobalb_Y); - const float x0 = hitGlobalb_X-axy*hitGlobalb_Y; - const float z0 = hitGlobalb_Z-azy*hitGlobalb_Y; - //Error X position - float errX = 0.05 + .03 * (*itC).size(); - addHitInZone( code, LHCb::LHCbID( (*itC).channelID() ), x0, z0, dxDy, dzDy, yMin, yMax, errX , zone, lay ); - if ( msgLevel(MSG::DEBUG) ) info() << " .. hit " << (*itC).channelID() << " zone " << zone << " x " << x0 << endmsg; + auto endPoints = mat->endPoints(id, fraction); + double invdy = 1./(endPoints.first.y()-endPoints.second.y()); + double dxdy = (endPoints.first.x()-endPoints.second.x())*invdy; + double dzdy = (endPoints.first.z()-endPoints.second.z())*invdy; + float x0 = endPoints.first.x()-dxdy*endPoints.first.y(); + float z0 = endPoints.first.z()-dzdy*endPoints.first.y(); + float yMin = std::min(endPoints.first.y(), endPoints.second.y()); + float yMax = std::max(endPoints.first.y(), endPoints.second.y()); + + float errX = 0.170; // TODO: this should be ~80 micron; get this from a tool + addHitInZone( code, LHCb::LHCbID( id ), x0, z0, dxdy, dzdy, yMin, yMax, errX , zone, lay ); + if ( msgLevel(MSG::VERBOSE) ) + verbose() << " .. hit " << id << " code=" << code << " x=" << x0 << " z=" << z0 << endmsg; } + for ( unsigned int lay = 0; nbZones() > lay ; ++lay ) { std::sort( getIterator_Begin(lay), getIterator_End(lay), PrHit::LowerByX0() ); } diff --git a/Pr/PrAlgorithms/src/PrFTHitManager.h b/Pr/PrAlgorithms/src/PrFTHitManager.h index 8e9ad8b1021342d8779842333bd54874d0d1c62a..7a84e409b0bede5ad13e3ee98c5e7cf12dc8d038 100644 --- a/Pr/PrAlgorithms/src/PrFTHitManager.h +++ b/Pr/PrAlgorithms/src/PrFTHitManager.h @@ -13,10 +13,6 @@ static const InterfaceID IID_PrFTHitManager ( "PrFTHitManager", 1, 0 ); * Tool that transforms clusters into 'hits' (spatial positions) which are then used by the pattern * recognition algorithms involving the FT. * - * Parameters: - * - XSmearing: Amount of gaussian smearing in x - * - ZSmearing: Switch on displacement in z (not used at the moment) - * * @author Olivier Callot * @date 2012-03-13 */ @@ -38,14 +34,9 @@ public: */ void decodeData() override; - protected: private: - DeFTDetector* m_ftDet; - bool m_doSmearing = false; - Rndm::Numbers m_gauss; ///< Random number generator for gaussian smearing - float m_xSmearing =-1.; ///< Amount of smearing in x - float m_zSmearing =-1.; + DeFTDetector* m_ftDet = nullptr; }; #endif // PRFTHITMANAGER_H diff --git a/Pr/PrMCTools/CMakeLists.txt b/Pr/PrMCTools/CMakeLists.txt index 82046778da0c73b4ae67e26e141ca32a76e53161..a682e63eeef664f02075552780a0e383abee2d1f 100644 --- a/Pr/PrMCTools/CMakeLists.txt +++ b/Pr/PrMCTools/CMakeLists.txt @@ -3,7 +3,8 @@ ################################################################################ gaudi_subdir(PrMCTools v2r12) -gaudi_depends_on_subdirs(Det/VPDet +gaudi_depends_on_subdirs(Det/FTDet + Det/VPDet Det/STDet Event/FTEvent Event/LinkerEvent @@ -19,5 +20,5 @@ include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}) gaudi_add_module(PrMCTools src/*.cpp INCLUDE_DIRS Event/FTEvent Pr/PrPixel - LINK_LIBRARIES LinkerEvent MCEvent GaudiAlgLib PrKernel VPDetLib STDetLib LoKiMCLib) + LINK_LIBRARIES LinkerEvent MCEvent GaudiAlgLib PrKernel VPDetLib STDetLib FTDetLib LoKiMCLib) diff --git a/Pr/PrMCTools/doc/release.notes b/Pr/PrMCTools/doc/release.notes index 93dd0ad667a418f32e0318397a01519152855a37..429d61f8032e048c22592422af169d59e162264e 100644 --- a/Pr/PrMCTools/doc/release.notes +++ b/Pr/PrMCTools/doc/release.notes @@ -4,6 +4,10 @@ ! Purpose : MC tools for pattern recognition for the Upgrade !----------------------------------------------------------------------------- +! 2016-12-14 - Jeroen van Tilburg + - Adapt to new FT geometry and numbering scheme + - Changed location of Linker tables. + !========================= PrMCTools v2r12 2016-01-28 ========================= ! 2016-01-21 Marco Cattaneo - Use std::abs instead of fabs for int argument, to silence clang warning diff --git a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp index 061c6a2bd81460305487cf02ec5da0bc4f1624d1..7d1d9575b79d905ad0695e19df4041575020c326 100644 --- a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp +++ b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp @@ -4,7 +4,7 @@ #include "Event/Track.h" #include "Event/StateParameters.h" #include "PrCheatedSciFiTracking.h" -#include "Event/FTCluster.h" +#include "Event/FTLiteCluster.h" #include "Linker/LinkedTo.h" #include "Event/MCTrackInfo.h" @@ -93,7 +93,7 @@ StatusCode PrCheatedSciFiTracking::finalize() { //========================================================================= void PrCheatedSciFiTracking::makeLHCbTracks ( LHCb::Tracks* result ) { - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default ); MCTrackInfo trackInfo( evtSvc(), msgSvc() ); diff --git a/Pr/PrMCTools/src/PrClustersResidual.cpp b/Pr/PrMCTools/src/PrClustersResidual.cpp index 56c3a23fcd15a828fa3b5d527ed5b1de0668102a..a1b5fc54bc7f401d47161f2c95e9fbbf909d9e94 100644 --- a/Pr/PrMCTools/src/PrClustersResidual.cpp +++ b/Pr/PrMCTools/src/PrClustersResidual.cpp @@ -38,16 +38,6 @@ PrClustersResidual::PrClustersResidual( const std::string& name, : GaudiTupleAlg ( name , pSvcLocator ), m_ftHitManager(nullptr), m_zone(24){ - declareProperty("MCHitsLocation",m_mcHitLocation = "/Event/MC/FT/Hits"); - declareProperty("HitManagerName",m_hitManagerName = "PrFTHitManager"); - declareProperty("DebugTracking",m_debugTracking = false); - declareProperty("DoClusterResidual",m_doClusterResidual = false); //Cluster residual tuple - declareProperty("DoTrackStudy",m_doTrackStudy = false); // Produce tuple for track studies - declareProperty("DecodeData",m_decodeData = false); //ask to decode data ( is False if other algorithms are runned - declareProperty("OnlyHasT",m_onlyHasT=false); //(keep only MCParticles with hasT = true - declareProperty("RemoveClones",m_removeClones = true); //Observed clones MCHits when processed ( not understood ) - declareProperty("DumpAllHits",m_dumpAllHits = false); - declareProperty("Occupancy",m_Occupancy = true); } //============================================================================= // Destructor @@ -83,8 +73,17 @@ StatusCode PrClustersResidual::initialize() { //============================================================================= // Main execution //============================================================================= -StatusCode PrClustersResidual::execute() { - +StatusCode PrClustersResidual::execute() { + m_ftDet = getDet<DeFTDetector>(DeFTDetectorLocation::Default); + if( UNLIKELY( !m_ftDet) ){ + error() << " Cannot load the FTDetector from default location" << endmsg; + return StatusCode::FAILURE; + }else if( m_ftDet->version() < 61){ + error()<< "This code requires FTDet v6.1 or higher "<<endmsg; + return StatusCode::FAILURE; + } + + if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg; m_nEvents++; if(m_doClusterResidual){ @@ -124,8 +123,8 @@ void PrClustersResidual::Occupancy(){ //"/Event/Prev/MC/FT/Hits" LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" ); - LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); - LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits"); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), m_ftClusterToParticleLinkerLocation ); + LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), m_ftClusterToHitLinkerLocation); char layerName[100]; char Title[100]; std::vector<int> nHits(12,0); @@ -149,7 +148,7 @@ void PrClustersResidual::Occupancy(){ int recoblewanted = 0; LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default ); LHCb::MCParticle* mcPart = nullptr; - for( LHCb::MCParticles::const_iterator iPart = mcParts->begin(); iPart != mcParts->end(); ++iPart){ + for( auto iPart = mcParts->begin(); iPart != mcParts->end(); ++iPart){ mcPart = (*iPart); bool isRecoble = trackInfo.hasT(mcPart); bool recobleWanted = isRecoble && (mcPart->originVertex()->isPrimary() || mcPart->originVertex()->isDecay()); @@ -174,19 +173,30 @@ void PrClustersResidual::Occupancy(){ for(auto itH = m_ftHitManager->getIterator_Begin(zoneI); m_ftHitManager->getIterator_End(zoneI) != itH;++itH){ hit = &(*itH); nHits[layer]++; - LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id()); + auto liteCluster = getLiteCluster(hit->id()); + + //---- new numbering scheme (v61 sci-fi) plots for occupancy! pseudochannel conversion needed! + /* + It requires to have the m_ftDet as private member of the algorithm and also modify the CMakeLists to compile + PrMCTools linking the Det/FTDet package. + */ + int quarter = liteCluster.channelID().quarter(); + const auto chanID = liteCluster.channelID(); + const DeFTModule* module = m_ftDet->findModule( chanID ); + if( module == nullptr){ + error()<<"Cannot find Module associated to FTChannelID : "<<chanID<<endmsg; + error()<<"Abort Occupancy plot checking"<<endmsg; + return; + } + int pseudochannel = module->pseudoChannel( chanID ); + //--- runs from [0 to 12,287] (96 SiPm *128 channels/SiPm) + //--- 4 quartes [0,1,2,3] , 128 channels squeezed to 1 value. (Occupancy in terms of Hits/SiPM array [128 channels] + //--- ToDo : make it depenent on xml geometry information through m_ftDet + int sipm = std::floor( pseudochannel/128); + int index = quarter * 128 + sipm ; + // old way [ < v61 ] : int index = liteCluster.channelID().uniqueSiPM() & 511; + - int Module = liteCluster.channelID().module(); - int Quarter = liteCluster.channelID().quarter(); - int SipmID = liteCluster.channelID().sipmId(); - int layer = liteCluster.channelID().layer(); - int module2 = Module; - if(Module == 11 && Module > 4 ) module2 = 1; - if(Module != 11 && Module > 4 ) module2 = Module-3; - if(Module < 5 ) module2 = 6-Module; - if(Module == 10 ) module2 = 1; - //---- VALID UNTIL v5x9 SciFi geometry. From v6x1, change in numbering scheme leads to different indexing. - int index = Quarter*128 + 16*module2+ std::floor(SipmID); sprintf(layerName,"Layer%i/AllContributions",layer); sprintf(Title,"Sipm Index Layer %i;SiPm Id;Clusters/Sipm/Event",layer); plot1D(index,layerName,Title,-0.5,511.5,512); @@ -308,8 +318,8 @@ void PrClustersResidual::TrackStudy(){ LHCb::ODIN* odin=getIfExists<LHCb::ODIN>(LHCb::ODINLocation::Default); //Linker for the MCHit LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" ); - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); - LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits"); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); + LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits"); //Get the linker for the generation of the track to study (the one after the linker) //Get the McParticles LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default ); @@ -430,6 +440,8 @@ void PrClustersResidual::TrackStudy(){ debug()<<"For the processed MCParticle i found "<<track.size()<<"PrHits associated to "<<assoc_MCHit.size()<<" MCHits"<<endmsg; tuple->farray("FiredLayers_Counter",firedLayers,"FiredLayers",24); tuple->column("CheatedSeeding_NHits",totHits); + std::vector<double> MCHit_EntryX,MCHit_EntryY,MCHit_EntryZ; + std::vector<double> MCHit_ExitX,MCHit_ExitY,MCHit_ExitZ; std::vector<double> MCHit_X,MCHit_Y,MCHit_Z; std::vector<double> MCHit_tx,MCHit_ty,MCHit_p; std::vector<double> MCHit_pathlength; @@ -443,12 +455,12 @@ void PrClustersResidual::TrackStudy(){ //ChannelID of track hits std::vector<double> ChID_Fraction; - std::vector<double> ChID_Charge; std::vector<double> ChID_SipmCell; + std::vector<double> ChID_Mat; std::vector<double> ChID_Module; std::vector<double> ChID_Layer; std::vector<double> ChID_Quarter; - std::vector<double> ChID_Mat; + std::vector<double> ChID_Station; std::vector<double> ChID_SipmID; std::vector<double> ChID_ID; @@ -477,6 +489,13 @@ void PrClustersResidual::TrackStudy(){ MCHit_ty.push_back(mcHit->dydz()); MCHit_p.push_back(mcHit->p()); MCHit_pathlength.push_back(mcHit->pathLength()); + + MCHit_EntryX.push_back(mcHit->entry().X()); + MCHit_EntryY.push_back(mcHit->entry().Y()); + MCHit_EntryZ.push_back(mcHit->entry().Z()); + MCHit_ExitX.push_back(mcHit->exit().X()); + MCHit_ExitY.push_back(mcHit->exit().Y()); + MCHit_ExitZ.push_back(mcHit->exit().Z()); MCHit_X.push_back(mcHit->midPoint().X()); MCHit_Y.push_back(mcHit->midPoint().Y()); MCHit_Z.push_back(mcHit->midPoint().Z()); @@ -501,14 +520,14 @@ void PrClustersResidual::TrackStudy(){ } LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id()); ChID_Fraction.push_back(liteCluster.fraction()); - ChID_SipmCell.push_back(liteCluster.channelID().sipmCell()); + ChID_SipmCell.push_back(liteCluster.channelID().channel()); ChID_Module.push_back(liteCluster.channelID().module()); + ChID_Mat.push_back(liteCluster.channelID().mat()); ChID_Layer.push_back(liteCluster.channelID().layer()); ChID_Quarter.push_back(liteCluster.channelID().quarter()); - ChID_Mat.push_back(liteCluster.channelID().mat()); - ChID_SipmID.push_back(liteCluster.channelID().sipmId()); + ChID_Station.push_back(liteCluster.channelID().station()); + ChID_SipmID.push_back(liteCluster.channelID().sipm()); ChID_ID.push_back(liteCluster.channelID()); - ChID_Charge.push_back(liteCluster.charge()); PrHit_LHCBID.push_back(hit->id().ftID()); PrHit_Xat0.push_back(hit->x()); PrHit_Zat0.push_back(hit->z()); @@ -530,6 +549,12 @@ void PrClustersResidual::TrackStudy(){ MCHit_time.push_back(-999999.); MCHit_Particle_key.push_back(-999999.); MCHit_pathlength.push_back(-999999.); + MCHit_EntryX.push_back(-999999.); + MCHit_EntryY.push_back(-999999.); + MCHit_EntryZ.push_back(-999999.); + MCHit_ExitX.push_back(-999999.); + MCHit_ExitY.push_back(-999999.); + MCHit_ExitZ.push_back(-999999.); MCHit_X.push_back(-999999.); MCHit_Y.push_back(-999999.); MCHit_Z.push_back(-999999.); @@ -552,6 +577,13 @@ void PrClustersResidual::TrackStudy(){ tuple->farray("MCHit_tx",MCHit_tx,"MC_ass",100); tuple->farray("MCHit_p",MCHit_p,"MC_ass",100); tuple->farray("MCHit_pathlength",MCHit_pathlength,"MC_ass",100); + tuple->farray("MCHit_Assoc_EntryX",MCHit_EntryX ,"MC_ass",100); + tuple->farray("MCHit_Assoc_EntryY",MCHit_EntryY ,"MC_ass",100); + tuple->farray("MCHit_Assoc_EntryZ",MCHit_EntryZ ,"MC_ass",100); + tuple->farray("MCHit_Assoc_ExitX",MCHit_ExitX ,"MC_ass",100); + tuple->farray("MCHit_Assoc_ExitY",MCHit_ExitY ,"MC_ass",100); + tuple->farray("MCHit_Assoc_ExitZ",MCHit_ExitZ ,"MC_ass",100); + tuple->farray("MCHit_Assoc_X",MCHit_X ,"MC_ass",100); tuple->farray("MCHit_Assoc_Y",MCHit_Y ,"MC_ass",100); tuple->farray("MCHit_Assoc_Z",MCHit_Z ,"MC_ass",100); @@ -571,13 +603,13 @@ void PrClustersResidual::TrackStudy(){ tuple->farray("PrHit_dzDy",PrHit_dzDy,"PrHit",100); tuple->farray("ChID_Fraction",ChID_Fraction,"ChID",100); - tuple->farray("ChID_SipmCell",ChID_SipmCell,"ChID",100); - tuple->farray("ChID_Module",ChID_Module,"ChID",100); - tuple->farray("ChID_Layer",ChID_Layer,"ChID",100); - tuple->farray("ChID_Quarter",ChID_Quarter,"ChID",100); - tuple->farray("ChID_Charge",ChID_Charge,"ChID",100); - tuple->farray("ChID_Mat",ChID_Mat,"ChID",100); + tuple->farray("ChID_SipmCell",ChID_SipmCell,"ChID",100); tuple->farray("ChID_SipmID",ChID_SipmID,"ChID",100); + tuple->farray("ChID_Mat",ChID_Mat,"ChID",100); + tuple->farray("ChID_Module",ChID_Module,"ChID",100); + tuple->farray("ChID_Quarter",ChID_Quarter,"ChID",100); + tuple->farray("ChID_Layer",ChID_Layer,"ChID",100); + tuple->farray("ChID_Station",ChID_Station,"ChID",100); tuple->farray("ChID_ID",ChID_ID,"ChID",100); tuple->column("N_PrHit_Assoc",track.size()); @@ -958,10 +990,10 @@ void PrClustersResidual::ClusterResidual(){ //LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward); //LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed); if(msgLevel(MSG::DEBUG)) debug()<<" Get ClusterLinker "<<endmsg; - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); if(msgLevel(MSG::DEBUG)) debug()<<" Get Cluster to MCHit Linker "<<endmsg; - - LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits"); + + LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits"); MCTrackInfo trackInfo( evtSvc(), msgSvc() ); if(msgLevel(MSG::DEBUG)) debug()<<" Add MCHitAndTrackStudy "<<endmsg; Tuples::Tuple tuple = GaudiTupleAlg::nTuple("ClusterMCHitAndTrackStudy","Events"); @@ -1032,6 +1064,8 @@ void PrClustersResidual::ClusterResidual(){ if(msgLevel(MSG::DEBUG)) debug()<<"Loaded"<<endmsg; //took first MCHit associated to Hit channelID std::vector<double> MCHit_X,MCHit_Y,MCHit_Z; + std::vector<double> MCHit_EntryX, MCHit_EntryY, MCHit_EntryZ; + std::vector<double> MCHit_ExitX, MCHit_ExitY, MCHit_ExitZ; std::vector<double> Residual_X,Residual_Z; std::vector<double> SlopeX; std::vector<double> SlopeY; @@ -1047,15 +1081,25 @@ void PrClustersResidual::ClusterResidual(){ MCHit_X.push_back(-9999.); MCHit_Y.push_back(-9999.); MCHit_Z.push_back(-9999.); - SlopeX.push_back(-999.); - SlopeY.push_back(-999.); - pathlength.push_back(-999.); - energy.push_back( -999.); + MCHit_EntryX.push_back(-9999.); + MCHit_EntryY.push_back(-9999.); + MCHit_EntryZ.push_back(-9999.); + + MCHit_ExitX.push_back(-9999.); + MCHit_ExitY.push_back(-9999.); + MCHit_ExitZ.push_back(-9999.); + + SlopeX.push_back(-9999.); + SlopeY.push_back(-9999.); + pathlength.push_back(-9999.); + energy.push_back( -9999.); hit_Physics.push_back(false); } while(mcHit != nullptr){ nMCHitToCluster++; - Gaudi::XYZPoint pMid = mcHit->midPoint(); + auto pMid = mcHit->midPoint(); + auto pEntry = mcHit->entry(); + auto pExit = mcHit->exit(); Residual_X.push_back( (*iHit).x(pMid.y()) - pMid.x()); //x at y - y at z Residual_Z.push_back( (*iHit).z(pMid.y()) - pMid.z()); SlopeX.push_back( mcHit->dxdz()); @@ -1065,6 +1109,13 @@ void PrClustersResidual::ClusterResidual(){ MCHit_X.push_back( pMid.x()); MCHit_Y.push_back( pMid.y()); MCHit_Z.push_back( pMid.z()); + + MCHit_ExitX.push_back( pExit.x()); + MCHit_ExitY.push_back( pExit.y()); + MCHit_ExitZ.push_back( pExit.z()); + MCHit_EntryX.push_back( pEntry.x()); + MCHit_EntryY.push_back( pEntry.y()); + MCHit_EntryZ.push_back( pEntry.z()); //mcHit = myFTCluster2MCHitLink.next(); // const LHCb::MCParticle *mcPart = mcHit->mcParticle(); const LHCb::MCParticle *mcPart = mcHit->mcParticle(); @@ -1074,6 +1125,7 @@ void PrClustersResidual::ClusterResidual(){ if(mcPart!= nullptr){ if(msgLevel(MSG::DEBUG)) debug() <<"Getting Vertex"<<endmsg; const LHCb::MCVertex* vertex = mcPart->originVertex(); + //track of physical interest if the associated vertex is a Primary one or from Decay ( i.e. veto the interaction one ) if( vertex!=nullptr && (vertex->isPrimary() || vertex->isDecay())){ hit_Physics.push_back( true ); }else{ @@ -1084,7 +1136,6 @@ void PrClustersResidual::ClusterResidual(){ } if(msgLevel(MSG::DEBUG)) debug()<<"Storing vectors"<<endmsg; tuple->column("isNoiseCluster",isNoise); - tuple->farray("MCHit_Physics", hit_Physics.begin(),hit_Physics.end(),"N_MC",100); tuple->farray("CosSlopeX",SlopeX.begin(),SlopeX.end(),"N_MC",100); tuple->farray("CosSlopeY",SlopeY.begin(),SlopeY.end(),"N_MC",100); @@ -1095,17 +1146,26 @@ void PrClustersResidual::ClusterResidual(){ tuple->farray("MCHit_X",MCHit_X.begin(),MCHit_X.end(),"N_MC",100); tuple->farray("MCHit_Z",MCHit_Z.begin(),MCHit_Z.end(),"N_MC",100); tuple->farray("MCHit_Y",MCHit_Y.begin(),MCHit_Y.end(),"N_MC",100); + + tuple->farray("MCHit_EntryX",MCHit_EntryX.begin(),MCHit_EntryX.end(),"N_MC",100); + tuple->farray("MCHit_EntryZ",MCHit_EntryZ.begin(),MCHit_EntryZ.end(),"N_MC",100); + tuple->farray("MCHit_EntryY",MCHit_EntryY.begin(),MCHit_EntryY.end(),"N_MC",100); + tuple->farray("MCHit_ExitX",MCHit_ExitX.begin(),MCHit_ExitX.end(),"N_MC",100); + tuple->farray("MCHit_ExitZ",MCHit_ExitZ.begin(),MCHit_ExitZ.end(),"N_MC",100); + tuple->farray("MCHit_ExitY",MCHit_ExitY.begin(),MCHit_ExitY.end(),"N_MC",100); + + //here fill the tuple with watever you want //get Cluster associated to the Hit if(msgLevel(MSG::DEBUG)) debug()<<"Getting LiteCluster from PrHit"<<endmsg; LHCb::FTLiteCluster litecluster = getLiteCluster( (*iHit).id()); - tuple->column("ClusterCharge",litecluster.charge()); - tuple->column("ClusterSize",litecluster.size()); + tuple->column("ClusterSize",litecluster.isLarge()); tuple->column("ClusterFraction",litecluster.fraction()); tuple->column("ClusterChannelID",litecluster.channelID()); - tuple->column("ClusterChannelIDSipmCell",litecluster.channelID().sipmCell()); - tuple->column("ClusterChannelIDSipmID",litecluster.channelID().sipmId()); + tuple->column("ClusterChannelIDSipmCell",litecluster.channelID().channel()); + tuple->column("ClusterChannelIDSipmID",litecluster.channelID().sipm()); tuple->column("ClusterChannelIDMat",litecluster.channelID().mat()); + tuple->column("ClusterChannelIDStation",litecluster.channelID().station()); tuple->column("ClusterChannelIDModule",litecluster.channelID().module()); tuple->column("ClusterChannelLayer",litecluster.channelID().layer()); tuple->column("ClusterChannelQuarter",litecluster.channelID().quarter()); @@ -1152,8 +1212,8 @@ void PrClustersResidual::ClusterResidual(){ tuple->column("Hit_Zone",(*iHit).zone()); tuple->column("Hit_dzDy_manually",((*iHit).z(1000.)-(*iHit).z(0.))/1000); //Get the list of tracks generated by the Hit i am lookign to - if(m_debugTracking) - { + if(m_debugTracking){ + //get linker Forward tracks and Seeding ones to see hit content! LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward); LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed); @@ -1311,3 +1371,4 @@ LHCb::FTLiteCluster PrClustersResidual::getLiteCluster(const LHCb::LHCbID id) return idTracks; } + diff --git a/Pr/PrMCTools/src/PrClustersResidual.h b/Pr/PrMCTools/src/PrClustersResidual.h index 0c6c12778d4656d1cc176c85f6927e87be051a29..6c48866a672273afd88ab593cacfb9d8dda1680a 100644 --- a/Pr/PrMCTools/src/PrClustersResidual.h +++ b/Pr/PrMCTools/src/PrClustersResidual.h @@ -23,7 +23,7 @@ #include "Event/Track.h" #include "Linker/LinkerTool.h" #include "PrKernel/IPrCounter.h" - +#include "FTDet/DeFTDetector.h" /** @class PrClustersResidual PrClustersResidual.h * Make nTuples of Cluster vs MCHits and Hit content on track + Occupancy study possible ( also in PrPlotFTHits ) * -MCHitsLocation : Where to find the MCHits ( /Event/MC/FT/FTHits by default ) @@ -103,37 +103,28 @@ private: std::vector<const LHCb::Track*> getTrack(const LHCb::LHCbID id, const std::string location); LHCb::FTLiteCluster getLiteCluster(const LHCb::LHCbID id); - std::string m_mcHitLocation; - // unsigned int m_zone; - bool m_debugTracking; + Gaudi::Property<std::string> m_mcHitLocation {this, "MCHitsLocation", "/Event/MC/FT/Hits"}; + Gaudi::Property<std::string> m_ftClusterToParticleLinkerLocation { this, "FTClusterToParticleLinkerLocation", LHCb::FTLiteClusterLocation::Default + "WithSpillover"}; + Gaudi::Property<std::string> m_ftClusterToHitLinkerLocation { this, "FTClusterToHitLinkerLocation", LHCb::FTLiteClusterLocation::Default + "2MCHitsWithSpillover"}; + Gaudi::Property<std::string> m_hitManagerName { this, "HitManagerName", "PrFTHitManager"}; + Gaudi::Property<bool> m_debugTracking { this, "DebugTracking", false}; + Gaudi::Property<bool> m_doClusterResidual { this, "DoClusterResidual", false}; //Cluster residual tuple + Gaudi::Property<bool> m_doTrackStudy { this, "DoTrackStudy", false}; // Produce tuple for track studies + Gaudi::Property<bool> m_decodeData { this, "DecodeData", false}; //ask to decode data ( is False if other algorithms are runned + Gaudi::Property<bool> m_onlyHasT { this, "OnlyHasT",false}; //(keep only MCParticles with hasT = true + Gaudi::Property<bool> m_removeClones { this, "RemoveClones", true}; //Observed clones MCHits when processed ( not understood ) + Gaudi::Property<bool> m_dumpAllHits { this, "DumpAllHits", false }; //Dump in tuple for a given event all the hits + Gaudi::Property<bool> m_Occupancy { this, "Occupancy", true };// Produce Occupancy plot int m_nEvents; - double m_NClusters ; double m_NMCHit ; double m_NMCHit_inClusters; - bool m_onlyHasT; PrHitManager* m_ftHitManager; unsigned int m_zone; - bool m_decodeData; - bool m_removeClones; - bool m_dumpAllHits; - bool m_Occupancy; - std::string m_hitManagerName; - bool m_doClusterResidual; - bool m_doTrackStudy; - - //from PrCounter - // typedef LinkerTool<LHCb::Track, LHCb::MCParticle> MyAsct; - // typedef MyAsct::DirectType Table; - // typedef MyAsct::DirectType::Range Range; - // typedef Table::iterator iterator; - // typedef MyAsct::InverseType InvTable; - // typedef InvTable::Range InvRange; - // typedef InvTable::iterator InvIterator; - // MyAsct* m_link; - // const InvTable* m_invTable; + //need the ftDet to produce the pseudochannel conversion numbering in occupancy plot for v61 + DeFTDetector * m_ftDet = nullptr; }; #endif // PRCLUSTERSRESIDUAL_H diff --git a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp index d7bf33cc1ec4d9e84208f6d7cc8764b3f5c948d3..9b979457bf71e5e91bee30e32dc005cf613ee132 100644 --- a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp +++ b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp @@ -218,7 +218,7 @@ StatusCode PrDebugTrackingLosses::execute() { } else if ( (*itId).isFT() ) { LHCb::FTChannelID ftID = (*itId).ftID(); info() << format( " FT St%2d La%2d Pm%2d Cel%4d ", - ftID.layer()/4, ftID.layer()&3, ftID.sipmId(), ftID.sipmCell() ); + ftID.station(), ftID.layer(), ftID.sipm(), ftID.channel() ); } else if ( (*itId).isOT() ) { LHCb::OTChannelID otID = (*itId).otID(); info() << format( " OT St%2d La%2d mo%2d Str%4d ", @@ -285,3 +285,4 @@ StatusCode PrDebugTrackingLosses::finalize() { } //============================================================================= + diff --git a/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp b/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp index ef1d3220cfbd53cca48e54eaec2a99b0552e548b..a7dc95b7e0c07b86a3b9c67f14803f8979b5d8da 100755 --- a/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp +++ b/Pr/PrMCTools/src/PrLHCbID2MCParticle.cpp @@ -8,7 +8,6 @@ #include "Event/VPCluster.h" #include "Event/STCluster.h" #include "Event/OTTime.h" -#include "Event/FTCluster.h" #include "Event/FTLiteCluster.h" // Range V3 @@ -137,7 +136,7 @@ StatusCode PrLHCbID2MCParticle::execute() { //== FT const auto* ft_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int>>( LHCb::FTLiteClusterLocation::Default ); if ( ft_clusters ) { - LinkedTo<LHCb::MCParticle> ftLink( evtSvc(), msgSvc(),LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCParticle> ftLink( evtSvc(), msgSvc(),LHCb::FTLiteClusterLocation::Default ); if ( !ftLink.notFound() ) { for(const auto& cluster : *ft_clusters) { linkAll( ftLink, lhcbLink, cluster.channelID(), {cluster.channelID()}); diff --git a/Pr/PrMCTools/src/PrPlotFTHits.cpp b/Pr/PrMCTools/src/PrPlotFTHits.cpp index faff722e190cdc4ae22a4ad3a1ac33a07b90b9f0..974d433549740622a9f8d2e9f9145caa6abf54e2 100644 --- a/Pr/PrMCTools/src/PrPlotFTHits.cpp +++ b/Pr/PrMCTools/src/PrPlotFTHits.cpp @@ -114,17 +114,17 @@ void PrPlotFTHits::plotOccupancy(){ noSpill = true; } //"/Event/Prev/MC/FT/Hits" - LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" ); - - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); - LinkedTo<LHCb::MCHit, LHCb::FTCluster> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits"); - //"/Event/Prev/MC/FT/Hits" - char layerName[100]; - char Title[100]; - std::vector<int> nHits(12,0); - MCTrackInfo trackInfo( evtSvc(), msgSvc() ); - - LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default ); + LinkedFrom<LHCb::MCHit,LHCb::MCParticle> myMCHitLink( evtSvc(), msgSvc(), LHCb::MCParticleLocation::Default + "2MC" + "FT" + "Hits" ); + + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); + LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits"); + //"/Event/Prev/MC/FT/Hits" + char layerName[100]; + char Title[100]; + std::vector<int> nHits(12,0); + MCTrackInfo trackInfo( evtSvc(), msgSvc() ); + + LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default ); if( mcVert == nullptr ){ error() << "Could not find MCVertices at " << LHCb::MCParticleLocation::Default << endmsg; } @@ -164,19 +164,7 @@ void PrPlotFTHits::plotOccupancy(){ hit = &*itH; nHits[layer]++; LHCb::FTLiteCluster liteCluster = getLiteCluster(hit->id()); - int Module = liteCluster.channelID().module(); - int Quarter = liteCluster.channelID().quarter(); - int SipmID = liteCluster.channelID().sipmId(); - int layer = liteCluster.channelID().layer(); - int module2 = Module; - if(Module == 11 && Module > 4 ) module2 = 1; - if(Module != 11 && Module > 4 ) module2 = Module-3; - if(Module < 5 ) module2 = 6-Module; - if(Module == 10 ) module2 = 1; - int index = Quarter*128 + 16*module2+ std::floor(SipmID); - // if(index<0 || index>){ - // always()<<"Problem in Sipm Index; Correct for that"<<endmsg; - // } + int index = liteCluster.channelID().uniqueSiPM() & 511; sprintf(layerName,"Layer%i/AllContributions",layer); sprintf(Title,"Sipm Index Layer %i;SiPm Id;Clusters/Sipm/Event",layer); m_histoTool->plot1D(index,layerName,Title,0.5,512.5,512); @@ -352,7 +340,7 @@ void PrPlotFTHits::plotHitEfficiency(){ // -- Link different track containers to MCParticles LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward); LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed); - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); MCTrackInfo trackInfo( evtSvc(), msgSvc() ); // ----------------------------------------------------------------------------------- @@ -739,8 +727,8 @@ void PrPlotFTHits::plotMCHits () { // -- MC linking - LinkedTo<LHCb::MCHit, LHCb::FTCluster> myMCHitLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default + "2MCHits"); - LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCHit> myMCHitLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits"); + LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); MCTrackInfo trackInfo( evtSvc(), msgSvc() ); diff --git a/Pr/PrMCTools/src/PrTStationDebugTool.cpp b/Pr/PrMCTools/src/PrTStationDebugTool.cpp index 304fbea5b1c195731ea92e7fe7f7c34cb29660df..67a93429c5be0e6b59f122e5f4e7a043016765e3 100644 --- a/Pr/PrMCTools/src/PrTStationDebugTool.cpp +++ b/Pr/PrMCTools/src/PrTStationDebugTool.cpp @@ -5,7 +5,7 @@ #include "Linker/LinkedTo.h" #include "Kernel/LHCbID.h" #include "Event/OTTime.h" -#include "Event/FTCluster.h" +#include "Event/FTLiteCluster.h" #include "Event/MCParticle.h" // local @@ -50,8 +50,8 @@ bool PrTStationDebugTool::matchKey ( LHCb::LHCbID id, int key ) { if ( key == part->key() ) return true; part = oLink.next(); } - } else if ( id.isFT() ) { - LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + } else if ( id.isFT() ) { + LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); LHCb::FTChannelID idO = id.ftID(); LHCb::MCParticle* part = fLink.first( idO ); @@ -77,7 +77,7 @@ void PrTStationDebugTool::printKey( MsgStream& msg, LHCb::LHCbID id ) { part = oLink.next(); } } else if ( id.isFT() ) { - LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default ); + LinkedTo<LHCb::MCParticle> fLink( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default ); LHCb::FTChannelID idO = id.ftID(); LHCb::MCParticle* part = fLink.first( idO ); diff --git a/Rec/GlobalReco/root/RichCorrelations.C b/Rec/GlobalReco/root/RichCorrelations.C new file mode 100644 index 0000000000000000000000000000000000000000..96d0e6757597b3b07360a42106c1b858302f2851 --- /dev/null +++ b/Rec/GlobalReco/root/RichCorrelations.C @@ -0,0 +1,66 @@ + +#include <memory> + +void RichCorrelations() +{ + + auto c = std::make_unique<TCanvas>( "RichCor", "RichCor", 1400, 1000 ); + + TTree tree; + + tree.AddFriend("hlt2 = ChargedProtoTuple/protoPtuple", + "/usera/jonesc/LHCbCMake/Brunel/output/future/Quartic100R1R2-AllHypos/protoparticles.tuples.root"); + + //tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/FastQuartic50R1R2-AllHypos-NoUnAmBigPhot-NoBeamPipeCheck-1It/protoparticles.tuples.root"); + //tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/CKEsti50R1R2-KaPi/protoparticles.tuples.root"); + tree.AddFriend("hlt1 = ChargedProtoTuple/protoPtuple", "/usera/jonesc/LHCbCMake/Brunel/output/future/CKEsti50R1R2-AllHypos/protoparticles.tuples.root"); + + tree.Draw("hlt2.TrackP:hlt1.TrackP","hlt2.TrackP<100000","zcol"); + c->SaveAs( "P-Correlations.pdf" ); + + auto sel = TCut("hlt2.RichAbovePiThres>0 && hlt1.RichAbovePiThres>0 && hlt2.TrackP<100000 && hlt2.TrackP>3000"); + + auto trueK = TCut("abs(hlt2.MCParticleType) == 321"); + + auto truePi = TCut("abs(hlt2.MCParticleType) == 211"); + + tree.Draw( "hlt2.RichDLLk:hlt1.RichDLLk",sel && trueK,"zcol"); + c->SaveAs( "TrueK-RichDLLk-2dCorrelations.pdf" ); + + tree.Draw( "hlt2.RichDLLk-hlt1.RichDLLk",sel && trueK ); + c->SaveAs( "TrueK-RichDLLk-1dCorrelations.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel ); + c->SaveAs( "All-RichDLLk-HLT1.pdf" ); + + tree.Draw( "hlt2.RichDLLk", sel ); + c->SaveAs( "All-RichDLLk-HLT2.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK ); + c->SaveAs( "TrueK-RichDLLk-HLT1.pdf" ); + + tree.Draw( "hlt2.RichDLLk", sel && trueK ); + c->SaveAs( "TrueK-RichDLLk-HLT2.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && truePi ); + c->SaveAs( "TruePi-RichDLLk-HLT1.pdf" ); + + tree.Draw( "hlt2.RichDLLk", sel && truePi ); + c->SaveAs( "TruePi-RichDLLk-HLT2.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>-2") ); + c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt-2.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>0") ); + c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt0.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>2") ); + c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt2.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>4") ); + c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt4.pdf" ); + + tree.Draw( "hlt1.RichDLLk", sel && trueK && TCut("hlt2.RichDLLk>10") ); + c->SaveAs( "TrueK-RichDLLk-HLT1-HLT2gt10.pdf" ); + +} diff --git a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C index 0d62054025282aa02a5d089761c48f3c21f8dbfd..abd6ca3a8a1b9a5e0e1e27f70a5d0edcd8443fcc 100755 --- a/Rec/GlobalReco/root/RichKaonIDCompareFiles.C +++ b/Rec/GlobalReco/root/RichKaonIDCompareFiles.C @@ -51,14 +51,14 @@ void RichKaonIDCompareFiles() //const std::string dir = "/Users/chris/LHCb/RootFiles/Run2"; //const std::string dir = "/usera/jonesc/LHCbCMake/lhcb-head/BrunelDevNightly/output/Run2"; - const std::string dir = "/usera/jonesc/LHCbCMake/Brunel/output/future"; + const std::string dir = "/home/chris/LHCb/future"; typedef std::vector< std::tuple<std::string,std::string,Color_t> > PlotData; const PlotData plotdata = { - std::make_tuple ( "rich", "RICH Current", kBlack ), - std::make_tuple ( "richfuture", "RICH Future", kRed-6 ) + std::make_tuple ( "old-quartic", "Old Quartic", kBlack ), + std::make_tuple ( "new-quartic", "New NR Quartic", kRed-6 ) }; // const PlotData plotdata = diff --git a/RecSys/CMakeLists.txt b/RecSys/CMakeLists.txt index 939db4b274ec0abcb564fb5bc7974532298c4d12..fce23855701a083474113c53348326665416e13c 100644 --- a/RecSys/CMakeLists.txt +++ b/RecSys/CMakeLists.txt @@ -60,6 +60,7 @@ gaudi_depends_on_subdirs(Calo/CaloMoniDst Rich/RichFutureRecSys Rich/RichFutureRecMonitors Rich/RichFutureRecCheckers + Rich/RichRecTests Tf/FastPV Tf/FastVelo Tf/PatAlgorithms diff --git a/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py b/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py index a051222bc35619bb9125deb68eb58812e1aa55a5..af38d34c09b81b0ed0d2d77f80439029c7e3e7ed 100644 --- a/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py +++ b/Rich/RichAlignment/scripts/RefractAndHPDJobs/CommonOptions.py @@ -31,7 +31,7 @@ importOptions("$L0TCK/L0DUConfig.opts") # For 2016 data importOptions("$APPCONFIGOPTS/Brunel/DataType-2016.py") LHCbApp().DDDBtag = "dddb-20150724" -LHCbApp().CondDBtag = "cond-20160517" +LHCbApp().CondDBtag = "cond-20170120-1" # Aerogel Sub Tiles #CondDB().LocalTags["LHCBCOND"] = ["rich1-20110624"] diff --git a/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py b/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py index 516bc597b7f174d8105309b277f845068c01b294..ef223eb1d1e8c9a3b2b8e52863cff120317466aa 100755 --- a/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py +++ b/Rich/RichAlignment/scripts/RefractAndHPDJobs/GetLFNsByRun.py @@ -42,6 +42,8 @@ elif 2013 == year : ConfigV = ['Collision13','Protonion13','Ionproton13','Ionsmog'] elif 2015 == year : ConfigV = ['Collision15','Collision15em','Protonhelium15','Protonneon15','Protonargon15','Lead15'] +elif 2016 == year : + ConfigV = ['Collision16'] else: print 'Unknown year', year DIRAC.exit(2) @@ -51,6 +53,10 @@ allOK = True # Dictionary for LFNs by Run RunLFNs = { } +polarity = "MagUp" +#polarity = "MagDown" +#polarity = "All" + startDate = getDate(year,startmonth,startday) endDate = getDate(year,lastmonth,lastday) dict = { 'StartDate' : startDate, @@ -81,80 +87,106 @@ else: for run in runs : if maxrun == -1 or run <= maxrun : - + nRun += 1 - for config in ConfigV : - - type = 91000000 # EXPRESS Stream - if year == 2009 : type = 90000000 # Use full stream for 2009 - if run > 71473 and run < 72332 : # Use FULL for first 14nb-1 in 2010 - type = 90000000 # FULL Stream - if run > 77595 and run < 77624 : # Express disappeared for unknown reasons - type = 90000000 # FULL Stream - if run > 100256 and run < 102177 : # Express turned off by accident after 09/2011 TS - type = 90000000 # FULL Stream - if config == 'Collision11_25' : # No express for 2011 25ns tests - type = 90000000 # FULL Stream - if config == 'Protonion12' : # (Currently) no express stream for pA data - type = 90000000 # FULL Stream - if config == 'Protonion13' or config == 'Ionproton13' : - type = 90000000 # FULL Stream - #if run > 77595 and run < 77624 : # Express disappeared for unknown reasons - if run > 157596 : # Express stream turned off in Run II - type = 90000000 # FULL Stream - - typeS = "EXPRESS" - if type == 90000000 : typeS = "FULL" - - # Processing pass. For 2011 express stream data used a Merged location - procpass = 'Real Data' - if year == 2011 and type != 90000000 : - procpass = 'Real Data/Merging' - - # Raw files - bkDict = { 'ConfigName' : 'LHCb', - 'ConfigVersion' : config, - 'ProcessingPass' : procpass, - 'FileType' : 'ALL', - 'StartRun' : run, - 'EndRun' : run, - 'EventType' : type } - resultB = database.getFilesWithGivenDataSets(bkDict) - - if not resultB['OK']: - print "Error querying BK with", bkDict, resultB['Message'] - allOK = False - else: - tmpLFNList = [ ] - for lfn in resultB['Value']: - filetype = lfn.split('.')[1] - if filetype == 'dst' or filetype == 'raw': - tmpLFNList += ["LFN:"+lfn.encode("ascii")] - tmpLFNList.sort() - if len(tmpLFNList) > 0 : - if run not in RunLFNs.keys() : RunLFNs[run] = [ ] - RunLFNs[run] += tmpLFNList - - if not allOK : break + # Polarity ? + res = database.getRunInformations(int(run)) + if polarity == "All" or polarity in res['Value']['DataTakingDescription'] : + + for config in ConfigV : + + type = 91000000 # EXPRESS Stream + if year == 2009 : type = 90000000 # Use full stream for 2009 + if run > 71473 and run < 72332 : # Use FULL for first 14nb-1 in 2010 + type = 90000000 # FULL Stream + if run > 77595 and run < 77624 : # Express disappeared for unknown reasons + type = 90000000 # FULL Stream + if run > 100256 and run < 102177 : # Express turned off by accident after 09/2011 TS + type = 90000000 # FULL Stream + if config == 'Collision11_25' : # No express for 2011 25ns tests + type = 90000000 # FULL Stream + if config == 'Protonion12' : # (Currently) no express stream for pA data + type = 90000000 # FULL Stream + if config == 'Protonion13' or config == 'Ionproton13' : + type = 90000000 # FULL Stream + #if run > 77595 and run < 77624 : # Express disappeared for unknown reasons + if run > 157596 : # Express stream turned off in Run II + type = 90000000 # FULL Stream + if config == 'Collision16' : + type = 90000000 # FULL Stream + + typeS = "EXPRESS" + if type == 90000000 : typeS = "FULL" + + # Processing pass. For 2011 express stream data used a Merged location + procpass = 'Real Data' + if year == 2011 and type != 90000000 : + procpass = 'Real Data/Merging' + + # Raw files + bkDict = { 'ConfigName' : 'LHCb', + 'ConfigVersion' : config, + 'ProcessingPass' : procpass, + 'FileType' : 'ALL', + 'StartRun' : run, + 'EndRun' : run, + 'EventType' : type } + resultB = database.getFilesWithGivenDataSets(bkDict) + + if not resultB['OK']: + + print "Error querying BK with", bkDict, resultB['Message'] + allOK = False + + else: + + tmpLFNList = [ ] - nLFNs = 0 - if run in RunLFNs.keys() : nLFNs = len(RunLFNs[run]) - if nLFNs > 0 : - print "(", nRun, "of", len(runs), ") Selected", nLFNs, "LFNs for run", run, ConfigV - print RunLFNs[run] - else: - print "(", nRun, "of", len(runs), ") No data selected for run", run, ConfigV + lfns = resultB['Value'] + tmplfns = list( set(lfns) ) + if len(tmplfns) != len(lfns) : + print "WARNING : run", run, "has duplicate LFNs" + print bkDict + print resultB['Value'] + lfns = tmplfns + + for lfn in lfns: + filetype = lfn.split('.')[1] + if filetype == 'dst' or filetype == 'raw': + tmpLFNList += ["LFN:"+lfn.encode("ascii")] + + tmpLFNList.sort() + if len(tmpLFNList) > 0 : + if run not in RunLFNs.keys() : RunLFNs[run] = [ ] + RunLFNs[run] += tmpLFNList + + if not allOK : break + + nLFNs = 0 + if run in RunLFNs.keys() : nLFNs = len(RunLFNs[run]) + if nLFNs > 0 : + print "(", nRun, "of", len(runs), ") Selected", nLFNs, "LFNs for run", run, ConfigV + #print RunLFNs[run] + else: + print "(", nRun, "of", len(runs), ") No data selected for run", run, ConfigV - else: + else : + + print "(", nRun, "of", len(runs), ") Run", run, "has wrong polarity" + else: + print "Skipping run", run if allOK : # Pickle the data if len(RunLFNs.keys()) > 0: - filename = "RunData/RunLFNs_"+startDate+"_"+endDate+".pck.bz2" + if polarity == "All" : + filename = "RunData/RunLFNs_"+startDate+"_"+endDate+".pck.bz2" + else : + filename = "RunData/RunLFNs_"+polarity+"_"+startDate+"_"+endDate+".pck.bz2" file = bz2.BZ2File(filename,"w") pickle.dump(RunLFNs,file) file.close() diff --git a/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h b/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h new file mode 100755 index 0000000000000000000000000000000000000000..d28f230cf03c26afe94d0b2a15302dd4491553c1 --- /dev/null +++ b/Rich/RichFutureGlobalPID/src/RichGlobalPIDNamespaces.h @@ -0,0 +1,20 @@ + +//----------------------------------------------------------------------------- +/** @namespace Rich::Future::Rec::GlobalPID + * + * Namespace for RICH Global PID software + * + * @author Chris Jones Christopher.Rob.Jones@cern.ch + * @date 04/12/2006 + */ +//----------------------------------------------------------------------------- + +//----------------------------------------------------------------------------- +/** @namespace Rich::Future::Rec::GlobalPID::MC + * + * Namespace for RICH Global PID MC related software + * + * @author Chris Jones Christopher.Rob.Jones@cern.ch + * @date 05/12/2006 + */ +//----------------------------------------------------------------------------- diff --git a/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h b/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h index c7472cd147a6058b29182d4282d99fc7b9d91dc8..8443cac3c582fd1d01833b22cc7ce3059a0fd789 100644 --- a/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h +++ b/Rich/RichFutureRecBase/RichFutureRecBase/RichRecNamespaces.h @@ -5,7 +5,7 @@ * General namespace for all RICH reconstruction software * * @author Chris Jones Christopher.Rob.Jones@cern.ch - * @date 08/07/2004 + * @date 01/02/2017 */ //----------------------------------------------------------------------------- @@ -15,6 +15,16 @@ * General namespace for RICH reconstruction MC related software * * @author Chris Jones Christopher.Rob.Jones@cern.ch - * @date 05/12/2006 + * @date 01/02/2017 + */ +//----------------------------------------------------------------------------- + +//----------------------------------------------------------------------------- +/** @namespace Rich::Future::Rec::Moni + * + * General namespace for RICH reconstruction monitoring utilities + * + * @author Chris Jones Christopher.Rob.Jones@cern.ch + * @date 01/02/2017 */ //----------------------------------------------------------------------------- diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h index 45351b5eed6973d666964a669290f0abf8d9f578..d7391683b51c73dc3398edb6b5c3b9f741c4fc8c 100644 --- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h +++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichRecPhotonPredictedPixelSignals.h @@ -24,6 +24,7 @@ namespace Rich /// Type for flags to say if a photon is active or not. i.e. has any signal. using PhotonFlags = LHCb::STL::Vector<bool>; + /// photon active flags TES locations namespace PhotonActiveFlagsLocation { /// Location in TES for the active photon flags diff --git a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h index 0c5854a24002c19c9d34e5bf6355306f563f6ce9..2c322afa970636a85e300d04a42c2efb5f72200b 100644 --- a/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h +++ b/Rich/RichFutureRecEvent/RichFutureRecEvent/RichSummaryEventData.h @@ -24,6 +24,14 @@ namespace Rich { // =========================================================================== + /** @namespace Summary + * + * General namespace for RICH reconstruction summary information + * + * @author Chris Jones Christopher.Rob.Jones@cern.ch + * @date 01/02/2017 + */ + namespace Summary { @@ -125,6 +133,14 @@ namespace Rich // =========================================================================== + /** @class Pixel RichSummaryEventData.h + * + * Summary of the immutable pixell data. + * + * @author Chris Jones + * @date 2016-09-30 + */ + class Pixel final { @@ -155,7 +171,7 @@ namespace Rich // =========================================================================== - /// TES locations + /// Summary information TES locations namespace TESLocations { /// Tracks diff --git a/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt b/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt index 5ccf26d32ef7c847a349b75e505bf3d92e17bd4d..85678715d1a9163c5a86e75d796e0601aea34a63 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt +++ b/Rich/RichFutureRecPhotonAlgorithms/CMakeLists.txt @@ -5,7 +5,8 @@ gaudi_subdir(RichFutureRecPhotonAlgorithms v1r0) gaudi_depends_on_subdirs(Rich/RichFutureRecBase Rich/RichFutureRecEvent - Rich/RichUtils) + Rich/RichUtils + Kernel/VectorClass) find_package(Boost) find_package(GSL) @@ -16,10 +17,23 @@ find_package(ROOT) include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}) -# For testing high vectorisation and optimisation levels -#set_property(SOURCE src/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -msse4 -O3 " ) - gaudi_add_module(RichFutureRecPhotonAlgorithms - src/*.cpp - INCLUDE_DIRS Boost GSL VDT Eigen Rich/RichFutureRecBase Rich/RichFutureRecEvent Rich/RichUtils + src/*.cpp src/generic/*.cpp src/sse4/*.cpp src/avx/*.cpp src/avx2/*.cpp + INCLUDE_DIRS Boost GSL VDT Eigen Kernel/VectorClass Rich/RichFutureRecBase Rich/RichFutureRecEvent Rich/RichUtils LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichFutureRecEvent RichUtils) + +set_property(SOURCE src/generic/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes " ) + +set_property(SOURCE src/sse4/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -msse4.2 " ) + +set_property(SOURCE src/avx/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx " ) + +# Until such a time that all build targets support AVX2+FMA, fallback to just AVX here +#set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx2 -mfma " ) +set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -O3 -mavx " ) + +if(LCG_COMP STREQUAL "gcc" OR + (BINARY_TAG_COMP_NAME STREQUAL "gcc" AND BINARY_TAG_COMP_VERSION VERSION_LESS "5")) + set_property(SOURCE src/avx/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/avx2/RichQuarticPhotonReco.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) +endif() diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h index ff0522a36413bd0cf966d4cc76528d816e0c6b94..09fadbbfe24f72ed61239dc1587e190d6cc7aeeb 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichBasePhotonReco.h @@ -166,18 +166,18 @@ namespace Rich } /// Check the final Cherenkov angles - template < typename HYPODATA > + template < typename HYPODATA, typename FTYPE > inline bool checkAngles( const Rich::RadiatorType rad, const HYPODATA& tkCkAngles, const HYPODATA& tkCkRes, - const float ckTheta, - const float ckPhi ) const noexcept + const FTYPE ckTheta, + const FTYPE ckPhi ) const noexcept { return ( // First the basic checks - ckPhi > 0 && - ckTheta > absMinCKTheta(rad) && ckTheta < absMaxCKTheta(rad) && + ckTheta > absMinCKTheta(rad) && + ckPhi > 0 && // Now check each hypo std::any_of( activeParticles().begin(), activeParticles().end(), diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp index 595911486692d9fc8fdac5ed3215862d017de142..14f6b40195c42be4b526f039f9ff8c0c58e7dc09 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.cpp @@ -47,6 +47,18 @@ StatusCode QuarticPhotonReco::initialize() auto sc = MultiTransformer::initialize(); if ( !sc ) return sc; + // Get the CPU capabilities and set dispatch method + if ( m_detectCPU ) + { + enum CPUID { GENERIC = 0, SSE4 = 6, AVX = 7, AVX2 = 8 }; + const auto cpuLevel = instrset_detect(); + _ri_debug << "Instruction set level = " << cpuLevel << endmsg; + if ( cpuLevel >= AVX2 ) { m_run = &QuarticPhotonReco::run_avx2; } + else if ( cpuLevel >= AVX ) { m_run = &QuarticPhotonReco::run_avx; } + else if ( cpuLevel >= SSE4 ) { m_run = &QuarticPhotonReco::run_sse4; } + else { m_run = &QuarticPhotonReco::run_generic; } + } + // get the detector elements m_rich[Rich::Rich1] = getDet<DeRich>( DeRichLocations::Rich1 ); m_rich[Rich::Rich2] = getDet<DeRich>( DeRichLocations::Rich2 ); @@ -112,678 +124,6 @@ StatusCode QuarticPhotonReco::initialize() //============================================================================= -OutData -QuarticPhotonReco::operator()( const LHCb::RichTrackSegment::Vector& segments, - const CherenkovAngles::Vector& ckAngles, - const CherenkovResolutions::Vector& ckResolutions, - const SegmentPanelSpacePoints::Vector& trHitPntsLoc, - const SegmentPhotonFlags::Vector& segPhotFlags, - const Rich::PDPixelCluster::Vector& clusters, - const SpacePointVector& globalPoints, - const SpacePointVector& localPoints, - const Relations::TrackToSegments::Vector& tkToSegRels ) const -{ - // use ranges v3 library - using namespace ranges::v3; - - // make the output data - OutData outData; - - // Shortcut to the photons and relations - auto & photons = std::get<0>(outData); - auto & relations = std::get<1>(outData); - // guess at reserve size - const auto resSize = segments.size() * globalPoints.size() / 20; - photons.reserve( resSize ); - relations.reserve( resSize ); - - // View range for pixel data - const auto allPixels = Ranges::ConstZip( clusters, globalPoints, localPoints ); - - // Make 'ranges' for RICH pixels - const auto RPixs = richRanges( allPixels ); - const auto & RichRanges = RPixs[0]; // Complete RICH ranges - const auto & Rich1Ranges = RPixs[1]; // ranges for RICH1 [top,bottom] - const auto & Rich2Ranges = RPixs[2]; // ranges for RICH2 [left,right] - - // local position corrector - // longer term need to remove this - const Rich::Rec::RadPositionCorrector corrector; - - // global photon index - int photonIndex(-1); - - // Loop over the track->segment relations - for ( const auto & inTkRel : tkToSegRels ) - { - // loop over segments for this track - for ( const auto & segIndex : inTkRel.segmentIndices ) - { - // Get the data from the segment containers - const auto & segment = segments[segIndex]; - const auto & tkCkAngles = ckAngles[segIndex]; - const auto & tkCkRes = ckResolutions[segIndex]; - const auto & tkLocPtn = trHitPntsLoc[segIndex]; - const auto & segFlags = segPhotFlags[segIndex]; - - // Is this segment above threshold. Should be implicit due to above segment loop - //if ( ! ( tkCkAngles[lightestActiveHypo()] > 0 ) ) continue; - - // which RICH and radiator - const auto rich = segment.rich(); - const auto rad = segment.radiator(); - - // This implementiation does not support Aerogel - if ( UNLIKELY( Rich::Aerogel == rad ) ) { Exception("Aerogel not supported"); } - - // get the best pixel range for this segment, based on where hits are expected - const auto& pixR = ( segFlags.inBothPanels() ? RichRanges[rich] : - Rich::Rich1 == rich ? - ( segFlags.inPanel(Rich::top) ? - Rich1Ranges[Rich::top] : Rich1Ranges[Rich::bottom] ) : - ( segFlags.inPanel(Rich::left) ? - Rich2Ranges[Rich::left] : Rich2Ranges[Rich::right] ) ); - - // Loop over pixel information - int pixIndex = pixR.first - allPixels.begin() - 1; // (start index - 1) for this range - for ( auto pix_it = pixR.first; pix_it != pixR.second; ++pix_it ) - { - // load the data from the tuple - const auto & cluster = std::get<0>(*pix_it); - const auto & gloPos = std::get<1>(*pix_it); - const auto & locPos = std::get<2>(*pix_it); - // increment the pixel index - ++pixIndex; - - //_ri_verbo << "Trying Segment P=" << segment.bestMomentum() << " Pixel " - // << cluster << endmsg; - - // Apply pre-selection - - // Pixel position, in local HPD coords corrected for average radiator distortion - // Might need to pre-compute and cache this as we repeat it here each segment... - const auto pixP = corrector.correct(locPos,rad); - - // Track local hit point on the same panel as the hit - /// @todo CRJ Need to figure out why I was not correcting this in the current stack... - /// Probably correct but should check. - //const auto segPanelPnt = corrector.correct(tkLocPtn.point(cluster.panel()),rad); - const auto& segPanelPnt = tkLocPtn.point(cluster.panel()); - - // compute the seperation squared - const auto sep2 = ( std::pow( (pixP.x()-segPanelPnt.x()), 2 ) + - std::pow( (pixP.y()-segPanelPnt.y()), 2 ) ); - - // presel. default to rejected - bool SelOK = false; - - // Check overall boundaries - //_ri_verbo << " -> sep2 = " << sep2 << " " << m_maxROI2PreSel[rad] - // << " " << m_minROI2PreSel[rad] << endmsg; - if ( sep2 < m_maxROI2PreSel[rad] && sep2 > m_minROI2PreSel[rad] ) - { - - // estimated CK theta - const auto ckThetaEsti = std::sqrt(sep2)*m_scalePreSel[rad]; - //_ri_verbo << " -> CK theta Esti " << ckThetaEsti << endmsg; - - // Is the hit close to any mass hypo in local coordinate space - SelOK = any_of( activeParticles(), - [sigma=m_nSigmaPreSel[rad],&ckThetaEsti,&tkCkAngles,&tkCkRes] - ( const auto hypo ) - { return fabs(tkCkAngles[hypo]-ckThetaEsti) < sigma*tkCkRes[hypo]; } ); - - } // boundary check - - // Did we pass the pre-sel - if ( !SelOK ) - { - //_ri_verbo << " -> FAILED pre-selection -> reject" << endmsg; - continue; - } - - // Move on to the full Quartic reconstruction - - //std::cout << "Starting photon reco" << std::endl; - - // make a photon object to work on - photons.emplace_back(); - auto & gPhoton = photons.back(); - DataGuard<decltype(photons)> photGuard(photons); - - // pixel side - const auto side = cluster.panel(); - - // Emission point to use for photon reconstruction - // operate directly on photon data member for efficiency - //auto & emissionPoint = gPhoton.emissionPoint(); - //emissionPoint = segment.bestPoint(); - auto emissionPoint = segment.bestPoint(); - - // Photon direction at emission point - // Again, operator directly on data member - //auto & photonDirection = gPhoton.emissionDir(); - Gaudi::XYZVector photonDirection; - - // Final reflection points on sec and spherical mirrors - // operate directly on photon data - //auto & sphReflPoint = gPhoton.sphMirReflectionPoint(); - //auto & secReflPoint = gPhoton.flatMirReflectionPoint(); - Gaudi::XYZPoint sphReflPoint, secReflPoint; - - // fraction of segment path length accessible to the photon - float fraction(1); - - // Pointers to best sec and spherical mirror segments - const DeRichSphMirror * sphSegment = nullptr; - const DeRichSphMirror * secSegment = nullptr; - - // flag to say if this photon candidate is un-ambiguous - default to false - bool unambigPhoton = false; - - // find the reflection of the detection point in the sec mirror - // (virtual detection point) starting with nominal values - // At this we are assuming a flat nominal mirror common to all segments - auto distance = m_rich[rich]->nominalPlane(side).Distance(gloPos); - auto virtDetPoint = gloPos - 2.0 * distance * m_rich[rich]->nominalNormal(side); - - // -------------------------------------------------------------------------------------- - // For gas radiators, try start and end points to see if photon is unambiguous - // NOTE : For this we are using the virtual detection point determined using - // the noimnal flat secondary mirror plane. Now the secondary mirrors are actually - // spherical this may be introducing some additional uncertainties. - // -------------------------------------------------------------------------------------- - if ( m_testForUnambigPhots[rad] ) - { - if ( UNLIKELY( rad == Rich::Aerogel ) ) - { - // use default emission point and assume unambiguous since path length is so short.. - unambigPhoton = true; - } - else // gas radiators - { - - // flag for beam pipe check. Default is OK - bool beamTestOK(true); - - // ------------------------------------------------------------------------------- - // First emission point, at start of track segment - const auto emissionPoint1 = segment.bestPoint(0.01); - // Find mirror segments for this emission point - const DeRichSphMirror* sphSegment1 = nullptr; - const DeRichSphMirror* secSegment1 = nullptr; - Gaudi::XYZPoint sphReflPoint1, secReflPoint1; - if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint1, - sphSegment1, secSegment1, sphReflPoint1, secReflPoint1 ) ) - { - //_ri_debug << rad << " : Failed to reconstruct photon for start of segment" << endmsg; - continue; - } - if ( m_checkBeamPipe[rad] ) - { - beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint1, sphReflPoint1 ); - } - // ------------------------------------------------------------------------------- - - // ------------------------------------------------------------------------------- - // now do it again for emission point #2, at end of segment - const auto emissionPoint2 = segment.bestPoint(0.99); - // Find mirror segments for this emission point - const DeRichSphMirror* sphSegment2 = nullptr; - const DeRichSphMirror* secSegment2 = nullptr; - Gaudi::XYZPoint sphReflPoint2, secReflPoint2; - if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint2, - sphSegment2, secSegment2, sphReflPoint2, secReflPoint2 ) ) - { - //_ri_debug << rad << " : Failed to reconstruct photon for end of segment" << endmsg; - continue; - } - if ( !beamTestOK && m_checkBeamPipe[rad] ) - { - beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint2, sphReflPoint2 ); - } - // ------------------------------------------------------------------------------- - - // ------------------------------------------------------------------------------- - if ( !beamTestOK ) - { - // both start and end points failed beam pipe test -> reject - //_ri_debug << rad << " : Failed ambiguous photon beampipe intersection checks" << endmsg; - continue; - } - // ------------------------------------------------------------------------------- - - // ------------------------------------------------------------------------------- - // Test to see if photon direction hits the real physical mirror segments - // ------------------------------------------------------------------------------- - if ( UNLIKELY(m_checkPrimMirrSegs[rad]) ) - { - // primary mirrors - const auto ok = ( sphSegment1->intersects( emissionPoint1, - sphReflPoint1-emissionPoint1 ) || - sphSegment2->intersects( emissionPoint2, - sphReflPoint2-emissionPoint2 ) ); - if ( !ok ) - { - //_ri_debug << rad << " : Failed mirror segment intersection checks" << endmsg; - continue; - } - } - // ------------------------------------------------------------------------------- - - // ------------------------------------------------------------------------------- - // Get the best gas emission point - if ( !getBestGasEmissionPoint( rad, - sphReflPoint1, - sphReflPoint2, - gloPos, - segment, - emissionPoint, - fraction ) ) - { - //_ri_debug << rad << " : Failed to compute best gas emission point" << endmsg; - continue; - } - // ------------------------------------------------------------------------------- - - // ------------------------------------------------------------------------------- - // Is this an unambiguous photon - I.e. only one possible mirror combination - if ( ( sphSegment1 == sphSegment2 ) && ( secSegment1 == secSegment2 ) ) - { - // Set pointers to the mirror detector objects - sphSegment = sphSegment1; - secSegment = secSegment1; - // rough guesses at reflection points (improved later on) - sphReflPoint = sphReflPoint1 + 0.5*(sphReflPoint2-sphReflPoint1); - secReflPoint = secReflPoint1 + 0.5*(secReflPoint2-secReflPoint1); - // photon is not unambiguous - unambigPhoton = true; - } - // ------------------------------------------------------------------------------- - - } // end radiator type if - - // if configured to do so reject ambiguous photons - if ( UNLIKELY( m_rejectAmbigPhots[rad] && !unambigPhoton ) ) - { - //_ri_debug << rad << " : Failed ambiguous photon test" << endmsg; - continue; - } - - } // end unambiguous photon check - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // Active segment fraction cut - // -------------------------------------------------------------------------------------- - if ( UNLIKELY( fraction < m_minActiveFrac[rad] ) ) - { - //_ri_debug << rad << " : Failed active segment fraction cut" << endmsg; - continue; - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // if ambiguous gas photon or if ambiguous photon check above has been skipped, try again - // using best emission point and nominal mirror geometries to get the spherical and sec - // mirrors. Also, force this reconstruction if the above unambiguous test was skipped - // -------------------------------------------------------------------------------------- - if ( UNLIKELY( !m_testForUnambigPhots[rad] || !unambigPhoton ) ) - { - if ( !findMirrorData(rich,side, - virtDetPoint,emissionPoint, - sphSegment,secSegment,sphReflPoint,secReflPoint ) ) - { - //_ri_debug << rad << " : Failed backup photon reconstruction" << endmsg; - continue; - } - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // Finally reconstruct the photon using best emission point and the best mirror segments - // -------------------------------------------------------------------------------------- - if ( m_useAlignedMirrSegs[rad] ) - { - - // If iterations are disabled for this radiator, update the virtual detection point - // using the selected mirror segment - if ( UNLIKELY( 0 == m_nMaxQits[rad] ) ) - { - - // Form the virtual detection point - distance = secSegment->centreNormalPlane().Distance(gloPos); - virtDetPoint = gloPos - 2.0 * distance * secSegment->centreNormal(); - - // solve the quartic using the new data - m_quarticSolver.solve<double>( emissionPoint, - sphSegment->centreOfCurvature(), - virtDetPoint, - sphSegment->radius(), - sphReflPoint ); - - // (re)find the spherical mirror segment - // CRJ - Is this needed ? - sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); - - } - else - { - - // Iterate to final solution, improving the secondary mirror info - int iIt(0); - Gaudi::XYZPoint last_mirror_point(0,0,0); - bool loopOK = true; - while ( iIt < m_nMaxQits[rad] ) - { - - // Get secondary mirror reflection point, - // using the best actual secondary mirror segment at this point - const auto dir = virtDetPoint - sphReflPoint; - const auto sc = ( m_treatSecMirrsFlat[rich] ? - intersectPlane( sphReflPoint, dir, - secSegment->centreNormalPlane(), - secReflPoint ) : - intersectSpherical( sphReflPoint, dir, - secSegment->centreOfCurvature(), - secSegment->radius(), - secReflPoint ) ); - if ( UNLIKELY(!sc) ) - { - //_ri_debug << rad << " : Failed to intersect nominal secondary mirror plane" << endmsg; - loopOK = false; - break; - } - - // (re)find the secondary mirror - secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint ); - - // Construct plane tangential to secondary mirror passing through reflection point - const Gaudi::Plane3D plane( secSegment->centreOfCurvature()-secReflPoint, secReflPoint ); - - // re-find the reflection of the detection point in the sec mirror - // (virtual detection point) with this mirror plane - distance = plane.Distance(gloPos); - virtDetPoint = gloPos - 2.0 * distance * plane.Normal(); - - // solve the quartic using the new data - m_quarticSolver.solve<double>( emissionPoint, - sphSegment->centreOfCurvature(), - virtDetPoint, - sphSegment->radius(), - sphReflPoint ); - - // (re)find the spherical mirror segment - sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); - - // for iterations after the first, see if we are still moving - // If not, abort iterations early - if ( iIt > m_nMinQits[rad] && - (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolIt[rad] ) { break; } - - // store last sec mirror point - last_mirror_point = secReflPoint; - - // increment iteration counter - ++iIt; - } - // if above while loop failed, abort this photon - if ( UNLIKELY(!loopOK) ) { continue; } - - } - - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // check that spherical mirror reflection point is on the same side as detection point - // and (if configured to do so) photon does not cross between detector sides - // -------------------------------------------------------------------------------------- - if ( UNLIKELY( !sameSide(rad,sphReflPoint,virtDetPoint) ) ) - { - //_ri_debug << rad << " : Reflection point on wrong side" << endmsg; - continue; - } - if ( UNLIKELY( m_checkPhotCrossSides[rad] && - !sameSide(rad,sphReflPoint,emissionPoint) ) ) - { - //_ri_debug << rad << " : Photon crosses between detector sides" << endmsg; - continue; - } - //else { continue; } // uncomment to select ONLY crossing photons - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // For as radiators if ambiguous photon checks are disabled (since this is - // already done for these photons during those checks), check if the photon would have - // intersected with the beampipe - // -------------------------------------------------------------------------------------- - if ( UNLIKELY( m_checkBeamPipe[rad] && !m_testForUnambigPhots[rad] ) ) - { - if ( UNLIKELY( deBeam(rich)->testForIntersection( emissionPoint, sphReflPoint ) ) ) - { - //_ri_debug << rad << " : Failed final beampipe intersection checks" << endmsg; - continue; - } - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // If using aligned mirror segments, get the final sec mirror reflection - // point using the best mirror segments available at this point - // For RICH2, use the spherical nature of the scondary mirrors - // For RICH1, where they are much flatter, assume complete flatness - // -------------------------------------------------------------------------------------- - if ( m_useAlignedMirrSegs[rad] ) - { - const auto dir = virtDetPoint - sphReflPoint; - const auto sc = ( m_treatSecMirrsFlat[rich] ? - intersectPlane( sphReflPoint, dir, - secSegment->centreNormalPlane(), - secReflPoint ) : - intersectSpherical( sphReflPoint, dir, - secSegment->centreOfCurvature(), - secSegment->radius(), - secReflPoint ) ); - if ( !sc ) - { - //_ri_debug << rad << " : Failed final secondary mirror plane intersection" << endmsg; - continue; - } - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // calculate the cherenkov angles using the photon and track vectors - // -------------------------------------------------------------------------------------- - photonDirection = (sphReflPoint-emissionPoint).Unit(); - float thetaCerenkov(0), phiCerenkov(0); - segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov ); - // -------------------------------------------------------------------------------------- - - //--------------------------------------------------------------------------------------- - // Apply fudge factor correction for small biases in CK theta - //--------------------------------------------------------------------------------------- - thetaCerenkov += ckThetaCorrection(rad); - //--------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // Final checks on the Cherenkov angles - // -------------------------------------------------------------------------------------- - if ( UNLIKELY( !checkAngles( rad, tkCkAngles, tkCkRes, - thetaCerenkov, phiCerenkov ) ) ) - { - //_ri_verbo << " -> photon FAILED checkAngleInRange test" << endmsg; - continue; - } - // -------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------- - // Set (remaining) photon parameters - // -------------------------------------------------------------------------------------- - gPhoton.setCherenkovTheta ( thetaCerenkov ); - gPhoton.setCherenkovPhi ( phiCerenkov ); - gPhoton.setActiveSegmentFraction ( fraction ); - //gPhoton.setDetectionPoint ( gloPos ); - gPhoton.setSmartID ( cluster.primaryID() ); - gPhoton.setUnambiguousPhoton ( unambigPhoton ); - //gPhoton.setPrimaryMirror ( sphSegment ); - //gPhoton.setSecondaryMirror ( secSegment ); - //_ri_verbo << "Created photon " << gPhoton << endmsg; - // -------------------------------------------------------------------------------------- - - // If we get here, keep the just made photon - photGuard.setOK(); - // Save relations - relations.emplace_back( ++photonIndex, pixIndex, segIndex, inTkRel.tkIndex ); - //info() << relations.back() << endmsg; - - } - - } - - } - - //_ri_debug << "Created " << photons.size() << " photons" << endmsg; - return outData; -} - -//========================================================================= -// Find mirror segments and reflection points for given data -//========================================================================= -bool -QuarticPhotonReco:: -findMirrorData( const Rich::DetectorType rich, - const Rich::Side side, - const Gaudi::XYZPoint& virtDetPoint, - const Gaudi::XYZPoint& emissionPoint, - const DeRichSphMirror*& sphSegment, - const DeRichSphMirror*& secSegment, - Gaudi::XYZPoint& sphReflPoint, - Gaudi::XYZPoint& secReflPoint ) const -{ - // solve quartic equation with nominal values and find spherical mirror reflection point - m_quarticSolver.solve<double>( emissionPoint, - m_rich[rich]->nominalCentreOfCurvature(side), - virtDetPoint, - m_rich[rich]->sphMirrorRadius(), - sphReflPoint ); - // find the spherical mirror segment - sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); - // Search for the secondary segment - // Direction vector from primary mirror point to virtual detection point - const auto dir ( virtDetPoint - sphReflPoint ); - // find the sec mirror intersction point and secondary mirror segment - const bool OK = intersectPlane( sphReflPoint, dir, - m_rich[rich]->nominalPlane(side), - secReflPoint ); - if ( OK ) - { - // find the secondary mirror - secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint ); - // Re-find the secondary mirror reflection point using the new mirror info - // OK = ( !m_treatSecMirrsFlat[rich] ? - // intersectSpherical( sphReflPoint, dir, - // secSegment->centreOfCurvature(), - // secSegment->radius(), - // secReflPoint ) : - // intersectPlane( sphReflPoint, dir, - // secSegment->centreNormalPlane(), - // secReflPoint ) ); - } - // return - return OK; -} - -//========================================================================= -// Compute the best emission point for the gas radiators using -// the given spherical mirror reflection points -//========================================================================= -bool -QuarticPhotonReco:: -getBestGasEmissionPoint( const Rich::RadiatorType radiator, - const Gaudi::XYZPoint& sphReflPoint1, - const Gaudi::XYZPoint& sphReflPoint2, - const Gaudi::XYZPoint& detectionPoint, - const LHCb::RichTrackSegment & segment, - Gaudi::XYZPoint & emissionPoint, - float & fraction ) const -{ - double alongTkFrac = 0.5; - - if ( radiator == Rich::Rich1Gas ) - { - // First reflection and hit point on same y side ? - const bool sameSide1 = ( sphReflPoint1.y() * detectionPoint.y() > 0 ); - const bool sameSide2 = ( sphReflPoint2.y() * detectionPoint.y() > 0 ); - if ( sameSide1 && sameSide2 ) - { - emissionPoint = segment.bestPoint(); - } - else if ( sameSide1 ) - { - fraction = (float)(std::fabs(sphReflPoint1.y()/(sphReflPoint1.y()-sphReflPoint2.y()))); - alongTkFrac = fraction/2.0; - emissionPoint = segment.bestPoint(alongTkFrac); - } - else if ( sameSide2 ) - { - fraction = (float)(std::fabs(sphReflPoint2.y()/(sphReflPoint1.y()-sphReflPoint2.y()))); - alongTkFrac = 1.0-fraction/2.0; - emissionPoint = segment.bestPoint(alongTkFrac); - } - else - { - //Warning( "Rich1Gas : Both RICH spherical mirror hits opposite side to detection point" ); - return false; - } - - } - else if ( radiator == Rich::Rich2Gas ) - { - // First sphReflPoint and hit point on same x side ? - const bool sameSide1 = ( sphReflPoint1.x() * detectionPoint.x() > 0 ); - const bool sameSide2 = ( sphReflPoint2.x() * detectionPoint.x() > 0 ); - if ( sameSide1 && sameSide2 ) - { - emissionPoint = segment.bestPoint(); - } - else if ( sameSide1 ) - { - fraction = (float)(std::fabs(sphReflPoint1.x()/(sphReflPoint1.x()-sphReflPoint2.x()))); - alongTkFrac = fraction/2.0; - emissionPoint = segment.bestPoint(alongTkFrac); - } - else if ( sameSide2 ) - { - fraction = (float)(std::fabs(sphReflPoint2.x()/(sphReflPoint1.x()-sphReflPoint2.x()))); - alongTkFrac = 1.0-fraction/2.0; - emissionPoint = segment.bestPoint(alongTkFrac); - } - else - { - //Warning( "Rich2Gas : Both RICH spherical mirror hits opposite side to detection point" ); - return false; - } - } - else { Error( "::getBestGasEmissionPoint() called for Aerogel segment !!" ); } - - // if ( msgLevel(MSG::VERBOSE) ) - // { - // verbose() << radiator << " best emission point correction :- " << endmsg - // << " -> Photon detection point = " << detectionPoint << endmsg - // << " -> Sph. Mirror ptns = " << sphReflPoint1 << " " << sphReflPoint2 << endmsg - // << " -> Segment entry/exit = " << trSeg.entryPoint() << " " << trSeg.exitPoint() << endmsg - // << " -> Segment fraction = " << fraction << endmsg - // << " -> Emm. Ptn. Along traj = " << alongTkFrac << endmsg - // << " -> Best Emission point = " << emissionPoint << endmsg; - // } - - return true; -} - -//============================================================================= - // Declaration of the Algorithm Factory DECLARE_ALGORITHM_FACTORY( QuarticPhotonReco ) diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h index 7f5d5fe2a010003617ada530115cbc3fbc371a40..a5a3b1e54d3de1285bc3c6fa9624e3dc8f4559f5 100644 --- a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.h @@ -36,11 +36,17 @@ #include "RichDet/DeRichSphMirror.h" // Quartic Solver -#include "RichRecUtils/QuarticSolver.h" +#include "RichRecUtils/QuarticSolverNewton.h" // interfaces #include "RichInterfaces/IRichMirrorSegFinderLookUpTable.h" +// Vector Class +#include "VectorClass/instrset.h" + +// VDT +#include "vdt/asin.h" + namespace Rich { namespace Future @@ -51,6 +57,9 @@ namespace Rich // Use the functional framework using namespace Gaudi::Functional; + // pull in methods from Rich::RayTracingUtils + using namespace Rich::RayTracingUtils; + namespace { /// Shortcut to the output data type @@ -98,8 +107,73 @@ namespace Rich const Rich::PDPixelCluster::Vector& clusters, const SpacePointVector& globalPoints, const SpacePointVector& localPoints, - const Relations::TrackToSegments::Vector& tkToSegRels ) const override; + const Relations::TrackToSegments::Vector& tkToSegRels ) const override + { + return (this->*m_run)( segments, ckAngles, ckResolutions, trHitPntsLoc, + segPhotFlags, clusters, globalPoints, localPoints, tkToSegRels ); + } + private: // dispatch methods + + /** Pointer to dispatch function for the operator call + * Default to lowest common denominator. + * Reset in initialize() according to runtime support */ + OutData (QuarticPhotonReco::*m_run) + ( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const + = &QuarticPhotonReco::run_generic; + + /// Generic method + OutData run_generic( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const; + + /// SSE4 method + OutData run_sse4( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const; + + /// AVX method + OutData run_avx( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const; + + /// AVX2 method + OutData run_avx2( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const; + private: // methods /// Access the DeRich beam pipe objects @@ -115,15 +189,39 @@ namespace Rich inline float myacos( const float x ) const noexcept { return vdt::fast_acosf(x); } /// Find mirror segments and reflection points for given data - bool findMirrorData( const Rich::DetectorType rich, - const Rich::Side side, - const Gaudi::XYZPoint & virtDetPoint, - const Gaudi::XYZPoint & emissionPoint, - const DeRichSphMirror *& sphSegment, - const DeRichSphMirror *& secSegment, - Gaudi::XYZPoint & sphReflPoint, - Gaudi::XYZPoint & secReflPoint ) const; - + inline bool findMirrorData( const Rich::DetectorType rich, + const Rich::Side side, + const Gaudi::XYZPoint & virtDetPoint, + const Gaudi::XYZPoint & emissionPoint, + const DeRichSphMirror *& sphSegment, + const DeRichSphMirror *& secSegment, + Gaudi::XYZPoint & sphReflPoint, + Gaudi::XYZPoint & secReflPoint ) const + { + // solve quartic equation with nominal values and find spherical mirror reflection point + m_quarticSolver.solve<double,3,3>( emissionPoint, + m_rich[rich]->nominalCentreOfCurvature(side), + virtDetPoint, + m_rich[rich]->sphMirrorRadius(), + sphReflPoint ); + // find the spherical mirror segment + sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); + // Search for the secondary segment + // Direction vector from primary mirror point to virtual detection point + const auto dir ( virtDetPoint - sphReflPoint ); + // find the sec mirror intersction point and secondary mirror segment + const bool OK = intersectPlane( sphReflPoint, dir, + m_rich[rich]->nominalPlane(side), + secReflPoint ); + if ( OK ) + { + // find the secondary mirror + secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint ); + } + // return + return OK; + } + /// Compute the best emission point for the gas radiators using /// the given spherical mirror reflection points bool getBestGasEmissionPoint( const Rich::RadiatorType radiator, @@ -136,8 +234,8 @@ namespace Rich private: - /// Quartic Solver - Rich::Rec::QuarticSolver m_quarticSolver; + /// Newton Quartic Solver + Rich::Rec::QuarticSolverNewton m_quarticSolver; /// Rich1 and Rich2 detector elements DetectorArray<const DeRich *> m_rich = {{}}; @@ -224,8 +322,98 @@ namespace Rich ToolHandle<const IMirrorSegFinderLookUpTable> m_mirrorSegFinder { "Rich::Future::MirrorSegFinderLookUpTable/MirrorFinder:PUBLIC", this }; + /// Flag to control if CPU detection should be performed or not + Gaudi::Property<bool> m_detectCPU { this, "DetectCPUCapabilites", true }; + }; + //========================================================================= + // Compute the best emission point for the gas radiators using + // the given spherical mirror reflection points + //========================================================================= + inline bool + QuarticPhotonReco:: + getBestGasEmissionPoint( const Rich::RadiatorType radiator, + const Gaudi::XYZPoint& sphReflPoint1, + const Gaudi::XYZPoint& sphReflPoint2, + const Gaudi::XYZPoint& detectionPoint, + const LHCb::RichTrackSegment & segment, + Gaudi::XYZPoint & emissionPoint, + float & fraction ) const + { + double alongTkFrac = 0.5; + + if ( radiator == Rich::Rich1Gas ) + { + // First reflection and hit point on same y side ? + const bool sameSide1 = ( sphReflPoint1.y() * detectionPoint.y() > 0 ); + const bool sameSide2 = ( sphReflPoint2.y() * detectionPoint.y() > 0 ); + if ( sameSide1 && sameSide2 ) + { + emissionPoint = segment.bestPoint(); + } + else if ( sameSide1 ) + { + fraction = (float)(std::fabs(sphReflPoint1.y()/(sphReflPoint1.y()-sphReflPoint2.y()))); + alongTkFrac = fraction/2.0; + emissionPoint = segment.bestPoint(alongTkFrac); + } + else if ( sameSide2 ) + { + fraction = (float)(std::fabs(sphReflPoint2.y()/(sphReflPoint1.y()-sphReflPoint2.y()))); + alongTkFrac = 1.0-fraction/2.0; + emissionPoint = segment.bestPoint(alongTkFrac); + } + else + { + //Warning( "Rich1Gas : Both RICH spherical mirror hits opposite side to detection point" ); + return false; + } + + } + else if ( radiator == Rich::Rich2Gas ) + { + // First sphReflPoint and hit point on same x side ? + const bool sameSide1 = ( sphReflPoint1.x() * detectionPoint.x() > 0 ); + const bool sameSide2 = ( sphReflPoint2.x() * detectionPoint.x() > 0 ); + if ( sameSide1 && sameSide2 ) + { + emissionPoint = segment.bestPoint(); + } + else if ( sameSide1 ) + { + fraction = (float)(std::fabs(sphReflPoint1.x()/(sphReflPoint1.x()-sphReflPoint2.x()))); + alongTkFrac = fraction/2.0; + emissionPoint = segment.bestPoint(alongTkFrac); + } + else if ( sameSide2 ) + { + fraction = (float)(std::fabs(sphReflPoint2.x()/(sphReflPoint1.x()-sphReflPoint2.x()))); + alongTkFrac = 1.0-fraction/2.0; + emissionPoint = segment.bestPoint(alongTkFrac); + } + else + { + //Warning( "Rich2Gas : Both RICH spherical mirror hits opposite side to detection point" ); + return false; + } + } + else { Error( "::getBestGasEmissionPoint() called for Aerogel segment !!" ); } + + // if ( msgLevel(MSG::VERBOSE) ) + // { + // verbose() << radiator << " best emission point correction :- " << endmsg + // << " -> Photon detection point = " << detectionPoint << endmsg + // << " -> Sph. Mirror ptns = " << sphReflPoint1 << " " << sphReflPoint2 << endmsg + // << " -> Segment entry/exit = " << trSeg.entryPoint() << " " << trSeg.exitPoint() << endmsg + // << " -> Segment fraction = " << fraction << endmsg + // << " -> Emm. Ptn. Along traj = " << alongTkFrac << endmsg + // << " -> Best Emission point = " << emissionPoint << endmsg; + // } + + return true; + } + } } } diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp new file mode 100644 index 0000000000000000000000000000000000000000..63eb84b8e528f4df75650fbd2bd017b21506af6e --- /dev/null +++ b/Rich/RichFutureRecPhotonAlgorithms/src/RichQuarticPhotonReco.icpp @@ -0,0 +1,534 @@ + +{ + // use ranges v3 library + using namespace ranges::v3; + + // make the output data + OutData outData; + + // Shortcut to the photons and relations + auto & photons = std::get<0>(outData); + auto & relations = std::get<1>(outData); + // guess at reserve size + const auto resSize = segments.size() * globalPoints.size() / 20; + photons.reserve( resSize ); + relations.reserve( resSize ); + + // View range for pixel data + const auto allPixels = Ranges::ConstZip( clusters, globalPoints, localPoints ); + + // Make 'ranges' for RICH pixels + const auto RPixs = richRanges( allPixels ); + const auto & RichRanges = RPixs[0]; // Complete RICH ranges + const auto & Rich1Ranges = RPixs[1]; // ranges for RICH1 [top,bottom] + const auto & Rich2Ranges = RPixs[2]; // ranges for RICH2 [left,right] + + // local position corrector + // longer term need to remove this + const Rich::Rec::RadPositionCorrector corrector; + + // global photon index + int photonIndex(-1); + + // Loop over the track->segment relations + for ( const auto & inTkRel : tkToSegRels ) + { + // loop over segments for this track + for ( const auto & segIndex : inTkRel.segmentIndices ) + { + // Get the data from the segment containers + const auto & segment = segments[segIndex]; + const auto & tkCkAngles = ckAngles[segIndex]; + const auto & tkCkRes = ckResolutions[segIndex]; + const auto & tkLocPtn = trHitPntsLoc[segIndex]; + const auto & segFlags = segPhotFlags[segIndex]; + + // Is this segment above threshold. Should be implicit due to above segment loop + //if ( ! ( tkCkAngles[lightestActiveHypo()] > 0 ) ) continue; + + // which RICH and radiator + const auto rich = segment.rich(); + const auto rad = segment.radiator(); + + // This implementiation does not support Aerogel + if ( UNLIKELY( Rich::Aerogel == rad ) ) { Exception("Aerogel not supported"); } + + // get the best pixel range for this segment, based on where hits are expected + const auto& pixR = ( segFlags.inBothPanels() ? RichRanges[rich] : + Rich::Rich1 == rich ? + ( segFlags.inPanel(Rich::top) ? + Rich1Ranges[Rich::top] : Rich1Ranges[Rich::bottom] ) : + ( segFlags.inPanel(Rich::left) ? + Rich2Ranges[Rich::left] : Rich2Ranges[Rich::right] ) ); + + // Loop over pixel information + int pixIndex = pixR.first - allPixels.begin() - 1; // (start index - 1) for this range + for ( auto pix_it = pixR.first; pix_it != pixR.second; ++pix_it ) + { + // load the data from the tuple + const auto & cluster = std::get<0>(*pix_it); + const auto & gloPos = std::get<1>(*pix_it); + const auto & locPos = std::get<2>(*pix_it); + // increment the pixel index + ++pixIndex; + + //_ri_verbo << "Trying Segment P=" << segment.bestMomentum() << " Pixel " + // << cluster << endmsg; + + // Apply pre-selection + + // Pixel position, in local HPD coords corrected for average radiator distortion + // Might need to pre-compute and cache this as we repeat it here each segment... + const auto pixP = corrector.correct(locPos,rad); + + // Track local hit point on the same panel as the hit + /// @todo CRJ Need to figure out why I was not correcting this in the current stack... + /// Probably correct but should check. + //const auto segPanelPnt = corrector.correct(tkLocPtn.point(cluster.panel()),rad); + const auto& segPanelPnt = tkLocPtn.point(cluster.panel()); + + // compute the seperation squared + const auto sep2 = ( std::pow( (pixP.x()-segPanelPnt.x()), 2 ) + + std::pow( (pixP.y()-segPanelPnt.y()), 2 ) ); + + // presel. default to rejected + bool SelOK = false; + + // Check overall boundaries + //_ri_verbo << " -> sep2 = " << sep2 << " " << m_maxROI2PreSel[rad] + // << " " << m_minROI2PreSel[rad] << endmsg; + if ( sep2 < m_maxROI2PreSel[rad] && sep2 > m_minROI2PreSel[rad] ) + { + + // estimated CK theta + const auto ckThetaEsti = std::sqrt(sep2)*m_scalePreSel[rad]; + //_ri_verbo << " -> CK theta Esti " << ckThetaEsti << endmsg; + + // Is the hit close to any mass hypo in local coordinate space + SelOK = any_of( activeParticles(), + [sigma=m_nSigmaPreSel[rad],&ckThetaEsti,&tkCkAngles,&tkCkRes] + ( const auto hypo ) + { return fabs(tkCkAngles[hypo]-ckThetaEsti) < sigma*tkCkRes[hypo]; } ); + + } // boundary check + + // Did we pass the pre-sel + if ( !SelOK ) + { + //_ri_verbo << " -> FAILED pre-selection -> reject" << endmsg; + continue; + } + + // Move on to the full Quartic reconstruction + + //std::cout << "Starting photon reco" << std::endl; + + // make a photon object to work on + photons.emplace_back(); + auto & gPhoton = photons.back(); + DataGuard<decltype(photons)> photGuard(photons); + + // pixel side + const auto side = cluster.panel(); + + // Emission point to use for photon reconstruction + // operate directly on photon data member for efficiency + //auto & emissionPoint = gPhoton.emissionPoint(); + //emissionPoint = segment.bestPoint(); + auto emissionPoint = segment.bestPoint(); + + // Photon direction at emission point + // Again, operator directly on data member + //auto & photonDirection = gPhoton.emissionDir(); + Gaudi::XYZVector photonDirection; + + // Final reflection points on sec and spherical mirrors + // operate directly on photon data + //auto & sphReflPoint = gPhoton.sphMirReflectionPoint(); + //auto & secReflPoint = gPhoton.flatMirReflectionPoint(); + Gaudi::XYZPoint sphReflPoint, secReflPoint; + + // fraction of segment path length accessible to the photon + float fraction(1); + + // Pointers to best sec and spherical mirror segments + const DeRichSphMirror * sphSegment = nullptr; + const DeRichSphMirror * secSegment = nullptr; + + // flag to say if this photon candidate is un-ambiguous - default to false + bool unambigPhoton = false; + + // find the reflection of the detection point in the sec mirror + // (virtual detection point) starting with nominal values + // At this we are assuming a flat nominal mirror common to all segments + const auto nom_dist = m_rich[rich]->nominalPlane(side).Distance(gloPos); + auto virtDetPoint = gloPos - 2.0 * nom_dist * m_rich[rich]->nominalNormal(side); + + // -------------------------------------------------------------------------------------- + // For gas radiators, try start and end points to see if photon is unambiguous + // NOTE : For this we are using the virtual detection point determined using + // the noimnal flat secondary mirror plane. Now the secondary mirrors are actually + // spherical this may be introducing some additional uncertainties. + // -------------------------------------------------------------------------------------- + if ( m_testForUnambigPhots[rad] ) + { + if ( UNLIKELY( rad == Rich::Aerogel ) ) + { + // use default emission point and assume unambiguous since path length is so short.. + unambigPhoton = true; + } + else // gas radiators + { + + // flag for beam pipe check. Default is OK + bool beamTestOK(true); + + // ------------------------------------------------------------------------------- + // First emission point, at start of track segment + const auto emissionPoint1 = segment.bestPoint(0.01); + // Find mirror segments for this emission point + const DeRichSphMirror* sphSegment1 = nullptr; + const DeRichSphMirror* secSegment1 = nullptr; + Gaudi::XYZPoint sphReflPoint1, secReflPoint1; + if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint1, + sphSegment1, secSegment1, sphReflPoint1, secReflPoint1 ) ) + { + //_ri_debug << rad << " : Failed to reconstruct photon for start of segment" << endmsg; + continue; + } + if ( m_checkBeamPipe[rad] ) + { + beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint1, sphReflPoint1 ); + } + // ------------------------------------------------------------------------------- + + // ------------------------------------------------------------------------------- + // now do it again for emission point #2, at end of segment + const auto emissionPoint2 = segment.bestPoint(0.99); + // Find mirror segments for this emission point + const DeRichSphMirror* sphSegment2 = nullptr; + const DeRichSphMirror* secSegment2 = nullptr; + Gaudi::XYZPoint sphReflPoint2, secReflPoint2; + if ( !findMirrorData( rich, side, virtDetPoint, emissionPoint2, + sphSegment2, secSegment2, sphReflPoint2, secReflPoint2 ) ) + { + //_ri_debug << rad << " : Failed to reconstruct photon for end of segment" << endmsg; + continue; + } + if ( !beamTestOK && m_checkBeamPipe[rad] ) + { + beamTestOK = !deBeam(rich)->testForIntersection( emissionPoint2, sphReflPoint2 ); + } + // ------------------------------------------------------------------------------- + + // ------------------------------------------------------------------------------- + if ( !beamTestOK ) + { + // both start and end points failed beam pipe test -> reject + //_ri_debug << rad << " : Failed ambiguous photon beampipe intersection checks" << endmsg; + continue; + } + // ------------------------------------------------------------------------------- + + // ------------------------------------------------------------------------------- + // Test to see if photon direction hits the real physical mirror segments + // ------------------------------------------------------------------------------- + if ( UNLIKELY(m_checkPrimMirrSegs[rad]) ) + { + // primary mirrors + const auto ok = ( sphSegment1->intersects( emissionPoint1, + sphReflPoint1-emissionPoint1 ) || + sphSegment2->intersects( emissionPoint2, + sphReflPoint2-emissionPoint2 ) ); + if ( !ok ) + { + //_ri_debug << rad << " : Failed mirror segment intersection checks" << endmsg; + continue; + } + } + // ------------------------------------------------------------------------------- + + // ------------------------------------------------------------------------------- + // Get the best gas emission point + if ( !getBestGasEmissionPoint( rad, + sphReflPoint1, + sphReflPoint2, + gloPos, + segment, + emissionPoint, + fraction ) ) + { + //_ri_debug << rad << " : Failed to compute best gas emission point" << endmsg; + continue; + } + // ------------------------------------------------------------------------------- + + // ------------------------------------------------------------------------------- + // Is this an unambiguous photon - I.e. only one possible mirror combination + if ( ( sphSegment1 == sphSegment2 ) && ( secSegment1 == secSegment2 ) ) + { + // Set pointers to the mirror detector objects + sphSegment = sphSegment1; + secSegment = secSegment1; + // rough guesses at reflection points (improved later on) + sphReflPoint = sphReflPoint1 + 0.5*(sphReflPoint2-sphReflPoint1); + secReflPoint = secReflPoint1 + 0.5*(secReflPoint2-secReflPoint1); + // photon is not unambiguous + unambigPhoton = true; + } + // ------------------------------------------------------------------------------- + + } // end radiator type if + + // if configured to do so reject ambiguous photons + if ( UNLIKELY( m_rejectAmbigPhots[rad] && !unambigPhoton ) ) + { + //_ri_debug << rad << " : Failed ambiguous photon test" << endmsg; + continue; + } + + } // end unambiguous photon check + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // Active segment fraction cut + // -------------------------------------------------------------------------------------- + if ( UNLIKELY( fraction < m_minActiveFrac[rad] ) ) + { + //_ri_debug << rad << " : Failed active segment fraction cut" << endmsg; + continue; + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // if ambiguous gas photon or if ambiguous photon check above has been skipped, try again + // using best emission point and nominal mirror geometries to get the spherical and sec + // mirrors. Also, force this reconstruction if the above unambiguous test was skipped + // -------------------------------------------------------------------------------------- + if ( UNLIKELY( !m_testForUnambigPhots[rad] || !unambigPhoton ) ) + { + if ( !findMirrorData(rich,side, + virtDetPoint,emissionPoint, + sphSegment,secSegment,sphReflPoint,secReflPoint ) ) + { + //_ri_debug << rad << " : Failed backup photon reconstruction" << endmsg; + continue; + } + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // Finally reconstruct the photon using best emission point and the best mirror segments + // -------------------------------------------------------------------------------------- + if ( m_useAlignedMirrSegs[rad] ) + { + + // If iterations are disabled for this radiator, update the virtual detection point + // using the selected mirror segment + if ( UNLIKELY( 0 == m_nMaxQits[rad] ) ) + { + + // Form the virtual detection point + const auto distance = secSegment->centreNormalPlane().Distance(gloPos); + virtDetPoint = gloPos - 2.0 * distance * secSegment->centreNormal(); + + // solve the quartic using the new data + m_quarticSolver.solve<double,3,4>( emissionPoint, + sphSegment->centreOfCurvature(), + virtDetPoint, + sphSegment->radius(), + sphReflPoint ); + + // (re)find the spherical mirror segment + // CRJ - Is this needed ? + sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); + + } + else + { + + // Iterate to final solution, improving the secondary mirror info + int iIt(0); + Gaudi::XYZPoint last_mirror_point(0,0,0); + bool loopOK = true; + while ( iIt < m_nMaxQits[rad] ) + { + + // Get secondary mirror reflection point, + // using the best actual secondary mirror segment at this point + const auto dir = virtDetPoint - sphReflPoint; + const auto sc = ( m_treatSecMirrsFlat[rich] ? + intersectPlane( sphReflPoint, dir, + secSegment->centreNormalPlane(), + secReflPoint ) : + intersectSpherical( sphReflPoint, dir, + secSegment->centreOfCurvature(), + secSegment->radius(), + secReflPoint ) ); + if ( UNLIKELY(!sc) ) + { + //_ri_debug << rad << " : Failed to intersect nominal secondary mirror plane" << endmsg; + loopOK = false; + break; + } + + // (re)find the secondary mirror + secSegment = m_mirrorSegFinder.get()->findSecMirror( rich, side, secReflPoint ); + + // Compute the virtual reflection point + + // Construct plane tangential to secondary mirror passing through reflection point + //const Gaudi::Plane3D plane( secSegment->centreOfCurvature()-secReflPoint, secReflPoint ); + // re-find the reflection of the detection point in the sec mirror + // (virtual detection point) with this mirror plane + //const auto distance = plane.Distance(gloPos); + //virtDetPoint = gloPos - 2.0 * distance * plane.Normal(); + + // Same as above, just written out by hand and simplified a bit + const auto normV = secSegment->centreOfCurvature() - secReflPoint; + const auto distance = ( ( normV.X() * ( gloPos.X() - secReflPoint.X() ) ) + + ( normV.Y() * ( gloPos.Y() - secReflPoint.Y() ) ) + + ( normV.Z() * ( gloPos.Z() - secReflPoint.Z() ) ) ); + virtDetPoint = gloPos - ( normV * ( 2.0 * distance / normV.Mag2() ) ); + + // solve the quartic using the new data + m_quarticSolver.solve<double,3,4>( emissionPoint, + sphSegment->centreOfCurvature(), + virtDetPoint, + sphSegment->radius(), + sphReflPoint ); + + // (re)find the spherical mirror segment + sphSegment = m_mirrorSegFinder.get()->findSphMirror( rich, side, sphReflPoint ); + + // for iterations after the first, see if we are still moving + // If not, abort iterations early + if ( iIt > m_nMinQits[rad] && + (last_mirror_point-secReflPoint).Mag2() < m_minSphMirrTolIt[rad] ) { break; } + + // store last sec mirror point + last_mirror_point = secReflPoint; + + // increment iteration counter + ++iIt; + } + // if above while loop failed, abort this photon + if ( UNLIKELY(!loopOK) ) { continue; } + + } + + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // check that spherical mirror reflection point is on the same side as detection point + // and (if configured to do so) photon does not cross between detector sides + // -------------------------------------------------------------------------------------- + if ( UNLIKELY( !sameSide(rad,sphReflPoint,virtDetPoint) ) ) + { + //_ri_debug << rad << " : Reflection point on wrong side" << endmsg; + continue; + } + if ( UNLIKELY( m_checkPhotCrossSides[rad] && + !sameSide(rad,sphReflPoint,emissionPoint) ) ) + { + //_ri_debug << rad << " : Photon crosses between detector sides" << endmsg; + continue; + } + //else { continue; } // uncomment to select ONLY crossing photons + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // For as radiators if ambiguous photon checks are disabled (since this is + // already done for these photons during those checks), check if the photon would have + // intersected with the beampipe + // -------------------------------------------------------------------------------------- + if ( UNLIKELY( m_checkBeamPipe[rad] && !m_testForUnambigPhots[rad] ) ) + { + if ( UNLIKELY( deBeam(rich)->testForIntersection( emissionPoint, sphReflPoint ) ) ) + { + //_ri_debug << rad << " : Failed final beampipe intersection checks" << endmsg; + continue; + } + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // If using aligned mirror segments, get the final sec mirror reflection + // point using the best mirror segments available at this point + // For RICH2, use the spherical nature of the scondary mirrors + // For RICH1, where they are much flatter, assume complete flatness + // -------------------------------------------------------------------------------------- + if ( m_useAlignedMirrSegs[rad] ) + { + const auto dir = virtDetPoint - sphReflPoint; + const auto sc = ( m_treatSecMirrsFlat[rich] ? + intersectPlane( sphReflPoint, dir, + secSegment->centreNormalPlane(), + secReflPoint ) : + intersectSpherical( sphReflPoint, dir, + secSegment->centreOfCurvature(), + secSegment->radius(), + secReflPoint ) ); + if ( !sc ) + { + //_ri_debug << rad << " : Failed final secondary mirror plane intersection" << endmsg; + continue; + } + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // calculate the cherenkov angles using the photon and track vectors + // -------------------------------------------------------------------------------------- + photonDirection = (sphReflPoint-emissionPoint).Unit(); + float thetaCerenkov(0), phiCerenkov(0); + segment.angleToDirection( photonDirection, thetaCerenkov, phiCerenkov ); + // -------------------------------------------------------------------------------------- + + //--------------------------------------------------------------------------------------- + // Apply fudge factor correction for small biases in CK theta + //--------------------------------------------------------------------------------------- + thetaCerenkov += ckThetaCorrection(rad); + //--------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // Final checks on the Cherenkov angles + // -------------------------------------------------------------------------------------- + if ( UNLIKELY( !checkAngles( rad, tkCkAngles, tkCkRes, + thetaCerenkov, phiCerenkov ) ) ) + { + //_ri_verbo << " -> photon FAILED checkAngleInRange test" << endmsg; + continue; + } + // -------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------- + // Set (remaining) photon parameters + // -------------------------------------------------------------------------------------- + gPhoton.setCherenkovTheta ( thetaCerenkov ); + gPhoton.setCherenkovPhi ( phiCerenkov ); + gPhoton.setActiveSegmentFraction ( fraction ); + //gPhoton.setDetectionPoint ( gloPos ); + gPhoton.setSmartID ( cluster.primaryID() ); + gPhoton.setUnambiguousPhoton ( unambigPhoton ); + //gPhoton.setPrimaryMirror ( sphSegment ); + //gPhoton.setSecondaryMirror ( secSegment ); + //_ri_verbo << "Created photon " << gPhoton << endmsg; + // -------------------------------------------------------------------------------------- + + // If we get here, keep the just made photon + photGuard.setOK(); + // Save relations + relations.emplace_back( ++photonIndex, pixIndex, segIndex, inTkRel.tkIndex ); + //info() << relations.back() << endmsg; + + } + + } + + } + + //_ri_debug << "Created " << photons.size() << " photons" << endmsg; + return outData; +} diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1840557ea2ddd9b2e1403f05c293205a77c93823 --- /dev/null +++ b/Rich/RichFutureRecPhotonAlgorithms/src/avx/RichQuarticPhotonReco.cpp @@ -0,0 +1,30 @@ + +// local +#include "../RichQuarticPhotonReco.h" + +// All code is in general Rich reconstruction namespace +using namespace Rich::Future; +using namespace Rich::Future::Rec; + +struct avx_guard { + // see Agner Fog's optimization guide, 12.1 about mixing AVX and non-AVX code, + // (http://www.agner.org/optimize/optimizing_cpp.pdf) + // and preserving the YMM register state. + // Invoking __mm256_zeroupper seems to reduce the overhead when switching. + ~avx_guard() { _mm256_zeroupper(); } +}; + +OutData +QuarticPhotonReco::run_avx( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const +{ + avx_guard guard{}; +#include "../RichQuarticPhotonReco.icpp" +} diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp new file mode 100644 index 0000000000000000000000000000000000000000..672aa639912aeaa8181868a4cf4af04f680a54d5 --- /dev/null +++ b/Rich/RichFutureRecPhotonAlgorithms/src/avx2/RichQuarticPhotonReco.cpp @@ -0,0 +1,30 @@ + +// local +#include "../RichQuarticPhotonReco.h" + +// All code is in general Rich reconstruction namespace +using namespace Rich::Future; +using namespace Rich::Future::Rec; + +struct avx_guard { + // see Agner Fog's optimization guide, 12.1 about mixing AVX and non-AVX code, + // (http://www.agner.org/optimize/optimizing_cpp.pdf) + // and preserving the YMM register state. + // Invoking __mm256_zeroupper seems to reduce the overhead when switching. + ~avx_guard() { _mm256_zeroupper(); } +}; + +OutData +QuarticPhotonReco::run_avx2( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const +{ + avx_guard guard{}; +#include "../RichQuarticPhotonReco.icpp" +} diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d5db58c9144febf65e7adeec47cd82df50ae8939 --- /dev/null +++ b/Rich/RichFutureRecPhotonAlgorithms/src/generic/RichQuarticPhotonReco.cpp @@ -0,0 +1,21 @@ + +// local +#include "../RichQuarticPhotonReco.h" + +// All code is in general Rich reconstruction namespace +using namespace Rich::Future; +using namespace Rich::Future::Rec; + +OutData +QuarticPhotonReco::run_generic( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const +{ +#include "../RichQuarticPhotonReco.icpp" +} diff --git a/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp b/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp new file mode 100644 index 0000000000000000000000000000000000000000..82131fc1330fee441861b78481e6720765de93bf --- /dev/null +++ b/Rich/RichFutureRecPhotonAlgorithms/src/sse4/RichQuarticPhotonReco.cpp @@ -0,0 +1,21 @@ + +// local +#include "../RichQuarticPhotonReco.h" + +// All code is in general Rich reconstruction namespace +using namespace Rich::Future; +using namespace Rich::Future::Rec; + +OutData +QuarticPhotonReco::run_sse4( const LHCb::RichTrackSegment::Vector& segments, + const CherenkovAngles::Vector& ckAngles, + const CherenkovResolutions::Vector& ckResolutions, + const SegmentPanelSpacePoints::Vector& trHitPntsLoc, + const SegmentPhotonFlags::Vector& segPhotFlags, + const Rich::PDPixelCluster::Vector& clusters, + const SpacePointVector& globalPoints, + const SpacePointVector& localPoints, + const Relations::TrackToSegments::Vector& tkToSegRels ) const +{ +#include "../RichQuarticPhotonReco.icpp" +} diff --git a/Rich/RichFutureRecSys/examples/Brunel.py b/Rich/RichFutureRecSys/examples/Brunel.py index a522b3d5058f1198d5b1860c74c74c1417ff9ae3..c6b69a4c3240d416ad4b05067c5500275a2511cb 100644 --- a/Rich/RichFutureRecSys/examples/Brunel.py +++ b/Rich/RichFutureRecSys/examples/Brunel.py @@ -45,15 +45,15 @@ Brunel().MCCheckSequence = richs+["PROTO"] ############################################################################# # Timestamps in messages -LHCbApp().TimeStamp = True +#LHCbApp().TimeStamp = True Brunel().OutputType = 'None' importOptions("$APPCONFIGOPTS/Persistency/Compression-ZLIB-1.py") -Brunel().EvtMax = 5000 +Brunel().EvtMax = 10000 Brunel().PrintFreq = 500 -#Brunel().Histograms = "Expert" +Brunel().Histograms = "Expert" #ApplicationMgr().ExtSvc += [ "AuditorSvc" ] diff --git a/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py b/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py index 93f5018ebebef6b0bf6c1a834b1a45a11fd0f93e..c56359d4d38a201539cafcc9bb41cc0392b84854 100644 --- a/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py +++ b/Rich/RichFutureRecSys/examples/Col16RawDataFiles.py @@ -21,3 +21,6 @@ from Configurables import Brunel, LHCbApp Brunel().DataType = "2016" Brunel().Simulation = False Brunel().WithMC = False + +#Brunel().DDDBtag = "dddb-20150724" +#Brunel().CondDBtag = "cond-20161011" diff --git a/Rich/RichFutureRecSys/examples/RichFuture.py b/Rich/RichFutureRecSys/examples/RichFuture.py index 0186612021d870948756bb5d8a4cb2263f1e10a5..43d45a745e14b777f2824863308e2766ab7ce2ae 100644 --- a/Rich/RichFutureRecSys/examples/RichFuture.py +++ b/Rich/RichFutureRecSys/examples/RichFuture.py @@ -94,7 +94,7 @@ RichCheck = RichRecCheckers( inputTrackLocations = tkLocs, # The final sequence to run all = GaudiSequencer( "All", Members = [ fetcher, odinDecode, richDecode, tkUnpack, tkFilt, RichRec ] ) -all.Members += [ RichMoni ] +#all.Members += [ RichMoni ] #all.Members += [ mcPUnpack, RichCheck ] all.MeasureTime = True diff --git a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py index d4124f78e8cec8fe1fc426ab08a4d493ee41944d..c21660c283ff5cfc0226b1423bff3ed8e7fcff21 100644 --- a/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py +++ b/Rich/RichFutureRecSys/python/RichFutureRecSys/ConfiguredRichReco.py @@ -137,7 +137,7 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group. preloadGeometry = False, # Number points for mass hypothesis ring ray tracing - NRingPoints = ( 100, 100, 100 ), + NRingPoints = ( 96, 96, 96 ), # Track Extrapolator type trackExtrapolator = "TrackRungeKuttaExtrapolator", @@ -243,9 +243,6 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group. if algprops["Detectors"][0] : pixname += "RICH1" if algprops["Detectors"][1] : pixname += "RICH2" - # The complete sequence - all = GaudiSequencer( "RichRecoSeq" + GroupName, MeasureTime = MeasureTime ) - # ================================================================== # RICH Pixel Treatment. Common to all track types and instanciations # of the RICH sequences. @@ -291,8 +288,9 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group. pixels = GaudiSequencer( "RichPixels"+GroupName, MeasureTime = MeasureTime, Members = [ richCluster, richGloPoints, richLocPoints ] ) - # Add to final sequence - all.Members += [ pixels ] + # The complete sequence + all = GaudiSequencer( "RichRecoSeq" + GroupName, MeasureTime = MeasureTime ) + all.Members = [pixels] # Loop over tracks for tktype,trackLocation in inputTrackLocations.iteritems() : @@ -308,6 +306,7 @@ def RichRecoSequence( GroupName = "", # Optional name given to this group. # Sequencer for this track type tkSeq = GaudiSequencer( "Rich"+name+"Reco", MeasureTime = MeasureTime ) + tkSeq.Members = [ ] # make sure empty all.Members += [ tkSeq ] # ================================================================== diff --git a/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt b/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt index d0eb77b7e8e02d7f9c925bea2cd98d97a0dab547..c4c1f41f2891f34a711319d20fec457c0315a3c2 100644 --- a/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt +++ b/Rich/RichFutureRecTrackAlgorithms/CMakeLists.txt @@ -15,10 +15,9 @@ find_package(Boost) find_package(GSL) find_package(VDT) -include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}) +include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${Vc_INCLUDE_DIR}) gaudi_add_module(RichFutureRecTrackAlgorithms src/*.cpp INCLUDE_DIRS Boost GSL VDT Rich/RichFutureRecBase Rich/RichUtils Rich/RichFutureUtils Tr/TrackInterfaces Tr/TrackKernel - LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichUtils RichFutureUtils) - + LINK_LIBRARIES Boost GSL VDT RichFutureRecBase RichUtils RichFutureUtils ) diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp index 14b9006693a24c30e7e2661373284d7326338442..65b1025f2b73c4393c208e6ce77a3a2affa4ab10 100755 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.cpp @@ -40,7 +40,7 @@ StatusCode BaseTrSegMaker::initialize() if ( !usedRads(Rich::Aerogel) ) { - debug() << "Track segments for Aerogel are disabled" << endmsg; + _ri_debug << "Track segments for Aerogel are disabled" << endmsg; //Warning("Track segments for Aerogel are disabled",StatusCode::SUCCESS).ignore(); //made into debug as it is still too verbose as a warning } diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h index 2aa35f224f33551ac6af02b9eb4c83a1f515cbc7..f4dc7eaec57324037845f00f1ce7938400eb58a6 100755 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichBaseTrSegMaker.h @@ -49,7 +49,7 @@ namespace Rich { //--------------------------------------------------------------------------------- - /** @class BaseTrSegMakerFromRecoTracks RichBaseTrSegMakerFromRecoTracks.h + /** @class BaseTrSegMaker RichBaseTrSegMaker.h * * Base class to tools that create RichTrackSegments from Tracks. * diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h index 66423a3cf4a201d53b5b3b235a472c7bb596070e..67e050fba025351d76cc7855a5483e4243793740 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichGeomEffCKMassRings.h @@ -42,7 +42,7 @@ namespace Rich SegmentPhotonFlags::Vector >; } - /** @class RichGeomEffCKMassRings RichRichGeomEffCKMassRings.h + /** @class GeomEffCKMassRings RichRichGeomEffCKMassRings.h * * Computes the geometrical efficiencies for a set of track segments. * diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp index b39a6afa73fec8fa9a164047f61a2b48c9500f78..64cc4a39e0ff7e59e01987f3cc502ec83dbcf6ca 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.cpp @@ -110,46 +110,106 @@ RayTraceCherenkovCones::operator()( const LHCb::RichTrackSegment::Vector& segmen // reserve size in the points container rings[id].reserve( m_nPoints[rad] ); - // loop around the ring - unsigned int nOK(0), nPhot(0); - for ( auto iP = m_cosSinPhiV[rad].begin(); iP != m_cosSinPhiV[rad].end(); ++iP, ++nPhot ) + // Container for photon directions + std::vector<Gaudi::XYZVector> vects; + vects.reserve( m_nPoints[rad] ); + + // loop around the ring to create the directions + for ( const auto& P : m_cosSinPhiV[rad] ) { // Photon direction around loop - const auto photDir = - segment.vectorAtCosSinThetaPhi( cosTheta, sinTheta, - (*iP).cosPhi, (*iP).sinPhi ); - - // do the ray tracing - // really need to speed this up by removing virtual function call and need for photon. - const auto result = - m_rayTrace.get()->traceToDetector( rich, emissionPoint, photDir, photon, + vects.emplace_back( segment.vectorAtCosSinThetaPhi( cosTheta, sinTheta, + P.cosPhi, P.sinPhi ) ); + } + + // Count the number of good photons + unsigned int nOK(0); + + // Which ray tracing to run + if ( m_useVectorised ) + { + + // The vectorised ray tracing + const auto results = + m_rayTrace.get()->traceToDetector( emissionPoint, vects, segment, m_traceModeRad[rad] ); + + // loop over the results and fill + unsigned int nPhot(0); + for ( const auto && data : Ranges::ConstZip(results,m_cosSinPhiV[rad]) ) + { + const auto & res = std::get<0>(data); + const auto & cosphi = std::get<1>(data); + + // count photons + ++nPhot; + + // Add a new point + const auto & gP = res.detectionPoint; + rings[id].emplace_back ( gP, + corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad), + res.smartID, + (RayTracedCKRingPoint::Acceptance)(res.result), + res.primaryMirror, + res.secondaryMirror, + cosphi.phi ); + + // count raytraces that are in HPD panel + if ( res.result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; } + + // bailout check + if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break; + + } + + } + else + { + + // Loop over the directions and ray trace them using the scalar method + unsigned int nPhot(0); + for ( const auto && data : Ranges::ConstZip(vects,m_cosSinPhiV[rad]) ) + { + const auto & photDir = std::get<0>(data); + const auto & cosphi = std::get<1>(data); + + // count photons + ++nPhot; + + // do the ray tracing + // really need to speed this up by removing virtual function call and need for photon. + const auto result = + m_rayTrace.get()->traceToDetector( rich, emissionPoint, photDir, photon, segment, m_traceModeRad[rad], Rich::top ); // Note forced side is not used.. - // Add a new point - const auto & gP = photon.detectionPoint(); - rings[id].emplace_back ( gP, - corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad), - photon.smartID(), - (RayTracedCKRingPoint::Acceptance)(result), - photon.primaryMirror(), - photon.secondaryMirror(), - (*iP).phi ); - - // count raytraces that are in HPD panel - if ( result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; } - - // bailout check - if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break; - + // Add a new point + const auto & gP = photon.detectionPoint(); + rings[id].emplace_back ( gP, + corrector.correct(m_idTool.get()->globalToPDPanel(gP),rad), + photon.smartID(), + (RayTracedCKRingPoint::Acceptance)(result), + photon.primaryMirror(), + photon.secondaryMirror(), + cosphi.phi ); + + // count raytraces that are in HPD panel + if ( result >= LHCb::RichTraceMode::InHPDPanel ) { ++nOK; } + + // bailout check + if ( 0 == nOK && nPhot >= m_nBailout[rad] ) break; + } + } + // if no good hits empty the container + //if ( 0 == nOK ) { rings[id].clear(); } + } - + } - + } - + return ringsV; } diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h index cb615571789f928efd44f83d4560fb66ece4c2ee..26dd924732da68f56166a598f8a1d9925ae5990f 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichRayTraceCherenkovCones.h @@ -100,7 +100,7 @@ namespace Rich /// Number of points to ray trace on each ring, for each radiator Gaudi::Property< RadiatorArray<unsigned int> > m_nPoints - { this, "NRingPoints", { 100u, 100u, 100u } }; + { this, "NRingPoints", { 96u, 96u, 96u } }; /// Bailout number RadiatorArray<unsigned int> m_nBailout = {{}}; @@ -122,6 +122,10 @@ namespace Rich Gaudi::Property< RadiatorArray<float> > m_bailoutFrac { this, "BailoutFraction", { 0.75, 0.75, 0.75 } }; + /// Use the scalar or vector ray tracing + Gaudi::Property<bool> m_useVectorised + { this, "UseVectorisedRayTracing", true }; + /// Ray tracing tool ToolHandle<const IRayTracing> m_rayTrace { "Rich::Future::RayTracing/RayTracing", this }; diff --git a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h index c8fc19ba0fa986d5687978f700ed08e362e521fc..3072c9e1147d043888cc3b189f2e28feb9f0282c 100644 --- a/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h +++ b/Rich/RichFutureRecTrackAlgorithms/src/RichTrackFunctionalCherenkovResolutions.h @@ -53,7 +53,7 @@ namespace Rich namespace Rec { - /** @class RichTrackFunctionalCherenkovResolutions RichRichTrackFunctionalCherenkovResolutions.h + /** @class TrackFunctionalCherenkovResolutions RichTrackFunctionalCherenkovResolutions.h * * Computes the expected Cherenkov resolutions for the given track segments * diff --git a/Rich/RichRecPhotonTools/CMakeLists.txt b/Rich/RichRecPhotonTools/CMakeLists.txt index c551830b1f52218da53fd8590664a146f8ecf67d..231f8c58a3e2acc3727de49df1b393c585490394 100644 --- a/Rich/RichRecPhotonTools/CMakeLists.txt +++ b/Rich/RichRecPhotonTools/CMakeLists.txt @@ -24,9 +24,3 @@ gaudi_add_module(RichRecPhotonTools INCLUDE_DIRS Eigen Boost GSL VectorClass src LINK_LIBRARIES Boost GSL RichDetLib GaudiKernel LHCbKernel LHCbMathLib RichKernelLib RichRecBase) -if(GAUDI_BUILD_TESTS) - gaudi_add_executable(RichPhotonRecoTest - src/application/PhotonReco/*.cpp - INCLUDE_DIRS ROOT Eigen GSL VectorClass src - LINK_LIBRARIES LHCbMathLib GSL GaudiKernel ROOT) -endif() diff --git a/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp b/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp deleted file mode 100644 index 20a2eb747de69fa4c1f156fe1aecec8d63ada915..0000000000000000000000000000000000000000 --- a/Rich/RichRecPhotonTools/src/application/PhotonReco/main.cpp +++ /dev/null @@ -1,81 +0,0 @@ - -// Quartic Solver -#include "RichRecUtils/QuarticSolver.h" - -// STL -#include <random> -#include <vector> -#include <iostream> -#include <string> -#include <typeinfo> - -// Make an instance of the quartic solver -Rich::Rec::QuarticSolver qSolver; - -Gaudi::XYZPoint sphReflPoint; - -class Data -{ -public: - typedef std::vector<Data> Vector; -public: - Gaudi::XYZPoint emissPnt; - Gaudi::XYZPoint centOfCurv; - Gaudi::XYZPoint virtDetPoint; - double radius; -public: - Data() - { - // randomn generator - static std::default_random_engine gen; - // Distributions for each member - static std::uniform_real_distribution<double> r_emiss_x(-800,800), r_emiss_y(-600,600), r_emiss_z(10000,10500); - static std::uniform_real_distribution<double> r_coc_x(-3000,3000), r_coc_y(-20,20), r_coc_z(3300,3400); - static std::uniform_real_distribution<double> r_vdp_x(-3000,3000), r_vdp_y(-200,200), r_vdp_z(8100,8200); - static std::uniform_real_distribution<float> r_rad(8500,8600); - // setup data - emissPnt = Gaudi::XYZPoint( r_emiss_x(gen), r_emiss_y(gen), r_emiss_z(gen) ); - centOfCurv = Gaudi::XYZPoint( r_coc_x(gen), r_coc_y(gen), r_coc_z(gen) ); - virtDetPoint = Gaudi::XYZPoint( r_vdp_x(gen), r_vdp_y(gen), r_vdp_z(gen) ); - radius = r_rad(gen); - } -}; - -template< class TYPE > -void solve( const Data& data ) -{ - qSolver.solve<TYPE>( data.emissPnt, - data.centOfCurv, - data.virtDetPoint, - data.radius, - sphReflPoint ); -} - -template< class TYPE > -void solve( const Data::Vector & dataV ) -{ - std::cout << "Solving Quartic Equation for " - << typeid(TYPE).name() - << " Photons ..." << std::endl; - // iterate over the data and solve it... - for ( const auto& data : dataV ) { solve<TYPE>(data); } -} - -int main ( int /*argc*/, char** /*argv*/ ) -{ - const unsigned int nPhotons = 1e6; - - Data::Vector dataV; - dataV.reserve( nPhotons ); - // Construct the data to work on - std::cout << "Creating " << nPhotons << " random photons ..." << std::endl; - for ( unsigned int i = 0; i < nPhotons; ++i ) { dataV.push_back( Data() ); } - - // run the solver for floats - solve<float>( dataV ); - - // run the solver for doubles - solve<double>( dataV ); - - return 0; -} diff --git a/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h b/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h index b1618bab6d96e3435962dfab59578025bbc03127..6331185cce1bc015c04b0f4840ae1eb8d5a998b8 100755 --- a/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h +++ b/Rich/RichRecPhotonTools/src/component/RichPhotonRecoUsingQuarticSoln.h @@ -56,7 +56,7 @@ #include "gsl/gsl_poly.h" // Quartic Solver -#include "RichRecUtils/QuarticSolver.h" +#include "RichRecUtils/QuarticSolverNewton.h" namespace Rich { @@ -175,7 +175,7 @@ namespace Rich mutable const ISnellsLawRefraction * m_snellsLaw = nullptr; /// Quartic Solver - QuarticSolver m_quarticSolver; + QuarticSolverNewton m_quarticSolver; /** @brief Flag to indicate if the unambiguous photon test should be performed * for each radiator diff --git a/Rich/RichRecTests/CMakeLists.txt b/Rich/RichRecTests/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..5a43b8ec44fc9b7ec830b5ab2b67558f9359979b --- /dev/null +++ b/Rich/RichRecTests/CMakeLists.txt @@ -0,0 +1,109 @@ +################################################################################ +# Package: RichRecTests +################################################################################ +gaudi_subdir(RichRecTests v1r0) + +gaudi_depends_on_subdirs(Det/RichDet + GaudiKernel + Kernel/LHCbKernel + Kernel/LHCbMath + Kernel/VectorClass + Kernel/TemplatedGenVector + Kernel/Vc + Rich/RichRecUtils) + +#find_package(vectorclass) +find_package(Boost) +find_package(GSL) +find_package(Eigen) +find_package(Vc) +find_package(ROOT COMPONENTS MathCore GenVector ) + +message(STATUS "Using Vc ${Vc_INCLUDE_DIR} ${Vc_LIB_DIR}") +#message(STATUS "Vc DEFINITIONS ${Vc_DEFINITIONS}") +#message(STATUS "Vc COMPILE flags ${Vc_COMPILE_FLAGS}") +#message(STATUS "Vc ARCHITECTURE flags ${Vc_ARCHITECTURE_FLAGS}") + +include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} + ${EIGEN_INCLUDE_DIRS} ${Vc_INCLUDE_DIR}) + +if(GAUDI_BUILD_TESTS) + + gaudi_add_executable(RichPhotonRecoTestSSE4 + src/PhotonReco/main_sse4.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichPhotonRecoTestSSE4 -lrt ) + set_property(SOURCE src/PhotonReco/main_sse4.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -msse4.2 " ) + + gaudi_add_executable(RichPhotonRecoTestAVX + src/PhotonReco/main_avx.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichPhotonRecoTestAVX -lrt ) + set_property(SOURCE src/PhotonReco/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " ) + + gaudi_add_executable(RichPhotonRecoTestAVX2 + src/PhotonReco/main_avx2.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichPhotonRecoTestAVX2 -lrt ) + # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here + #set_property(SOURCE src/PhotonReco/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx2 " ) + set_property(SOURCE src/PhotonReco/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " ) + + gaudi_add_executable(RichPhotonRecoTestAVX2FMA + src/PhotonReco/main_avx2fma.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichPhotonRecoTestAVX2FMA -lrt ) + # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here + #set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mfma -mavx2 " ) + set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -mavx " ) + + + + gaudi_add_executable(RichRayTracingTestSSE4 + src/RayTracing/main_sse4.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichRayTracingTestSSE4 "-lrt ${Vc_LIB_DIR}/libVc.a" ) + set_property(SOURCE src/RayTracing/main_sse4.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -msse4.2 " ) + + gaudi_add_executable(RichRayTracingTestAVX + src/RayTracing/main_avx.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichRayTracingTestAVX "-lrt ${Vc_LIB_DIR}/libVc.a" ) + set_property(SOURCE src/RayTracing/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " ) + + gaudi_add_executable(RichRayTracingTestAVX2 + src/RayTracing/main_avx2.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichRayTracingTestAVX2 "-lrt ${Vc_LIB_DIR}/libVc.a" ) + # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here + #set_property(SOURCE src/RayTracing/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx2 " ) + set_property(SOURCE src/RayTracing/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " ) + + gaudi_add_executable(RichRayTracingTestAVX2FMA + src/RayTracing/main_avx2fma.cpp + INCLUDE_DIRS Kernel/TemplatedGenVector ROOT Vc Eigen GSL VectorClass src + LINK_LIBRARIES LHCbMathLib GSL GaudiKernel TemplatedGenVectorLib ) + target_link_libraries( RichRayTracingTestAVX2FMA "-lrt ${Vc_LIB_DIR}/libVc.a" ) + # Until such a time that all build targets support AVX2+FMA, fallback to just AVX here + #set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mfma -mavx2 " ) + set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -Wno-ignored-attributes -mavx " ) + + if(LCG_COMP STREQUAL "gcc" OR + (BINARY_TAG_COMP_NAME STREQUAL "gcc" AND BINARY_TAG_COMP_VERSION VERSION_LESS "5")) + set_property(SOURCE src/PhotonReco/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/PhotonReco/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/PhotonReco/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/RayTracing/main_avx.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/RayTracing/main_avx2.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + set_property(SOURCE src/RayTracing/main_avx2fma.cpp APPEND_STRING PROPERTY COMPILE_FLAGS " -fabi-version=0 " ) + endif() + + +endif() diff --git a/Rich/RichRecTests/src/PhotonReco/main.icpp b/Rich/RichRecTests/src/PhotonReco/main.icpp new file mode 100644 index 0000000000000000000000000000000000000000..dbc074555b940527d629f5e3882824aa66bbfd69 --- /dev/null +++ b/Rich/RichRecTests/src/PhotonReco/main.icpp @@ -0,0 +1,124 @@ + +// Gaudi +//#include "GaudiKernel/Transform3DTypes.h" +#include "GaudiKernel/Point3DTypes.h" +#include "GaudiKernel/Vector3DTypes.h" + +// Quartic Solver +#include "RichRecUtils/QuarticSolverNewton.h" + +// STL +#include <random> +#include <vector> +#include <iostream> +#include <string> +#include <typeinfo> + +// Make an instance of the quartic solver +Rich::Rec::QuarticSolverNewton qSolver; + +ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<float> > sphReflPointF; +ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<double> > sphReflPointD; + +unsigned long long int time_diff( struct timespec *start, + struct timespec *stop ) +{ + if ((stop->tv_nsec - start->tv_nsec) < 0) { + return 1000000000ULL * (stop->tv_sec - start->tv_sec - 1) + + stop->tv_nsec - start->tv_nsec + 1000000000ULL; + } else { + return 1000000000ULL * (stop->tv_sec - start->tv_sec) + + stop->tv_nsec - start->tv_nsec; + } +} + +template< typename TYPE > +class Data +{ +public: + typedef std::vector<Data> Vector; +public: + ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > emissPnt; + ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > centOfCurv; + ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D<TYPE> > virtDetPoint; + TYPE radius; +public: + Data() + { + // randomn generator + static std::default_random_engine gen; + // Distributions for each member + static std::uniform_real_distribution<TYPE> r_emiss_x(-800,800), r_emiss_y(-600,600), r_emiss_z(10000,10500); + static std::uniform_real_distribution<TYPE> r_coc_x(-3000,3000), r_coc_y(-20,20), r_coc_z(3300,3400); + static std::uniform_real_distribution<TYPE> r_vdp_x(-3000,3000), r_vdp_y(-200,200), r_vdp_z(8100,8200); + static std::uniform_real_distribution<TYPE> r_rad(8500,8600); + // setup data + emissPnt = { r_emiss_x(gen), r_emiss_y(gen), r_emiss_z(gen) }; + centOfCurv = { r_coc_x(gen), r_coc_y(gen), r_coc_z(gen) }; + virtDetPoint = { r_vdp_x(gen), r_vdp_y(gen), r_vdp_z(gen) }; + radius = r_rad(gen); + } +}; + +template< class TYPE, class POINT > +inline void solve( const Data<TYPE>& data, POINT & sphReflPoint ) +{ + qSolver.solve<TYPE,3,2>( data.emissPnt, + data.centOfCurv, + data.virtDetPoint, + data.radius, + sphReflPoint ); +} + +template< class TYPE, class POINT > +unsigned long long int __attribute__((noinline)) +solve( const typename Data<TYPE>::Vector & dataV, POINT & sphReflPoint ) +{ + unsigned long long int best_dur{ 99999999999999999 }; + + timespec start, end; + + const unsigned int nTests = 1000; + for ( unsigned int i = 0; i < nTests; ++i ) + { + clock_gettime(CLOCK_MONOTONIC_RAW,&start); + // iterate over the data and solve it... + for ( const auto& data : dataV ) { solve<TYPE>(data,sphReflPoint); } + clock_gettime(CLOCK_MONOTONIC_RAW,&end); + const auto duration = time_diff(&start,&end); + if ( duration < best_dur ) { best_dur = duration; } + } + + return best_dur ; +} + +int main ( int /*argc*/, char** /*argv*/ ) +{ + const unsigned int nPhotons = 1e5; + + Data<double>::Vector dataVD; + Data<float>::Vector dataVF; + dataVD.reserve( nPhotons ); + dataVF.reserve( nPhotons ); + // Construct the data to work on + std::cout << "Creating " << nPhotons << " random photons ..." << std::endl; + for ( unsigned int i = 0; i < nPhotons; ++i ) + { + dataVD.emplace_back( ); + dataVF.emplace_back( ); + } + + // run the solver for doubles + auto dTime = solve<double>( dataVD, sphReflPointD ); + std::cout << "Double " << dTime << std::endl; + + // run the solver for floats + auto fTime = solve<float>( dataVF, sphReflPointF ); + std::cout << "Float " << fTime << std::endl; + + // make sure we are not optimized away + asm volatile ("" : "+x"(fTime)); + asm volatile ("" : "+x"(dTime)); + + return 0; +} diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/PhotonReco/main_avx.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/PhotonReco/main_avx2.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp b/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/PhotonReco/main_avx2fma.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp b/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/PhotonReco/main_sse4.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/RayTracing/main.icpp b/Rich/RichRecTests/src/RayTracing/main.icpp new file mode 100644 index 0000000000000000000000000000000000000000..8fb4bbf6928ceeb611dc2373fc4bdb2a26c3195e --- /dev/null +++ b/Rich/RichRecTests/src/RayTracing/main.icpp @@ -0,0 +1,270 @@ + +// STL +#include <random> +#include <vector> +#include <iostream> +#include <string> +#include <typeinfo> + +// LHCbMath +#include "LHCbMath/EigenTypes.h" +//#include "LHCbMath/VectorClassTypes.h" + +// LHCb TemplatedGenVector +#include "TemplatedGenVector/DisplacementVector3D.h" +#include "TemplatedGenVector/PositionVector3D.h" +#include "TemplatedGenVector/Plane3D.h" + +// Rich Utils +#include "RichUtils/RichRayTracingUtils.h" +#include "RichUtils/ZipRange.h" + +#include <Vc/Vc> + +unsigned long long int time_diff( struct timespec *start, + struct timespec *stop ) +{ + if ((stop->tv_nsec - start->tv_nsec) < 0) { + return 1000000000ULL * (stop->tv_sec - start->tv_sec - 1) + + stop->tv_nsec - start->tv_nsec + 1000000000ULL; + } else { + return 1000000000ULL * (stop->tv_sec - start->tv_sec) + + stop->tv_nsec - start->tv_nsec; + } +} + +// randomn generator +static std::default_random_engine gen; +// Distributions for each member +static std::uniform_real_distribution<double> p_x(-800,800), p_y(-600,600), p_z(10000,10500); +static std::uniform_real_distribution<double> d_x(-0.2,0.2), d_y(-0.1,0.1), d_z(0.95,0.99); +static std::uniform_real_distribution<double> c_x(3100,3200), c_y(10,15), c_z(3200,3300); +static std::uniform_real_distribution<double> r_rad(8500,8600); +static std::uniform_real_distribution<double> p0(-0.002,0.002), p1(-0.2,0.2), p2(0.97,0.99), p3(-1300,1300); + +template < typename POINT, typename VECTOR, typename PLANE, typename FTYPE > +class Data +{ +public: + typedef std::vector< Data, Vc::Allocator<Data> > Vector; +public: + POINT position; + VECTOR direction; + POINT CoC; + PLANE plane; + FTYPE radius{0}; +public: + template< typename INDATA > + Data( const INDATA & ind ) + : position(ind.position), + direction(ind.direction), + CoC(ind.CoC), + plane(ind.plane), + radius(ind.radius) { } + Data( ) + : position ( p_x(gen), p_y(gen), p_z(gen) ), + direction ( d_x(gen), d_y(gen), d_z(gen) ), + CoC ( c_x(gen), c_y(gen), c_z(gen) ), + plane ( p0(gen), p1(gen), p2(gen), p3(gen) ), + radius ( r_rad(gen) ) { } +}; + +template < typename DATA > +inline void trace( DATA & dataV ) +{ + using namespace Rich::RayTracingUtils; + for ( auto & data : dataV ) + { + reflectSpherical( data.position, data.direction, data.CoC, data.radius ); + reflectPlane ( data.position, data.direction, data.plane ); + //std::cout << data.position << std::endl; + } +} + + +template < typename DATA > +unsigned long long int __attribute__((noinline)) rayTrace( const DATA & in_dataV ) +{ + timespec start, end; + + unsigned long long int best_dur{ 99999999999999999 }; + + const unsigned int nTests = 10000; + //std::cout << "Running " << nTests << " tests on " << in_dataV.size() << " data objects" << std::endl; + DATA dataV; // working data copy; + for ( unsigned int i = 0; i < nTests; ++i ) + { + dataV = in_dataV; // copy data, as it will be modified + clock_gettime(CLOCK_MONOTONIC_RAW,&start); + trace( dataV ); + clock_gettime(CLOCK_MONOTONIC_RAW,&end); + const auto duration = time_diff(&start,&end); + if ( duration < best_dur ) { best_dur = duration; } + } + + return best_dur ; +} + +int main ( int /*argc*/, char** /*argv*/ ) +{ + + const unsigned int nPhotons = 96 * 1e2; + //const unsigned int nPhotons = 8; + std::cout << "Creating " << nPhotons << " random photons ..." << std::endl; + + // Templated GenVector + using LHCbXYZPointD = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<double>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using LHCbXYZVectorD = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<double>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using LHCbXYZPointF = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<float>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using LHCbXYZVectorF = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<float>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using LHCbPlane3DD = LHCbROOT::Math::Impl::Plane3D<double>; + using LHCbPlane3DF = LHCbROOT::Math::Impl::Plane3D<float>; + + // Vc + using VcXYZPointD = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<Vc::double_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using VcXYZVectorD = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<Vc::double_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using VcXYZPointF = LHCbROOT::Math::PositionVector3D<LHCbROOT::Math::Cartesian3D<Vc::float_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using VcXYZVectorF = LHCbROOT::Math::DisplacementVector3D<LHCbROOT::Math::Cartesian3D<Vc::float_v>,LHCbROOT::Math::DefaultCoordinateSystemTag>; + using VcPlane3DD = LHCbROOT::Math::Impl::Plane3D<Vc::double_v>; + using VcPlane3DF = LHCbROOT::Math::Impl::Plane3D<Vc::float_v>; + + //if ( false ) + { + + // Data in 'ROOT' SMatrix double format + Data<LHCbXYZPointD,LHCbXYZVectorD,LHCbPlane3DD,LHCbXYZPointD::Scalar>::Vector root_dataV( nPhotons ); + auto resROOT = rayTrace( root_dataV ); + std::cout << "ROOT double " << resROOT << std::endl; + + // Same data in Eigen double format + Data<LHCb::Math::Eigen::XYZPoint<double>,LHCb::Math::Eigen::XYZVector<double>,LHCb::Math::Eigen::Plane3D<double>,double>::Vector eigen_dataV( nPhotons ); + auto resEIGEN = rayTrace( eigen_dataV ); + std::cout << "Eigen double " << resEIGEN << " speedup " << (float)resROOT/(float)resEIGEN << std::endl; + + // Same data in VectorClass double format + //Data<LHCb::Math::VectorClass::XYZPoint<double>,LHCb::Math::VectorClass::XYZVector<double>,double>::Vector vclass_dataV( nPhotons ); + //auto resVClass = rayTrace( vclass_dataV ); + //std::cout << "VectorClass double " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl; + + // Vc float + Data<VcXYZPointD,VcXYZVectorD,VcPlane3DD,Vc::double_v>::Vector VC_dataV( nPhotons / Vc::double_v::Size ); + auto resVC = rayTrace( VC_dataV ); + std::cout << "Vc double " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl; + + // make sure we are not optimized away + asm volatile ("" : "+x"(resROOT)); + asm volatile ("" : "+x"(resEIGEN)); + //asm volatile ("" : "+x"(resVClass)); + asm volatile ("" : "+x"(resVC)); + + } + + //if ( false ) + { + + // Data in 'ROOT' SMatrix float format + Data<LHCbXYZPointF,LHCbXYZVectorF,LHCbPlane3DF,LHCbXYZPointF::Scalar>::Vector root_dataV( nPhotons ); + auto resROOT = rayTrace( root_dataV ); + std::cout << "ROOT float " << resROOT << std::endl; + + // Same data in Eigen float format + Data<LHCb::Math::Eigen::XYZPoint<float>,LHCb::Math::Eigen::XYZVector<float>,LHCb::Math::Eigen::Plane3D<float>,float>::Vector eigen_dataV( nPhotons ); + auto resEIGEN = rayTrace( eigen_dataV ); + std::cout << "Eigen float " << resEIGEN << " speedup " << (float)resROOT/(float)resEIGEN << std::endl; + + // Same data in VectorClass float format + //Data<LHCb::Math::VectorClass::XYZPoint<float>,LHCb::Math::VectorClass::XYZVector<float>,float>::Vector vclass_dataV( nPhotons ); + //auto resVClass = rayTrace( vclass_dataV ); + //std::cout << "VectorClass float " << resVClass << " speedup " << (float)resROOT/(float)resVClass << std::endl; + + // Vc float + Data<VcXYZPointF,VcXYZVectorF,VcPlane3DF,Vc::float_v>::Vector VC_dataV( nPhotons / Vc::float_v::Size ); + auto resVC = rayTrace( VC_dataV ); + std::cout << "Vc float " << resVC << " speedup " << (float)resROOT/(float)resVC << std::endl; + + // make sure we are not optimized away + asm volatile ("" : "+x"(resROOT)); + asm volatile ("" : "+x"(resEIGEN)); + //asm volatile ("" : "+x"(resVClass)); + asm volatile ("" : "+x"(resVC)); + + } + + // Basic checks + if ( false ) + { + + std::vector< Vc::float_v > a = { 1, 2, 3, 4, 5, 6, 7, 8 }; + std::vector< Vc::float_v > b = { 1, 2, 3, 4, 5, 6, 7, 8 }; + + std::cout << "Size " << a.size() << " " << Vc::float_v::Size << std::endl; + + for ( unsigned int i = 0; i < (a.size()/Vc::float_v::Size); ++i ) + { + a[i] += b[i]; + std::cout << "Filling " << i << " " << a[i] << " " << b[i] << std::endl; + } + + for ( unsigned int i = 0; i < a.size(); ++i ) + { + std::cout << "Reading " << i << " " << a[i] << " " << b[i] << std::endl; + } + + } + + + // VC example + if ( false ) + { + + using Vc::float_v; + //! [includes] + + //! [memory allocation] + + // allocate memory for our initial x and y coordinates. Note that you can also put it into a + // normal float C-array but that you then must ensure alignment to Vc::VectorAlignment! + Vc::Memory<float_v, 100> x_mem; + Vc::Memory<float_v, 100> y_mem; + Vc::Memory<float_v, 100> r_mem; + Vc::Memory<float_v, 100> phi_mem; + //! [memory allocation] + + //! [random init] + // fill the memory with values from -1.f to 1.f + std::cout << "vectorsCount = " << x_mem.vectorsCount() << std::endl; + for (size_t i = 0; i < x_mem.vectorsCount(); ++i) + { + x_mem.vector(i) = float_v::Random() * 2.f - 1.f; + y_mem.vector(i) = float_v::Random() * 2.f - 1.f; + } + //! [random init] + + //! [conversion] + // calculate the polar coordinates for all coordinates and overwrite the euclidian coordinates + // with the result + for (size_t i = 0; i < x_mem.vectorsCount(); ++i) + { + const float_v x = x_mem.vector(i); + const float_v y = y_mem.vector(i); + + r_mem.vector(i) = Vc::sqrt(x * x + y * y); + float_v phi = Vc::atan2(y, x) * 57.295780181884765625f; // 180/pi + phi(phi < 0.f) += 360.f; + phi_mem.vector(i) = phi; + } + //! [conversion] + + //! [output] + // print the results + std::cout << "entriesCount = " << x_mem.entriesCount() << std::endl; + for (size_t i = 0; i < x_mem.entriesCount(); ++i) { + std::cout << std::setw(3) << i << ": "; + std::cout << std::setw(10) << x_mem[i] << ", " << std::setw(10) << y_mem[i] << " -> "; + std::cout << std::setw(10) << r_mem[i] << ", " << std::setw(10) << phi_mem[i] << '\n'; + } + + } + + return 0; +} diff --git a/Rich/RichRecTests/src/RayTracing/main_avx.cpp b/Rich/RichRecTests/src/RayTracing/main_avx.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/RayTracing/main_avx.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/RayTracing/main_avx2.cpp b/Rich/RichRecTests/src/RayTracing/main_avx2.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/RayTracing/main_avx2.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp b/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/RayTracing/main_avx2fma.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecTests/src/RayTracing/main_sse4.cpp b/Rich/RichRecTests/src/RayTracing/main_sse4.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c6736566335e5e6e0bbcf64768dbf6d083437d1e --- /dev/null +++ b/Rich/RichRecTests/src/RayTracing/main_sse4.cpp @@ -0,0 +1,2 @@ + +#include "main.icpp" diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h index 4fbabe905e080b2efc52b75a53859659f7167d8f..067801e2169127a263cd49a9eed6e7d105459503 100755 --- a/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h +++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolver.h @@ -7,8 +7,7 @@ */ //---------------------------------------------------------------------- -#ifndef RICHRECPHOTONTOOLS_QuarticSolver_H -#define RICHRECPHOTONTOOLS_QuarticSolver_H 1 +#pragma once // Gaudi #include "GaudiKernel/Kernel.h" @@ -52,100 +51,6 @@ namespace Rich class QuarticSolver final { - public: - - // Use eigen types - typedef LHCb::Math::Eigen::XYZPoint Point; ///< Point type - typedef LHCb::Math::Eigen::XYZVector Vector; ///< vector type - - public: - - /** Solves the characteristic quartic equation for the RICH optical system. - * - * See note LHCB/98-040 RICH section 3 for more details - * - * @param emissionPoint Assumed photon emission point on track - * @param CoC Spherical mirror centre of curvature - * @param virtDetPoint Virtual detection point - * @param radius Spherical mirror radius of curvature - * @param sphReflPoint The reconstructed reflection pont on the spherical mirror - * - * @return boolean indicating status of the quartic solution - * @retval true Calculation was successful. sphReflPoint is valid. - * @retval false Calculation failed. sphReflPoint is not valid. - */ - template< class TYPE > - inline void solve( const Gaudi::XYZPoint& emissionPoint, - const Gaudi::XYZPoint& CoC, - const Gaudi::XYZPoint& virtDetPoint, - const TYPE radius, - Gaudi::XYZPoint& sphReflPoint ) const - { - - // typedefs vectorised types - typedef Eigen::Matrix< TYPE , 3 , 1 > Eigen3Vector; - using Vec4x = - typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type; - - // vector from mirror centre of curvature to assumed emission point - const Vector evec( emissionPoint - CoC ); - const TYPE e2 = evec.dot(evec); - - // vector from mirror centre of curvature to virtual detection point - const Vector dvec( virtDetPoint - CoC ); - const TYPE d2 = dvec.dot(dvec); - - // various quantities needed to create quartic equation - // see LHCB/98-040 section 3, equation 3 - const auto ed2 = e2 * d2; - const TYPE cosgamma2 = ( ed2 > 0 ? std::pow(evec.dot(dvec),2)/ed2 : 1.0 ); - // vectorise 4 square roots into 1 - const auto tmp_sqrt = sqrt( Vec4x( e2, d2, - cosgamma2 < 1.0 ? 1.0-cosgamma2 : 0.0, - cosgamma2 ) ); - const auto e = tmp_sqrt[0]; - const auto d = tmp_sqrt[1]; - const auto singamma = tmp_sqrt[2]; - const auto cosgamma = tmp_sqrt[3]; - - const auto dx = d * cosgamma; - const auto dy = d * singamma; - const auto r2 = radius * radius; - const auto dy2 = dy * dy; - const auto edx = e + dx; - - // Fill array for quartic equation - const auto a0 = 4.0 * ed2; - const auto inv_a0 = ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() ); - const auto dyrad2 = 2.0 * dy * radius; - const auto a1 = - ( 2.0 * dyrad2 * e2 ) * inv_a0; - const auto a2 = ( (dy2 * r2) + ( edx * edx * r2 ) - a0 ) * inv_a0; - const auto a3 = ( dyrad2 * e * (e-dx) ) * inv_a0; - const auto a4 = ( ( e2 - r2 ) * dy2 ) * inv_a0; - - // use simplified RICH version of quartic solver - const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a - a2, // b - a3, // c - a4 // d - ); - - // (normalised) normal vector to reflection plane - auto n = evec.cross3(dvec); - n /= std::sqrt( n.dot(n) ); - - // construct rotation transformation - // Set vector magnitude to radius - // rotate vector and update reflection point - const Eigen::AngleAxis<TYPE> angleaxis( vdt::fast_asinf(sinbeta), - Eigen3Vector(n[0],n[1],n[2]) ); - sphReflPoint = ( CoC + - Gaudi::XYZVector( angleaxis * - Eigen3Vector(evec[0],evec[1],evec[2]) * - ( radius / e ) ) ); - - } - private: /// The cube root implementaton to use @@ -169,10 +74,10 @@ namespace Rich */ //---------------------------------------------------------------------- template < class TYPE > - inline TYPE solve_quartic_RICH ( const TYPE& a, - const TYPE& b, - const TYPE& c, - const TYPE& d ) const + inline TYPE solve_quartic_RICH ( const TYPE a, + const TYPE b, + const TYPE c, + const TYPE d ) const { constexpr const TYPE r4 = 1.0 / 4.0; @@ -237,9 +142,97 @@ namespace Rich res ); } + public: + + /** Solves the characteristic quartic equation for the RICH optical system. + * + * See note LHCB/98-040 RICH section 3 for more details + * + * @param emissionPoint Assumed photon emission point on track + * @param CoC Spherical mirror centre of curvature + * @param virtDetPoint Virtual detection point + * @param radius Spherical mirror radius of curvature + * @param sphReflPoint The reconstructed reflection pont on the spherical mirror + * + * @return boolean indicating status of the quartic solution + * @retval true Calculation was successful. sphReflPoint is valid. + * @retval false Calculation failed. sphReflPoint is not valid. + */ + template< class TYPE > + inline void solve( const Gaudi::XYZPoint& emissionPoint, + const Gaudi::XYZPoint& CoC, + const Gaudi::XYZPoint& virtDetPoint, + const TYPE radius, + Gaudi::XYZPoint& sphReflPoint ) const + { + + // typedefs vectorised types + using Eigen3Vector = Eigen::Matrix< TYPE , 3 , 1 >; + using Vector = LHCb::Math::Eigen::XYZVector<TYPE>; + using Point = LHCb::Math::Eigen::XYZPoint<TYPE>; + using Vec4x = + typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type; + + // vector from mirror centre of curvature to assumed emission point + const Vector evec( emissionPoint - CoC ); + const TYPE e2 = evec.dot(evec); + + // vector from mirror centre of curvature to virtual detection point + const Vector dvec( virtDetPoint - CoC ); + const TYPE d2 = dvec.dot(dvec); + + // various quantities needed to create quartic equation + // see LHCB/98-040 section 3, equation 3 + const auto ed2 = e2 * d2; + const TYPE cosgamma2 = ( ed2 > 0 ? std::pow(evec.dot(dvec),2)/ed2 : 1.0 ); + // vectorise 4 square roots into 1 + const auto tmp_sqrt = sqrt( Vec4x( e2, d2, + cosgamma2 < 1.0 ? 1.0-cosgamma2 : 0.0, + cosgamma2 ) ); + const auto e = tmp_sqrt[0]; + const auto d = tmp_sqrt[1]; + const auto singamma = tmp_sqrt[2]; + const auto cosgamma = tmp_sqrt[3]; + + const auto dx = d * cosgamma; + const auto dy = d * singamma; + const auto r2 = radius * radius; + const auto dy2 = dy * dy; + const auto edx = e + dx; + + // Fill array for quartic equation + const auto a0 = 4.0 * ed2; + const auto inv_a0 = ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() ); + const auto dyrad2 = 2.0 * dy * radius; + const auto a1 = - ( 2.0 * dyrad2 * e2 ) * inv_a0; + const auto a2 = ( (dy2 * r2) + ( edx * edx * r2 ) - a0 ) * inv_a0; + const auto a3 = ( dyrad2 * e * (e-dx) ) * inv_a0; + const auto a4 = ( ( e2 - r2 ) * dy2 ) * inv_a0; + + // use simplified RICH version of quartic solver + const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a + a2, // b + a3, // c + a4 // d + ); + + // (normalised) normal vector to reflection plane + auto n = evec.cross3(dvec); + n /= std::sqrt( n.dot(n) ); + + // construct rotation transformation + // Set vector magnitude to radius + // rotate vector and update reflection point + const Eigen::AngleAxis<TYPE> angleaxis( vdt::fast_asinf(sinbeta), + Eigen3Vector(n[0],n[1],n[2]) ); + sphReflPoint = ( CoC + + Gaudi::XYZVector( angleaxis * + Eigen3Vector(evec[0],evec[1],evec[2]) * + ( radius / e ) ) ); + + } + }; } } - -#endif // RICHRECPHOTONTOOLS_QuarticSolver_H diff --git a/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h new file mode 100755 index 0000000000000000000000000000000000000000..dfc5da3a9f3e2e8f388a6d1de300b0ad1acb6481 --- /dev/null +++ b/Rich/RichRecUtils/RichRecUtils/QuarticSolverNewton.h @@ -0,0 +1,324 @@ + +//---------------------------------------------------------------------- +/** @file QuarticSolverNewton.h + * + * @author Christina Quast, Rainer Schwemmer + * @date 2017-02-03 + */ +//---------------------------------------------------------------------- + +#pragma once + +// Gaudi +#include "GaudiKernel/Kernel.h" + +// VectorClass +#include "VectorClass/vectorclass.h" + +// STL +#include <math.h> +#include <type_traits> +#include <array> + +namespace Rich +{ + namespace Rec + { + + //----------------------------------------------------------------------------- + /** @class QuarticSolverNewton + * + * Utility class that implements the solving of the Quartic equation for the RICH + * Based on original code by Chris Jones + * + * @author Christina Quast, Rainer Schwemmer + * @date 2017-02-03 + */ + //----------------------------------------------------------------------------- + class QuarticSolverNewton + { + + public: + + /** Solves the characteristic quartic equation for the RICH optical system. + * + * See note LHCB/98-040 RICH section 3 for more details + * + * @param emissionPoint Assumed photon emission point on track + * @param CoC Spherical mirror centre of curvature + * @param virtDetPoint Virtual detection point + * @param radius Spherical mirror radius of curvature + * @param sphReflPoint The reconstructed reflection pont on the spherical mirror + * + * @return boolean indicating status of the quartic solution + * @retval true Calculation was successful. sphReflPoint is valid. + * @retval false Calculation failed. sphReflPoint is not valid. + */ + template< class TYPE, + std::size_t BISECTITS = 3, std::size_t NEWTONITS = 4, + class POINT = Gaudi::XYZPoint > + inline void solve( const POINT& emissionPoint, + const POINT& CoC, + const POINT& virtDetPoint, + const TYPE radius, + POINT& sphReflPoint ) const + { + using namespace std; + + // vector from mirror centre of curvature to assumed emission point + const auto evec( emissionPoint - CoC ); + const TYPE e2 = evec.Dot(evec); + + // vector from mirror centre of curvature to virtual detection point + const auto dvec( virtDetPoint - CoC ); + const TYPE d2 = dvec.Dot(dvec); + + // various quantities needed to create quartic equation + // see LHCB/98-040 section 3, equation 3 + const TYPE ed2 = e2 * d2; + const TYPE cosgamma2 = ( ed2 > TYPE(0.0) ? pow(evec.Dot(dvec),2)/ed2 : TYPE(1.0) ); + // vectorise 4 square roots into 1 + using Vec4x = typename std::conditional<std::is_same<TYPE,float>::value,Vec4f,Vec4d>::type; + const auto tmp_sqrt = sqrt( Vec4x( e2, d2, + cosgamma2 < TYPE(1.0) ? TYPE(1.0)-cosgamma2 : TYPE(0.0), + cosgamma2 ) ); + const TYPE e = tmp_sqrt[0]; + const TYPE d = tmp_sqrt[1]; + const TYPE singamma = tmp_sqrt[2]; + const TYPE cosgamma = tmp_sqrt[3]; + // scalar versions + //const TYPE e = std::sqrt(e2); + //const TYPE d = std::sqrt(d2); + //const TYPE singamma = std::sqrt( cosgamma2 < TYPE(1.0) ? TYPE(1.0)-cosgamma2 : TYPE(0.0) ); + //const TYPE cosgamma = std::sqrt(cosgamma2); + + const TYPE dx = d * cosgamma; + const TYPE dy = d * singamma; + const TYPE r2 = radius * radius; + const TYPE dy2 = dy * dy; + const TYPE edx = e + dx; + + // Fill array for quartic equation + const TYPE a0 = TYPE(4.0) * ed2; + //Newton solver doesn't care about a0 being not 1.0. Remove costly division and several multiplies. + //This has some downsides though. The a-values are hovering around a numerical value of 10^15. single + //precision float max is 10^37. A single square and some multiplies will push it over the limit of what + //single precision float can handle. It's ok for the newton method, but Halley or higher order Housholder + //will fail without this normalization. + //const auto inv_a0 = ( a0 > 0 ? 1.0 / a0 : std::numeric_limits<TYPE>::max() ); + const TYPE dyrad2 = TYPE(2.0) * dy * radius; + const TYPE a1 = - ( TYPE(2.0) * dyrad2 * e2 );// * inv_a0; + const TYPE a2 = ( (dy2 * r2) + ( edx * edx * r2 ) - a0 );// * inv_a0; + const TYPE a3 = ( dyrad2 * e * (e-dx) );// * inv_a0; + const TYPE a4 = ( ( e2 - r2 ) * dy2 );// * inv_a0; + + // use simplified RICH version of quartic solver + // const auto sinbeta = solve_quartic_RICH<TYPE>( a1, // a + // a2, // b + // a3, // c + // a4 // d + // ); + + //Use optimized newton solver on quartic equation. + const auto sinbeta = solve_quartic_newton_RICH<TYPE,BISECTITS,NEWTONITS>( a0, a1, a2, a3, a4 ); + + //TODO: This method should be better but has problems still for some reasons + //const auto sinbeta = solve_quartic_housholder_RICH<TYPE, 3>(a1, a2, a3, a4); + + // construct rotation transformation + // Set vector magnitude to radius + // rotate vector and update reflection point + //rotation matrix uses sin(beta) and cos(beta) to perform rotation + //even fast_asinf (which is only single precision and defeats the purpose + //of this class being templatable to double btw) is still too slow + //plus there is a cos and sin call inside AngleAxis ... + //We can do much better by just using the cos(beta) we already have to calculate + //sin(beta) and do our own rotation. On top of that we rotate non-normalized and save several + //Divisions by normalizing only once at the very end + //Again, care has to be taken since we are close to float_max here without immediate normalization. + //As far as we have tried with extreme values in the rich coordinate systems this is fine. + const TYPE nx = (evec.y()*dvec.z()) - (evec.z()*dvec.y()); + const TYPE ny = (evec.z()*dvec.x()) - (evec.x()*dvec.z()); + const TYPE nz = (evec.x()*dvec.y()) - (evec.y()*dvec.x()); + + const TYPE nx2 = nx*nx; + const TYPE ny2 = ny*ny; + const TYPE nz2 = nz*nz; + + const TYPE norm = nx2 + ny2 + nz2; + const TYPE norm_sqrt = sqrt(norm); + + const TYPE a = sinbeta * norm_sqrt; + const TYPE sinbeta2 = sinbeta * sinbeta; + const TYPE b = ( sinbeta2 < TYPE(1.0) ? ( TYPE(1.0) - sqrt( TYPE(1.0) - sinbeta2 ) ) : TYPE(1.0) ); + const TYPE enorm = radius/(e*norm); + + const TYPE bnxny = b*nx*ny; + const TYPE bnxnz = b*nx*nz; + const TYPE bnynz = b*ny*nz; + + //Perform non-normalized rotation + const std::array<TYPE,9> M = + { norm - b*(nz2+ny2), a*nz+bnxny, -a*ny+bnxnz, + -a*nz+bnxny, norm - b*(nx2+nz2), a*nx+bnynz, + a*ny+bnxnz, -a*nx+bnynz, norm-b*(ny2+nx2) }; + + //re-normalize rotation and scale to radius in one step + const TYPE ex = enorm*(evec.x()*M[0]+evec.y()*M[3]+evec.z()*M[6]); + const TYPE ey = enorm*(evec.x()*M[1]+evec.y()*M[4]+evec.z()*M[7]); + const TYPE ez = enorm*(evec.x()*M[2]+evec.y()*M[5]+evec.z()*M[8]); + + //sphReflPoint = ( CoC + Gaudi::XYZVector( ex, ey, ez ) ); + sphReflPoint = { CoC.x() + ex, CoC.y() + ey, CoC.z() + ez }; + } + + private: + + //A newton iteration solver for the Rich quartic equation + //Since the polynomial that is evaluated here is extremely constrained + //(root is in small interval, one root guaranteed), we can use a much more + //efficient approximation (which still has the same precision) instead of the + //full blown mathematically absolute correct method and still end up with + //usable results + + template < class TYPE > + inline TYPE f4( const TYPE &a0, const TYPE &a1, const TYPE &a2, + const TYPE &a3, const TYPE &a4, TYPE &x ) const + { + //return (a0 * x*x*x*x + a1 * x*x*x + a2 * x*x + a3 * x + a4); + //A bit more FMA friendly + return ( ( ( ( ( ( ( a0 * x ) + a1 ) * x ) + a2 ) * x ) + a3 ) * x ) + a4; + } + + template < class TYPE > + inline TYPE df4( const TYPE &a0, const TYPE &a1, const TYPE &a2, const TYPE &a3, TYPE &x ) const + { + //return (4.0f*a0 * x*x*x + 3.0f*a1 * x*x + 2.0f*a2 * x + a3); + return ( ( ( ( TYPE(4.0) * a0*x ) + TYPE(3.0) * a1 ) * x + ( TYPE(2.0) * a2 ) ) * x ) + a3; + } + + /** Horner's method to evaluate the polynomial and its derivatives with as little math operations as + * possible. We use a template here to allow the compiler to unroll the for loops and produce code + * that is free from branches and optimized for the grade of polynomial and derivatives as necessary. + */ + template < class TYPE, std::size_t ORDER = 4, std::size_t DIFFGRADE = 3 > + inline void evalPolyHorner( const TYPE (&a)[ORDER+1], TYPE (&res)[DIFFGRADE+1], const TYPE x ) const + { + for ( std::size_t i = 0; i <= DIFFGRADE; ++i ) + { + res[i] = a[0]; + } + + for ( std::size_t j = 1; j <= ORDER; ++j ) + { + res[0] = res[0] * x + a[j]; + const int l = (ORDER - j) > DIFFGRADE ? DIFFGRADE : ORDER - j; + for ( int i = 1; i <= l ; ++i ) + { + res[i] = res[i] * x + res[i-1]; + } + } + + //TODO: Check assembly to see if this is optimized away if DIFFGRADE is <= 2 + TYPE l = 1.0; + for ( std::size_t i = 2; i <= DIFFGRADE; ++i ) + { + l *= i; + res[i] = res[i] * l; + } + + } + + /** 3rd grade Housholder's method for iteratively finding the root of a function. Tests have shown that we seem to be + * already too constrained in our polynomial and input data that we are not gaining anything over newton. In fact it seems + * to make performance worse. + * TODO: find out why performance of this is so bad. We might be missing something. It should converge much faster than newton. + */ + template < class TYPE, std::size_t ITER = 3 > + inline TYPE householder(const TYPE (&a)[5], TYPE x0) const + { + TYPE res[4]; + for ( std::size_t i = 0; i < ITER; ++i ) + { + evalPolyHorner<TYPE, 4, 3>(a, res, x0); + x0 -= ((TYPE(6.0) * res[0]*res[1]*res[1] - TYPE(3.0) * res[0]*res[0]*res[2])/ + (TYPE(6.0) * res[1]*res[1]*res[1] - TYPE(6.0) * res[0]*res[1]*res[2] + res[0]*res[0]*res[3])); + } + return x0; + } + + /** Use Housholder's method to solve the rich quartic equation. Important! If this function is used, the coefficients of the + * polynomial have to be normalized. They are typically around 10^15 and floating point overflows will occur in here if these + * large values are used. + */ + template < class TYPE, std::size_t ITER = 3 > + inline TYPE solve_quartic_housholder_RICH( const TYPE& a0, + const TYPE& a1, + const TYPE& a2, + const TYPE& a3) const + { + //Starting value for housholder iteration. Empirically values seem to be + //between 0 0and 0.4 so we chose the value in the middle as starting point + TYPE x0 = TYPE(0.2); + TYPE a[5] = { TYPE(1.0), a0, a1, a2, a3 }; + return householder<TYPE, ITER>( a, x0 ); + } + + /** Newton-Rhapson method for calculating the root of the rich polynomial. It uses the bisection method in the beginning + * to get close enough to the root to allow the second stage newton method to converge faster. After 4 iterations of newton + * precision is as good as single precision floating point will get you. We have introduced a few tuning parameters like the + * newton gain factor and a slightly skewed bisection division, which in this particular case help to speed things up. + * TODO: Once we are happy with the number of newton and bisection iterations, this function should be templated to the number of + * these iterations to allow loop unrolling and elimination of unnecessary branching + * TODO: These tuning parameters have been found by low effort experimentation on random input data. A more detailed + * study should be done with real data to find the best values + */ + template < class TYPE, std::size_t BISECTITS = 3, std::size_t NEWTONITS = 4 > + inline TYPE solve_quartic_newton_RICH( const TYPE& a0, + const TYPE& a1, + const TYPE& a2, + const TYPE& a3, + const TYPE& a4 ) const + { + //Use N steps of bisection method to find starting point for newton + TYPE l(0.0); + TYPE u(0.5); + TYPE m(0.2); //We start a bit off center since the distribution of roots tends to be more to the left side + const TYPE a[5] = {a0, a1, a2, a3, a4}; + TYPE res[2]; + for ( std::size_t i = 0; i < BISECTITS; ++i ) + { + const auto oppositeSign = std::signbit(f4(a0,a1,a2,a3,a4,m) * f4(a0,a1,a2,a3,a4,l)); + + l = oppositeSign ? l : m; + u = oppositeSign ? m : u; + //0.4 instead of 0.5 to speed up convergence. Most roots seem to be closer to 0 than to the extreme end + m = ( u + l ) * TYPE(0.4); + } + + //Newton for the rest + // CRJ : why copy the value. original value is not needed later on... + //TYPE x = m; + + //Most of the times we are approaching the root of the polynomial from one side + //and fall short by a certain fraction. This fraction seems to be around 1.04 of the + //quotient which is subtracted from x. By scaling it up, we take bigger steps towards + //the root and thus converge faster. + //TODO: study this factor more closely it's pure guesswork right now. We might get + //away with 3 iterations if we can find an exact value + const TYPE gain = 1.04; + + for ( std::size_t i = 0; i < NEWTONITS; ++i ) + { + evalPolyHorner<TYPE, 4, 1>( a, res, m ); + //epsilon = f4(a0,a1,a2,a3,a4,x) / df4(a0,a1,a2,a3,x); + m -= gain * ( res[0] / res[1] ); + } + + return m; + } + + }; + + } +} diff --git a/Tr/TrackFitEvent/doc/release.notes b/Tr/TrackFitEvent/doc/release.notes index ebe90661e0e49d37af3d3b2cc72eaedb6d3a43a0..861836c30abee94a39ae7d826356a33652c15d5e 100755 --- a/Tr/TrackFitEvent/doc/release.notes +++ b/Tr/TrackFitEvent/doc/release.notes @@ -4,6 +4,9 @@ ! Purpose : Package for the tracking fitting-related classes !----------------------------------------------------------------------------- +! 2016-12-14 - Jeroen van Tilburg + - Adapt to new FT geometry and numbering scheme + !========================= TrackFitEvent v6r6 2016-01-28 ========================= ! 2015-12-05 - Gerhard Raven - reduce mention of 'static' and 'std::auto_ptr' diff --git a/Tr/TrackFitEvent/src/FTMeasurement.cpp b/Tr/TrackFitEvent/src/FTMeasurement.cpp index d18b80557421aa60fff41c250d7f996eaf8583f5..bb21204db05b325ee4078aae0a1adb8e43868d01 100644 --- a/Tr/TrackFitEvent/src/FTMeasurement.cpp +++ b/Tr/TrackFitEvent/src/FTMeasurement.cpp @@ -5,7 +5,7 @@ // from Event #include "Event/FTCluster.h" -#include "Kernel/LineTraj.h" + // local #include "Event/FTMeasurement.h" @@ -18,17 +18,13 @@ using namespace LHCb; //----------------------------------------------------------------------------- void FTMeasurement::init( const DeFTDetector& geom ) { - const DeFTFibreMat* ftMat = geom.findFibreMat( m_cluster.channelID() ); + const DeFTMat* ftMat = geom.findMat( m_cluster.channelID() ); m_detectorElement = ftMat; - DetectorSegment mySeg = ftMat->createDetSegment( m_cluster.channelID(), m_cluster.fraction() ); - Gaudi::XYZPoint begPoint( mySeg.x( mySeg.yMin() ), mySeg.yMin(), mySeg.z( mySeg.yMin() ) ); - Gaudi::XYZPoint endPoint( mySeg.x( mySeg.yMax() ), mySeg.yMax(), mySeg.z( mySeg.yMax() ) ); - - m_size = m_cluster.size(); - m_errMeasure = 0.050 + .030 * m_size; - m_trajectory.reset( new LHCb::LineTraj( begPoint, endPoint ) ); - m_z = ftMat->layerCenterZ(); - m_measure = mySeg.x( 0. ); + m_size = m_cluster.isLarge(); + m_errMeasure = 0.080; // + .030 * m_size; + m_trajectory = ftMat->trajectory(m_cluster.channelID(), m_cluster.fraction()); + m_z = ftMat->globalZ(); + m_measure = 0. ; // JvT: I believe this is purely historical. Should remove it? } diff --git a/Tr/TrackFitEvent/xml/FTMeasurement.xml b/Tr/TrackFitEvent/xml/FTMeasurement.xml index 3a2bed3ab7646c0f5406febbf97f68cab9022c41..e1f3a31fb0be8f40130356b642ee553ed5a21a30 100644 --- a/Tr/TrackFitEvent/xml/FTMeasurement.xml +++ b/Tr/TrackFitEvent/xml/FTMeasurement.xml @@ -44,7 +44,7 @@ </code> </method> - <method + <!--method type = 'double' name = 'charge' desc = 'Retrieve the cluster charge' @@ -52,7 +52,7 @@ <code> return m_cluster.charge(); </code> - </method> + </method--> <constructor diff --git a/Tr/TrackMCTools/doc/release.notes b/Tr/TrackMCTools/doc/release.notes index 4e14d1a8882064704099f58b1fa67c42d37fc2cd..d631e91644827b7335d29a0ce4f805710ef5db23 100755 --- a/Tr/TrackMCTools/doc/release.notes +++ b/Tr/TrackMCTools/doc/release.notes @@ -4,6 +4,10 @@ ! Purpose : package for tracking tools accessing MC information !----------------------------------------------------------------------------- +! 2016-12-14 - Jeroen van Tilburg + - Adapt to new FT geometry and numbering scheme + - Changed location of Linker tables. + !========================= TrackMCTools v3r13 2016-03-18 ========================= ! 2016-02-08 - Adam Davis - Add more debugging tools for PatLongLivedTracking diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp b/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp index 4bcfc474c3ce058df99a5b5753135aafe35bcfe9..c8f0f80999675197740e78a2a3dcd79a4b6b546e 100755 --- a/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp +++ b/Tr/TrackMCTools/src/LHCbIDsToMCHits.cpp @@ -10,7 +10,7 @@ #include "Event/VPCluster.h" #include "Event/OTTime.h" #include "Event/MuonCoord.h" -#include "Event/FTCluster.h" +#include "Event/FTLiteCluster.h" DECLARE_TOOL_FACTORY( LHCbIDsToMCHits ) @@ -219,7 +219,7 @@ void LHCbIDsToMCHits::linkFT(const LHCbID& lhcbid, LinkMap& output) const { if (!m_configuredFT){ m_configuredFT = true; - m_ftLinks = FTLinks( evtSvc(), msgSvc(),LHCb::FTClusterLocation::Default+m_endString); + m_ftLinks = FTLinks( evtSvc(), msgSvc(),LHCb::FTLiteClusterLocation::Default+m_endString); if (m_ftLinks.notFound()) { throw GaudiException("no FTLinker", "LHCbIDsToMCHits" , StatusCode::FAILURE); diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCHits.h b/Tr/TrackMCTools/src/LHCbIDsToMCHits.h index 75ce6ead745e648276c48b80f5514f2b2cd41540..506c1bc64da4eb345dad335f5ec2b8fcc24b9cc0 100755 --- a/Tr/TrackMCTools/src/LHCbIDsToMCHits.h +++ b/Tr/TrackMCTools/src/LHCbIDsToMCHits.h @@ -24,8 +24,7 @@ namespace LHCb{ class VeloCluster; class MuonCoord; class VPCluster; - class UTCluster; - class FTCluster; + class UTCluster; } class LHCbIDsToMCHits: public GaudiTool, virtual public ILHCbIDsToMCHits, @@ -90,9 +89,7 @@ private: typedef LinkedTo<LHCb::MCHit,LHCb::VeloCluster> VeloLinks; typedef LinkedTo<LHCb::MCHit,LHCb::VPCluster> VPLinks; typedef LinkedTo<LHCb::MCHit,LHCb::MuonCoord> MuonLinks; - typedef LinkedTo<LHCb::MCHit,LHCb::FTCluster> FTLinks; - - + typedef LinkedTo<LHCb::MCHit> FTLinks; template<typename ID, typename LINKER> void linkToDetTruth(const ID& id, LINKER& aLinker, LinkMap& output ) const; diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp index 8a5a2eb761f10a28d7105a2ca006c7d6efddc730..d65fce775b6076e4f00814e1f33f9440b5f018b6 100755 --- a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp +++ b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.cpp @@ -8,7 +8,7 @@ #include "Event/STCluster.h" #include "Event/VeloCluster.h" #include "Event/VPCluster.h" -#include "Event/FTCluster.h" +#include "Event/FTLiteCluster.h" #include "Event/OTTime.h" #include "Event/MuonCoord.h" @@ -237,7 +237,7 @@ StatusCode LHCbIDsToMCParticles::linkFT(const LHCbID& lhcbid, LinkMap& output) c if (!m_configuredFT){ m_configuredFT = true; m_ftLinks = FTLinks(evtSvc(), msgSvc(), - LHCb::FTClusterLocation::Default); + LHCb::FTLiteClusterLocation::Default); if (m_ftLinks.notFound()) { return Error("no FTLinker",StatusCode::FAILURE,10); } diff --git a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h index 80718e92543e4875ff2915707640634fbd936dd7..8a6397df1fcfa7a54aa99369f2a5357e17ad7717 100755 --- a/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h +++ b/Tr/TrackMCTools/src/LHCbIDsToMCParticles.h @@ -24,7 +24,6 @@ namespace LHCb{ class VeloCluster; class MuonCoord; class VPCluster; - class FTCluster; } class LHCbIDsToMCParticles: public GaudiTool, virtual public ILHCbIDsToMCParticles, @@ -90,7 +89,7 @@ private: typedef LinkedTo<LHCb::MCParticle,LHCb::MuonCoord> MuonLinks; // -- upgrade typedef LinkedTo<LHCb::MCParticle,LHCb::VPCluster> VPLinks; - typedef LinkedTo<LHCb::MCParticle,LHCb::FTCluster> FTLinks; + typedef LinkedTo<LHCb::MCParticle> FTLinks; template<typename ID, typename LINKER> void linkToDetTruth(const ID& id, LINKER& aLinker, LinkMap& output ) const; diff --git a/Tr/TrackTools/doc/release.notes b/Tr/TrackTools/doc/release.notes index 403a8e9a81dd3b4b817ccd368d497a8b9b0e74f1..944ec15f975577d41695b037faab8c235a7ee258 100755 --- a/Tr/TrackTools/doc/release.notes +++ b/Tr/TrackTools/doc/release.notes @@ -4,6 +4,10 @@ ! Purpose : Package for tracking tools not accessing MC information !----------------------------------------------------------------------------- +! 2016-12-14 - Jeroen van Tilburg + - Adapt to new FT geometry and numbering scheme + - Changed location of Linker tables. + ! 2016-06-22 - David Hutchcroft - Added code to set the pT of the VELO only tracks with a sign based on the first strip in the first cluster rather than the track key (which was random). Also an option to reverse the sign diff --git a/Tr/TrackTools/src/FTHitExpectation.cpp b/Tr/TrackTools/src/FTHitExpectation.cpp index 18cd99f1e4fb23d443e571f0129b13500b31607a..671f3aef6f346d3929f04c070adf33ae2359b70c 100644 --- a/Tr/TrackTools/src/FTHitExpectation.cpp +++ b/Tr/TrackTools/src/FTHitExpectation.cpp @@ -5,18 +5,13 @@ #include "Event/Track.h" #include "Event/State.h" -//#include "Kernel/ISTChannelIDSelector.h" - // Tsa #include "TsaKernel/Line.h" #include "TsaKernel/Parabola.h" #include "TsaKernel/Line3D.h" -//#include "TsaKernel/IITExpectedHits.h" #include "LHCbMath/GeomFun.h" #include "FTDet/DeFTDetector.h" -#include "FTDet/DeFTLayer.h" -#include "FTDet/DeFTFibreMat.h" // local #include "FTHitExpectation.h" @@ -36,30 +31,14 @@ DECLARE_TOOL_FACTORY( FTHitExpectation ) //============================================================================= StatusCode FTHitExpectation::initialize() { - StatusCode sc = THitExpectation::initialize(); - if (sc.isFailure()){ - return Error("Failed to initialize", sc); - } + if (sc.isFailure()) return Error("Failed to initialize", sc); - //m_expectedITHits = tool<Tf::Tsa::IITExpectedHits>("Tf::Tsa::ITExpectedHits"); - - // need to know the z of T stations.... + // Get all layers (need to know the z-plane of each layer) m_ftDet = getDet<DeFTDetector>(DeFTDetectorLocation::Default); - if ( m_ftDet -> version() < 20 ){ - return Error("FTDetector version not implemented (different than v2 or v5)",StatusCode::FAILURE); - } - - m_layers = m_ftDet -> layers(); - m_fibremats = m_ftDet -> fibremats(); + if ( m_ftDet -> version() < 61 ) + return Error("This version requires FTDet v6.1 or higher", StatusCode::FAILURE); - /* - // (selector - if (m_selectorName != "ALL"){ - info() << "Selector name " << m_selectorName << endmsg; - m_selector = tool< ISTChannelIDSelector >( m_selectorType,m_selectorName ); - } - */ return StatusCode::SUCCESS; } @@ -69,135 +48,100 @@ IHitExpectation::Info FTHitExpectation::expectation(const LHCb::Track& aTrack) c IHitExpectation::Info info; info.likelihood = 0.0; info.nFound = 0; - info.nExpected = 0;; - + info.nExpected = 0; const auto& ids = aTrack.lhcbIDs(); std::vector<LHCb::LHCbID> ftHits; ftHits.reserve(ids.size()); std::copy_if( ids.begin(), ids.end(), std::back_inserter(ftHits), [](const LHCb::LHCbID& id) { return id.isFT(); } ); - for ( auto iterL = m_layers.begin(); iterL != m_layers.end() ; ++iterL ) - { - LHCb::FTChannelID elemID( (unsigned int)std::distance( m_layers.begin(), iterL ), 0u, 0u, 0u, 0u ); - double layerZ = ((*iterL)->geometry()->toGlobal( Gaudi::XYZPoint(0. ,0., 0.) )).z(); - - Tf::Tsa::Line line = yLine ( aTrack, layerZ ); - Tf::Tsa::Parabola aParab = xParabola( aTrack, layerZ ); - Tf::Tsa::Line tanLine = aParab.tangent( layerZ ); - Tf::Tsa::Line3D aLine3D = Tf::Tsa::createLine3D( tanLine, line, layerZ ); - std::vector<std::pair<LHCb::FTChannelID, double> > vectFTPairs; - - collectHits( aLine3D, elemID, vectFTPairs ); - - int old = -1; - for ( const auto& p : vectFTPairs ) - { - if ( int(p.second) != old) { - old = int(p.second); - if ( p.first ) { - ++(info.nExpected); - auto itIter = std::find( ftHits.begin(), ftHits.end(), p.first); - if (itIter != ftHits.end() ) ++info.nFound; - } - } - } // pairs - }//layers - return info; -} + // Loop over all layers + for( auto station : m_ftDet->stations() ) { + for( auto layer : station->layers() ) { + + const DeFTMat* mat; + Gaudi::XYZPoint intersectionPoint; + bool expectHit = findIntersectingMat(layer, aTrack, mat, intersectionPoint); + + if( expectHit ) { + ++(info.nExpected); -void FTHitExpectation::collectHits( Tf::Tsa::Line3D aLine3D, - LHCb::FTChannelID chan, - std::vector<std::pair<LHCb::FTChannelID, double> >& vectFTPairs - ) const{ - - for ( const auto& f : m_fibremats) - { - if( (unsigned int)(f -> layer()) != chan.layer() ) continue; - // propagate to z of sector - Gaudi::Plane3D entryplane, exitplane; - EntryExitPlanes( *f, entryplane, exitplane ); - - Gaudi::XYZPoint globalEntry = intersection( aLine3D, entryplane ); - Gaudi::XYZPoint globalExit = intersection( aLine3D, exitplane ); - - LHCb::MCHit hit; - hit.setEntry( globalEntry ); - hit.setDisplacement ( globalExit-globalEntry ); - hit.setEnergy( 1.0 ); - if( f -> isInside( hit.midPoint() ) ){ - f -> calculateHits( hit, vectFTPairs ).ignore(); + // Check of this layer is present in the hits + auto itIter = std::find_if( ftHits.begin(), ftHits.end(), [&mat] (const LHCb::LHCbID& id) + {return id.ftID().uniqueLayer() == mat->elementID().uniqueLayer(); } ); + if (itIter != ftHits.end() ) ++info.nFound; } } -} - + } -void FTHitExpectation::EntryExitPlanes( DeFTFibreMat& mat, Gaudi::Plane3D& entry, Gaudi::Plane3D& exit ) const -{ - double xmin = mat.fibreMatMinX(); - double ymin = mat.fibreMatMinY(); - double xmax = mat.fibreMatMaxX(); - double ymax = mat.fibreMatMaxY(); - - double xLower = -0.5*(xmax-xmin); - double xUpper = 0.5*(xmax-xmin); - double yLower = -0.5*(ymax-ymin); - double yUpper = 0.5*(ymax-ymin); - - const Gaudi::XYZPoint g1 = mat.geometry() -> toGlobal( Gaudi::XYZPoint( xLower, yLower, 0. ) ); - const Gaudi::XYZPoint g2 = mat.geometry() -> toGlobal( Gaudi::XYZPoint( xLower, yUpper, 0. ) ); - // trajectory of middle - //const Gaudi::XYZPoint g3 = mat.geometry()->toGlobal(xLower, 0., 0.); - const Gaudi::XYZPoint g4 = mat.geometry()->toGlobal(Gaudi::XYZPoint( xUpper, 0., 0.) ); - // plane - Gaudi::Plane3D plane(g1,g2,g4); - - double thickness = 1.2; - entry = Gaudi::Plane3D(plane.Normal(), mat.geometry()->toGlobal( Gaudi::XYZPoint(0.,0.,-0.5*thickness))); - exit = Gaudi::Plane3D(plane.Normal(), mat.geometry()->toGlobal( Gaudi::XYZPoint(0.,0., 0.5*thickness))); + return info; } -void FTHitExpectation::collect(const LHCb::Track& aTrack ,std::vector<LHCb::LHCbID>& ids) const{ +bool FTHitExpectation::findIntersectingMat(const DeFTLayer* layer, const LHCb::Track& aTrack, + const DeFTMat*& mat, Gaudi::XYZPoint& intersectionPoint) const { + + // make a Line3D from the track + double layerZ = layer->globalZ(); + Tf::Tsa::Line line = yLine ( aTrack, layerZ ); + Tf::Tsa::Parabola aParab = xParabola( aTrack, layerZ ); + Tf::Tsa::Line tanLine = aParab.tangent( layerZ ); + Tf::Tsa::Line3D aLine3D = Tf::Tsa::createLine3D( tanLine, line, layerZ ); - for ( auto iterL = m_layers.begin(); iterL != m_layers.end() ; ++iterL ) - { - LHCb::FTChannelID elemID( (unsigned int)std::distance( m_layers.begin(), iterL ), 0u, 0u, 0u, 0u ); - double layerZ = ((*iterL)->geometry()->toGlobal( Gaudi::XYZPoint(0. ,0., 0.) )).z(); + // find intersection point of track and plane of layer + const Gaudi::Plane3D layerPlane = layer->plane(); + intersectionPoint = intersection(aLine3D, layerPlane); - Tf::Tsa::Line line = yLine ( aTrack, layerZ ); - Tf::Tsa::Parabola aParab = xParabola( aTrack, layerZ ); - Tf::Tsa::Line tanLine = aParab.tangent( layerZ ); - Tf::Tsa::Line3D aLine3D = Tf::Tsa::createLine3D( tanLine, line, layerZ ); - std::vector<std::pair<LHCb::FTChannelID, double> > vectFTPairs; + // find the module that is hit + const DeFTModule* module = layer->findModule( intersectionPoint ); + if( module == nullptr ) return false; - collectHits( aLine3D, elemID, vectFTPairs ); + // find intersection point of track and plane of module + // (to account for misalignments between module and layer) + tanLine = aParab.tangent( intersectionPoint.z() ); + aLine3D = Tf::Tsa::createLine3D( tanLine, line, intersectionPoint.z() ); + const Gaudi::Plane3D modulePlane = module->plane(); + intersectionPoint = intersection(aLine3D, modulePlane); + // Find the corresponding mat + mat = module->findMat(intersectionPoint); - int old = -1; - ids.reserve(vectFTPairs.size()); - for ( std::vector<std::pair<LHCb::FTChannelID, double> >::iterator iter = vectFTPairs.begin(); - iter != vectFTPairs.end(); - ++iter ) - { - if ( int(iter->second) != old ) { - old = int(iter->second); - if ( iter->first) ids.push_back( LHCb::LHCbID(iter->first) ); - } - } // pairs - }//layers + // check if intersection point is inside the fibres of the module + return (mat == nullptr ) ? false : mat->isInside( intersectionPoint ); } +void FTHitExpectation::collect(const LHCb::Track& aTrack, + std::vector<LHCb::LHCbID>& ids) const{ + // Loop over all layers + for( auto station : m_ftDet->stations() ) { + for( auto layer : station->layers() ) { + + const DeFTMat* mat; + Gaudi::XYZPoint intersectionPoint; + bool expectHit = findIntersectingMat(layer, aTrack, mat, intersectionPoint); + + if( expectHit ) { + // Find the channel that is closest + Gaudi::XYZPoint localP = mat->geometry()->toLocal( intersectionPoint ); + double frac; + LHCb::FTChannelID ftChan = mat->calculateChannelAndFrac(localP.x(), frac); + + // Add the channels + // JvT: Without the fraction the bare FTChannelID is pretty useless... + if( std::abs(frac) <= 0.5 ) ids.push_back( LHCb::LHCbID(ftChan) ); + } + } + } +} unsigned int FTHitExpectation::nExpected(const LHCb::Track& aTrack) const{ - return expectation(aTrack).nExpected; } + Gaudi::XYZPoint FTHitExpectation::intersection(const Tf::Tsa::Line3D& line, const Gaudi::Plane3D& aPlane) const{ - // make a plane Gaudi::XYZPoint inter; double mu = 0; Gaudi::Math::intersection( line, aPlane, inter, mu ); diff --git a/Tr/TrackTools/src/FTHitExpectation.h b/Tr/TrackTools/src/FTHitExpectation.h index 010f9fcad1e23f050d12961d2ff3c972f2b09bf9..2fe8ba3de7034c709f48f3f4344a0fb0e7d42521 100644 --- a/Tr/TrackTools/src/FTHitExpectation.h +++ b/Tr/TrackTools/src/FTHitExpectation.h @@ -47,7 +47,6 @@ public: */ unsigned int nExpected ( const LHCb::Track& aTrack ) const override; - /** Returns number of hits expected * * @param aTrack Reference to the Track to test @@ -67,27 +66,13 @@ public: private: + bool findIntersectingMat(const DeFTLayer* layer, const LHCb::Track& aTrack, + const DeFTMat*& mat, Gaudi::XYZPoint& intersectionPoint) const ; - //bool select(const LHCb::FTChannelID& chan) const; - - void EntryExitPlanes( DeFTFibreMat& mat, Gaudi::Plane3D& entry, Gaudi::Plane3D& exit ) const; - void collectHits( Tf::Tsa::Line3D aLine3D, - LHCb::FTChannelID chan, - std::vector<std::pair<LHCb::FTChannelID, double> >& vectFTPairs - ) const; Gaudi::XYZPoint intersection(const Tf::Tsa::Line3D& line, const Gaudi::Plane3D& aPlane) const; - DeFTDetector* m_ftDet = nullptr; - std::vector<DeFTLayer*> m_layers; - std::vector<DeFTFibreMat*> m_fibremats; - + }; -/* -inline bool ITHitExpectation::select(const LHCb::STChannelID& chan) const{ - return !m_selector || m_selector->select(chan); -} -*/ - #endif diff --git a/Tr/TrackTools/src/FTMeasurementProvider.cpp b/Tr/TrackTools/src/FTMeasurementProvider.cpp index f3beec1cd9c7a05f9d6d9ab07be503a8575afc67..12978959fa04794a10a0cb2de9ec338c246a6a28 100644 --- a/Tr/TrackTools/src/FTMeasurementProvider.cpp +++ b/Tr/TrackTools/src/FTMeasurementProvider.cpp @@ -31,8 +31,11 @@ public: using base_class::base_class; StatusCode initialize() override; - LHCb::Measurement* measurement( const LHCb::LHCbID& id, bool localY=false ) const override; - LHCb::Measurement* measurement( const LHCb::LHCbID& id, const LHCb::ZTrajectory& refvector, bool localY=false) const override; + LHCb::Measurement* measurement( const LHCb::LHCbID& id, + bool localY=false ) const override; + LHCb::Measurement* measurement( const LHCb::LHCbID& id, + const LHCb::ZTrajectory& refvector, + bool localY=false) const override; inline LHCb::FTMeasurement* ftmeasurement( const LHCb::LHCbID& id ) const ; void handle ( const Incident& incident ) override; @@ -91,13 +94,10 @@ const FastClusterContainer<LHCb::FTLiteCluster,int>* FTMeasurementProvider::clus { /// If there is a new event, get the clusters and sort them according to their channel ID if( UNLIKELY(!m_clusters) ){ - m_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int> >( LHCb::FTLiteClusterLocation::Default ); - if(m_clusters){ //TODO: sort before putting clusters on the TES! - std::sort(m_clusters->begin(), m_clusters->end(), - [](const LHCb::FTLiteCluster& lhs, const LHCb::FTLiteCluster& rhs){ return lhs.channelID() < rhs.channelID(); }); - }else{ - error() << "Could not find FTLiteClusters at: " << LHCb::FTLiteClusterLocation::Default << endmsg; - } + m_clusters = getIfExists<FastClusterContainer<LHCb::FTLiteCluster,int> > + ( LHCb::FTLiteClusterLocation::Default ); + if(!m_clusters) error() << "Could not find FTLiteClusters at: " + << LHCb::FTLiteClusterLocation::Default << endmsg; } return m_clusters ; } @@ -111,8 +111,8 @@ LHCb::FTMeasurement* FTMeasurementProvider::ftmeasurement( const LHCb::LHCbID& i error() << "Not an FT measurement" << endmsg ; return nullptr; } - /// The clusters are sorted, so we can use a binary search (lower bound search) to find the element - /// corresponding to the channel ID + /// The clusters are sorted, so we can use a binary search (lower bound search) + /// to find the element corresponding to the channel ID const auto& c = *clusters(); auto itH = std::lower_bound( c.begin(), c.end(), id.ftID(), [](const LHCb::FTLiteCluster clus, const LHCb::FTChannelID id){ @@ -124,7 +124,8 @@ LHCb::FTMeasurement* FTMeasurementProvider::ftmeasurement( const LHCb::LHCbID& i //----------------------------------------------------------------------------- /// Return the measurement //----------------------------------------------------------------------------- -LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, bool /*localY*/ ) const { +LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, + bool /*localY*/ ) const { return ftmeasurement(id) ; } @@ -133,7 +134,8 @@ LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, b //----------------------------------------------------------------------------- LHCb::Measurement* FTMeasurementProvider::measurement( const LHCb::LHCbID& id, - const LHCb::ZTrajectory& /* reftraj */, bool /*localY*/ ) const { + const LHCb::ZTrajectory& /* reftraj */, + bool /*localY*/ ) const { // default implementation return ftmeasurement(id) ; }