diff --git a/Configuration/ProcessModifiers/python/stripNtupletFit_cff.py b/Configuration/ProcessModifiers/python/stripNtupletFit_cff.py new file mode 100644 index 0000000000000..9089f3574b987 --- /dev/null +++ b/Configuration/ProcessModifiers/python/stripNtupletFit_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier can be used in combination with pixelNtupletFit to include strip hits in the tracks + +stripNtupletFit = cms.Modifier() diff --git a/Configuration/PyReleaseValidation/python/relval_gpu.py b/Configuration/PyReleaseValidation/python/relval_gpu.py index 71e474303bf87..13c3488630aad 100644 --- a/Configuration/PyReleaseValidation/python/relval_gpu.py +++ b/Configuration/PyReleaseValidation/python/relval_gpu.py @@ -54,6 +54,7 @@ 12850.402, 12850.403, 12850.404, 12450.406, 12450.407, 12450.408, 12861.402, + 16834.409, # 2024 with PU, Alpaka-based 13034.402, 13034.403, 13034.404, @@ -66,6 +67,7 @@ 13050.402, 13050.403, 13050.404, 13050.406, 13050.407, 13050.408, 13061.402, + 17034.0,17034.409, # Run4, Alpaka-based noPU 29634.402, 29634.403, 29634.404, 29634.406, 29634.704, diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 8275172ecfbda..1fc588ebc8a2f 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -1402,6 +1402,20 @@ def setup_(self, step, stepName, stepDict, k, properties): offset = 0.408, ) +upgradeWFs['PatatrackPixelPlusStripAlpaka'] = PatatrackWorkflow( + digi = { + '--procModifiers': 'alpaka', + '--customise': 'HLTrigger/Configuration/customizeHLTforAlpakaStripNoDoubletRecovery.customizeHLTforAlpakaStripNoDoubletRecovery' + }, + reco = { + '-s': 'RAW2DIGI,RECO', + '--procModifiers': 'alpaka,stripNtupletFit' + }, + harvest = None, + suffix = 'Patatrack_PixelPlusStripAlpaka', + offset = 0.409, +) + # end of Patatrack workflows ############################################################################################################### diff --git a/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareRecHitsSoA.cc b/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareRecHitsSoA.cc index 6e2a908b59b38..3b8170f1d7198 100644 --- a/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareRecHitsSoA.cc +++ b/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareRecHitsSoA.cc @@ -27,7 +27,7 @@ #include "Geometry/CommonTopologies/interface/PixelTopology.h" #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" - +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" template class SiPixelCompareRecHitsSoA : public DQMEDAnalyzer { public: @@ -246,9 +246,12 @@ void SiPixelCompareRecHitsSoA::fillDescriptions(edm::ConfigurationDescription } using SiPixelPhase1CompareRecHitsSoA = SiPixelCompareRecHitsSoA; +using SiPixelPhase1StripCompareRecHitsSoA = SiPixelCompareRecHitsSoA; using SiPixelPhase2CompareRecHitsSoA = SiPixelCompareRecHitsSoA; using SiPixelHIonPhase1CompareRecHitsSoA = SiPixelCompareRecHitsSoA; +#include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(SiPixelPhase1CompareRecHitsSoA); +DEFINE_FWK_MODULE(SiPixelPhase1StripCompareRecHitsSoA); DEFINE_FWK_MODULE(SiPixelPhase2CompareRecHitsSoA); DEFINE_FWK_MODULE(SiPixelHIonPhase1CompareRecHitsSoA); diff --git a/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareTrackSoA.cc b/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareTrackSoA.cc index f3635d6df45da..456645df65576 100644 --- a/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareTrackSoA.cc +++ b/DQM/SiPixelHeterogeneous/plugins/SiPixelCompareTrackSoA.cc @@ -22,6 +22,7 @@ #include "DQMServices/Core/interface/DQMStore.h" #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousHost.h" #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousDevice.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" // for string manipulations #include @@ -362,7 +363,9 @@ void SiPixelCompareTrackSoA::fillDescriptions(edm::ConfigurationDescriptions& using SiPixelPhase1CompareTrackSoA = SiPixelCompareTrackSoA; using SiPixelPhase2CompareTrackSoA = SiPixelCompareTrackSoA; using SiPixelHIonPhase1CompareTrackSoA = SiPixelCompareTrackSoA; +using SiPixelPhase1StripCompareTrackSoA = SiPixelCompareTrackSoA; DEFINE_FWK_MODULE(SiPixelPhase1CompareTrackSoA); +DEFINE_FWK_MODULE(SiPixelPhase1StripCompareTrackSoA); DEFINE_FWK_MODULE(SiPixelPhase2CompareTrackSoA); DEFINE_FWK_MODULE(SiPixelHIonPhase1CompareTrackSoA); diff --git a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorRecHitsSoAAlpaka.cc b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorRecHitsSoAAlpaka.cc index f4c8968fafb16..f206d6cd9bb83 100644 --- a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorRecHitsSoAAlpaka.cc +++ b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorRecHitsSoAAlpaka.cc @@ -189,10 +189,12 @@ void SiPixelMonitorRecHitsSoAAlpaka::fillDescriptions(edm::ConfigurationDescr } using SiPixelPhase1MonitorRecHitsSoAAlpaka = SiPixelMonitorRecHitsSoAAlpaka; +using SiPixelPhase1StripMonitorRecHitsSoAAlpaka = SiPixelMonitorRecHitsSoAAlpaka; using SiPixelPhase2MonitorRecHitsSoAAlpaka = SiPixelMonitorRecHitsSoAAlpaka; using SiPixelHIonPhase1MonitorRecHitsSoAAlpaka = SiPixelMonitorRecHitsSoAAlpaka; #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(SiPixelPhase1MonitorRecHitsSoAAlpaka); +DEFINE_FWK_MODULE(SiPixelPhase1StripMonitorRecHitsSoAAlpaka); DEFINE_FWK_MODULE(SiPixelPhase2MonitorRecHitsSoAAlpaka); DEFINE_FWK_MODULE(SiPixelHIonPhase1MonitorRecHitsSoAAlpaka); diff --git a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoA.cc b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoA.cc index f3ccb74bc3fea..df68a82ecf905 100644 --- a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoA.cc +++ b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoA.cc @@ -22,6 +22,7 @@ #include "DQMServices/Core/interface/DQMStore.h" #include "CUDADataFormats/Track/interface/PixelTrackUtilities.h" #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousHost.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" // for string manipulations #include @@ -195,7 +196,9 @@ void SiPixelMonitorTrackSoA::fillDescriptions(edm::ConfigurationDescriptions& using SiPixelPhase1MonitorTrackSoA = SiPixelMonitorTrackSoA; using SiPixelPhase2MonitorTrackSoA = SiPixelMonitorTrackSoA; using SiPixelHIonPhase1MonitorTrackSoA = SiPixelMonitorTrackSoA; +using SiPixelPhase1StripMonitorTrackSoA = SiPixelMonitorTrackSoA; DEFINE_FWK_MODULE(SiPixelPhase1MonitorTrackSoA); DEFINE_FWK_MODULE(SiPixelPhase2MonitorTrackSoA); DEFINE_FWK_MODULE(SiPixelHIonPhase1MonitorTrackSoA); +DEFINE_FWK_MODULE(SiPixelPhase1StripMonitorTrackSoA); diff --git a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoAAlpaka.cc b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoAAlpaka.cc index 8bd1cdfa2e429..478f11ce82906 100644 --- a/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoAAlpaka.cc +++ b/DQM/SiPixelHeterogeneous/plugins/SiPixelMonitorTrackSoAAlpaka.cc @@ -25,6 +25,8 @@ #include "DQMServices/Core/interface/DQMStore.h" #include "DataFormats/TrackSoA/interface/TracksHost.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" + template class SiPixelMonitorTrackSoAAlpaka : public DQMEDAnalyzer { public: @@ -195,7 +197,9 @@ void SiPixelMonitorTrackSoAAlpaka::fillDescriptions(edm::ConfigurationDescrip using SiPixelPhase1MonitorTrackSoAAlpaka = SiPixelMonitorTrackSoAAlpaka; using SiPixelPhase2MonitorTrackSoAAlpaka = SiPixelMonitorTrackSoAAlpaka; using SiPixelHIonPhase1MonitorTrackSoAAlpaka = SiPixelMonitorTrackSoAAlpaka; +using SiPixelPhase1StripMonitorTrackSoAAlpaka = SiPixelMonitorTrackSoAAlpaka; DEFINE_FWK_MODULE(SiPixelPhase1MonitorTrackSoAAlpaka); DEFINE_FWK_MODULE(SiPixelPhase2MonitorTrackSoAAlpaka); DEFINE_FWK_MODULE(SiPixelHIonPhase1MonitorTrackSoAAlpaka); +DEFINE_FWK_MODULE(SiPixelPhase1StripMonitorTrackSoAAlpaka); diff --git a/DQM/SiPixelHeterogeneous/python/SiPixelHeterogenousDQM_FirstStep_cff.py b/DQM/SiPixelHeterogeneous/python/SiPixelHeterogenousDQM_FirstStep_cff.py index a08b3df2eea33..b2da30f4ab5df 100644 --- a/DQM/SiPixelHeterogeneous/python/SiPixelHeterogenousDQM_FirstStep_cff.py +++ b/DQM/SiPixelHeterogeneous/python/SiPixelHeterogenousDQM_FirstStep_cff.py @@ -7,6 +7,7 @@ from DQM.SiPixelHeterogeneous.siPixelPhase2MonitorTrackSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelHIonPhase1MonitorTrackSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelMonitorVertexSoA_cfi import * +from DQM.SiPixelHeterogeneous.siPixelPhase1StripMonitorTrackSoA_cfi import * # Alpaka Modules from Configuration.ProcessModifiers.alpaka_cff import alpaka from DQM.SiPixelHeterogeneous.siPixelPhase1MonitorRecHitsSoAAlpaka_cfi import * @@ -17,6 +18,8 @@ from DQM.SiPixelHeterogeneous.siPixelHIonPhase1MonitorTrackSoAAlpaka_cfi import * from DQM.SiPixelHeterogeneous.siPixelMonitorVertexSoAAlpaka_cfi import * + + # Run-3 sequence monitorpixelSoASource = cms.Sequence(siPixelPhase1MonitorRecHitsSoA * siPixelPhase1MonitorTrackSoA * siPixelMonitorVertexSoA) # Run-3 Alpaka sequence @@ -45,6 +48,7 @@ from DQM.SiPixelHeterogeneous.siPixelPhase1CompareTrackSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelPhase2CompareTrackSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelHIonPhase1CompareTrackSoA_cfi import * +from DQM.SiPixelHeterogeneous.siPixelPhase1StripCompareTrackSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelCompareVertexSoA_cfi import * from DQM.SiPixelHeterogeneous.siPixelPhase1RawDataErrorComparator_cfi import * from DQM.SiPixelPhase1Common.SiPixelPhase1RawData_cfi import * @@ -211,6 +215,33 @@ topFolderName = cms.string('SiPixelHeterogeneous/PixelVertexDevice') ) +monitorPixelTracksAlpaka = cms.Sequence(siPixelTrackSoAMonitorSerial * + siPixelTrackSoAMonitorDevice * + siPixelPhase1CompareTrackSoA) + +# PixelTracks: monitor of CPUSerial product (Alpaka backend: 'serial_sync') +siPixelTrackSoAMonitorSerialStrip = siPixelPhase1StripMonitorTrackSoA.clone( + pixelTrackSrc = cms.InputTag('pixelTracksAlpakaSerial'), + topFolderName = cms.string('SiPixelHeterogeneous/PixelTrackSerial') +) + +# PixelTracks: monitor of CPUSerial product (Alpaka backend: 'serial_sync') +siPixelTrackSoAMonitorDeviceStrip = siPixelPhase1StripMonitorTrackSoA.clone( + pixelTrackSrc = cms.InputTag('pixelTracksAlpaka'), + topFolderName = cms.string('SiPixelHeterogeneous/PixelTrackDevice') +) + +monitorPixelTracksAlpakaStrip = cms.Sequence( siPixelTrackSoAMonitorSerialStrip * + siPixelTrackSoAMonitorDeviceStrip * + siPixelPhase1StripCompareTrackSoA) + + +from Configuration.ProcessModifiers.stripNtupletFit_cff import stripNtupletFit +stripNtupletFit.toReplaceWith(monitorPixelTracksAlpaka, monitorPixelTracksAlpakaStrip) +stripNtupletFit.toReplaceWith(siPixelPhase1CompareTrackSoA, siPixelPhase1StripCompareTrackSoA) +stripNtupletFit.toReplaceWith(siPixelTrackSoAMonitorSerial, siPixelTrackSoAMonitorSerial) +stripNtupletFit.toReplaceWith(siPixelTrackSoAMonitorDevice, siPixelTrackSoAMonitorDevice) + # Run-3 sequence monitorpixelSoACompareSource = cms.Sequence(siPixelPhase1MonitorRawDataACPU * siPixelPhase1MonitorRawDataAGPU * diff --git a/DataFormats/TrackSoA/interface/TracksHost.h b/DataFormats/TrackSoA/interface/TracksHost.h index 69b500b9672cd..37087b54ff49c 100644 --- a/DataFormats/TrackSoA/interface/TracksHost.h +++ b/DataFormats/TrackSoA/interface/TracksHost.h @@ -4,12 +4,12 @@ #include #include - #include "DataFormats/Common/interface/Uninitialized.h" #include "DataFormats/Portable/interface/PortableHostCollection.h" #include "DataFormats/TrackSoA/interface/TrackDefinitions.h" #include "DataFormats/TrackSoA/interface/TracksSoA.h" #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" // TODO: The class is created via inheritance of the PortableHostCollection. // This is generally discouraged, and should be done via composition. @@ -42,6 +42,7 @@ namespace pixelTrack { using TracksHostPhase1 = TracksHost; using TracksHostPhase2 = TracksHost; using TracksHostHIonPhase1 = TracksHost; + using TracksHostPhase1Strip = TracksHost; } // namespace pixelTrack diff --git a/DataFormats/TrackSoA/interface/TracksSoA.h b/DataFormats/TrackSoA/interface/TracksSoA.h index ed4ef2e5a4c93..e6cc332434f5b 100644 --- a/DataFormats/TrackSoA/interface/TracksSoA.h +++ b/DataFormats/TrackSoA/interface/TracksSoA.h @@ -7,6 +7,7 @@ #include "HeterogeneousCore/AlpakaInterface/interface/OneToManyAssoc.h" #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "DataFormats/SoATemplate/interface/SoALayout.h" #include "DataFormats/TrackSoA/interface/TrackDefinitions.h" @@ -66,33 +67,39 @@ namespace reco { struct IsTrackSoAConstView> : std::true_type {}; template <> struct IsTrackSoAConstView> : std::true_type {}; + template <> + struct IsTrackSoAConstView> : std::true_type {}; + template <> + struct IsTrackSoAConstView> : std::true_type {}; template constexpr bool isTrackSoAConstView = IsTrackSoAConstView::value; - template >> + // enable_if should be used when there is another implementation, + // please use static_assert to report invalid template arguments + template //, typename = std::enable_if_t>> ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float charge(ConstView const& tracks, int32_t i) { //was: std::copysign(1.f, tracks[i].state()(2)). Will be constexpr with C++23 float v = tracks[i].state()(2); return float((0.0f < v) - (v < 0.0f)); } - template >> + template //, typename = std::enable_if_t>> ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float phi(ConstView const& tracks, int32_t i) { return tracks[i].state()(0); } - template >> + template //, typename = std::enable_if_t>> ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float tip(ConstView const& tracks, int32_t i) { return tracks[i].state()(1); } - template >> + template //, typename = std::enable_if_t>> ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float zip(ConstView const& tracks, int32_t i) { return tracks[i].state()(4); } - template >> + template //, typename = std::enable_if_t>> ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr bool isTriplet(ConstView const& tracks, int32_t i) { return tracks[i].nLayers() == 3; } diff --git a/DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h b/DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h index f50756f3ddbca..610af5933608b 100644 --- a/DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h +++ b/DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h @@ -140,6 +140,7 @@ namespace pixelTrack { float pt = std::min(tracks.pt(it), chi2MaxPt); float chi2Cut = chi2Scale * (chi2Coeff[0] + roughLog(pt) * chi2Coeff[1]); + //chi2Cut = chi2Scale; if (tracks.chi2(it) >= chi2Cut) { #ifdef NTUPLE_FIT_DEBUG printf("Bad chi2 %d pt %f eta %f chi2 %f\n", it, tracks.pt(it), tracks.eta(it), tracks.chi2(it)); @@ -173,6 +174,7 @@ namespace pixelTrack { // TODO: Should those be placed in the ALPAKA_ACCELERATOR_NAMESPACE template struct TracksUtilities; +template struct TracksUtilities; template struct TracksUtilities; #endif // DataFormats_TrackSoA_interface_alpaka_TrackUtilities_h diff --git a/DataFormats/TrackSoA/interface/alpaka/TracksSoACollection.h b/DataFormats/TrackSoA/interface/alpaka/TracksSoACollection.h index c9294d693d4c4..08690a9e8bbb1 100644 --- a/DataFormats/TrackSoA/interface/alpaka/TracksSoACollection.h +++ b/DataFormats/TrackSoA/interface/alpaka/TracksSoACollection.h @@ -8,7 +8,7 @@ #include "DataFormats/Portable/interface/alpaka/PortableCollection.h" #include "DataFormats/TrackSoA/interface/TracksDevice.h" #include "DataFormats/TrackSoA/interface/TracksHost.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/AssertDeviceMatchesHostCollection.h" #include "HeterogeneousCore/AlpakaInterface/interface/CopyToHost.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" @@ -29,6 +29,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using TracksSoACollectionPhase1 = TracksSoACollection; using TracksSoACollectionPhase2 = TracksSoACollection; using TracksSoACollectionHIonPhase1 = TracksSoACollection; + using TracksSoACollectionPhase1Strip = TracksSoACollection; } // namespace pixelTrack } // namespace ALPAKA_ACCELERATOR_NAMESPACE @@ -50,5 +51,6 @@ namespace cms::alpakatools { ASSERT_DEVICE_MATCHES_HOST_COLLECTION(pixelTrack::TracksSoACollectionPhase1, pixelTrack::TracksHostPhase1); ASSERT_DEVICE_MATCHES_HOST_COLLECTION(pixelTrack::TracksSoACollectionPhase2, pixelTrack::TracksHostPhase2); ASSERT_DEVICE_MATCHES_HOST_COLLECTION(pixelTrack::TracksSoACollectionHIonPhase1, pixelTrack::TracksHostHIonPhase1); +ASSERT_DEVICE_MATCHES_HOST_COLLECTION(pixelTrack::TracksSoACollectionPhase1Strip, pixelTrack::TracksHostPhase1Strip); #endif // DataFormats_TrackSoA_interface_alpaka_TracksSoACollection_h diff --git a/DataFormats/TrackSoA/src/alpaka/classes_cuda_def.xml b/DataFormats/TrackSoA/src/alpaka/classes_cuda_def.xml index 08d7eb9724b34..d5fdc4f378e21 100644 --- a/DataFormats/TrackSoA/src/alpaka/classes_cuda_def.xml +++ b/DataFormats/TrackSoA/src/alpaka/classes_cuda_def.xml @@ -13,4 +13,9 @@ + + + + + diff --git a/DataFormats/TrackSoA/src/alpaka/classes_rocm_def.xml b/DataFormats/TrackSoA/src/alpaka/classes_rocm_def.xml index 0cd8cf9bb49d7..fd0fed505783c 100644 --- a/DataFormats/TrackSoA/src/alpaka/classes_rocm_def.xml +++ b/DataFormats/TrackSoA/src/alpaka/classes_rocm_def.xml @@ -13,4 +13,9 @@ + + + + + diff --git a/DataFormats/TrackSoA/src/classes.cc b/DataFormats/TrackSoA/src/classes.cc index 97e00cc5b5638..3aaeaf00d4b93 100644 --- a/DataFormats/TrackSoA/src/classes.cc +++ b/DataFormats/TrackSoA/src/classes.cc @@ -1,9 +1,10 @@ #include "DataFormats/Portable/interface/PortableHostCollectionReadRules.h" #include "DataFormats/TrackSoA/interface/TracksSoA.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" using namespace reco; SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); // SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); //TODO: For the moment we live without HIons +SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); diff --git a/DataFormats/TrackSoA/src/classes_def.xml b/DataFormats/TrackSoA/src/classes_def.xml index dcbe554334b93..f43436e85137a 100644 --- a/DataFormats/TrackSoA/src/classes_def.xml +++ b/DataFormats/TrackSoA/src/classes_def.xml @@ -25,4 +25,14 @@ + + + + + + + + + + diff --git a/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsDevice.h b/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsDevice.h index 885aba8f106a5..83de285ba9fb6 100644 --- a/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsDevice.h +++ b/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsDevice.h @@ -9,6 +9,7 @@ #include "DataFormats/Portable/interface/PortableDeviceCollection.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" template class TrackingRecHitDevice : public PortableDeviceCollection, TDev> { diff --git a/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h b/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h index 1480236f9517b..32499a8b8b97a 100644 --- a/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h +++ b/DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h @@ -9,6 +9,7 @@ #include "DataFormats/Portable/interface/PortableHostCollection.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" template class TrackingRecHitHost : public PortableHostCollection> { @@ -50,5 +51,6 @@ class TrackingRecHitHost : public PortableHostCollection; using TrackingRecHitHostPhase2 = TrackingRecHitHost; using TrackingRecHitHostHIonPhase1 = TrackingRecHitHost; +using TrackingRecHitHostPhase1Strip = TrackingRecHitHost; #endif // DataFormats_TrackingRecHitSoA_interface_TrackingRecHitsHost_h diff --git a/DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h b/DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h index 7c3fef745c669..cc31ddee5425c 100644 --- a/DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h +++ b/DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h @@ -10,7 +10,7 @@ #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsDevice.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsHost.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/CopyToHost.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" @@ -25,6 +25,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using TrackingRecHitSoAPhase1 = TrackingRecHitsSoACollection; using TrackingRecHitSoAPhase2 = TrackingRecHitsSoACollection; using TrackingRecHitSoAHIonPhase1 = TrackingRecHitsSoACollection; + using TrackingRecHitSoAPhase1Strip = TrackingRecHitsSoACollection; } // namespace ALPAKA_ACCELERATOR_NAMESPACE @@ -74,5 +75,5 @@ namespace cms::alpakatools { ASSERT_DEVICE_MATCHES_HOST_COLLECTION(TrackingRecHitSoAPhase1, TrackingRecHitHostPhase1); ASSERT_DEVICE_MATCHES_HOST_COLLECTION(TrackingRecHitSoAPhase2, TrackingRecHitHostPhase2); ASSERT_DEVICE_MATCHES_HOST_COLLECTION(TrackingRecHitSoAHIonPhase1, TrackingRecHitHostHIonPhase1); - +ASSERT_DEVICE_MATCHES_HOST_COLLECTION(TrackingRecHitSoAPhase1Strip, TrackingRecHitHostPhase1Strip); #endif // DataFormats_TrackingRecHitSoA_interface_alpaka_TrackingRecHitsSoACollection_h diff --git a/DataFormats/TrackingRecHitSoA/src/alpaka/classes_cuda_def.xml b/DataFormats/TrackingRecHitSoA/src/alpaka/classes_cuda_def.xml index 80c267b57d585..7f88c04299571 100644 --- a/DataFormats/TrackingRecHitSoA/src/alpaka/classes_cuda_def.xml +++ b/DataFormats/TrackingRecHitSoA/src/alpaka/classes_cuda_def.xml @@ -13,4 +13,9 @@ + + + + + diff --git a/DataFormats/TrackingRecHitSoA/src/alpaka/classes_rocm_def.xml b/DataFormats/TrackingRecHitSoA/src/alpaka/classes_rocm_def.xml index bc4c969137121..76f2fcf60f77b 100644 --- a/DataFormats/TrackingRecHitSoA/src/alpaka/classes_rocm_def.xml +++ b/DataFormats/TrackingRecHitSoA/src/alpaka/classes_rocm_def.xml @@ -13,5 +13,10 @@ + + + + + diff --git a/DataFormats/TrackingRecHitSoA/src/classes.cc b/DataFormats/TrackingRecHitSoA/src/classes.cc index bbcc923b04373..dcdd2c8378d76 100644 --- a/DataFormats/TrackingRecHitSoA/src/classes.cc +++ b/DataFormats/TrackingRecHitSoA/src/classes.cc @@ -1,7 +1,8 @@ #include "DataFormats/Portable/interface/PortableHostCollectionReadRules.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); +SET_PORTABLEHOSTCOLLECTION_READ_RULES(PortableHostCollection>); diff --git a/DataFormats/TrackingRecHitSoA/src/classes_def.xml b/DataFormats/TrackingRecHitSoA/src/classes_def.xml index 54c0a3d30a365..c7b7f09b19c1f 100644 --- a/DataFormats/TrackingRecHitSoA/src/classes_def.xml +++ b/DataFormats/TrackingRecHitSoA/src/classes_def.xml @@ -25,4 +25,13 @@ + + + + + + + + + diff --git a/Geometry/CommonTopologies/interface/SimplePixelTopology.h b/Geometry/CommonTopologies/interface/SimplePixelTopology.h index 98ea51c18ed01..efa9222dd35a5 100644 --- a/Geometry/CommonTopologies/interface/SimplePixelTopology.h +++ b/Geometry/CommonTopologies/interface/SimplePixelTopology.h @@ -1,8 +1,10 @@ #ifndef Geometry_CommonTopologies_SimplePixelTopology_h #define Geometry_CommonTopologies_SimplePixelTopology_h +#include #include #include +#include #include #include "FWCore/Utilities/interface/HostDeviceConstant.h" @@ -27,6 +29,7 @@ namespace pixelTopology { constexpr int16_t phi0p06 = 626; // round(625.82270...) = phi2short(0.06); constexpr int16_t phi0p07 = 730; // round(730.12648...) = phi2short(0.07); constexpr int16_t phi0p09 = 900; + constexpr int16_t phi5deg = 1820; constexpr uint16_t last_barrel_layer = 3; // this is common between all the topologies @@ -131,8 +134,9 @@ namespace phase1PixelTopology { using pixelTopology::phi0p06; using pixelTopology::phi0p07; - constexpr uint32_t numberOfLayers = 28; + constexpr uint32_t numberOfLayers = 10; constexpr int nPairs = 13 + 2 + 4; + ; // without jump + jumping barrel + jumping forward constexpr uint16_t numberOfModules = 1856; constexpr uint32_t maxNumClustersPerModules = 1024; @@ -230,35 +234,36 @@ namespace phase2PixelTopology { 4, 6, 5, 7, 6, 8, 7, 9, 8, 10, 9, 11, 10, 12, // POS Jump (48) 16, 18, 17, 19, 18, 20, 19, 21, 20, 22, 21, 23, 22, 24, // NEG Jump (55) }; - HOST_DEVICE_CONSTANT uint32_t layerStart[numberOfLayers + 1] = {0, - 108, - 324, - 504, // Barrel - 756, - 864, - 972, - 1080, - 1188, - 1296, - 1404, - 1512, - 1620, - 1796, - 1972, - 2148, // Fp - 2324, - 2432, - 2540, - 2648, - 2756, - 2864, - 2972, - 3080, - 3188, - 3364, - 3540, - 3716, // Np - numberOfModules}; + + static constexpr uint32_t layerStart[numberOfLayers + 1] = {0, + 108, + 324, + 504, // Barrel + 756, + 864, + 972, + 1080, + 1188, + 1296, + 1404, + 1512, + 1620, + 1796, + 1972, + 2148, // Fp + 2324, + 2432, + 2540, + 2648, + 2756, + 2864, + 2972, + 3080, + 3188, + 3364, + 3540, + 3716, // Np + numberOfModules}; HOST_DEVICE_CONSTANT int16_t phicuts[nPairs]{ phi0p05, phi0p05, phi0p05, phi0p06, phi0p07, phi0p07, phi0p06, phi0p07, phi0p07, phi0p05, phi0p05, @@ -337,7 +342,10 @@ namespace pixelTopology { static constexpr uint32_t maxNumOfActiveDoublets = maxNumberOfDoublets / 8; static constexpr uint32_t maxNumberOfQuadruplets = maxNumberOfTuples; static constexpr uint32_t maxDepth = 12; - static constexpr uint32_t numberOfLayers = 28; + + static constexpr uint32_t numberOfPixelLayers = 28; + static constexpr uint32_t numberOfStripLayers = 0; + static constexpr uint32_t numberOfLayers = numberOfStripLayers + numberOfPixelLayers; static constexpr uint32_t maxSizeCluster = 2047; @@ -378,7 +386,9 @@ namespace pixelTopology { static constexpr int maxDYsize = 10; static constexpr int maxDYPred = 20; - static constexpr uint16_t numberOfModules = phase2PixelTopology::numberOfModules; + static constexpr uint16_t numberOfModules = 3892; + static constexpr uint16_t numberOfPixelModules = 3892; + static constexpr uint16_t numberOfStripModules = 0; // 1000 bins < 1024 bins (10 bits) must be: // - < 32*32 (warpSize*warpSize for block prefix scan for CUDA) @@ -413,12 +423,14 @@ namespace pixelTopology { static constexpr inline uint16_t localX(uint16_t px) { return px; } static constexpr inline uint16_t localY(uint16_t py) { return py; } + + static constexpr int mapIndex(int i) { return i; } }; struct Phase1 { // types using hindex_type = uint32_t; // FIXME from siPixelRecHitsHeterogeneousProduct - using tindex_type = uint16_t; // for tuples + using tindex_type = uint32_t; // for tuples using cindex_type = uint32_t; // for cells static constexpr uint32_t maxCellNeighbors = 36; @@ -434,7 +446,10 @@ namespace pixelTopology { static constexpr uint32_t maxNumOfActiveDoublets = maxNumberOfDoublets / 8; static constexpr uint32_t maxNumberOfQuadruplets = maxNumberOfTuples; static constexpr uint32_t maxDepth = 6; - static constexpr uint32_t numberOfLayers = 10; + + static constexpr uint32_t numberOfPixelLayers = 10; + static constexpr uint32_t numberOfStripLayers = 0; + static constexpr uint32_t numberOfLayers = numberOfPixelLayers + numberOfStripLayers; static constexpr uint32_t maxSizeCluster = 1023; @@ -464,8 +479,8 @@ namespace pixelTopology { static constexpr float dzdrFact = 8 * 0.0285 / 0.015; // from dz/dr to "DY" - static constexpr int minYsizeB1 = 36; - static constexpr int minYsizeB2 = 28; + static constexpr int minYsizeB1 = 1; + static constexpr int minYsizeB2 = 1; static constexpr int nPairsForQuadruplets = 13; // quadruplets require hits in all layers static constexpr int nPairsForTriplets = nPairsForQuadruplets + 2; // include barrel "jumping" layer pairs @@ -475,7 +490,9 @@ namespace pixelTopology { static constexpr int maxDYsize = 20; static constexpr int maxDYPred = 20; - static constexpr uint16_t numberOfModules = phase1PixelTopology::numberOfModules; + static constexpr uint16_t numberOfModules = 1856; + static constexpr uint16_t numberOfPixelModules = phase1PixelTopology::layerStart[numberOfPixelLayers]; + static constexpr uint16_t numberOfStripModules = numberOfModules - numberOfPixelModules; static constexpr uint16_t numRowsInRoc = 80; static constexpr uint16_t numColsInRoc = 52; @@ -544,6 +561,8 @@ namespace pixelTopology { shift += 1; return py + shift; } + + static constexpr int mapIndex(int i) { return i; } }; struct HIonPhase1 : public Phase1 { @@ -577,6 +596,9 @@ namespace pixelTopology { template using isPhase2Topology = typename std::enable_if::value>::type; + template + using base_traits_t = std::conditional_t, Phase1, Phase2>; + } // namespace pixelTopology #endif // Geometry_CommonTopologies_SimplePixelTopology_h diff --git a/Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h b/Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h new file mode 100644 index 0000000000000..9566ad2662253 --- /dev/null +++ b/Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h @@ -0,0 +1,293 @@ +#ifndef Geometry_CommonTopologies_SimpleSeedingLayersTopology_h +#define Geometry_CommonTopologies_SimpleSeedingLayersTopology_h +#include +#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +namespace phase1PixelStripTopology { + + struct LayerData { + uint32_t start; + uint32_t end; + bool isStrip2D = false; // if true then map every two module indices to the same sequential ID + }; + + enum Layer : uint8_t; + + struct LayerPairData { + Layer inner; + Layer outer; + int16_t phicut; + float minz; + float maxz; + float maxr; + }; + + enum Layer : uint8_t { + BPIX1 = 0, + BPIX2, + BPIX3, + BPIX4, + FPIX1Pos, + FPIX2Pos, + FPIX3Pos, + FPIX1Neg, + FPIX2Neg, + FPIX3Neg, + TIB1, + TIB2, + TID1Pos2D, + TID2Pos2D, + TID3Pos2D, + TID1Neg2D, + TID2Neg2D, + TID3Neg2D, + nLayers + }; + + constexpr LayerData layerData[nLayers] = { + {0, 96}, // BPIX1 + {96, 320}, // BPIX2 + {320, 672}, // BPIX3 + {672, 1184}, // BPIX4 + {1184, 1296}, // FPIX1Pos + {1296, 1408}, // FPIX2Pos + {1408, 1520}, // FPIX3Pos + {1520, 1632}, // FPIX1Neg + {1632, 1744}, // FPIX2Neg + {1744, 1856}, // FPIX3Neg + {1856, 2528, true}, // TIB1 + {2528, 3392, true}, // TIB2 + {4580, 4676, true}, // TID1Pos2D + {4716, 4812, true}, // TID2Pos2D + {4852, 4948, true}, // TID3Pos2D + {4988, 5084, true}, // TID1Neg2D + {5124, 5220, true}, // TID2Neg2D + {5260, 5356, true}, // TID3Neg2D + }; + + using pixelTopology::phi0p05; + using pixelTopology::phi0p06; + using pixelTopology::phi0p07; + using pixelTopology::phi0p09; + using pixelTopology::phi5deg; + //constexpr int16_t phi0p05 = 522; // round(521.52189...) = phi2short(0.05); + //constexpr int16_t phi0p06 = 626; // round(625.82270...) = phi2short(0.06); + //constexpr int16_t phi0p07 = 730; // round(730.12648...) = phi2short(0.07); + //constexpr int16_t phi0p09 = 900; + //constexpr int16_t phi5deg = 1820; + constexpr LayerPairData layerPairData[] = { + {BPIX1, BPIX2, phi0p05, -20., 20., 20.}, // 0 + {BPIX1, FPIX1Pos, phi0p07, 0., 30., 9.}, // 1 + {BPIX1, FPIX1Neg, phi0p07, -30., 0., 9.}, // 2 + {BPIX2, BPIX3, phi0p05, -22., 22., 20.}, // 3 + {BPIX2, FPIX1Pos, phi0p07, 10., 30., 7.}, // 4 + {BPIX2, FPIX1Neg, phi0p06, -30., -10., 7.}, // 5 + {FPIX1Pos, FPIX2Pos, phi0p06, -70., 70., 5.}, // 6 + {FPIX1Neg, FPIX2Neg, phi0p05, -70., 70., 5.}, // 7 + {BPIX1, BPIX3, phi0p05, -20., 20., 20.}, // 8 + {BPIX2, BPIX4, phi0p05, -22., 22., 20.}, // 9 + {BPIX1, FPIX2Pos, phi0p06, 0., 30., 9.}, // 10 + {BPIX1, FPIX2Neg, phi0p05, -30., 0., 9.}, // 11 + {FPIX1Pos, TIB1, 1200, -70., 70., 1000.}, // 12 + {FPIX1Neg, TIB1, 1200, -70., 70., 1000.}, // 13 + {BPIX3, BPIX4, phi0p06, -22., 22., 20.}, // 14 + {BPIX3, FPIX1Pos, phi0p07, 15., 30., 6.}, // 15 + {BPIX3, FPIX1Neg, phi0p06, -30, -15., 6.}, // 16 + {FPIX2Pos, FPIX3Pos, phi0p06, -70., 70., 5.}, // 17 + {FPIX2Neg, FPIX3Neg, phi0p05, -70., 70., 5.}, // 18 + {BPIX3, TIB1, phi5deg, -22., 22., 1000.}, // 19 + {BPIX4, TIB1, phi5deg, -22., 22., 1000.}, // 20 + {BPIX4, TIB2, phi5deg, -22., 22., 1000.}, // 21 + {TIB1, TIB2, phi5deg, -55., 55., 1000.}, // 22 + {FPIX2Neg, TIB1, phi5deg, -70., 70., 1000.}, // 23 + {FPIX3Neg, TIB1, phi5deg, -70., 70., 1000.}, // 24 + {TIB1, TID1Neg2D, phi5deg, -55., 0., 1000.}, // 25 + {TIB2, TID1Neg2D, phi5deg, -55., 0., 1000.}, // 26 + {BPIX2, TIB1, phi5deg, -22., 0., 1000.}, // 27 + {BPIX2, TIB2, phi5deg, -22., 0., 1000.}, // 28 + {BPIX1, TIB1, phi5deg, -22., 0., 1000.}, // 29 + {BPIX3, TIB2, phi5deg, -22., 22., 1000.}, // 30 + {BPIX4, TID1Neg2D, phi5deg, -55., 0., 1000.}, // 31 + {FPIX1Pos, FPIX3Pos, phi0p06, -70., 70., 9.}, // 32 + {FPIX1Neg, FPIX3Neg, phi0p05, -70., 70., 9.}, // 33 + {FPIX1Neg, TIB2, phi5deg, -70., 70., 1000.}, // 34 + {FPIX2Neg, TIB2, phi5deg, -70., 70., 1000.}, // 35 + {FPIX3Neg, TIB2, phi5deg, -70., 70., 1000.}, // 36 + {BPIX2, FPIX2Neg, phi5deg, -30., 0., 1000.} // 37 + }; + + constexpr uint32_t maxNumClustersPerModules = 1024; + constexpr auto numberOfLayers = nLayers; + constexpr auto nPairs = std::size(layerPairData); + + constexpr auto makeLayerStart() { + std::array layerStart = {{0}}; + for (auto i = 0u; i < numberOfLayers; ++i) + if (layerData[i].isStrip2D) + layerStart[i + 1] = layerStart[i] + (layerData[i].end - layerData[i].start) / 2; + else + layerStart[i + 1] = layerStart[i] + layerData[i].end - layerData[i].start; + return layerStart; + } + + HOST_DEVICE_CONSTANT std::array layerStart = makeLayerStart(); + + constexpr uint16_t numberOfModules = layerStart[numberOfLayers]; + + constexpr auto makeLayerPairs() { + std::array layerPairs = {{0}}; + for (auto i = 0u; i < nPairs; ++i) { + layerPairs[2 * i] = layerPairData[i].inner; + layerPairs[2 * i + 1] = layerPairData[i].outer; + } + return layerPairs; + } + + HOST_DEVICE_CONSTANT std::array layerPairs = makeLayerPairs(); + + template + constexpr auto makePairCutsArray(F&& f) { + std::array result{}; + for (auto i = 0u; i < nPairs; ++i) + result[i] = f(layerPairData[i]); + return result; + } + + HOST_DEVICE_CONSTANT std::array phicuts = + makePairCutsArray([](const auto& x) { return x.phicut; }); + HOST_DEVICE_CONSTANT std::array minz = makePairCutsArray([](const auto& x) { return x.minz; }); + HOST_DEVICE_CONSTANT std::array maxz = makePairCutsArray([](const auto& x) { return x.maxz; }); + HOST_DEVICE_CONSTANT std::array maxr = makePairCutsArray([](const auto& x) { return x.maxr; }); + + using IndexMap = std::array; + + constexpr IndexMap makeIndexMap() { + IndexMap indexMap = {{0}}; + for (auto& i : indexMap) + i = numberOfModules; + int newIndex = 0; + for (const auto& layer : layerData) { + for (auto i = layer.start; i < layer.end; ++i) { + indexMap[i] = newIndex; + if (!layer.isStrip2D || i % 2 == 1) + ++newIndex; + } + } + return indexMap; + } + + constexpr IndexMap indexMap = makeIndexMap(); + + // constexpr uint32_t numberOfLayers = 12; + // constexpr int nPairs = 21 + 4 + 10 + 1; // without jump + jumping barrel + jumping forward + // constexpr uint16_t numberOfModules = 3392; + + // HOST_DEVICE_CONSTANT uint8_t layerPairs[2 * nPairs] = { + // 0, 1, 0, 4, 0, 7, // BPIX1 (3) + // 1, 2, 1, 4, 1, 7, // BPIX2 (6) + // 4, 5, 7, 8, // FPIX1 (8) + // 2, 3, 2, 4, 2, 7, 5, 6, 8, 9, // BPIX3 & FPIX2 (13) + // 0, 2, 1, 3, // Jumping Barrel (15) + // 0, 5, 0, 8, // Jumping Forward (BPIX1,FPIX2) + // 4, 6, 7, 9, // Jumping Forward (19) + // 3, 10, // BPIX4 (20) + // 4, 10, 5, 10, 6, 10, // Pixel Positive Endcap (23) + // 7, 10, 8, 10, 9, 10, // Pixel Negative Endcap (26) + // 10, 11, // TIB1 (27) + // 1, 10, 2, 10, 3, 11, // Jumping from Pixel Barrel (30) + // 4, 11, 5, 11, 6, 11, // Jumping from Pixel Positive Endcap (33) + // 7, 11, 8, 11, 9, 11 // Jumping from Pixel Negative Endcap (36) + // }; + + // HOST_DEVICE_CONSTANT int16_t phicuts[nPairs]{phi0p05, phi0p07, phi0p07, phi0p05, phi0p06, + // phi0p06, phi0p05, phi0p05, phi0p06, phi0p06, + // phi0p06, phi0p05, phi0p05, phi0p05, phi0p05, + + // phi0p05, phi0p05, phi0p05, phi0p05, phi5deg, + // phi5deg, phi5deg, phi5deg, phi0p09, phi0p09, + // phi0p09, phi0p09, phi0p09, phi0p09, phi0p09, + + // phi0p09, phi0p09, phi0p09, phi0p09, phi0p09, + // phi0p09 + // }; + + // HOST_DEVICE_CONSTANT float minz[nPairs] = { + // -20., 0., -30., -22., 10., -30., -70., -70., -22., 15., -30, -70., -70., -20., -22., 0, -30., -70., -70., + // -22.,-70.,-70.,-70.,-70.,-70.,-70.,-80.,-22.,-22.,-70.,-70.,-70.,-70.,-70.,-70.,-70.}; + // HOST_DEVICE_CONSTANT float maxz[nPairs] = { + // 20., 30., 0., 22., 30., -10., 70., 70., 22., 30., -15., 70., 70., 20., 22., 30., 0., 70., 70., + // 22.,70.,70.,70.,70.,70.,70.,80.,22.,22.,70.,70., 70.,70.,70.,70.,70.}; + // HOST_DEVICE_CONSTANT float maxr[nPairs] = { + // 20., 9., 9., 20., 7., 7., 5., 5., 20., 6., 6., 5., 5., 20., 20., 9., 9., 9., 9., + // 10000.,10000.,10000.,10000.,10000.,10000.,10000.,10000.,10000.,10000.,10000.,10000., 10000.,10000.,10000.,10000.,10000.}; + + // static constexpr uint32_t layerStart[numberOfLayers + 1] = {0, + // 96, + // 320, + // 672, // barrel + // 1184, + // 1296, + // 1408, // positive endcap + // 1520, + // 1632, + // 1744, // negative endcap + // 1856, + // 2528, + // numberOfModules}; + +} // namespace phase1PixelStripTopology +namespace pixelTopology { + + struct Phase1Strip : public Phase1 { + typedef Phase1 PixelBase; //Could be removed using based class + static constexpr uint32_t maxNumClustersPerModules = phase1PixelStripTopology::maxNumClustersPerModules; + static constexpr uint32_t maxHitsInModule = phase1PixelStripTopology::maxNumClustersPerModules; + static constexpr uint32_t maxCellNeighbors = 64; + static constexpr uint32_t maxCellTracks = 90; + static constexpr uint32_t maxHitsOnTrack = 15; + static constexpr uint32_t maxHitsOnTrackForFullFit = 6; + static constexpr uint32_t avgHitsPerTrack = 7; + static constexpr uint32_t maxCellsPerHit = 256; + static constexpr uint32_t avgTracksPerHit = 10; + static constexpr uint32_t maxNumberOfTuples = 32 * 1024 * 4; + //this is well above thanks to maxNumberOfTuples + static constexpr uint32_t maxHitsForContainers = avgHitsPerTrack * maxNumberOfTuples; + static constexpr uint32_t maxNumberOfDoublets = 10 * 256 * 1024; + static constexpr uint32_t maxNumOfActiveDoublets = maxNumberOfDoublets / 8; + static constexpr uint32_t maxNumberOfQuadruplets = maxNumberOfTuples; + static constexpr uint32_t maxDepth = 12; + static constexpr int minYsizeB1 = 1; + static constexpr int minYsizeB2 = 1; + static constexpr uint32_t const* layerStart = phase1PixelStripTopology::layerStart.data(); + + static constexpr float const* minz = phase1PixelStripTopology::minz.data(); + static constexpr float const* maxz = phase1PixelStripTopology::maxz.data(); + static constexpr float const* maxr = phase1PixelStripTopology::maxr.data(); + + static constexpr uint8_t const* layerPairs = phase1PixelStripTopology::layerPairs.data(); + static constexpr int16_t const* phicuts = phase1PixelStripTopology::phicuts.data(); + + static constexpr uint32_t numberOfLayers = phase1PixelStripTopology::numberOfLayers; + static constexpr uint32_t numberOfStripLayers = numberOfLayers - numberOfPixelLayers; + + static constexpr uint16_t numberOfModules = phase1PixelStripTopology::numberOfModules; + static constexpr uint16_t numberOfPixelModules = phase1PixelStripTopology::layerStart[numberOfPixelLayers]; + static constexpr uint16_t numberOfStripModules = numberOfModules - numberOfPixelModules; + + static constexpr int nPairsForQuadruplets = + phase1PixelStripTopology::nPairs; // quadruplets require hits in all layers + static constexpr int nPairsForTriplets = nPairsForQuadruplets; // include barrel "jumping" layer pairs + static constexpr int nPairs = nPairsForTriplets; // include forward "jumping" layer pairs + + static constexpr char const* nameModifier = "Phase1Strip"; + + static constexpr int mapIndex(int i) { + if (i > 0 && i < (int)phase1PixelStripTopology::indexMap.size()) + return phase1PixelStripTopology::indexMap[i]; + else + return i; + } + }; + +} // namespace pixelTopology + +#endif // Geometry_CommonTopologies_SimpleSeedingLayersTopology_h diff --git a/HLTrigger/Configuration/python/customizeHLTforAlpakaStripNoDoubletRecovery.py b/HLTrigger/Configuration/python/customizeHLTforAlpakaStripNoDoubletRecovery.py new file mode 100644 index 0000000000000..98086b3c0c810 --- /dev/null +++ b/HLTrigger/Configuration/python/customizeHLTforAlpakaStripNoDoubletRecovery.py @@ -0,0 +1,363 @@ +import FWCore.ParameterSet.Config as cms +import re +import itertools + +from FWCore.ParameterSet.MassReplace import massReplaceInputTag +from HeterogeneousCore.AlpakaCore.functions import * +from HLTrigger.Configuration.common import * + +def customizeHLTforAlpakaPixelRecoLocal(process): + '''Customisation to introduce the Local Pixel and Strip Reconstruction in Alpaka + ''' + + if not hasattr(process, 'HLTDoLocalPixelSequence'): + return process + + process.hltSiStripRawToClustersFacility = cms.EDProducer("SiStripClusterizerFromRaw", + Algorithms = cms.PSet( + CommonModeNoiseSubtractionMode = cms.string('Median'), + PedestalSubtractionFedMode = cms.bool(True), + SiStripFedZeroSuppressionMode = cms.uint32(4), + TruncateInSuppressor = cms.bool(True), + Use10bitsTruncation = cms.bool(False), + doAPVRestore = cms.bool(False), + useCMMeanMap = cms.bool(False) + ), + Clusterizer = cms.PSet( + Algorithm = cms.string('ThreeThresholdAlgorithm'), + ChannelThreshold = cms.double(2.0), + ClusterThreshold = cms.double(5.0), + ConditionsLabel = cms.string(''), + MaxAdjacentBad = cms.uint32(0), + MaxClusterSize = cms.uint32(8), + MaxSequentialBad = cms.uint32(1), + MaxSequentialHoles = cms.uint32(0), + RemoveApvShots = cms.bool(True), + SeedThreshold = cms.double(3.0), + clusterChargeCut = cms.PSet( + refToPSet_ = cms.string('HLTSiStripClusterChargeCutNone') + ), + setDetId = cms.bool(True) + ), + ConditionsLabel = cms.string(''), + DoAPVEmulatorCheck = cms.bool(False), + HybridZeroSuppressed = cms.bool(False), + LegacyUnpacker = cms.bool(False), + ProductLabel = cms.InputTag("rawDataCollector"), + onDemand = cms.bool(False) + ) + + + process.hltSiStripMatchedRecHitsFull = cms.EDProducer( "SiStripRecHitConverter", + ClusterProducer = cms.InputTag( "hltSiStripRawToClustersFacility" ), + rphiRecHits = cms.string( "rphiRecHit" ), + stereoRecHits = cms.string( "stereoRecHit" ), + matchedRecHits = cms.string( "matchedRecHit" ), + useSiStripQuality = cms.bool( False ), + MaskBadAPVFibers = cms.bool( False ), + doMatching = cms.bool( True ), + StripCPE = cms.ESInputTag( "hltESPStripCPEfromTrackAngle","hltESPStripCPEfromTrackAngle" ), + Matcher = cms.ESInputTag( "SiStripRecHitMatcherESProducer","StandardMatcher" ), + siStripQualityLabel = cms.ESInputTag( "","" ) + ) + + process.hltSiPixelOnlyRecHitsSoA = cms.EDProducer('SiPixelRecHitAlpakaPhase1@alpaka', + beamSpot = cms.InputTag('hltOnlineBeamSpotDevice'), + src = cms.InputTag('hltSiPixelClustersSoA'), + CPE = cms.string('PixelCPEFastParams'), + mightGet = cms.optional.untracked.vstring, + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('') + ) + ) + + process.hltSiPixelRecHitsSoA = cms.EDProducer('SiStripRecHitSoAPhase1@alpaka', + stripRecHitSource = cms.InputTag('hltSiStripMatchedRecHitsFull', 'matchedRecHit'), + beamSpot = cms.InputTag('hltOnlineBeamSpot'), + pixelRecHitSoASource = cms.InputTag('hltSiPixelOnlyRecHitsSoA'), + mightGet = cms.optional.untracked.vstring, + + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('') + ) + ) + process.hltSiPixelRecHits = cms.EDProducer('SiPixelRecHitFromSoAAlpakaPhase1', + pixelRecHitSrc = cms.InputTag('hltSiPixelOnlyRecHitsSoA'), + src = cms.InputTag('hltSiPixelClusters'), + ) + + ### + ### Sequence: Addition of strip modules in Pixel Local Reconstruction + ### + + process.HLTDoLocalPixelSequence.insert(process.HLTDoLocalPixelSequence.index(process.hltSiPixelDigiErrors)+1,process.hltSiStripRawToClustersFacility+process.hltSiStripMatchedRecHitsFull+process.hltSiPixelOnlyRecHitsSoA) + + ### + ### SerialSync version of Pixel Local Reconstruction + ### + + process.hltSiPixelOnlyRecHitsSoACPUSerial = makeSerialClone(process.hltSiPixelOnlyRecHitsSoA, + beamSpot = 'hltOnlineBeamSpotDeviceSerialSync', + src = 'hltSiPixelClustersSoASerialSync' + ) + process.hltSiPixelRecHitsSoASerialSync = makeSerialClone(process.hltSiPixelRecHitsSoA, + pixelRecHitSoASource = cms.InputTag('hltSiPixelOnlyRecHitsSoACPUSerial'), + ) + + process.hltSiPixelRecHitsSerialSync = process.hltSiPixelRecHits.clone( + pixelRecHitSrc = 'hltSiPixelOnlyRecHitsSoACPUSerial', + src = 'hltSiPixelClustersSerialSync', + ) + + if(not hasattr(process,'hltSiPixelOnlyRecHitsSoA')): + return process + if not hasattr(process, 'HLTDoLocalPixelSequenceSerialSync'): + return process + process.HLTDoLocalPixelSequenceSerialSync.insert(process.HLTDoLocalPixelSequenceSerialSync.index(process.hltSiPixelDigiErrorsSerialSync)+1,process.hltSiStripRawToClustersFacility+process.hltSiStripMatchedRecHitsFull+process.hltSiPixelOnlyRecHitsSoACPUSerial) + return process + + +def customizeHLTforAlpakaPixelRecoTracking(process): + '''Customisation to introduce the Pixel-Track Reconstruction in Alpaka + ''' + + if not hasattr(process, 'HLTRecoPixelTracksSequence'): + return process + + for producer in producers_by_type(process, "TrackListMerger"): + current_producers = producer.TrackProducers + if ( + 'hltIter0PFlowTrackSelectionHighPurity' in current_producers and + 'hltDoubletRecoveryPFlowTrackSelectionHighPurity' in current_producers + ): + setattr(producer, "TrackProducers",cms.VInputTag('hltIter0PFlowTrackSelectionHighPurity')) + setattr(producer,"hasSelector",cms.vint32( 0)) + setattr(producer,"indivShareFrac",cms.vdouble( 1.0)) + setattr(producer, "selectedTrackQuals", cms.VInputTag( 'hltIter0PFlowTrackSelectionHighPurity')) + setattr(producer,"setsToMerge",cms.VPSet( cms.PSet( pQual = cms.bool( False ), tLists = cms.vint32( 0)))) + + if hasattr(process, "hltDoubletRecoveryPFlowTrackSelectionHighPurity"): + del process.hltDoubletRecoveryPFlowTrackSelectionHighPurity + + for producer in producers_by_type(process, "TrackListMerger"): + current_producers = producer.TrackProducers + if ( + 'hltIter0PFlowTrackSelectionHighPuritySerialSync' in current_producers and + 'hltDoubletRecoveryPFlowTrackSelectionHighPuritySerialSync' in current_producers + ): + setattr(producer, "TrackProducers",cms.VInputTag('hltIter0PFlowTrackSelectionHighPuritySerialSync')) + setattr(producer,"hasSelector",cms.vint32( 0)) + setattr(producer,"indivShareFrac",cms.vdouble( 1.0)) + setattr(producer, "selectedTrackQuals", cms.VInputTag( 'hltIter0PFlowTrackSelectionHighPuritySerialSync')) + setattr(producer,"setsToMerge",cms.VPSet( cms.PSet( pQual = cms.bool( False ), tLists = cms.vint32( 0)))) + if hasattr(process, "HLTIterativeTrackingDoubletRecovery"): + del process.HLTIterativeTrackingDoubletRecovery + if hasattr(process, "HLTIterativeTrackingDoubletRecoverySerialSync"): + del process.HLTIterativeTrackingDoubletRecoverySerialSync + if hasattr(process, "hltDoubletRecoveryPFlowTrackSelectionHighPuritySerialSync"): + del process.hltDoubletRecoveryPFlowTrackSelectionHighPuritySerialSync + + # alpaka EDProducer + # consumes + # - TrackingRecHitsSoACollection + # produces + # - TkSoADevice + + process.frameSoAESProducerPhase1Strip = cms.ESProducer('FrameSoAESProducerPhase1Strip@alpaka', + ComponentName = cms.string('FrameSoAPhase1Strip'), + appendToDataLabel = cms.string(''), + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('') + ) + ) + + process.hltPixelTracksSoA = cms.EDProducer('CAHitNtupletAlpakaPhase1Strip@alpaka', + pixelRecHitSrc = cms.InputTag('hltSiPixelRecHitsSoA'), + frameSoA = cms.string('FrameSoAPhase1Strip'), + ptmin = cms.double(0.9), + maxNumberOfDoublets = cms.uint32(8*256*1024), + CAThetaCutBarrel = cms.double(0.002), + CAThetaCutForward = cms.double(0.003), + hardCurvCut = cms.double(0.0328407225), + dcaCutInnerTriplet = cms.double(0.15), + dcaCutOuterTriplet = cms.double(0.25), + CAThetaCutBarrelPixelBarrelStrip = cms.double(0.002), + CAThetaCutBarrelPixelForwardStrip = cms.double(0.002), + CAThetaCutBarrelStripForwardStrip = cms.double(0.002), + CAThetaCutBarrelStrip = cms.double(0.002), + CAThetaCutDefault = cms.double(0.003), + dcaCutInnerTripletPixelStrip = cms.double(0.15), + dcaCutOuterTripletPixelStrip = cms.double(0.25), + dcaCutTripletStrip = cms.double(0.25), + dcaCutTripletDefault = cms.double(0.25), + earlyFishbone = cms.bool(True), + lateFishbone = cms.bool(False), + fillStatistics = cms.bool(False), + minHitsPerNtuplet = cms.uint32(3), + minHitsForSharingCut = cms.uint32(10), + fitNas4 = cms.bool(False), + doClusterCut = cms.bool(True), + doZ0Cut = cms.bool(True), + cellZ0Cut = cms.double(10.0), + cellPtCut = cms.double(0.5), + doPtCut = cms.bool(True), + useRemovers = cms.bool(True), + useRiemannFit = cms.bool(False), + doSharedHitCut = cms.bool(True), + dupPassThrough = cms.bool(False), + useSimpleTripletCleaner = cms.bool(True), + idealConditions = cms.bool(False), + includeJumpingForwardDoublets = cms.bool(True), + trackQualityCuts = cms.PSet( + chi2MaxPt = cms.double(10), + chi2Coeff = cms.vdouble( + 0.9, + 1.8 + ), + chi2Scale = cms.double(8), + tripletMinPt = cms.double(0.5), + tripletMaxTip = cms.double(0.3), + tripletMaxZip = cms.double(12), + quadrupletMinPt = cms.double(0.3), + quadrupletMaxTip = cms.double(0.5), + quadrupletMaxZip = cms.double(12) + ), + phiCuts = cms.vint32( + 522, 730, 730, 522, 730, 626, 626, 522, 522, 522, 626, 522, 1200, 1200, 626, 730, 626, 626, 522, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 626, 522, 2000, 2000, 2000, 2000), + # autoselect the alpaka backend + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('') + ) + ) + + if hasattr(process, 'hltL2TauTagNNProducer'): + process.hltL2TauTagNNProducer = cms.EDProducer("L2TauNNProducerAlpakaStrip", **process.hltL2TauTagNNProducer.parameters_()) + process.hltPixelTracksSoASerialSync = makeSerialClone(process.hltPixelTracksSoA, + pixelRecHitSrc = 'hltSiPixelRecHitsSoASerialSync' + ) + + process.hltPixelTracks = cms.EDProducer("PixelTrackProducerFromSoAAlpakaPhase1Strip", + beamSpot = cms.InputTag("hltOnlineBeamSpot"), + minNumberOfHits = cms.int32(0), + minQuality = cms.string('loose'), + pixelRecHitLegacySrc = cms.InputTag("hltSiPixelRecHits"), + trackSrc = cms.InputTag("hltPixelTracksSoA"), + useStripHits = cms.bool(True), + hitModuleStartSrc = cms.InputTag("hltSiPixelRecHitsSoA"), + stripRecHitLegacySrc = cms.InputTag('hltSiStripMatchedRecHitsFull', 'matchedRecHit'), + mightGet = cms.optional.untracked.vstring + ) + process.hltESPTTRHBuilderPixelOnly = cms.ESProducer("TkTransientTrackingRecHitBuilderESProducer", + ComponentName = cms.string('hltESPTTRHBuilderPixelOnly'), + ComputeCoarseLocalPositionFromDisk = cms.bool(False), + Matcher = cms.string('StandardMatcher'), + Phase2StripCPE = cms.string(''), + PixelCPE = cms.string('hltESPPixelCPEGeneric'), + StripCPE = cms.string('hltESPStripCPEfromTrackAngle'), + appendToDataLabel = cms.string('') + ) + process.hltPixelTracksSerialSync = process.hltPixelTracks.clone( + pixelRecHitLegacySrc = cms.InputTag("hltSiPixelRecHitsSerialSync"), + hitModuleStartSrc = cms.InputTag("hltSiPixelRecHitsSoASerialSync"), + trackSrc = cms.InputTag("hltPixelTracksSoASerialSync") + ) + process.HLTRecoPixelTracksTask = cms.ConditionalTask( + process.hltPixelTracksSoA, + process.hltPixelTracks, + ) + process.HLTRecoPixelTracksCPUSerialTask = cms.ConditionalTask( + process.hltPixelTracksSoASerialSync, + process.hltPixelTracksSerialSync, + ) + process.HLTRecoPixelTracksSequence = cms.Sequence(process.hltPixelTracksSoA+process.hltPixelTracks) + if hasattr(process, 'hltPixelTracksSerialSync'): + process.HLTRecoPixelTracksCPUSerialSequence = cms.Sequence(process.hltPixelTracksSoASerialSync+process.hltPixelTracksSerialSync) + + return process + +def customizeHLTforAlpakaPixelRecoVertexing(process): + '''Customisation to introduce the Pixel-Vertex Reconstruction in Alpaka + ''' + + # alpaka EDProducer + # consumes + # - TkSoADevice + # produces + # - ZVertexDevice + + process.hltPixelVerticesSoA = cms.EDProducer('PixelVertexProducerAlpakaPhase1Strip@alpaka', + oneKernel = cms.bool(True), + useDensity = cms.bool(True), + useDBSCAN = cms.bool(False), + useIterative = cms.bool(False), + minT = cms.int32(2), + eps = cms.double(0.07), + errmax = cms.double(0.01), + chi2max = cms.double(9), + PtMin = cms.double(0.5), + PtMax = cms.double(75), + pixelTrackSrc = cms.InputTag('hltPixelTracksSoA'), + # autoselect the alpaka backend + alpaka = cms.untracked.PSet( + backend = cms.untracked.string('') + ) + ) + + process.hltPixelVerticesSoASerialSync = makeSerialClone(process.hltPixelVerticesSoA, + pixelTrackSrc = 'hltPixelTracksSoASerialSync' + ) + + process.hltPixelVertices = cms.EDProducer("PixelVertexProducerFromSoAAlpaka", + TrackCollection = cms.InputTag("hltPixelTracks"), + beamSpot = cms.InputTag("hltOnlineBeamSpot"), + src = cms.InputTag("hltPixelVerticesSoA") + ) + + process.hltPixelVerticesSerialSync = process.hltPixelVertices.clone( + TrackCollection = cms.InputTag("hltPixelTracksSerialSync"), + src = cms.InputTag("hltPixelVerticesSoASerialSync") + ) + + if hasattr(process, 'hltPixelVerticesCPU'): + del process.hltPixelVerticesCPU + if hasattr(process, 'hltPixelVerticesCPUOnly'): + del process.hltPixelVerticesCPUOnly + if hasattr(process, 'hltPixelVerticesFromGPU'): + del process.hltPixelVerticesFromGPU + if hasattr(process, 'hltPixelVerticesGPU'): + del process.hltPixelVerticesGPU + + ## failsafe for fake menus + if hasattr(process, 'hltTrimmedPixelVertices'): + process.HLTRecopixelvertexingTask = cms.ConditionalTask( + process.HLTRecoPixelTracksTask, + process.hltPixelVerticesSoA, + process.hltPixelVertices, + process.hltTrimmedPixelVertices + ) + process.HLTRecopixelvertexingSequence = cms.Sequence( process.HLTRecopixelvertexingTask ) + if hasattr(process, 'hltTrimmedPixelVerticesSerialSync'): + process.HLTRecopixelvertexingCPUSerialTask = cms.ConditionalTask( + process.HLTRecoPixelTracksCPUSerialTask, + process.hltPixelVerticesSoASerialSync, + process.hltPixelVerticesSerialSync, + process.hltTrimmedPixelVerticesSerialSync + ) + process.HLTRecopixelvertexingSequenceSerialSync = cms.Sequence( process.HLTRecopixelvertexingCPUSerialTask ) + + return process + +def customizeHLTforAlpakaPixelReco(process): + '''Customisation to introduce the Pixel Local+Track+Vertex Reconstruction in Alpaka + ''' + process = customizeHLTforAlpakaPixelRecoLocal(process) + process = customizeHLTforAlpakaPixelRecoTracking(process) + process = customizeHLTforAlpakaPixelRecoVertexing(process) + + return process + + +def customizeHLTforAlpakaStripNoDoubletRecovery(process): + process = customizeHLTforAlpakaPixelReco(process) + + return process diff --git a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py index dd28b63d73517..56edea241d72f 100644 --- a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py +++ b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py @@ -126,12 +126,43 @@ def customizeHLTfor47079(process): delattr(prod,'HFStripFilter') return process +def configureFrameSoAESProducers(process): + """ + Configures the appropriate FrameSoAESProducer based on the pixel topology (Phase1, Phase2, etc.). + If the corresponding producer is not found, it will add it to the process. + """ + + # Define a mapping of pixel topology to corresponding ESProducer and component name + topology_to_es_producer = { + 'Phase1': ('frameSoAESProducerPhase1', 'FrameSoAPhase1', 'FrameSoAESProducerPhase1@alpaka'), + 'HIonPhase1': ('frameSoAESProducerHIonPhase1', 'FrameSoAPhase1HIonPhase1', 'FrameSoAESProducerHIonPhase1@alpaka'), + 'Phase2': ('frameSoAESProducerPhase2', 'FrameSoAPhase2', 'FrameSoAESProducerPhase2@alpaka'), + 'Phase1Strip': ('frameSoAESProducerPhase1Strip', 'FrameSoAPhase1Strip', 'FrameSoAESProducerPhase1Strip@alpaka'), + } + + has_alpaka_named_module = any("Alpaka" in name for name in process.__dict__.keys()) + + if not has_alpaka_named_module: + print("No modules with 'Alpaka' in their names found in the process. Skipping configuration of FrameSoAESProducers.") + return process + for pixel_topology, (es_name, component_name, producer_type) in topology_to_es_producer.items(): + if not hasattr(process, es_name): + # If the producer does not exist, create and add it + #print(f"Adding {es_name} with component name {component_name}") + setattr(process, es_name, cms.ESProducer(producer_type, + ComponentName=cms.string(component_name), + appendToDataLabel=cms.string(''))) + else: + print(f"{es_name} already configured.") + + return process # CMSSW version specific customizations def customizeHLTforCMSSW(process, menuType="GRun"): process = customiseForOffline(process) - + process = configureFrameSoAESProducers(process) + # add call to action function in proper order: newest last! # process = customiseFor12718(process) diff --git a/RecoLocalTracker/ClusterParameterEstimator/BuildFile.xml b/RecoLocalTracker/ClusterParameterEstimator/BuildFile.xml index 5ed83e90681fe..aeb932d9229ae 100644 --- a/RecoLocalTracker/ClusterParameterEstimator/BuildFile.xml +++ b/RecoLocalTracker/ClusterParameterEstimator/BuildFile.xml @@ -1,6 +1,13 @@ + + + + + + + diff --git a/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoADevice.h b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoADevice.h new file mode 100644 index 0000000000000..f513d0ef7746b --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoADevice.h @@ -0,0 +1,10 @@ +#ifndef RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoADevice_h +#define RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoADevice_h + +#include "DataFormats/Portable/interface/PortableDeviceCollection.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h" + +template +using FrameSoADevice = PortableDeviceCollection; + +#endif // RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoADevice_h diff --git a/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h new file mode 100644 index 0000000000000..02bc01cc33867 --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h @@ -0,0 +1,10 @@ +#ifndef RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoAHost_h +#define RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoAHost_h + +#include "DataFormats/Portable/interface/PortableHostCollection.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h" +#include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" + +using FrameSoAHost = PortableHostCollection; + +#endif // RecoLocalTracker_ClusterParameterEstimator_interface_FrameSoAHost_h diff --git a/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h new file mode 100644 index 0000000000000..26fc7baaa142a --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h @@ -0,0 +1,18 @@ +#ifndef RecoLocalTracker_ClusterParameterEstimator_interface_FrameLayout_h +#define RecoLocalTracker_ClusterParameterEstimator_interface_FrameLayout_h + +#include +#include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/memory.h" +#include "DataFormats/GeometrySurface/interface/SOARotation.h" + +#include "DataFormats/SoATemplate/interface/SoALayout.h" + +GENERATE_SOA_LAYOUT(FrameLayout, SOA_COLUMN(SOAFrame, detFrame)) + +using FrameSoALayout = FrameLayout<>; +using FrameSoAView = FrameSoALayout::View; +using FrameSoAConstView = FrameSoALayout::ConstView; + +#endif //RecoLocalTracker_ClusterParameterEstimator_interface_FrameLayout_h diff --git a/RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h b/RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h new file mode 100644 index 0000000000000..ef79473dd3a0e --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h @@ -0,0 +1,36 @@ +#ifndef RecoLocalTracker_ClusterParameterEstimator_interface_alpaka_FrameSoACollection_h +#define RecoLocalTracker_ClusterParameterEstimator_interface_alpaka_FrameSoACollection_h + +#include + +#include + +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoADevice.h" +#include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + +#ifdef ALPAKA_ACC_CPU_B_SEQ_T_SEQ_ENABLED + using FrameSoACollection = FrameSoAHost; +#else + using FrameSoACollection = FrameSoADevice; +#endif + +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +namespace cms::alpakatools { + template <> + struct CopyToDevice { + template + static auto copyAsync(TQueue& queue, FrameSoAHost const& srcData) { + using TDevice = typename alpaka::trait::DevType::type; + FrameSoADevice dstData(srcData->metadata().size(), queue); + alpaka::memcpy(queue, dstData.buffer(), srcData.buffer()); + return dstData; + } + }; +} // namespace cms::alpakatools + +#endif // RecoLocalTracker_ClusterParameterEstimator_interface_alpaka_FrameSoACollection_h diff --git a/RecoLocalTracker/ClusterParameterEstimator/plugins/BuildFile.xml b/RecoLocalTracker/ClusterParameterEstimator/plugins/BuildFile.xml new file mode 100644 index 0000000000000..e5e20cfbcf455 --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/plugins/BuildFile.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/RecoLocalTracker/ClusterParameterEstimator/plugins/alpaka/FrameSoAESProducer.cc b/RecoLocalTracker/ClusterParameterEstimator/plugins/alpaka/FrameSoAESProducer.cc new file mode 100644 index 0000000000000..426fabb143b64 --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/plugins/alpaka/FrameSoAESProducer.cc @@ -0,0 +1,118 @@ +#include +#include +#include +#include + +#include "DataFormats/GeometrySurface/interface/SOARotation.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "Geometry/CommonTopologies/interface/TrackerGeomDet.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" + +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" +#include "Geometry/CommonTopologies/interface/PixelTopology.h" +#include "Geometry/CommonTopologies/interface/Topology.h" + +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESProducer.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/EventSetup.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/ModuleFactory.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/memory.h" + +#include "RecoLocalTracker/Records/interface/FrameSoARecord.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + template + class FrameSoAESProducer : public ESProducer { + using Rotation = SOARotation; + using Frame = SOAFrame; + + public: + FrameSoAESProducer(edm::ParameterSet const& iConfig); + std::unique_ptr produce(const FrameSoARecord& iRecord); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + edm::ESGetToken geometry_; + edm::ESGetToken topology_; + }; + + using namespace edm; + + template + FrameSoAESProducer::FrameSoAESProducer(const edm::ParameterSet& p) : ESProducer(p) { + auto const& myname = p.getParameter("ComponentName"); + + auto cc = setWhatProduced(this, myname); + geometry_ = cc.consumes(); + topology_ = cc.consumes(); + } + + template + std::unique_ptr FrameSoAESProducer::produce(const FrameSoARecord& iRecord) { + const TrackerGeometry* geometry = &iRecord.get(geometry_); + const TrackerTopology* topology = &iRecord.get(topology_); + + auto const& detUnits = geometry->detUnits(); + + auto product = std::make_unique(TrackerTraits::numberOfModules, cms::alpakatools::host()); + + if constexpr (std::is_same_v) { + int i = 0; + for (auto layer : phase1PixelStripTopology::layerData) { + auto step = layer.isStrip2D ? 2 : 1; + for (auto j = layer.start; j != layer.end; j += step) { + auto& s = layer.isStrip2D ? geometry->idToDet(topology->glued(detUnits[i]->geographicalId()))->surface() + : detUnits[j]->surface(); + product->view()[i].detFrame() = Frame(s.position().x(), s.position().y(), s.position().z(), s.rotation()); + ++i; + } + } + } else { + constexpr auto n_detectors = + TrackerTraits::numberOfModules; // converting only up to the modules used in the CA topology + + assert(n_detectors < + detUnits.size()); //still there shouldn't be more modules than what we have from the TrackerGeometry + + for (unsigned i = 0; i != n_detectors; ++i) { + auto det = detUnits[i]; + auto vv = det->surface().position(); + auto rr = Rotation(det->surface().rotation()); + product->view()[i].detFrame() = Frame(vv.x(), vv.y(), vv.z(), rr); + } + } + + return product; + } + + template + void FrameSoAESProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + + std::string name = "FrameSoAPhase1"; + name += TrackerTraits::nameModifier; + desc.add("ComponentName", name); + + descriptions.addWithDefaultLabel(desc); + } + + using FrameSoAESProducerPhase1 = FrameSoAESProducer; + using FrameSoAESProducerPhase2 = FrameSoAESProducer; + using FrameSoAESProducerHIonPhase1 = FrameSoAESProducer; + using FrameSoAESProducerPhase1Strip = FrameSoAESProducer; +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(FrameSoAESProducerPhase1); +DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(FrameSoAESProducerPhase1Strip); +DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(FrameSoAESProducerPhase2); +DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(FrameSoAESProducerHIonPhase1); diff --git a/RecoLocalTracker/ClusterParameterEstimator/src/ES_FrameSoA.cc b/RecoLocalTracker/ClusterParameterEstimator/src/ES_FrameSoA.cc new file mode 100644 index 0000000000000..7ec90366be1e5 --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/src/ES_FrameSoA.cc @@ -0,0 +1,4 @@ +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoAHost.h" +#include "FWCore/Utilities/interface/typelookup.h" + +TYPELOOKUP_DATA_REG(FrameSoAHost); diff --git a/RecoLocalTracker/ClusterParameterEstimator/src/alpaka/ES_FrameSoA.cc b/RecoLocalTracker/ClusterParameterEstimator/src/alpaka/ES_FrameSoA.cc new file mode 100644 index 0000000000000..3293d432f258a --- /dev/null +++ b/RecoLocalTracker/ClusterParameterEstimator/src/alpaka/ES_FrameSoA.cc @@ -0,0 +1,4 @@ +#include "RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/typelookup.h" + +TYPELOOKUP_ALPAKA_DATA_REG(FrameSoACollection); diff --git a/RecoLocalTracker/Records/BuildFile.xml b/RecoLocalTracker/Records/BuildFile.xml index 629b1aa7a1ebc..fae441b2f4ac2 100644 --- a/RecoLocalTracker/Records/BuildFile.xml +++ b/RecoLocalTracker/Records/BuildFile.xml @@ -8,3 +8,4 @@ + diff --git a/RecoLocalTracker/Records/interface/FrameSoARecord.h b/RecoLocalTracker/Records/interface/FrameSoARecord.h new file mode 100644 index 0000000000000..5c149a55ed149 --- /dev/null +++ b/RecoLocalTracker/Records/interface/FrameSoARecord.h @@ -0,0 +1,13 @@ +#ifndef RecoLocalTracker_Records_FrameSoARecord_h +#define RecoLocalTracker_Records_FrameSoARecord_h +#include "FWCore/Framework/interface/EventSetupRecordImplementation.h" +#include "FWCore/Framework/interface/DependentRecordImplementation.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include "FWCore/Utilities/interface/mplVector.h" + +class FrameSoARecord : public edm::eventsetup::DependentRecordImplementation< + FrameSoARecord, + edm::mpl::Vector > {}; +#endif diff --git a/RecoLocalTracker/Records/src/FrameSoARecord.cc b/RecoLocalTracker/Records/src/FrameSoARecord.cc new file mode 100644 index 0000000000000..75f8e1e13469a --- /dev/null +++ b/RecoLocalTracker/Records/src/FrameSoARecord.cc @@ -0,0 +1,5 @@ +#include "RecoLocalTracker/Records/interface/FrameSoARecord.h" +#include "FWCore/Framework/interface/eventsetuprecord_registration_macro.h" +#include "FWCore/Utilities/interface/typelookup.h" + +EVENTSETUP_RECORD_REG(FrameSoARecord); diff --git a/RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h b/RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h index 1338bc457c1f7..b74ddab614ae1 100644 --- a/RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h +++ b/RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h @@ -8,6 +8,7 @@ #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFastParamsDevice.h" #include "DataFormats/Portable/interface/alpaka/PortableCollection.h" #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" // TODO: The class is created via inheritance of the PortableCollection. // This is generally discouraged, and should be done via composition. diff --git a/RecoLocalTracker/SiPixelRecHits/python/PixelCPEESProducers_cff.py b/RecoLocalTracker/SiPixelRecHits/python/PixelCPEESProducers_cff.py index d0532d33ff260..db899124bb8e3 100644 --- a/RecoLocalTracker/SiPixelRecHits/python/PixelCPEESProducers_cff.py +++ b/RecoLocalTracker/SiPixelRecHits/python/PixelCPEESProducers_cff.py @@ -17,10 +17,14 @@ from CalibTracker.SiPixelESProducers.SiPixelTemplateDBObjectESProducer_cfi import * from CalibTracker.SiPixelESProducers.SiPixel2DTemplateDBObjectESProducer_cfi import * +# Alpaka specic def _addProcessCPEsAlpaka(process): process.load("RecoLocalTracker.SiPixelRecHits.pixelCPEFastParamsESProducerAlpakaPhase1_cfi") process.load("RecoLocalTracker.SiPixelRecHits.pixelCPEFastParamsESProducerAlpakaPhase2_cfi") process.load("RecoLocalTracker.SiPixelRecHits.pixelCPEFastParamsESProducerAlpakaHIonPhase1_cfi") + process.load("RecoLocalTracker.ClusterParameterEstimator.frameSoAESProducerPhase1Strip_cfi") + process.load("RecoLocalTracker.ClusterParameterEstimator.frameSoAESProducerPhase1_cfi") + process.load("RecoLocalTracker.ClusterParameterEstimator.frameSoAESProducerPhase2_cfi") modifyConfigurationForAlpakaCPEs_ = alpaka.makeProcessModifier(_addProcessCPEsAlpaka) diff --git a/RecoLocalTracker/SiPixelRecHits/src/ES_PixelCPEFastParams.cc b/RecoLocalTracker/SiPixelRecHits/src/ES_PixelCPEFastParams.cc index 2515e71b44852..4aa24bd6f5171 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/ES_PixelCPEFastParams.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/ES_PixelCPEFastParams.cc @@ -1,6 +1,6 @@ #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFastParamsHost.h" #include "FWCore/Utilities/interface/typelookup.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" using PixelCPEFastParamsHostPhase1 = PixelCPEFastParamsHost; using PixelCPEFastParamsHostHIonPhase1 = PixelCPEFastParamsHost; diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc index 8ab0c265a866f..e15721cda2a91 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFast.cc @@ -4,7 +4,8 @@ #include "DataFormats/DetId/interface/DetId.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h" #include "MagneticField/Engine/interface/MagneticField.h" #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFast.h" @@ -529,3 +530,4 @@ void PixelCPEFast::fillPSetDescription(edm::ParameterSetDescripti template class PixelCPEFast; template class PixelCPEFast; template class PixelCPEFast; +template class PixelCPEFast; diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParams.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParams.cc index ea42202746b72..8337853ca386f 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParams.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParams.cc @@ -1,6 +1,6 @@ #include "FWCore/Utilities/interface/typelookup.h" #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFastParamsDevice.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" using PixelCPEFastParamsPhase1 = PixelCPEFastParamsDevice; using PixelCPEFastParamsHIonPhase1 = PixelCPEFastParamsDevice; diff --git a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParamsHost.cc b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParamsHost.cc index b290aabc194d1..33aeccfeb15cb 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParamsHost.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/PixelCPEFastParamsHost.cc @@ -4,10 +4,10 @@ #include "DataFormats/GeometrySurface/interface/SOARotation.h" #include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h" #include "DataFormats/TrackingRecHitSoA/interface/SiPixelHitStatus.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" #include "HeterogeneousCore/AlpakaInterface/interface/memory.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEFastParamsHost.h" //----------------------------------------------------------------------------- diff --git a/RecoLocalTracker/SiStripRecHitConverter/BuildFile.xml b/RecoLocalTracker/SiStripRecHitConverter/BuildFile.xml index 0bbf434606371..fae5ce8a8f5ce 100644 --- a/RecoLocalTracker/SiStripRecHitConverter/BuildFile.xml +++ b/RecoLocalTracker/SiStripRecHitConverter/BuildFile.xml @@ -3,6 +3,7 @@ + @@ -14,6 +15,7 @@ + diff --git a/RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h b/RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h index d88ea518d35ce..2ad7047ba4961 100644 --- a/RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h +++ b/RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h @@ -2,7 +2,8 @@ #define RecoLocalTracker_SiStripRecHitConverter_StripCPEfromTrackAngle_H #include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPE.h" - +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" class StripCPEfromTrackAngle : public StripCPE { private: using StripCPE::localParameters; diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/BuildFile.xml b/RecoLocalTracker/SiStripRecHitConverter/plugins/BuildFile.xml index df0f7979aa74f..5dd7d1d7921cf 100644 --- a/RecoLocalTracker/SiStripRecHitConverter/plugins/BuildFile.xml +++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/BuildFile.xml @@ -1,5 +1,14 @@ + + + + + + + + + diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoA.cc b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoA.cc new file mode 100644 index 0000000000000..4eb7c5a53ebcf --- /dev/null +++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoA.cc @@ -0,0 +1,258 @@ + +#include +#include +#include +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/BeamSpot/interface/BeamSpotPOD.h" +#include "DataFormats/BeamSpot/interface/alpaka/BeamSpotDevice.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/SiPixelDigiSoA/interface/SiPixelDigisDevice.h" +#include "DataFormats/SiPixelDigiSoA/interface/alpaka/SiPixelDigisSoACollection.h" +#include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" +#include "DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/Event.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/EventSetup.h" +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/stream/SynchronizingEDProducer.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/memory.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoLocalTracker/Records/interface/PixelCPEFastParamsRecord.h" + +#include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h" +#include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforDevice.h" +#include "RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h" + +// #include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h" +#include "DataFormats/Math/interface/approx_atan2.h" + +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/GluedGeomDet.h" + +#include "SiStripRecHitSoAKernel.h" +#include "alpaka/mem/view/Traits.hpp" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + + template + class SiStripRecHitSoA : public stream::SynchronizingEDProducer<> { + using PixelBase = typename TrackerTraits::PixelBase; + + using StripHits = TrackingRecHitsSoACollection; + using PixelHits = TrackingRecHitsSoACollection; + + using StripHitsHost = TrackingRecHitHost; + using PixelHitsHost = TrackingRecHitHost; + + using Algo = hitkernels::SiStripRecHitSoAKernel; + + public: + explicit SiStripRecHitSoA(const edm::ParameterSet& iConfig); + ~SiStripRecHitSoA() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + void acquire(device::Event const& iEvent, device::EventSetup const& iSetup) override {}; + void produce(device::Event& iEvent, device::EventSetup const& iSetup) override; + + const edm::ESGetToken geomToken_; + const edm::EDGetTokenT recHitToken_; + const edm::EDGetTokenT beamSpotToken_; + const edm::EDGetTokenT pixelRecHitSoAToken_; + + const device::EDPutToken stripSoA_; + const edm::EDPutTokenT> hmsToken_; + + const Algo Algo_; + }; + + template + SiStripRecHitSoA::SiStripRecHitSoA(const edm::ParameterSet& iConfig) + : geomToken_(esConsumes()), + recHitToken_{consumes(iConfig.getParameter("stripRecHitSource"))}, + //beamSpotToken(consumes(edm::InputTag("offlineBeamSpot"))), + beamSpotToken_(consumes(iConfig.getParameter("beamSpot"))), + pixelRecHitSoAToken_{consumes(iConfig.getParameter("pixelRecHitSoASource"))}, + stripSoA_{produces()}, + hmsToken_{produces()} {} + + template + void SiStripRecHitSoA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + + desc.add("stripRecHitSource", edm::InputTag("siStripMatchedRecHits", "matchedRecHit")); + desc.add("beamSpot", edm::InputTag("offlineBeamSpot")); + desc.add("pixelRecHitSoASource", edm::InputTag("siPixelRecHitsPreSplittingAlpaka")); + descriptions.addWithDefaultLabel(desc); + + // desc.setUnknown(); + // descriptions.addDefault(desc); + } + + template + void SiStripRecHitSoA::produce(device::Event& iEvent, device::EventSetup const& iSetup) { + // Get the objects that we need + const TrackerGeometry* trackerGeometry = &iSetup.getData(geomToken_); + auto const& stripHits = iEvent.get(recHitToken_); + auto const& pixelHitsHost = iEvent.get(pixelRecHitSoAToken_); + auto& bs = iEvent.get(beamSpotToken_); + + // Count strip hits + size_t nStripHits = 0; + //std::cout << "number of modules: " << TrackerTraits::numberOfModules << std::endl; + //std::cout << "stripHits size: " << stripHits.size() << std::endl; + for (const auto& detSet : stripHits) { + const GluedGeomDet* det = static_cast(trackerGeometry->idToDet(detSet.detId())); + //std::cout << "detSet.detId()" << detSet.detId() << std::endl; + //std::cout << "det->stereoDet()->index()" << det->stereoDet()->index() << std::endl; + if (TrackerTraits::mapIndex(det->stereoDet()->index()) < TrackerTraits::numberOfModules) + nStripHits += detSet.size(); + } + + size_t nPixelHits = pixelHitsHost.view().metadata().size(); + + //std::cout << "nStripHits = " << nStripHits << std::endl; + //std::cout << "nPixelHits = " << nPixelHits << std::endl; + + // HostView pixelHitsHostView(pixelHits, iEvent.queue()); + // PixelHitsHost& pixelHitsHost = pixelHitsHostView.get(); + // PixelHitsHost pixelHitsHost(nPixelHits, iEvent.queue()); + + // alpaka::memcpy(iEvent.queue(), pixelHitsHost.buffer(), pixelHits.buffer()); + + StripHitsHost allHitsHost(iEvent.queue(), nPixelHits + nStripHits); + + // Copy pixel data + std::copy(pixelHitsHost.view().xLocal(), pixelHitsHost.view().xLocal() + nPixelHits, allHitsHost.view().xLocal()); + std::copy(pixelHitsHost.view().yLocal(), pixelHitsHost.view().yLocal() + nPixelHits, allHitsHost.view().yLocal()); + std::copy(pixelHitsHost.view().xerrLocal(), + pixelHitsHost.view().xerrLocal() + nPixelHits, + allHitsHost.view().xerrLocal()); + std::copy(pixelHitsHost.view().yerrLocal(), + pixelHitsHost.view().yerrLocal() + nPixelHits, + allHitsHost.view().yerrLocal()); + std::copy( + pixelHitsHost.view().xGlobal(), pixelHitsHost.view().xGlobal() + nPixelHits, allHitsHost.view().xGlobal()); + std::copy( + pixelHitsHost.view().yGlobal(), pixelHitsHost.view().yGlobal() + nPixelHits, allHitsHost.view().yGlobal()); + std::copy( + pixelHitsHost.view().zGlobal(), pixelHitsHost.view().zGlobal() + nPixelHits, allHitsHost.view().zGlobal()); + std::copy( + pixelHitsHost.view().rGlobal(), pixelHitsHost.view().rGlobal() + nPixelHits, allHitsHost.view().rGlobal()); + std::copy(pixelHitsHost.view().iphi(), pixelHitsHost.view().iphi() + nPixelHits, allHitsHost.view().iphi()); + std::copy(pixelHitsHost.view().chargeAndStatus(), + pixelHitsHost.view().chargeAndStatus() + nPixelHits, + allHitsHost.view().chargeAndStatus()); + std::copy(pixelHitsHost.view().clusterSizeX(), + pixelHitsHost.view().clusterSizeX() + nPixelHits, + allHitsHost.view().clusterSizeX()); + std::copy(pixelHitsHost.view().clusterSizeY(), + pixelHitsHost.view().clusterSizeY() + nPixelHits, + allHitsHost.view().clusterSizeY()); + std::copy(pixelHitsHost.view().detectorIndex(), + pixelHitsHost.view().detectorIndex() + nPixelHits, + allHitsHost.view().detectorIndex()); + + std::copy(pixelHitsHost.view().phiBinnerStorage(), + pixelHitsHost.view().phiBinnerStorage() + nPixelHits, + allHitsHost.view().phiBinnerStorage()); + + allHitsHost.view().offsetBPIX2() = pixelHitsHost.view().offsetBPIX2(); + + auto& hitsModuleStart = allHitsHost.view().hitsModuleStart(); + + std::copy(pixelHitsHost.view().hitsModuleStart().begin(), + pixelHitsHost.view().hitsModuleStart().end(), + hitsModuleStart.begin()); + + std::map> mappedModuleHits; + + // Loop over strip RecHits + for (auto detSet : stripHits) { + const GluedGeomDet* det = static_cast(trackerGeometry->idToDet(detSet.detId())); + size_t index = TrackerTraits::mapIndex(det->stereoDet()->index()); + + if (index >= TrackerTraits::numberOfModules) + continue; + + mappedModuleHits[index] = detSet; + } + + size_t i = 0; + size_t lastIndex = TrackerTraits::numberOfPixelModules; + + for (auto& [index, detSet] : mappedModuleHits) { + const GluedGeomDet* det = static_cast(trackerGeometry->idToDet(detSet.detId())); + + // no hits since lastIndex: hitsModuleStart[lastIndex:index] = hitsModuleStart[lastIndex] + for (auto j = lastIndex + 1; j < index + 1; ++j) + hitsModuleStart[j] = hitsModuleStart[lastIndex]; + + hitsModuleStart[index + 1] = hitsModuleStart[index] + detSet.size(); + lastIndex = index + 1; + + for (const auto& recHit : detSet) { + allHitsHost.view()[nPixelHits + i].xLocal() = recHit.localPosition().x(); + allHitsHost.view()[nPixelHits + i].yLocal() = recHit.localPosition().y(); + allHitsHost.view()[nPixelHits + i].xerrLocal() = recHit.localPositionError().xx(); + allHitsHost.view()[nPixelHits + i].yerrLocal() = recHit.localPositionError().yy(); + auto globalPosition = det->toGlobal(recHit.localPosition()); + double gx = globalPosition.x() - bs.x0(); + double gy = globalPosition.y() - bs.y0(); + double gz = globalPosition.z() - bs.z0(); + allHitsHost.view()[nPixelHits + i].xGlobal() = gx; + allHitsHost.view()[nPixelHits + i].yGlobal() = gy; + allHitsHost.view()[nPixelHits + i].zGlobal() = gz; + allHitsHost.view()[nPixelHits + i].rGlobal() = sqrt(gx * gx + gy * gy); + allHitsHost.view()[nPixelHits + i].iphi() = unsafe_atan2s<7>(gy, gx); + // allHitsHost.view()[nPixelHits + i].chargeAndStatus().charge = ? + // allHitsHost.view()[nPixelHits + i].chargeAndStatus().status = ? + // allHitsHost.view()[nPixelHits + i].clusterSizeX() = ? + // allHitsHost.view()[nPixelHits + i].clusterSizeY() = ? + allHitsHost.view()[nPixelHits + i].detectorIndex() = index; + // ??? + ++i; + } + } + + for (auto j = lastIndex + 1; j < TrackerTraits::numberOfModules + 1; ++j) + hitsModuleStart[j] = hitsModuleStart[lastIndex]; + + for (auto layer = 0U; layer < TrackerTraits::numberOfLayers + 1; ++layer) { + allHitsHost.view().hitsLayerStart()[layer] = hitsModuleStart[TrackerTraits::layerStart[layer]]; + } + + iEvent.emplace(hmsToken_, std::vector(hitsModuleStart.begin(), hitsModuleStart.end())); + + iEvent.emplace(stripSoA_, Algo_.fillHitsAsync(allHitsHost, iEvent.queue())); + + //std::cout << "produce done" << std::endl; + } + using SiStripRecHitSoAPhase1 = SiStripRecHitSoA; +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h" +DEFINE_FWK_ALPAKA_MODULE(SiStripRecHitSoAPhase1); + +// using SiPixelRecHitSoAFromLegacyPhase2 = SiStripRecHitSoA; +// DEFINE_FWK_MODULE(SiPixelRecHitSoAFromLegacyPhase2); diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.dev.cc b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.dev.cc new file mode 100644 index 0000000000000..153dd1a28ca0c --- /dev/null +++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.dev.cc @@ -0,0 +1,71 @@ +// C++ headers +#include +#include + +// Alpaka headers +#include + +// CMSSW headers +#include "HeterogeneousCore/AlpakaInterface/interface/HistoContainer.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" + +#include "SiStripRecHitSoAKernel.h" + +//#define GPU_DEBUG + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + using namespace cms::alpakatools; + + namespace hitkernels { + + template + TrackingRecHitsSoACollection SiStripRecHitSoAKernel::fillHitsAsync( + StripHitsHost const& hits_h, Queue queue) const { + TrackingRecHitsSoACollection hits_d( + queue, hits_h.nHits(), hits_h.offsetBPIX2(), hits_h.hitsModuleStart()); + alpaka::memcpy(queue, hits_d.buffer(), hits_h.buffer()); + + // assuming full warp of threads is better than a smaller number... + if (hits_h.nHits()) { + constexpr auto nLayers = TrackerTraits::numberOfLayers; + + typename TrackingRecHitSoA::PhiBinnerView hrv_d; + hrv_d.assoc = &(hits_d.view().phiBinner()); + hrv_d.offSize = -1; + hrv_d.offStorage = nullptr; + hrv_d.contentSize = hits_h.nHits(); + hrv_d.contentStorage = hits_d.view().phiBinnerStorage(); + + // cms::alpakatools::fillManyFromVector(&(hits_d.view().phiBinner()), + // nLayers, + // hits_d.view().iphi(), + // hits_d.view().hitsLayerStart().data(), + // hits_h.nHits(), + // (uint32_t)256, + // queue); + cms::alpakatools::fillManyFromVector(&(hits_d.view().phiBinner()), + hrv_d, + nLayers, + hits_d.view().iphi(), + hits_d.view().hitsLayerStart().data(), + hits_h.nHits(), + (uint32_t)256, + queue); + +#ifdef GPU_DEBUG + alpaka::wait(queue); +#endif + } + +#ifdef GPU_DEBUG + alpaka::wait(queue); + std::cout << "SiStripRecHitSoAKernel -> DONE!" << std::endl; +#endif + + return hits_d; + } + + template class SiStripRecHitSoAKernel; + + } // namespace hitkernels +} // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.h b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.h new file mode 100644 index 0000000000000..4dcba16c59d1d --- /dev/null +++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/alpaka/SiStripRecHitSoAKernel.h @@ -0,0 +1,42 @@ +#ifndef RecoLocalTracker_SiPixelRecHits_SiStripRecHitSoAKernel_h +#define RecoLocalTracker_SiPixelRecHits_SiStripRecHitSoAKernel_h + +#include + +#include + +#include "DataFormats/BeamSpot/interface/BeamSpotPOD.h" +#include "DataFormats/SiPixelClusterSoA/interface/alpaka/SiPixelClustersSoACollection.h" +#include "DataFormats/SiPixelClusterSoA/interface/SiPixelClustersDevice.h" +#include "DataFormats/SiPixelDigiSoA/interface/SiPixelDigisDevice.h" +#include "DataFormats/SiPixelDigiSoA/interface/alpaka/SiPixelDigisSoACollection.h" +#include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" +#include "DataFormats/TrackingRecHitSoA/interface/alpaka/TrackingRecHitsSoACollection.h" +#include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" +#include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforDevice.h" + +namespace ALPAKA_ACCELERATOR_NAMESPACE { + namespace hitkernels { + using namespace cms::alpakatools; + + template + class SiStripRecHitSoAKernel { + using StripHits = TrackingRecHitsSoACollection; + using StripHitsHost = TrackingRecHitHost; + + public: + SiStripRecHitSoAKernel() = default; + ~SiStripRecHitSoAKernel() = default; + + SiStripRecHitSoAKernel(const SiStripRecHitSoAKernel&) = delete; + SiStripRecHitSoAKernel(SiStripRecHitSoAKernel&&) = delete; + SiStripRecHitSoAKernel& operator=(const SiStripRecHitSoAKernel&) = delete; + SiStripRecHitSoAKernel& operator=(SiStripRecHitSoAKernel&&) = delete; + + StripHits fillHitsAsync(StripHitsHost const& hits_h, Queue queue) const; + }; + } // namespace hitkernels +} // namespace ALPAKA_ACCELERATOR_NAMESPACE + +#endif // RecoLocalTracker_SiPixelRecHits_SiStripRecHitSoAKernel_h diff --git a/RecoTauTag/HLTProducers/src/L2TauTagNNProducerAlpakaStrip.cc b/RecoTauTag/HLTProducers/src/L2TauTagNNProducerAlpakaStrip.cc new file mode 100644 index 0000000000000..7e01585998c9d --- /dev/null +++ b/RecoTauTag/HLTProducers/src/L2TauTagNNProducerAlpakaStrip.cc @@ -0,0 +1,830 @@ +/* + * \class L2TauTagProducer + * + * L2Tau identification using Convolutional NN. + * + * \author Valeria D'Amante, Università di Siena and INFN Pisa + * Konstantin Androsov, EPFL and ETHZ +*/ +#include +#include +#include +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/isFinite.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloTopology/interface/HcalTopology.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "DataFormats/CaloRecHit/interface/CaloRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h" +#include "DataFormats/HcalDetId/interface/HcalDetId.h" +#include "DataFormats/HcalRecHit/interface/HBHERecHit.h" +#include "DataFormats/HcalRecHit/interface/HcalRecHitDefs.h" +#include "DataFormats/HcalRecHit/interface/HFRecHit.h" +#include "DataFormats/HcalRecHit/interface/HORecHit.h" +#include "DataFormats/HLTReco/interface/TriggerTypeDefs.h" +#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h" +#include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h" +#include "RecoTracker/PixelTrackFitting/interface/FitUtils.h" +#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h" +#include "DataFormats/TrackReco/interface/HitPattern.h" +#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h" +#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h" +#include "DataFormats/GeometrySurface/interface/Plane.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h" + +#include "DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h" +#include "DataFormats/TrackSoA/interface/TracksHost.h" +#include "DataFormats/VertexSoA/interface/ZVertexHost.h" + +namespace L2TauTagNNv1 { + constexpr int nCellEta = 5; + constexpr int nCellPhi = 5; + constexpr int nVars = 31; + constexpr float dR_max = 0.5; + enum class NNInputs { + nVertices = 0, + l1Tau_pt, + l1Tau_eta, + l1Tau_hwIso, + EcalEnergySum, + EcalSize, + EcalEnergyStdDev, + EcalDeltaEta, + EcalDeltaPhi, + EcalChi2, + EcalEnergySumForPositiveChi2, + EcalSizeForPositiveChi2, + HcalEnergySum, + HcalSize, + HcalEnergyStdDev, + HcalDeltaEta, + HcalDeltaPhi, + HcalChi2, + HcalEnergySumForPositiveChi2, + HcalSizeForPositiveChi2, + PatatrackPtSum, + PatatrackSize, + PatatrackSizeWithVertex, + PatatrackPtSumWithVertex, + PatatrackChargeSum, + PatatrackDeltaEta, + PatatrackDeltaPhi, + PatatrackChi2OverNdof, + PatatrackNdof, + PatatrackDxy, + PatatrackDz + }; + + const std::map varNameMap = { + {NNInputs::nVertices, "nVertices"}, + {NNInputs::l1Tau_pt, "l1Tau_pt"}, + {NNInputs::l1Tau_eta, "l1Tau_eta"}, + {NNInputs::l1Tau_hwIso, "l1Tau_hwIso"}, + {NNInputs::EcalEnergySum, "EcalEnergySum"}, + {NNInputs::EcalSize, "EcalSize"}, + {NNInputs::EcalEnergyStdDev, "EcalEnergyStdDev"}, + {NNInputs::EcalDeltaEta, "EcalDeltaEta"}, + {NNInputs::EcalDeltaPhi, "EcalDeltaPhi"}, + {NNInputs::EcalChi2, "EcalChi2"}, + {NNInputs::EcalEnergySumForPositiveChi2, "EcalEnergySumForPositiveChi2"}, + {NNInputs::EcalSizeForPositiveChi2, "EcalSizeForPositiveChi2"}, + {NNInputs::HcalEnergySum, "HcalEnergySum"}, + {NNInputs::HcalSize, "HcalSize"}, + {NNInputs::HcalEnergyStdDev, "HcalEnergyStdDev"}, + {NNInputs::HcalDeltaEta, "HcalDeltaEta"}, + {NNInputs::HcalDeltaPhi, "HcalDeltaPhi"}, + {NNInputs::HcalChi2, "HcalChi2"}, + {NNInputs::HcalEnergySumForPositiveChi2, "HcalEnergySumForPositiveChi2"}, + {NNInputs::HcalSizeForPositiveChi2, "HcalSizeForPositiveChi2"}, + {NNInputs::PatatrackPtSum, "PatatrackPtSum"}, + {NNInputs::PatatrackSize, "PatatrackSize"}, + {NNInputs::PatatrackSizeWithVertex, "PatatrackSizeWithVertex"}, + {NNInputs::PatatrackPtSumWithVertex, "PatatrackPtSumWithVertex"}, + {NNInputs::PatatrackChargeSum, "PatatrackChargeSum"}, + {NNInputs::PatatrackDeltaEta, "PatatrackDeltaEta"}, + {NNInputs::PatatrackDeltaPhi, "PatatrackDeltaPhi"}, + {NNInputs::PatatrackChi2OverNdof, "PatatrackChi2OverNdof"}, + {NNInputs::PatatrackNdof, "PatatrackNdof"}, + {NNInputs::PatatrackDxy, "PatatrackDxy"}, + {NNInputs::PatatrackDz, "PatatrackDz"}}; +} // namespace L2TauTagNNv1 +namespace { + inline float& getCellImpl( + tensorflow::Tensor& cellGridMatrix, int tau_idx, int phi_idx, int eta_idx, L2TauTagNNv1::NNInputs NNInput_idx) { + return cellGridMatrix.tensor()(tau_idx, phi_idx, eta_idx, static_cast(NNInput_idx)); + } +} // namespace +struct normDictElement { + float mean; + float std; + float min; + float max; +}; + +struct L2TauNNProducerAlpakaStripCacheData { + L2TauNNProducerAlpakaStripCacheData() : graphDef(nullptr), session(nullptr) {} + tensorflow::GraphDef* graphDef; + tensorflow::Session* session; + std::vector normVec; +}; + +class L2TauNNProducerAlpakaStrip + : public edm::stream::EDProducer> { +public: + using TracksHost = pixelTrack::TracksHostPhase1Strip; + + struct caloRecHitCollections { + const HBHERecHitCollection* hbhe; + const HORecHitCollection* ho; + const EcalRecHitCollection* eb; + const EcalRecHitCollection* ee; + const CaloGeometry* geometry; + }; + + struct InputDescTau { + std::string CollectionName; + edm::EDGetTokenT inputToken_; + }; + + static constexpr float dR2_max = L2TauTagNNv1::dR_max * L2TauTagNNv1::dR_max; + static constexpr float dEta_width = 2 * L2TauTagNNv1::dR_max / static_cast(L2TauTagNNv1::nCellEta); + static constexpr float dPhi_width = 2 * L2TauTagNNv1::dR_max / static_cast(L2TauTagNNv1::nCellPhi); + + explicit L2TauNNProducerAlpakaStrip(const edm::ParameterSet&, const L2TauNNProducerAlpakaStripCacheData*); + static void fillDescriptions(edm::ConfigurationDescriptions&); + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet&); + static void globalEndJob(L2TauNNProducerAlpakaStripCacheData*); + +private: + void checknan(tensorflow::Tensor& tensor, int debugLevel); + void standardizeTensor(tensorflow::Tensor& tensor); + std::vector getTauScore(const tensorflow::Tensor& cellGridMatrix); + void produce(edm::Event& event, const edm::EventSetup& eventsetup) override; + void fillL1TauVars(tensorflow::Tensor& cellGridMatrix, const std::vector& allTaus); + void fillCaloRecHits(tensorflow::Tensor& cellGridMatrix, + const std::vector& allTaus, + const caloRecHitCollections& caloRecHits); + void fillPatatracks(tensorflow::Tensor& cellGridMatrix, + const std::vector& allTaus, + const TracksHost& patatracks_tsoa, + const ZVertexHost& patavtx_soa, + const reco::BeamSpot& beamspot, + const MagneticField* magfi); + void selectGoodTracksAndVertices(const ZVertexHost& patavtx_soa, + const TracksHost& patatracks_tsoa, + std::vector& trkGood, + std::vector& vtxGood); + std::pair impactParameter(int it, + const TracksHost& patatracks_tsoa, + float patatrackPhi, + const reco::BeamSpot& beamspot, + const MagneticField* magfi); + template + std::tuple getEtaPhiIndices(const VPos& position, const LVec& tau_p4); + template + std::tuple getEtaPhiIndices(float eta, float phi, const LVec& tau_p4); + +private: + const int debugLevel_; + const edm::EDGetTokenT tauTriggerToken_; + std::vector L1TauDesc_; + const edm::EDGetTokenT hbheToken_; + const edm::EDGetTokenT hoToken_; + const edm::EDGetTokenT ebToken_; + const edm::EDGetTokenT eeToken_; + const edm::ESGetToken geometryToken_; + const edm::ESGetToken bFieldToken_; + const edm::EDGetTokenT pataVerticesToken_; + const edm::EDGetTokenT pataTracksToken_; + const edm::EDGetTokenT beamSpotToken_; + const unsigned int maxVtx_; + const float fractionSumPt2_; + const float minSumPt2_; + const float trackPtMin_; + const float trackPtMax_; + const float trackChi2Max_; + std::string inputTensorName_; + std::string outputTensorName_; + const L2TauNNProducerAlpakaStripCacheData* L2cacheData_; +}; + +std::unique_ptr L2TauNNProducerAlpakaStrip::initializeGlobalCache( + const edm::ParameterSet& cfg) { + std::unique_ptr cacheData = + std::make_unique(); + cacheData->normVec.reserve(L2TauTagNNv1::nVars); + + auto const graphPath = edm::FileInPath(cfg.getParameter("graphPath")).fullPath(); + + cacheData->graphDef = tensorflow::loadGraphDef(graphPath); + cacheData->session = tensorflow::createSession(cacheData->graphDef); + + boost::property_tree::ptree loadPtreeRoot; + auto const normalizationDict = edm::FileInPath(cfg.getParameter("normalizationDict")).fullPath(); + boost::property_tree::read_json(normalizationDict, loadPtreeRoot); + for (const auto& [key, val] : L2TauTagNNv1::varNameMap) { + boost::property_tree::ptree var = loadPtreeRoot.get_child(val); + normDictElement current_element; + current_element.mean = var.get_child("mean").get_value(); + current_element.std = var.get_child("std").get_value(); + current_element.min = var.get_child("min").get_value(); + current_element.max = var.get_child("max").get_value(); + cacheData->normVec.push_back(current_element); + } + return cacheData; +} +void L2TauNNProducerAlpakaStrip::globalEndJob(L2TauNNProducerAlpakaStripCacheData* cacheData) { + if (cacheData->graphDef != nullptr) { + delete cacheData->graphDef; + } + tensorflow::closeSession(cacheData->session); +} +void L2TauNNProducerAlpakaStrip::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("debugLevel", 0)->setComment("set debug level for printing out info"); + edm::ParameterSetDescription l1TausPset; + l1TausPset.add("L1CollectionName", "DoubleTau")->setComment("Name of collections"); + l1TausPset.add("L1TauTrigger", edm::InputTag("hltL1sDoubleTauBigOR")) + ->setComment("Which trigger should the L1 Taus collection pass"); + edm::ParameterSet l1TausPSetDefault; + l1TausPSetDefault.addParameter("L1CollectionName", "DoubleTau"); + l1TausPSetDefault.addParameter("L1TauTrigger", edm::InputTag("hltL1sDoubleTauBigOR")); + desc.addVPSet("L1Taus", l1TausPset, {l1TausPSetDefault}); + desc.add("hbheInput", edm::InputTag("hltHbhereco"))->setComment("HBHE recHit collection"); + desc.add("hoInput", edm::InputTag("hltHoreco"))->setComment("HO recHit Collection"); + desc.add("ebInput", edm::InputTag("hltEcalRecHit:EcalRecHitsEB"))->setComment("EB recHit Collection"); + desc.add("eeInput", edm::InputTag("hltEcalRecHit:EcalRecHitsEE"))->setComment("EE recHit Collection"); + desc.add("pataVertices", edm::InputTag("hltPixelVerticesSoA")) + ->setComment("patatrack vertices collection"); + desc.add("pataTracks", edm::InputTag("hltPixelTracksSoA"))->setComment("patatrack collection"); + desc.add("BeamSpot", edm::InputTag("hltOnlineBeamSpot"))->setComment("BeamSpot Collection"); + desc.add("maxVtx", 100)->setComment("max output collection size (number of accepted vertices)"); + desc.add("fractionSumPt2", 0.3)->setComment("threshold on sumPt2 fraction of the leading vertex"); + desc.add("minSumPt2", 0.)->setComment("min sumPt2"); + desc.add("track_pt_min", 1.0)->setComment("min track p_T"); + desc.add("track_pt_max", 10.0)->setComment("max track p_T"); + desc.add("track_chi2_max", 99999.)->setComment("max track chi2"); + desc.add("graphPath", "RecoTauTag/TrainingFiles/data/L2TauNNTag/L2TauTag_Run3v1.pb") + ->setComment("path to the saved CNN"); + desc.add("normalizationDict", "RecoTauTag/TrainingFiles/data/L2TauNNTag/NormalizationDict.json") + ->setComment("path to the dictionary for variable standardization"); + descriptions.addWithDefaultLabel(desc); +} + +L2TauNNProducerAlpakaStrip::L2TauNNProducerAlpakaStrip(const edm::ParameterSet& cfg, + const L2TauNNProducerAlpakaStripCacheData* cacheData) + : debugLevel_(cfg.getParameter("debugLevel")), + hbheToken_(consumes(cfg.getParameter("hbheInput"))), + hoToken_(consumes(cfg.getParameter("hoInput"))), + ebToken_(consumes(cfg.getParameter("ebInput"))), + eeToken_(consumes(cfg.getParameter("eeInput"))), + geometryToken_(esConsumes()), + bFieldToken_(esConsumes()), + pataVerticesToken_(consumes(cfg.getParameter("pataVertices"))), + pataTracksToken_(consumes(cfg.getParameter("pataTracks"))), + beamSpotToken_(consumes(cfg.getParameter("BeamSpot"))), + maxVtx_(cfg.getParameter("maxVtx")), + fractionSumPt2_(cfg.getParameter("fractionSumPt2")), + minSumPt2_(cfg.getParameter("minSumPt2")), + trackPtMin_(cfg.getParameter("track_pt_min")), + trackPtMax_(cfg.getParameter("track_pt_max")), + trackChi2Max_(cfg.getParameter("track_chi2_max")) { + if (cacheData->graphDef == nullptr) { + throw cms::Exception("InvalidCacheData") << "Invalid Cache Data."; + } + inputTensorName_ = cacheData->graphDef->node(0).name(); + outputTensorName_ = cacheData->graphDef->node(cacheData->graphDef->node_size() - 1).name(); + L2cacheData_ = cacheData; + std::vector L1TauCollections = cfg.getParameter>("L1Taus"); + L1TauDesc_.reserve(L1TauCollections.size()); + for (const auto& l1TauInput : L1TauCollections) { + InputDescTau toInsert; + toInsert.CollectionName = l1TauInput.getParameter("L1CollectionName"); + toInsert.inputToken_ = + consumes(l1TauInput.getParameter("L1TauTrigger")); + L1TauDesc_.push_back(toInsert); + } + for (const auto& desc : L1TauDesc_) + produces>(desc.CollectionName); +} + +void L2TauNNProducerAlpakaStrip::checknan(tensorflow::Tensor& tensor, int debugLevel) { + using NNInputs = L2TauTagNNv1::NNInputs; + std::vector tensor_shape(tensor.shape().dims()); + for (int d = 0; d < tensor.shape().dims(); d++) { + tensor_shape.at(d) = tensor.shape().dim_size(d); + } + if (tensor_shape.size() != 4) { + throw cms::Exception("InvalidTensor") << "Tensor shape does not have 4 dimensions!"; + } + for (int tau_idx = 0; tau_idx < tensor_shape.at(0); tau_idx++) { + for (int phi_idx = 0; phi_idx < tensor_shape.at(1); phi_idx++) { + for (int eta_idx = 0; eta_idx < tensor_shape.at(2); eta_idx++) { + for (int var_idx = 0; var_idx < tensor_shape.at(3); var_idx++) { + auto getCell = [&](NNInputs input) -> float& { + return getCellImpl(tensor, tau_idx, phi_idx, eta_idx, input); + }; + auto nonstd_var = getCell(static_cast(var_idx)); + if (edm::isNotFinite(nonstd_var)) { + edm::LogWarning("InputVar") << "var is nan \nvar name= " + << L2TauTagNNv1::varNameMap.at(static_cast(var_idx)) + << "\t var_idx = " << var_idx << "\t eta_idx = " << eta_idx + << "\t phi_idx = " << phi_idx << "\t tau_idx = " << tau_idx; + if (debugLevel > 2) { + edm::LogWarning("InputVar") << "other vars in same cell \n"; + if (var_idx + 1 < tensor_shape.at(3)) + edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast(var_idx + 1)) + << "\t = " << getCell(static_cast(var_idx + 1)); + if (var_idx + 2 < tensor_shape.at(3)) + edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast(var_idx + 2)) + << "\t = " << getCell(static_cast(var_idx + 2)); + if (var_idx + 3 < tensor_shape.at(3)) + edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast(var_idx + 3)) + << "\t = " << getCell(static_cast(var_idx + 3)); + if (var_idx + 4 < tensor_shape.at(3)) + edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast(var_idx + 4)) + << "\t = " << getCell(static_cast(var_idx + 4)); + } + } + } + } + } + } +} + +void L2TauNNProducerAlpakaStrip::standardizeTensor(tensorflow::Tensor& tensor) { + using NNInputs = L2TauTagNNv1::NNInputs; + std::vector tensor_shape(tensor.shape().dims()); + for (int d = 0; d < tensor.shape().dims(); d++) { + tensor_shape.at(d) = tensor.shape().dim_size(d); + } + if (tensor_shape.size() != 4) { + throw cms::Exception("InvalidTensor") << "Tensor shape does not have 4 dimensions!"; + } + for (int tau_idx = 0; tau_idx < tensor_shape.at(0); tau_idx++) { + for (int phi_idx = 0; phi_idx < tensor_shape.at(1); phi_idx++) { + for (int eta_idx = 0; eta_idx < tensor_shape.at(2); eta_idx++) { + for (int var_idx = 0; var_idx < tensor_shape.at(3); var_idx++) { + auto getCell = [&](NNInputs input) -> float& { + return getCellImpl(tensor, tau_idx, phi_idx, eta_idx, input); + }; + float mean = L2cacheData_->normVec.at(var_idx).mean; + float std = L2cacheData_->normVec.at(var_idx).std; + float min = L2cacheData_->normVec.at(var_idx).min; + float max = L2cacheData_->normVec.at(var_idx).max; + float nonstd_var = getCell(static_cast(var_idx)); + float std_var = static_cast((nonstd_var - mean) / std); + if (std_var > max) { + std_var = static_cast(max); + } else if (std_var < min) { + std_var = static_cast(min); + } + getCell(static_cast(var_idx)) = std_var; + } + } + } + } +} + +void L2TauNNProducerAlpakaStrip::fillL1TauVars(tensorflow::Tensor& cellGridMatrix, + const std::vector& allTaus) { + using NNInputs = L2TauTagNNv1::NNInputs; + + const int nTaus = allTaus.size(); + for (int tau_idx = 0; tau_idx < nTaus; tau_idx++) { + for (int eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) { + for (int phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) { + auto getCell = [&](NNInputs input) -> float& { + return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input); + }; + getCell(NNInputs::l1Tau_pt) = allTaus[tau_idx]->pt(); + getCell(NNInputs::l1Tau_eta) = allTaus[tau_idx]->eta(); + getCell(NNInputs::l1Tau_hwIso) = allTaus[tau_idx]->hwIso(); + } + } + } +} + +template +std::tuple L2TauNNProducerAlpakaStrip::getEtaPhiIndices(float eta, + float phi, + const LVec& tau_p4) { + const float deta = eta - tau_p4.eta(); + const float dphi = reco::deltaPhi(phi, tau_p4.phi()); + const int eta_idx = static_cast(floor((deta + L2TauTagNNv1::dR_max) / dEta_width)); + const int phi_idx = static_cast(floor((dphi + L2TauTagNNv1::dR_max) / dPhi_width)); + return std::make_tuple(deta, dphi, eta_idx, phi_idx); +} + +template +std::tuple L2TauNNProducerAlpakaStrip::getEtaPhiIndices(const VPos& position, + const LVec& tau_p4) { + return getEtaPhiIndices(position.eta(), position.phi(), tau_p4); +} + +void L2TauNNProducerAlpakaStrip::fillCaloRecHits(tensorflow::Tensor& cellGridMatrix, + const std::vector& allTaus, + const caloRecHitCollections& caloRecHits) { + using NNInputs = L2TauTagNNv1::NNInputs; + + const int nTaus = allTaus.size(); + float deta, dphi; + int eta_idx = 0; + int phi_idx = 0; + int tau_idx = 0; + + auto getCell = [&](NNInputs input) -> float& { + return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input); + }; + for (tau_idx = 0; tau_idx < nTaus; tau_idx++) { + // calorechit_EE + for (const auto& caloRecHit_ee : *caloRecHits.ee) { + if (caloRecHit_ee.energy() <= 0) + continue; + const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_ee.id())->getPosition(); + const float eeCalEn = caloRecHit_ee.energy(); + const float eeCalChi2 = caloRecHit_ee.chi2(); + if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) { + std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4()); + getCell(NNInputs::EcalEnergySum) += eeCalEn; + getCell(NNInputs::EcalSize) += 1.; + getCell(NNInputs::EcalEnergyStdDev) += eeCalEn * eeCalEn; + getCell(NNInputs::EcalDeltaEta) += deta * eeCalEn; + getCell(NNInputs::EcalDeltaPhi) += dphi * eeCalEn; + if (eeCalChi2 >= 0) { + getCell(NNInputs::EcalChi2) += eeCalChi2 * eeCalEn; + getCell(NNInputs::EcalEnergySumForPositiveChi2) += eeCalEn; + getCell(NNInputs::EcalSizeForPositiveChi2) += 1.; + } + } + } + + // calorechit_EB + for (const auto& caloRecHit_eb : *caloRecHits.eb) { + if (caloRecHit_eb.energy() <= 0) + continue; + const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_eb.id())->getPosition(); + const float ebCalEn = caloRecHit_eb.energy(); + const float ebCalChi2 = caloRecHit_eb.chi2(); + if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) { + std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4()); + getCell(NNInputs::EcalEnergySum) += ebCalEn; + getCell(NNInputs::EcalSize) += 1.; + getCell(NNInputs::EcalEnergyStdDev) += ebCalEn * ebCalEn; + getCell(NNInputs::EcalDeltaEta) += deta * ebCalEn; + getCell(NNInputs::EcalDeltaPhi) += dphi * ebCalEn; + if (ebCalChi2 >= 0) { + getCell(NNInputs::EcalChi2) += ebCalChi2 * ebCalEn; + getCell(NNInputs::EcalEnergySumForPositiveChi2) += ebCalEn; + getCell(NNInputs::EcalSizeForPositiveChi2) += 1.; + } + } + } + + // calorechit_HBHE + for (const auto& caloRecHit_hbhe : *caloRecHits.hbhe) { + if (caloRecHit_hbhe.energy() <= 0) + continue; + const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_hbhe.id())->getPosition(); + const float hbheCalEn = caloRecHit_hbhe.energy(); + const float hbheCalChi2 = caloRecHit_hbhe.chi2(); + if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) { + std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4()); + getCell(NNInputs::HcalEnergySum) += hbheCalEn; + getCell(NNInputs::HcalEnergyStdDev) += hbheCalEn * hbheCalEn; + getCell(NNInputs::HcalSize) += 1.; + getCell(NNInputs::HcalDeltaEta) += deta * hbheCalEn; + getCell(NNInputs::HcalDeltaPhi) += dphi * hbheCalEn; + if (hbheCalChi2 >= 0) { + getCell(NNInputs::HcalChi2) += hbheCalChi2 * hbheCalEn; + getCell(NNInputs::HcalEnergySumForPositiveChi2) += hbheCalEn; + getCell(NNInputs::HcalSizeForPositiveChi2) += 1.; + } + } + } + + // calorechit_HO + for (const auto& caloRecHit_ho : *caloRecHits.ho) { + if (caloRecHit_ho.energy() <= 0) + continue; + const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_ho.id())->getPosition(); + const float hoCalEn = caloRecHit_ho.energy(); + if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) { + std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4()); + getCell(NNInputs::HcalEnergySum) += hoCalEn; + getCell(NNInputs::HcalEnergyStdDev) += hoCalEn * hoCalEn; + getCell(NNInputs::HcalSize) += 1.; + getCell(NNInputs::HcalDeltaEta) += deta * hoCalEn; + getCell(NNInputs::HcalDeltaPhi) += dphi * hoCalEn; + } + } + + // normalize to sum and define stdDev + for (eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) { + for (phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) { + /* normalize eCal vars*/ + if (getCell(NNInputs::EcalEnergySum) > 0.) { + getCell(NNInputs::EcalDeltaEta) /= getCell(NNInputs::EcalEnergySum); + getCell(NNInputs::EcalDeltaPhi) /= getCell(NNInputs::EcalEnergySum); + } + if (getCell(NNInputs::EcalEnergySumForPositiveChi2) > 0.) { + getCell(NNInputs::EcalChi2) /= getCell(NNInputs::EcalEnergySumForPositiveChi2); + } + if (getCell(NNInputs::EcalSize) > 1.) { + // (stdDev - (enSum*enSum)/size) / (size-1) + getCell(NNInputs::EcalEnergyStdDev) = + (getCell(NNInputs::EcalEnergyStdDev) - + (getCell(NNInputs::EcalEnergySum) * getCell(NNInputs::EcalEnergySum)) / getCell(NNInputs::EcalSize)) / + (getCell(NNInputs::EcalSize) - 1); + } else { + getCell(NNInputs::EcalEnergyStdDev) = 0.; + } + /* normalize hCal Vars */ + if (getCell(NNInputs::HcalEnergySum) > 0.) { + getCell(NNInputs::HcalDeltaEta) /= getCell(NNInputs::HcalEnergySum); + getCell(NNInputs::HcalDeltaPhi) /= getCell(NNInputs::HcalEnergySum); + } + if (getCell(NNInputs::HcalEnergySumForPositiveChi2) > 0.) { + getCell(NNInputs::HcalChi2) /= getCell(NNInputs::HcalEnergySumForPositiveChi2); + } + if (getCell(NNInputs::HcalSize) > 1.) { + // (stdDev - (enSum*enSum)/size) / (size-1) + getCell(NNInputs::HcalEnergyStdDev) = + (getCell(NNInputs::HcalEnergyStdDev) - + (getCell(NNInputs::HcalEnergySum) * getCell(NNInputs::HcalEnergySum)) / getCell(NNInputs::HcalSize)) / + (getCell(NNInputs::HcalSize) - 1); + } else { + getCell(NNInputs::HcalEnergyStdDev) = 0.; + } + } + } + } +} + +void L2TauNNProducerAlpakaStrip::selectGoodTracksAndVertices(const ZVertexHost& patavtx_soa, + const TracksHost& patatracks_tsoa, + std::vector& trkGood, + std::vector& vtxGood) { + using patatrackHelpers = TracksUtilities; + const auto maxTracks = patatracks_tsoa.view().metadata().size(); + const int nv = patavtx_soa.view().nvFinal(); + trkGood.clear(); + trkGood.reserve(maxTracks); + vtxGood.clear(); + vtxGood.reserve(nv); + auto const* quality = patatracks_tsoa.view().quality(); + + // No need to sort either as the algorithms is just using the max (not even the location, just the max value of pt2sum). + std::vector pTSquaredSum(nv, 0); + std::vector nTrkAssociated(nv, 0); + + for (int32_t trk_idx = 0; trk_idx < maxTracks; ++trk_idx) { + auto nHits = patatrackHelpers::nHits(patatracks_tsoa.view(), trk_idx); + if (nHits == 0) { + break; + } + int vtx_ass_to_track = patavtx_soa.view()[trk_idx].idv(); + if (vtx_ass_to_track >= 0 && vtx_ass_to_track < nv) { + auto patatrackPt = patatracks_tsoa.view()[trk_idx].pt(); + ++nTrkAssociated[vtx_ass_to_track]; + if (patatrackPt >= trackPtMin_ && patatracks_tsoa.const_view()[trk_idx].chi2() <= trackChi2Max_) { + patatrackPt = std::min(patatrackPt, trackPtMax_); + pTSquaredSum[vtx_ass_to_track] += patatrackPt * patatrackPt; + } + } + if (nHits > 0 and quality[trk_idx] >= pixelTrack::Quality::loose) { + trkGood.push_back(trk_idx); + } + } + if (nv > 0) { + const auto minFOM_fromFrac = (*std::max_element(pTSquaredSum.begin(), pTSquaredSum.end())) * fractionSumPt2_; + for (int j = nv - 1; j >= 0 && vtxGood.size() < maxVtx_; --j) { + auto vtx_idx = patavtx_soa.view()[j].sortInd(); + assert(vtx_idx < nv); + if (nTrkAssociated[vtx_idx] >= 2 && pTSquaredSum[vtx_idx] >= minFOM_fromFrac && + pTSquaredSum[vtx_idx] > minSumPt2_) { + vtxGood.push_back(vtx_idx); + } + } + } +} + +std::pair L2TauNNProducerAlpakaStrip::impactParameter(int it, + const TracksHost& patatracks_tsoa, + float patatrackPhi, + const reco::BeamSpot& beamspot, + const MagneticField* magfi) { + /* dxy and dz */ + riemannFit::Vector5d ipar, opar; + riemannFit::Matrix5d icov, ocov; + TracksUtilities::copyToDense(patatracks_tsoa.view(), ipar, icov, it); + riemannFit::transformToPerigeePlane(ipar, icov, opar, ocov); + LocalTrajectoryParameters lpar(opar(0), opar(1), opar(2), opar(3), opar(4), 1.); + float sp = std::sin(patatrackPhi); + float cp = std::cos(patatrackPhi); + Surface::RotationType Rotation(sp, -cp, 0, 0, 0, -1.f, cp, sp, 0); + GlobalPoint BeamSpotPoint(beamspot.x0(), beamspot.y0(), beamspot.z0()); + Plane impPointPlane(BeamSpotPoint, Rotation); + GlobalTrajectoryParameters gp( + impPointPlane.toGlobal(lpar.position()), impPointPlane.toGlobal(lpar.momentum()), lpar.charge(), magfi); + GlobalPoint vv = gp.position(); + math::XYZPoint pos(vv.x(), vv.y(), vv.z()); + GlobalVector pp = gp.momentum(); + math::XYZVector mom(pp.x(), pp.y(), pp.z()); + auto lambda = M_PI_2 - pp.theta(); + auto phi = pp.phi(); + float patatrackDxy = -vv.x() * std::sin(phi) + vv.y() * std::cos(phi); + float patatrackDz = + (vv.z() * std::cos(lambda) - (vv.x() * std::cos(phi) + vv.y() * std::sin(phi)) * std::sin(lambda)) / + std::cos(lambda); + return std::make_pair(patatrackDxy, patatrackDz); +} + +void L2TauNNProducerAlpakaStrip::fillPatatracks(tensorflow::Tensor& cellGridMatrix, + const std::vector& allTaus, + const TracksHost& patatracks_tsoa, + const ZVertexHost& patavtx_soa, + const reco::BeamSpot& beamspot, + const MagneticField* magfi) { + using NNInputs = L2TauTagNNv1::NNInputs; + using patatrackHelpers = TracksUtilities; + float deta, dphi; + int eta_idx = 0; + int phi_idx = 0; + int tau_idx = 0; + + auto getCell = [&](NNInputs input) -> float& { + return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input); + }; + + std::vector trkGood; + std::vector vtxGood; + + selectGoodTracksAndVertices(patavtx_soa, patatracks_tsoa, trkGood, vtxGood); + + const int nTaus = allTaus.size(); + for (tau_idx = 0; tau_idx < nTaus; tau_idx++) { + const float tauEta = allTaus[tau_idx]->eta(); + const float tauPhi = allTaus[tau_idx]->phi(); + + for (const auto it : trkGood) { + const float patatrackPt = patatracks_tsoa.const_view()[it].pt(); + if (patatrackPt <= 0) + continue; + const float patatrackPhi = reco::phi(patatracks_tsoa.const_view(), it); + const float patatrackEta = patatracks_tsoa.const_view()[it].eta(); + const float patatrackCharge = reco::charge(patatracks_tsoa.const_view(), it); + const float patatrackChi2OverNdof = patatracks_tsoa.view()[it].chi2(); + const auto nHits = patatrackHelpers::nHits(patatracks_tsoa.const_view(), it); + if (nHits <= 0) + continue; + const int patatrackNdof = 2 * std::min(6, nHits) - 5; + + const int vtx_idx_assTrk = patavtx_soa.view()[it].idv(); + if (reco::deltaR2(patatrackEta, patatrackPhi, tauEta, tauPhi) < dR2_max) { + std::tie(deta, dphi, eta_idx, phi_idx) = + getEtaPhiIndices(patatrackEta, patatrackPhi, allTaus[tau_idx]->polarP4()); + getCell(NNInputs::PatatrackPtSum) += patatrackPt; + getCell(NNInputs::PatatrackSize) += 1.; + getCell(NNInputs::PatatrackChargeSum) += patatrackCharge; + getCell(NNInputs::PatatrackDeltaEta) += deta * patatrackPt; + getCell(NNInputs::PatatrackDeltaPhi) += dphi * patatrackPt; + getCell(NNInputs::PatatrackChi2OverNdof) += patatrackChi2OverNdof * patatrackPt; + getCell(NNInputs::PatatrackNdof) += patatrackNdof * patatrackPt; + std::pair impactParameters = impactParameter(it, patatracks_tsoa, patatrackPhi, beamspot, magfi); + getCell(NNInputs::PatatrackDxy) += impactParameters.first * patatrackPt; + getCell(NNInputs::PatatrackDz) += impactParameters.second * patatrackPt; + if ((std::find(vtxGood.begin(), vtxGood.end(), vtx_idx_assTrk) != vtxGood.end())) { + getCell(NNInputs::PatatrackPtSumWithVertex) += patatrackPt; + getCell(NNInputs::PatatrackSizeWithVertex) += 1.; + } + } + } + + // normalize to sum and define stdDev + for (eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) { + for (phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) { + getCell(NNInputs::nVertices) = vtxGood.size(); + if (getCell(NNInputs::PatatrackPtSum) > 0.) { + getCell(NNInputs::PatatrackDeltaEta) /= getCell(NNInputs::PatatrackPtSum); + getCell(NNInputs::PatatrackDeltaPhi) /= getCell(NNInputs::PatatrackPtSum); + getCell(NNInputs::PatatrackChi2OverNdof) /= getCell(NNInputs::PatatrackPtSum); + getCell(NNInputs::PatatrackNdof) /= getCell(NNInputs::PatatrackPtSum); + getCell(NNInputs::PatatrackDxy) /= getCell(NNInputs::PatatrackPtSum); + getCell(NNInputs::PatatrackDz) /= getCell(NNInputs::PatatrackPtSum); + } + } + } + } +} + +std::vector L2TauNNProducerAlpakaStrip::getTauScore(const tensorflow::Tensor& cellGridMatrix) { + const int nTau = cellGridMatrix.shape().dim_size(0); + if (nTau == 0) { + return std::vector(); + } else { + std::vector pred_tensor; + tensorflow::run(L2cacheData_->session, {{inputTensorName_, cellGridMatrix}}, {outputTensorName_}, &pred_tensor); + std::vector pred_vector(nTau); + for (int tau_idx = 0; tau_idx < nTau; ++tau_idx) { + pred_vector[tau_idx] = pred_tensor[0].matrix()(tau_idx, 0); + } + + return pred_vector; + } +} + +void L2TauNNProducerAlpakaStrip::produce(edm::Event& event, const edm::EventSetup& eventsetup) { + std::vector> TauCollectionMap(L1TauDesc_.size()); + l1t::TauVectorRef allTaus; + + for (size_t inp_idx = 0; inp_idx < L1TauDesc_.size(); inp_idx++) { + l1t::TauVectorRef l1Taus; + auto const& l1TriggeredTaus = event.get(L1TauDesc_[inp_idx].inputToken_); + l1TriggeredTaus.getObjects(trigger::TriggerL1Tau, l1Taus); + TauCollectionMap.at(inp_idx).resize(l1Taus.size()); + + for (size_t l1_idx = 0; l1_idx < l1Taus.size(); l1_idx++) { + size_t tau_idx; + const auto iter = std::find(allTaus.begin(), allTaus.end(), l1Taus[l1_idx]); + if (iter != allTaus.end()) { + tau_idx = std::distance(allTaus.begin(), iter); + } else { + allTaus.push_back(l1Taus[l1_idx]); + tau_idx = allTaus.size() - 1; + } + TauCollectionMap.at(inp_idx).at(l1_idx) = tau_idx; + } + } + const auto ebCal = event.getHandle(ebToken_); + const auto eeCal = event.getHandle(eeToken_); + const auto hbhe = event.getHandle(hbheToken_); + const auto ho = event.getHandle(hoToken_); + auto const& patatracks_SoA = event.get(pataTracksToken_); + auto const& vertices_SoA = event.get(pataVerticesToken_); + const auto bsHandle = event.getHandle(beamSpotToken_); + + auto const fieldESH = eventsetup.getHandle(bFieldToken_); + auto const geometry = eventsetup.getHandle(geometryToken_); + + caloRecHitCollections caloRecHits; + caloRecHits.hbhe = &*hbhe; + caloRecHits.ho = &*ho; + caloRecHits.eb = &*ebCal; + caloRecHits.ee = &*eeCal; + caloRecHits.geometry = &*geometry; + + const int nTaus = allTaus.size(); + tensorflow::Tensor cellGridMatrix(tensorflow::DT_FLOAT, + {nTaus, L2TauTagNNv1::nCellEta, L2TauTagNNv1::nCellPhi, L2TauTagNNv1::nVars}); + const int n_inputs = nTaus * L2TauTagNNv1::nCellEta * L2TauTagNNv1::nCellPhi * L2TauTagNNv1::nVars; + for (int input_idx = 0; input_idx < n_inputs; ++input_idx) { + cellGridMatrix.flat()(input_idx) = 0; + } + fillL1TauVars(cellGridMatrix, allTaus); + + fillCaloRecHits(cellGridMatrix, allTaus, caloRecHits); + + fillPatatracks(cellGridMatrix, allTaus, patatracks_SoA, vertices_SoA, *bsHandle, fieldESH.product()); + + standardizeTensor(cellGridMatrix); + + if (debugLevel_ > 0) { + checknan(cellGridMatrix, debugLevel_); + } + + std::vector tau_score = getTauScore(cellGridMatrix); + + for (size_t inp_idx = 0; inp_idx < L1TauDesc_.size(); inp_idx++) { + const size_t nTau = TauCollectionMap[inp_idx].size(); + auto tau_tags = std::make_unique>(nTau); + for (size_t tau_pos = 0; tau_pos < nTau; ++tau_pos) { + const auto tau_idx = TauCollectionMap[inp_idx][tau_pos]; + if (debugLevel_ > 0) { + edm::LogInfo("DebugInfo") << event.id().event() << " \t " << (allTaus[tau_idx])->pt() << " \t " + << tau_score.at(tau_idx) << std::endl; + } + (*tau_tags)[tau_pos] = tau_score.at(tau_idx); + } + event.put(std::move(tau_tags), L1TauDesc_[inp_idx].CollectionName); + } +} +//define this as a plug-in +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(L2TauNNProducerAlpakaStrip); diff --git a/RecoTracker/PixelSeeding/BuildFile.xml b/RecoTracker/PixelSeeding/BuildFile.xml index 7bc10578b4448..e3875d0fcfdb3 100644 --- a/RecoTracker/PixelSeeding/BuildFile.xml +++ b/RecoTracker/PixelSeeding/BuildFile.xml @@ -16,6 +16,7 @@ + diff --git a/RecoTracker/PixelSeeding/plugins/BuildFile.xml b/RecoTracker/PixelSeeding/plugins/BuildFile.xml index a387c35caa691..5afedd5ca7574 100644 --- a/RecoTracker/PixelSeeding/plugins/BuildFile.xml +++ b/RecoTracker/PixelSeeding/plugins/BuildFile.xml @@ -11,6 +11,9 @@ + + + diff --git a/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorKernels.h b/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorKernels.h index 250aef21c1d6a..44a19d653cd06 100644 --- a/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorKernels.h +++ b/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorKernels.h @@ -40,6 +40,15 @@ namespace caHitNtupletGenerator { const float hardCurvCut_; const float dcaCutInnerTriplet_; const float dcaCutOuterTriplet_; + const float CAThetaCutBarrelPixelBarrelStrip_; + const float CAThetaCutBarrelPixelForwardStrip_; + const float CAThetaCutBarrelStripForwardStrip_; + const float CAThetaCutBarrelStrip_; + const float CAThetaCutDefault_; + const float dcaCutInnerTripletPixelStrip_; + const float dcaCutOuterTripletPixelStrip_; + const float dcaCutTripletStrip_; + const float dcaCutTripletDefault_; }; template diff --git a/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorOnGPU.cc b/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorOnGPU.cc index 5100cf734142c..4492c7943a9a3 100644 --- a/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorOnGPU.cc +++ b/RecoTracker/PixelSeeding/plugins/CAHitNtupletGeneratorOnGPU.cc @@ -61,16 +61,24 @@ namespace { template struct topologyCuts> { static constexpr CAParamsT makeCACuts(edm::ParameterSet const& cfg) { - return CAParamsT{{ - cfg.getParameter("maxNumberOfDoublets"), - cfg.getParameter("minHitsPerNtuplet"), - static_cast(cfg.getParameter("ptmin")), - static_cast(cfg.getParameter("CAThetaCutBarrel")), - static_cast(cfg.getParameter("CAThetaCutForward")), - static_cast(cfg.getParameter("hardCurvCut")), - static_cast(cfg.getParameter("dcaCutInnerTriplet")), - static_cast(cfg.getParameter("dcaCutOuterTriplet")), - }}; + return CAParamsT{ + {cfg.getParameter("maxNumberOfDoublets"), + cfg.getParameter("minHitsPerNtuplet"), + static_cast(cfg.getParameter("ptmin")), + static_cast(cfg.getParameter("CAThetaCutBarrel")), + static_cast(cfg.getParameter("CAThetaCutForward")), + static_cast(cfg.getParameter("hardCurvCut")), + static_cast(cfg.getParameter("dcaCutInnerTriplet")), + static_cast(cfg.getParameter("dcaCutOuterTriplet")), + static_cast(cfg.getParameter("CAThetaCutBarrelPixelBarrelStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelPixelForwardStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelStripForwardStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelStrip")), + static_cast(cfg.getParameter("CAThetaCutDefault")), + static_cast(cfg.getParameter("dcaCutInnerTripletPixelStrip")), + static_cast(cfg.getParameter("dcaCutOuterTripletPixelStrip")), + static_cast(cfg.getParameter("dcaCutTripletStrip")), + static_cast(cfg.getParameter("dcaCutTripletDefault"))}}; }; static constexpr pixelTrack::QualityCutsT makeQualityCuts(edm::ParameterSet const& pset) { @@ -98,15 +106,25 @@ namespace { template struct topologyCuts> { static constexpr CAParamsT makeCACuts(edm::ParameterSet const& cfg) { - return CAParamsT{{cfg.getParameter("maxNumberOfDoublets"), - cfg.getParameter("minHitsPerNtuplet"), - static_cast(cfg.getParameter("ptmin")), - static_cast(cfg.getParameter("CAThetaCutBarrel")), - static_cast(cfg.getParameter("CAThetaCutForward")), - static_cast(cfg.getParameter("hardCurvCut")), - static_cast(cfg.getParameter("dcaCutInnerTriplet")), - static_cast(cfg.getParameter("dcaCutOuterTriplet"))}, - {(bool)cfg.getParameter("includeFarForwards")}}; + return CAParamsT{ + {cfg.getParameter("maxNumberOfDoublets"), + cfg.getParameter("minHitsPerNtuplet"), + static_cast(cfg.getParameter("ptmin")), + static_cast(cfg.getParameter("CAThetaCutBarrel")), + static_cast(cfg.getParameter("CAThetaCutForward")), + static_cast(cfg.getParameter("hardCurvCut")), + static_cast(cfg.getParameter("dcaCutInnerTriplet")), + static_cast(cfg.getParameter("dcaCutOuterTriplet")), + static_cast(cfg.getParameter("CAThetaCutBarrelPixelBarrelStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelPixelForwardStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelStripForwardStrip")), + static_cast(cfg.getParameter("CAThetaCutBarrelStrip")), + static_cast(cfg.getParameter("CAThetaCutDefault")), + static_cast(cfg.getParameter("dcaCutInnerTripletPixelStrip")), + static_cast(cfg.getParameter("dcaCutOuterTripletPixelStrip")), + static_cast(cfg.getParameter("dcaCutTripletStrip")), + static_cast(cfg.getParameter("dcaCutTripletDefault"))}, + {(bool)cfg.getParameter("includeFarForwards")}}; } static constexpr pixelTrack::QualityCutsT makeQualityCuts(edm::ParameterSet const& pset) { diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/BrokenLineFit.dev.cc b/RecoTracker/PixelSeeding/plugins/alpaka/BrokenLineFit.dev.cc index f8cadf653114d..502cf8eeddeca 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/BrokenLineFit.dev.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/BrokenLineFit.dev.cc @@ -7,9 +7,13 @@ #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" +#include "HeterogeneousCore/AlpakaInterface/interface/traits.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h" #include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforDevice.h" #include "RecoTracker/PixelTrackFitting/interface/alpaka/BrokenLine.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" + #include "HelixFit.h" template @@ -30,7 +34,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { Tuples const *__restrict__ foundNtuplets, TupleMultiplicity const *__restrict__ tupleMultiplicity, TrackingRecHitSoAConstView hh, - pixelCPEforDevice::ParamsOnDeviceT const *__restrict__ cpeParams, + // pixelCPEforDevice::ParamsOnDeviceT const *__restrict__ cpeParams, + FrameSoAConstView fr, typename TrackerTraits::tindex_type *__restrict__ ptkids, double *__restrict__ phits, float *__restrict__ phits_ge, @@ -112,15 +117,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { n += incr; auto hit = hitId[j]; float ge[6]; - + auto const &frame = fr.detFrame(hh.detectorIndex(hit)); #ifdef YERR_FROM_DC - auto const &dp = cpeParams->detParams(hh.detectorIndex(hit)); auto status = hh[hit].chargeAndStatus().status; int qbin = CPEFastParametrisation::kGenErrorQBins - 1 - status.qBin; ALPAKA_ASSERT_ACC(qbin >= 0 && qbin < 5); bool nok = (status.isBigY | status.isOneY); // compute cotanbeta and use it to recompute error - dp.frame.rotation().multiply(dx, dy, dz, ux, uy, uz); + frame.rotation().multiply(dx, dy, dz, ux, uy, uz); auto cb = std::abs(uy / uz); int bin = int(cb * (float(phase1PixelTopology::pixelThickess) / float(phase1PixelTopology::pixelPitchY)) * 8.f) - 4; @@ -133,9 +137,23 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { yerr *= yerr; yerr += dp.apeYY; yerr = nok ? hh[hit].yerrLocal() : yerr; - dp.frame.toGlobal(hh[hit].xerrLocal(), 0, yerr, ge); + frame.toGlobal(hh[hit].xerrLocal(), 0, yerr, ge); #else - cpeParams->detParams(hh[hit].detectorIndex()).frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); + frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); + // if (hh[hit].detectorIndex() <= TrackerTraits::numberOfPixelModules) + // cpeParams->detParams(hh[hit].detectorIndex()).frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); + // else + // { + // auto xe = hh[hit].xerrLocal(); + // auto ye = hh[hit].yerrLocal(); + + // ge[0] = xe; + // ge[1] = sqrt(xe*xe + ye*ye); + // ge[2] = ye; + // ge[3] = sqrt(xe*xe + ye*ye); + // ge[4] = sqrt(xe*xe + ye*ye); + // ge[5] = sqrt(xe*xe + ye*ye); + // } #endif #ifdef BL_DUMP_HITS @@ -244,7 +262,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template void HelixFit::launchBrokenLineKernels( const TrackingRecHitSoAConstView &hv, - pixelCPEforDevice::ParamsOnDeviceT const *cpeParams, + // pixelCPEforDevice::ParamsOnDeviceT const *cpeParams, + const FrameSoAConstView &fr, uint32_t hitsInFit, uint32_t maxNumberOfTuples, Queue &queue) { @@ -274,7 +293,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tuples_, tupleMultiplicity_, hv, - cpeParams, + fr, tkidDevice.data(), hitsDevice.data(), hits_geDevice.data(), @@ -298,7 +317,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // fit all as 4 riemannFit::rolling_fits<4, TrackerTraits::maxHitsOnTrack, 1>([this, &hv, - &cpeParams, + &fr, &tkidDevice, &hitsDevice, &hits_geDevice, @@ -312,7 +331,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tuples_, tupleMultiplicity_, hv, - cpeParams, + fr, tkidDevice.data(), hitsDevice.data(), hits_geDevice.data(), @@ -336,7 +355,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { } else { riemannFit::rolling_fits<4, TrackerTraits::maxHitsOnTrackForFullFit, 1>([this, &hv, - &cpeParams, + &fr, &tkidDevice, &hitsDevice, &hits_geDevice, @@ -350,7 +369,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tuples_, tupleMultiplicity_, hv, - cpeParams, + fr, tkidDevice.data(), hitsDevice.data(), hits_geDevice.data(), @@ -380,7 +399,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tuples_, tupleMultiplicity_, hv, - cpeParams, + fr, tkidDevice.data(), hitsDevice.data(), hits_geDevice.data(), @@ -407,5 +426,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class HelixFit; template class HelixFit; template class HelixFit; + template class HelixFit; } // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CACell.h b/RecoTracker/PixelSeeding/plugins/alpaka/CACell.h index 44b7bea5075fb..4c2c03079b3a4 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CACell.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CACell.h @@ -2,16 +2,14 @@ #define RecoTracker_PixelSeeding_plugins_alpaka_CACell_h // #define ONLY_TRIPLETS_IN_HOLE - #include #include - #include #include "DataFormats/TrackSoA/interface/TrackDefinitions.h" #include "DataFormats/TrackSoA/interface/TracksSoA.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/SimpleVector.h" #include "HeterogeneousCore/AlpakaInterface/interface/VecArray.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" @@ -43,7 +41,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using HitContainer = typename reco::TrackSoA::HitContainer; using Quality = ::pixelTrack::Quality; static constexpr auto bad = ::pixelTrack::Quality::bad; - enum class StatusBit : uint16_t { kUsed = 1, kInTrack = 2, kKilled = 1 << 15 }; CACellT() = default; @@ -164,10 +161,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const float caThetaCutBarrel, const float caThetaCutForward, const float dcaCutInnerTriplet, - const float dcaCutOuterTriplet) const { + const float dcaCutOuterTriplet, + const float caThetaCutBarrelPixelBarrelStrip, + const float caThetaCutBarrelPixelForwardStrip, + const float caThetaCutBarrelStripForwardStrip, + const float caThetaCutBarrelStrip, + const float caThetaCutDefault, + const float dcaCutInnerTripletPixelStrip, + const float dcaCutOuterTripletPixelStrip, + const float dcaCutTripletStrip, + const float dcaCutTripletDefault) const { // detIndex of the layerStart for the Phase1 Pixel Detector: // [BPX1, BPX2, BPX3, BPX4, FP1, FP2, FP3, FN1, FN2, FN3, LAST_VALID] - // [ 0, 96, 320, 672, 1184, 1296, 1408, 1520, 1632, 1744, 1856] + // [ 0, 96, 320, 672, 1184, 1296, 1408, 1520, 1632, 1744, 1856] , 3392 TIB2 auto ri = inner_r(hh); auto zi = inner_z(hh); @@ -176,14 +182,47 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { auto r1 = otherCell.inner_r(hh); auto z1 = otherCell.inner_z(hh); - auto isBarrel = otherCell.outer_detIndex(hh) < TrackerTraits::last_barrel_detIndex; // TODO tune CA cuts below (theta and dca) - bool aligned = areAlignedRZ(r1, z1, ri, zi, ro, zo, ptmin, isBarrel ? caThetaCutBarrel : caThetaCutForward); - return (aligned && dcaCut(hh, - otherCell, - otherCell.inner_detIndex(hh) < TrackerTraits::last_bpix1_detIndex ? dcaCutInnerTriplet - : dcaCutOuterTriplet, - hardCurvCut)); + // Distinguish caThetaCuts for different cases + float caThetaCut; + + auto isOuterBarrelPixel = otherCell.outer_detIndex(hh) < TrackerTraits::last_barrel_detIndex; + auto isInnerBarrelPixel = otherCell.inner_detIndex(hh) < TrackerTraits::last_barrel_detIndex; + auto isOuterForwardPixel = otherCell.outer_detIndex(hh) >= TrackerTraits::last_barrel_detIndex && + otherCell.outer_detIndex(hh) < TrackerTraits::numberOfPixelModules; + auto isInnerForwardPixel = otherCell.inner_detIndex(hh) >= TrackerTraits::last_barrel_detIndex && + otherCell.inner_detIndex(hh) < TrackerTraits::numberOfPixelModules; + auto isOuterBarrelStrip = + otherCell.outer_detIndex(hh) >= TrackerTraits::numberOfPixelModules && otherCell.outer_detIndex(hh) < 3392; + auto isInnerBarrelStrip = + otherCell.inner_detIndex(hh) >= TrackerTraits::numberOfPixelModules && otherCell.inner_detIndex(hh) < 3392; + auto isOuterForwardStrip = otherCell.outer_detIndex(hh) >= 3392; + auto isInnerForwardStrip = otherCell.inner_detIndex(hh) >= 3392; + caThetaCut = (isInnerBarrelPixel && isOuterBarrelPixel) ? caThetaCutBarrel + : (isInnerBarrelPixel && isOuterForwardPixel) ? caThetaCutForward + : (isInnerBarrelPixel && isOuterBarrelStrip) ? caThetaCutBarrelPixelBarrelStrip + : (isInnerBarrelPixel && isOuterForwardStrip) ? caThetaCutBarrelPixelForwardStrip + : (isInnerBarrelStrip && isOuterForwardStrip) ? caThetaCutBarrelStripForwardStrip + : (isInnerBarrelStrip && isOuterBarrelStrip) ? caThetaCutBarrelStrip + : caThetaCutDefault; + + auto isFirstInnerBarrelPixel = otherCell.inner_detIndex(hh) < TrackerTraits::last_bpix1_detIndex; + auto isBeyondFirstInnerBarrelPixel = otherCell.inner_detIndex(hh) > TrackerTraits::last_bpix1_detIndex && + otherCell.inner_detIndex(hh) < TrackerTraits::numberOfPixelModules; + float dcaCutTriplet; + + dcaCutTriplet = + (isFirstInnerBarrelPixel && (isOuterBarrelStrip || isOuterForwardStrip)) ? dcaCutInnerTripletPixelStrip + : (isBeyondFirstInnerBarrelPixel && (isOuterBarrelStrip || isOuterForwardStrip)) + ? dcaCutOuterTripletPixelStrip + : (isFirstInnerBarrelPixel && (isOuterBarrelPixel || isOuterForwardPixel)) ? dcaCutInnerTriplet + : (isBeyondFirstInnerBarrelPixel && (isOuterBarrelPixel || isOuterForwardPixel)) ? dcaCutOuterTriplet + : ((isInnerBarrelStrip || isInnerForwardStrip) && (isOuterBarrelStrip || isOuterForwardStrip)) + ? dcaCutTripletStrip + : dcaCutTripletDefault; + + bool aligned = areAlignedRZ(r1, z1, ri, zi, ro, zo, ptmin, caThetaCut); + return (aligned && dcaCut(hh, otherCell, dcaCutTriplet, hardCurvCut)); } ALPAKA_FN_ACC ALPAKA_FN_INLINE __attribute__((always_inline)) static bool areAlignedRZ( @@ -337,10 +376,13 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ALPAKA_ASSERT_ACC(nh < TrackerTraits::maxHitsOnTrack); hits[nh] = theOuterHitId; auto it = foundNtuplets.bulkFill(acc, apc, hits, nh + 1); + if (it >= 0) { // if negative is overflow.... for (auto c : tmpNtuplet) cells[c].addTrack(acc, it, cellTracks); quality[it] = bad; // initialize to bad + } else { + //printf("Going into overflow from bulkFill"); } } } diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtuplet.cc b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtuplet.cc index dc4639d40e0ee..f983a99ca0424 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtuplet.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtuplet.cc @@ -23,6 +23,9 @@ #include "RecoLocalTracker/Records/interface/PixelCPEFastParamsRecord.h" #include "RecoLocalTracker/SiPixelRecHits/interface/alpaka/PixelCPEFastParamsCollection.h" +#include "RecoLocalTracker/Records/interface/FrameSoARecord.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h" + #include "CAHitNtupletGenerator.h" namespace ALPAKA_ACCELERATOR_NAMESPACE { @@ -45,7 +48,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { private: const edm::ESGetToken tokenField_; - const device::ESGetToken, PixelCPEFastParamsRecord> cpeToken_; + //const device::ESGetToken>, PixelCPEFastParamsRecord> cpeToken_; + const device::ESGetToken frameToken_; const device::EDGetToken tokenHit_; const device::EDPutToken tokenTrack_; @@ -55,7 +59,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template CAHitNtupletAlpaka::CAHitNtupletAlpaka(const edm::ParameterSet& iConfig) : tokenField_(esConsumes()), - cpeToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter("CPE")))), + frameToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter("frameSoA")))), tokenHit_(consumes(iConfig.getParameter("pixelRecHitSrc"))), tokenTrack_(produces()), deviceAlgo_(iConfig) {} @@ -65,11 +69,14 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { edm::ParameterSetDescription desc; desc.add("pixelRecHitSrc", edm::InputTag("siPixelRecHitsPreSplittingAlpaka")); - - std::string cpe = "PixelCPEFastParams"; - cpe += TrackerTraits::nameModifier; + std::string cpe = "DUMMY"; + //cpe += TrackerTraits::nameModifier; desc.add("CPE", cpe); + std::string frame = "FrameSoAPhase1"; + frame += TrackerTraits::nameModifier; + desc.add("frameSoA", frame); + Algo::fillPSetDescription(desc); descriptions.addWithDefaultLabel(desc); } @@ -78,16 +85,19 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { void CAHitNtupletAlpaka::produce(device::Event& iEvent, const device::EventSetup& es) { auto bf = 1. / es.getData(tokenField_).inverseBzAtOriginInGeV(); - auto& fcpe = es.getData(cpeToken_); + auto const& frame = es.getData(frameToken_); auto const& hits = iEvent.get(tokenHit_); - iEvent.emplace(tokenTrack_, deviceAlgo_.makeTuplesAsync(hits, fcpe.const_buffer().data(), bf, iEvent.queue())); + // iEvent.emplace(tokenTrack_, deviceAlgo_.makeTuplesAsync(hits, fcpe.const_buffer().data(), bf, iEvent.queue())); + iEvent.emplace(tokenTrack_, deviceAlgo_.makeTuplesAsync(hits, frame, bf, iEvent.queue())); } using CAHitNtupletAlpakaPhase1 = CAHitNtupletAlpaka; using CAHitNtupletAlpakaHIonPhase1 = CAHitNtupletAlpaka; using CAHitNtupletAlpakaPhase2 = CAHitNtupletAlpaka; + using CAHitNtupletAlpakaPhase1Strip = CAHitNtupletAlpaka; + } // namespace ALPAKA_ACCELERATOR_NAMESPACE #include "HeterogeneousCore/AlpakaCore/interface/alpaka/MakerMacros.h" @@ -95,3 +105,4 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { DEFINE_FWK_ALPAKA_MODULE(CAHitNtupletAlpakaPhase1); DEFINE_FWK_ALPAKA_MODULE(CAHitNtupletAlpakaHIonPhase1); DEFINE_FWK_ALPAKA_MODULE(CAHitNtupletAlpakaPhase2); +DEFINE_FWK_ALPAKA_MODULE(CAHitNtupletAlpakaPhase1Strip); diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.cc b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.cc index a9abc99425a2e..126a2e7571be8 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.cc @@ -21,6 +21,8 @@ #include "CAPixelDoublets.h" #include "CAPixelDoubletsAlgos.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" + namespace ALPAKA_ACCELERATOR_NAMESPACE { namespace { @@ -46,6 +48,17 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ->setComment("Cut on minimum curvature, used in DCA ntuplet selection"); desc.add("dcaCutInnerTriplet", 0.15f)->setComment("Cut on origin radius when the inner hit is on BPix1"); desc.add("dcaCutOuterTriplet", 0.25f)->setComment("Cut on origin radius when the outer hit is on BPix1"); + desc.add("CAThetaCutBarrelPixelBarrelStrip", 0.002f)->setComment("Cut on RZ alignement for Barrel"); + desc.add("CAThetaCutBarrelPixelForwardStrip", 0.003f)->setComment("Cut on RZ alignment for Forward"); + desc.add("CAThetaCutBarrelStripForwardStrip", 0.003f)->setComment("Cut on RZ alignment for Forward"); + desc.add("CAThetaCutBarrelStrip", 0.002f)->setComment("Cut on RZ alignement for Barrel"); + desc.add("CAThetaCutDefault", 0.003f)->setComment("Cut on RZ alignment for Default"); + desc.add("dcaCutInnerTripletPixelStrip", 0.15f) + ->setComment("Cut on origin radius when the inner hit is on BPix1"); + desc.add("dcaCutOuterTripletPixelStrip", 0.25f) + ->setComment("Cut on origin radius when the outer hit is on BPix1"); + desc.add("dcaCutTripletStrip", 0.25f)->setComment("Cut on origin radius when the outer hit is on Strip"); + desc.add("dcaCutTripletDefault", 0.25f)->setComment("Cut on origin radius default"); desc.add("earlyFishbone", true); desc.add("lateFishbone", false); desc.add("fillStatistics", false); @@ -61,6 +74,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { desc.add("doSharedHitCut", true)->setComment("Sharing hit nTuples cleaning"); desc.add("dupPassThrough", false)->setComment("Do not reject duplicate"); desc.add("useSimpleTripletCleaner", true)->setComment("use alternate implementation"); + desc.add("useRemovers", true); } AlgoParams makeCommonParams(edm::ParameterSet const& cfg) { @@ -73,7 +87,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { cfg.getParameter("fillStatistics"), cfg.getParameter("doSharedHitCut"), cfg.getParameter("dupPassThrough"), - cfg.getParameter("useSimpleTripletCleaner")}); + cfg.getParameter("useSimpleTripletCleaner"), + cfg.getParameter("useRemovers")}); } //This is needed to have the partial specialization for isPhase1Topology/isPhase2Topology @@ -90,7 +105,16 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (float)cfg.getParameter("CAThetaCutForward"), (float)cfg.getParameter("hardCurvCut"), (float)cfg.getParameter("dcaCutInnerTriplet"), - (float)cfg.getParameter("dcaCutOuterTriplet")}}; + (float)cfg.getParameter("dcaCutOuterTriplet"), + (float)cfg.getParameter("CAThetaCutBarrelPixelBarrelStrip"), + (float)cfg.getParameter("CAThetaCutBarrelPixelForwardStrip"), + (float)cfg.getParameter("CAThetaCutBarrelStripForwardStrip"), + (float)cfg.getParameter("CAThetaCutBarrelStrip"), + (float)cfg.getParameter("CAThetaCutDefault"), + (float)cfg.getParameter("dcaCutInnerTripletPixelStrip"), + (float)cfg.getParameter("dcaCutOuterTripletPixelStrip"), + (float)cfg.getParameter("dcaCutTripletStrip"), + (float)cfg.getParameter("dcaCutTripletDefault")}}; }; static constexpr ::pixelTrack::QualityCutsT makeQualityCuts(edm::ParameterSet const& pset) { @@ -125,7 +149,16 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { (float)cfg.getParameter("CAThetaCutForward"), (float)cfg.getParameter("hardCurvCut"), (float)cfg.getParameter("dcaCutInnerTriplet"), - (float)cfg.getParameter("dcaCutOuterTriplet")}, + (float)cfg.getParameter("dcaCutOuterTriplet"), + (float)cfg.getParameter("CAThetaCutBarrelPixelBarrelStrip"), + (float)cfg.getParameter("CAThetaCutBarrelPixelForwardStrip"), + (float)cfg.getParameter("CAThetaCutBarrelStripForwardStrip"), + (float)cfg.getParameter("CAThetaCutBarrelStrip"), + (float)cfg.getParameter("CAThetaCutDefault"), + (float)cfg.getParameter("dcaCutInnerTripletPixelStrip"), + (float)cfg.getParameter("dcaCutOuterTripletPixelStrip"), + (float)cfg.getParameter("dcaCutTripletStrip"), + (float)cfg.getParameter("dcaCutTripletDefault")}, {(bool)cfg.getParameter("includeFarForwards")}}; } @@ -150,7 +183,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { cfg.getParameter("idealConditions"), (float)cfg.getParameter("cellZ0Cut"), (float)cfg.getParameter("cellPtCut"), - cfg.getParameter>("phiCuts")}; + cfg.getParameter>("phiCuts"), + cfg.getParameter>("minz"), + cfg.getParameter>("maxz"), + cfg.getParameter>("maxr")}; } } // namespace @@ -197,7 +233,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { fillDescriptionsCommon(desc); desc.add("maxNumberOfDoublets", pixelTopology::Phase1::maxNumberOfDoublets); - desc.add("idealConditions", true); + desc.add("idealConditions", false); desc.add("includeJumpingForwardDoublets", false); desc.add("cellZ0Cut", 12.0); desc.add("cellPtCut", 0.5); @@ -225,6 +261,63 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { "phiCuts", std::vector(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts))) ->setComment("Cuts in phi for cells"); + + desc.add>( + "minz", std::vector(std::begin(phase1PixelTopology::minz), std::end(phase1PixelTopology::minz))) + ->setComment("Cuts in minz for cells"); + desc.add>( + "maxz", std::vector(std::begin(phase1PixelTopology::maxz), std::end(phase1PixelTopology::maxz))) + ->setComment("Cuts in maxz for cells"); + desc.add>( + "maxr", std::vector(std::begin(phase1PixelTopology::maxr), std::end(phase1PixelTopology::maxr))) + ->setComment("Cuts in maxr for cells"); + } + + template <> + void CAHitNtupletGenerator::fillPSetDescription(edm::ParameterSetDescription& desc) { + fillDescriptionsCommon(desc); + + desc.add("maxNumberOfDoublets", pixelTopology::Phase1::maxNumberOfDoublets); + desc.add("idealConditions", false); + desc.add("includeJumpingForwardDoublets", false); + desc.add("cellZ0Cut", 12.0); + desc.add("cellPtCut", 0.5); + + edm::ParameterSetDescription trackQualityCuts; + trackQualityCuts.add("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut"); + trackQualityCuts.add>("chi2Coeff", {0.9, 1.8})->setComment("chi2 at 1GeV and at ptMax above"); + trackQualityCuts.add("chi2Scale", 8.) + ->setComment( + "Factor to multiply the pT-dependent chi2 cut (currently: 8 for the broken line fit, ?? for the Riemann " + "fit)"); + trackQualityCuts.add("tripletMinPt", 1.0)->setComment("Min pT for triplets, in GeV"); + trackQualityCuts.add("tripletMaxTip", 0.3)->setComment("Max |Tip| for triplets, in cm"); + trackQualityCuts.add("tripletMaxZip", 12.)->setComment("Max |Zip| for triplets, in cm"); + trackQualityCuts.add("quadrupletMinPt", 0.5)->setComment("Min pT for quadruplets, in GeV"); + trackQualityCuts.add("quadrupletMaxTip", 0.5)->setComment("Max |Tip| for quadruplets, in cm"); + trackQualityCuts.add("quadrupletMaxZip", 12.)->setComment("Max |Zip| for quadruplets, in cm"); + desc.add("trackQualityCuts", trackQualityCuts) + ->setComment( + "Quality cuts based on the results of the track fit:\n - apply a pT-dependent chi2 cut;\n - apply " + "\"region " + "cuts\" based on the fit results (pT, Tip, Zip)."); + + desc.add>("phiCuts", + std::vector(std::begin(phase1PixelStripTopology::phicuts), + std::end(phase1PixelStripTopology::phicuts))) + ->setComment("Cuts in phi for cells"); + desc.add>( + "minz", + std::vector(std::begin(phase1PixelStripTopology::minz), std::end(phase1PixelStripTopology::minz))) + ->setComment("Cuts in minz for cells"); + desc.add>( + "maxz", + std::vector(std::begin(phase1PixelStripTopology::maxz), std::end(phase1PixelStripTopology::maxz))) + ->setComment("Cuts in maxz for cells"); + desc.add>( + "maxr", + std::vector(std::begin(phase1PixelStripTopology::maxr), std::end(phase1PixelStripTopology::maxr))) + ->setComment("Cuts in maxr for cells"); } template <> @@ -261,6 +354,16 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { "phiCuts", std::vector(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts))) ->setComment("Cuts in phi for cells"); + + desc.add>( + "minz", std::vector(std::begin(phase1PixelTopology::minz), std::end(phase1PixelTopology::minz))) + ->setComment("Cuts in minz for cells"); + desc.add>( + "maxz", std::vector(std::begin(phase1PixelTopology::maxz), std::end(phase1PixelTopology::maxz))) + ->setComment("Cuts in maxz for cells"); + desc.add>( + "maxr", std::vector(std::begin(phase1PixelTopology::maxr), std::end(phase1PixelTopology::maxr))) + ->setComment("Cuts in maxr for cells"); } template <> @@ -288,11 +391,24 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { "phiCuts", std::vector(std::begin(phase2PixelTopology::phicuts), std::end(phase2PixelTopology::phicuts))) ->setComment("Cuts in phi for cells"); + desc.add>( + "minz", std::vector(std::begin(phase2PixelTopology::minz), std::end(phase2PixelTopology::minz))) + ->setComment("Cuts in minz for cells"); + desc.add>( + "maxz", std::vector(std::begin(phase2PixelTopology::maxz), std::end(phase2PixelTopology::maxz))) + ->setComment("Cuts in maxz for cells"); + desc.add>( + "maxr", std::vector(std::begin(phase2PixelTopology::maxr), std::end(phase2PixelTopology::maxr))) + ->setComment("Cuts in maxr for cells"); } template TracksSoACollection CAHitNtupletGenerator::makeTuplesAsync( - HitsOnDevice const& hits_d, ParamsOnDevice const* cpeParams, float bfield, Queue& queue) const { + // HitsOnDevice const& hits_d, ParamsOnDevice const* cpeParams, float bfield, Queue& queue) const { + HitsOnDevice const& hits_d, + FrameOnDevice const& frame, + float bfield, + Queue& queue) const { using HelixFit = HelixFit; using TrackSoA = TracksSoACollection; using GPUKernels = CAHitNtupletGeneratorKernels; @@ -315,10 +431,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { fitter.allocate(kernels.tupleMultiplicity(), tracks.view()); if (m_params.useRiemannFit_) { fitter.launchRiemannKernels( - hits_d.view(), cpeParams, hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue); + hits_d.view(), frame.view(), hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue); } else { fitter.launchBrokenLineKernels( - hits_d.view(), cpeParams, hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue); + hits_d.view(), frame.view(), hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue); } kernels.classifyTuples(hits_d.view(), tracks.view(), queue); #ifdef GPU_DEBUG @@ -332,4 +448,5 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class CAHitNtupletGenerator; template class CAHitNtupletGenerator; template class CAHitNtupletGenerator; + template class CAHitNtupletGenerator; } // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.h b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.h index ec3273a89dee6..91e5dc2a87fa3 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGenerator.h @@ -52,7 +52,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using Params = caHitNtupletGenerator::ParamsT; using Counters = caHitNtupletGenerator::Counters; - using ParamsOnDevice = pixelCPEforDevice::ParamsOnDeviceT; + //using ParamsOnDevice = pixelCPEforDevice::ParamsOnDeviceT>; + using FrameOnDevice = FrameSoACollection; public: CAHitNtupletGenerator(const edm::ParameterSet& cfg); @@ -68,7 +69,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // void endJob(); TkSoADevice makeTuplesAsync(HitsOnDevice const& hits_d, - ParamsOnDevice const* cpeParams, + FrameOnDevice const& frame_d, + //ParamsOnDevice const* cpeParams, float bfield, Queue& queue) const; diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.dev.cc b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.dev.cc index 2d778f5e6e9de..30af199382a23 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.dev.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.dev.cc @@ -1,4 +1,5 @@ // C++ headers + #ifdef DUMP_GPU_TK_TUPLES #include #endif @@ -19,6 +20,7 @@ //#define GPU_DEBUG //#define NTUPLE_DEBUG +//#define NTUPLE_DEBUGS namespace ALPAKA_ACCELERATOR_NAMESPACE { @@ -215,14 +217,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // remove duplicates (tracks that share a doublet) numberOfBlocks = cms::alpakatools::divide_up_by(3 * m_params.caParams_.maxNumberOfDoublets_ / 4, blockSize); workDiv1D = cms::alpakatools::make_workdiv(numberOfBlocks, blockSize); - - alpaka::exec(queue, - workDiv1D, - Kernel_earlyDuplicateRemover{}, - this->device_theCells_.data(), - this->device_nCells_.data(), - tracks_view, - this->m_params.dupPassThrough_); + if (this->m_params.useRemovers_) { + alpaka::exec(queue, + workDiv1D, + Kernel_earlyDuplicateRemover{}, + this->device_theCells_.data(), + this->device_nCells_.data(), + tracks_view, + this->m_params.dupPassThrough_); + } #ifdef GPU_DEBUG alpaka::wait(queue); #endif @@ -391,13 +394,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { // mark duplicates (tracks that share a doublet) numberOfBlocks = cms::alpakatools::divide_up_by(3 * m_params.caParams_.maxNumberOfDoublets_ / 4, blockSize); workDiv1D = cms::alpakatools::make_workdiv(numberOfBlocks, blockSize); - alpaka::exec(queue, - workDiv1D, - Kernel_fastDuplicateRemover{}, - this->device_theCells_.data(), - this->device_nCells_.data(), - tracks_view, - this->m_params.dupPassThrough_); + if (this->m_params.useRemovers_) { + alpaka::exec(queue, + workDiv1D, + Kernel_fastDuplicateRemover{}, + this->device_theCells_.data(), + this->device_nCells_.data(), + tracks_view, + this->m_params.dupPassThrough_); + } #ifdef GPU_DEBUG alpaka::wait(queue); #endif @@ -425,6 +430,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { numberOfBlocks = cms::alpakatools::divide_up_by(3 * TrackerTraits::maxNumberOfQuadruplets / 4, blockSize); // TODO: Check if correct workDiv1D = cms::alpakatools::make_workdiv(numberOfBlocks, blockSize); + alpaka::exec(queue, workDiv1D, Kernel_rejectDuplicate{}, @@ -464,6 +470,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { this->m_params.dupPassThrough_, this->device_hitToTuple_.data()); } + #ifdef GPU_DEBUG alpaka::wait(queue); #endif @@ -557,5 +564,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class CAHitNtupletGeneratorKernels; template class CAHitNtupletGeneratorKernels; template class CAHitNtupletGeneratorKernels; + template class CAHitNtupletGeneratorKernels; } // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.h b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.h index 6999169a46d73..cf18747bb19c8 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernels.h @@ -36,6 +36,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const bool doSharedHitCut_; const bool dupPassThrough_; const bool useSimpleTripletCleaner_; + const bool useRemovers_; }; //CAParams @@ -48,6 +49,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { const float hardCurvCut_; const float dcaCutInnerTriplet_; const float dcaCutOuterTriplet_; + const float CAThetaCutBarrelPixelBarrelStrip_; + const float CAThetaCutBarrelPixelForwardStrip_; + const float CAThetaCutBarrelStripForwardStrip_; + const float CAThetaCutBarrelStrip_; + const float CAThetaCutDefault_; + const float dcaCutInnerTripletPixelStrip_; + const float dcaCutOuterTripletPixelStrip_; + const float dcaCutTripletStrip_; + const float dcaCutTripletDefault_; }; template @@ -60,15 +70,25 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { struct CAParamsT> : public CACommon { /// Is is a starting layer pair? ALPAKA_FN_ACC ALPAKA_FN_INLINE bool startingLayerPair(int16_t pid) const { - return minHitsPerNtuplet_ > 3 ? pid < 3 : pid < 8 || pid > 12; + if constexpr (std::is_same_v) { + return (pid < 12 || pid == 32 || pid == 33 || pid == 29 || pid == 27 || pid == 18 || pid == 37); + } else { + return minHitsPerNtuplet_ > 3 ? pid < 3 : pid < 8 || pid > 12; + } } /// Is this a pair with inner == 0? ALPAKA_FN_ACC ALPAKA_FN_INLINE bool startAt0(int16_t pid) const { - ALPAKA_ASSERT_ACC( - (pixelTopology::Phase1::layerPairs[pid * 2] == 0) == - (pid < 3 || pid == 13 || pid == 15 || pid == 16)); // to be 100% sure it's working, may be removed - return pixelTopology::Phase1::layerPairs[pid * 2] == 0; + if constexpr (std::is_same_v) { + assert((pixelTopology::Phase1Strip::layerPairs[pid * 2] == 0) == + (pid < 3 || pid == 8 || pid == 10 || pid == 11 || + pid == 29)); // to be 100% sure it's working, may be removed + return pixelTopology::Phase1Strip::layerPairs[pid * 2] == 0; + } else { + assert((pixelTopology::Phase1::layerPairs[pid * 2] == 0) == + (pid < 3 || pid == 13 || pid == 15 || pid == 16)); // to be 100% sure it's working, may be removed + return pixelTopology::Phase1::layerPairs[pid * 2] == 0; + } } }; diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernelsImpl.h b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernelsImpl.h index e72d221f7e21c..6c2a855e325ba 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernelsImpl.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAHitNtupletGeneratorKernelsImpl.h @@ -124,6 +124,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { #endif if (cms::alpakatools::once_per_grid(acc)) { +#ifdef GPU_DEBUG if (apc->get().first >= TrackerTraits::maxNumberOfQuadruplets) printf("Tuples overflow\n"); if (*nCells >= maxNumberOfDoublets) @@ -134,7 +135,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { printf("cellTracks overflow\n"); if (int(hitToTuple->nOnes()) < nHits) printf("ERROR hitToTuple overflow %d %d\n", hitToTuple->nOnes(), nHits); -#ifdef GPU_DEBUG printf("size of cellNeighbors %d \n cellTracks %d \n hitToTuple %d \n", cellNeighbors->size(), cellTracks->size(), @@ -198,7 +198,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { TkSoAView tracks_view, bool dupPassThrough) const { // quality to mark rejected - constexpr auto reject = Quality::edup; /// cannot be loose + constexpr auto reject = pixelTrack::Quality::edup; /// cannot be loose ALPAKA_ASSERT_ACC(nCells); for (auto idx : cms::alpakatools::uniform_elements(acc, *nCells)) { auto const &thisCell = cells[idx]; @@ -350,29 +350,63 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { auto zi = thisCell.inner_z(hh); auto ro = thisCell.outer_r(hh); auto zo = thisCell.outer_z(hh); - auto isBarrel = thisCell.inner_detIndex(hh) < TrackerTraits::last_barrel_detIndex; - + float caThetaCut; + auto isOuterBarrelPixel = thisCell.outer_detIndex(hh) < TrackerTraits::last_barrel_detIndex; + auto isInnerBarrelPixel = thisCell.inner_detIndex(hh) < TrackerTraits::last_barrel_detIndex; + auto isOuterForwardPixel = thisCell.outer_detIndex(hh) >= TrackerTraits::last_barrel_detIndex && + thisCell.outer_detIndex(hh) < TrackerTraits::numberOfPixelModules; + auto isOuterBarrelStrip = + thisCell.outer_detIndex(hh) >= TrackerTraits::numberOfPixelModules && thisCell.outer_detIndex(hh) < 3392; + auto isInnerBarrelStrip = + thisCell.inner_detIndex(hh) >= TrackerTraits::numberOfPixelModules && thisCell.inner_detIndex(hh) < 3392; + auto isOuterForwardStrip = thisCell.outer_detIndex(hh) >= 3392; + caThetaCut = (isInnerBarrelPixel && isOuterBarrelPixel) ? params.CAThetaCutBarrel_ + : (isInnerBarrelPixel && isOuterForwardPixel) ? params.CAThetaCutForward_ + : (isInnerBarrelPixel && isOuterBarrelStrip) ? params.CAThetaCutBarrelPixelBarrelStrip_ + : (isInnerBarrelPixel && isOuterForwardStrip) ? params.CAThetaCutBarrelPixelForwardStrip_ + : (isInnerBarrelStrip && isOuterForwardStrip) ? params.CAThetaCutBarrelStripForwardStrip_ + : (isInnerBarrelStrip && isOuterBarrelStrip) ? params.CAThetaCutBarrelStrip_ + : params.CAThetaCutDefault_; // loop on inner cells for (uint32_t j : cms::alpakatools::independent_group_elements_x(acc, numberOfPossibleNeighbors)) { auto otherCell = (vi[j]); auto &oc = cells[otherCell]; auto r1 = oc.inner_r(hh); auto z1 = oc.inner_z(hh); - bool aligned = Cell::areAlignedRZ( - r1, - z1, - ri, - zi, - ro, - zo, - params.ptmin_, - isBarrel ? params.CAThetaCutBarrel_ : params.CAThetaCutForward_); // 2.f*thetaCut); // FIXME tune cuts - if (aligned && - thisCell.dcaCut(hh, - oc, - oc.inner_detIndex(hh) < TrackerTraits::last_bpix1_detIndex ? params.dcaCutInnerTriplet_ - : params.dcaCutOuterTriplet_, - params.hardCurvCut_)) { // FIXME tune cuts + bool aligned = Cell::areAlignedRZ(r1, + z1, + ri, + zi, + ro, + zo, + params.ptmin_, + caThetaCut); // 2.f*thetaCut); // FIXME tune cuts + auto isOuterBarrelPixel = oc.outer_detIndex(hh) < TrackerTraits::last_barrel_detIndex; + auto isOuterForwardPixel = oc.outer_detIndex(hh) >= TrackerTraits::last_barrel_detIndex && + oc.outer_detIndex(hh) < TrackerTraits::numberOfPixelModules; + auto isOuterBarrelStrip = + oc.outer_detIndex(hh) >= TrackerTraits::numberOfPixelModules && oc.outer_detIndex(hh) < 3392; + auto isInnerBarrelStrip = + oc.inner_detIndex(hh) >= TrackerTraits::numberOfPixelModules && oc.inner_detIndex(hh) < 3392; + auto isOuterForwardStrip = oc.outer_detIndex(hh) >= 3392; + auto isInnerForwardStrip = oc.inner_detIndex(hh) >= 3392; + auto isFirstInnerBarrelPixel = oc.inner_detIndex(hh) < TrackerTraits::last_bpix1_detIndex; + auto isBeyondFirstInnerBarrelPixel = oc.inner_detIndex(hh) > TrackerTraits::last_bpix1_detIndex && + oc.inner_detIndex(hh) < TrackerTraits::numberOfPixelModules; + float dcaCutTriplet; + dcaCutTriplet = (isFirstInnerBarrelPixel && (isOuterBarrelStrip || isOuterForwardStrip)) + ? params.dcaCutInnerTripletPixelStrip_ + : (isBeyondFirstInnerBarrelPixel && (isOuterBarrelStrip || isOuterForwardStrip)) + ? params.dcaCutOuterTripletPixelStrip_ + : (isFirstInnerBarrelPixel && (isOuterBarrelPixel || isOuterForwardPixel)) + ? params.dcaCutInnerTriplet_ + : (isBeyondFirstInnerBarrelPixel && (isOuterBarrelPixel || isOuterForwardPixel)) + ? params.dcaCutOuterTriplet_ + : ((isInnerBarrelStrip || isInnerForwardStrip) && (isOuterBarrelStrip || isOuterForwardStrip)) + ? params.dcaCutTripletStrip_ + : params.dcaCutTripletDefault_; + if (aligned && thisCell.dcaCut(hh, oc, dcaCutTriplet, + params.hardCurvCut_)) { // FIXME tune cuts oc.addOuterNeighbor(acc, cellIndex, *cellNeighbors); thisCell.setStatusBits(Cell::StatusBit::kUsed); oc.setStatusBits(Cell::StatusBit::kUsed); @@ -396,12 +430,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { // recursive: not obvious to widen using Cell = CACellT; - #ifdef GPU_DEBUG if (cms::alpakatools::once_per_grid(acc)) printf("starting producing ntuplets from %d cells \n", *nCells); #endif - for (auto idx : cms::alpakatools::uniform_elements(acc, (*nCells))) { auto const &thisCell = cells[idx]; @@ -927,7 +959,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caHitNtupletGeneratorKernels { // mark worse ambiguities for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { auto const it = *ip; - if (tracks_view[it].quality() > reject && reco::isTriplet(tracks_view, it) && it != im) + if (tracks_view[it].quality() > reject && tracks_view.hitIndices().size(it) < 4 && it != im) tracks_view[it].quality() = reject; //no race: simple assignment of the same constant } diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/CAPixelDoubletsAlgos.h b/RecoTracker/PixelSeeding/plugins/alpaka/CAPixelDoubletsAlgos.h index 11363c4e2d58a..01672cd29b6f3 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/CAPixelDoubletsAlgos.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/CAPixelDoubletsAlgos.h @@ -13,6 +13,7 @@ #include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/VecArray.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h" @@ -52,7 +53,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { const bool idealConditions, const float z0Cut, const float ptCut, - const std::vector& phiCutsV) + const std::vector& phiCutsV, + const std::vector& minzV, + const std::vector& maxzV, + const std::vector& maxrV) : doClusterCut_(doClusterCut), doZ0Cut_(doZ0Cut), doPtCut_(doPtCut), @@ -61,6 +65,12 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { ptCut_(ptCut) { assert(phiCutsV.size() == TrackerTraits::nPairs); std::copy(phiCutsV.begin(), phiCutsV.end(), &phiCuts[0]); + assert(minzV.size() == TrackerTraits::nPairs); + std::copy(minzV.begin(), minzV.end(), &minz[0]); + assert(maxzV.size() == TrackerTraits::nPairs); + std::copy(maxzV.begin(), maxzV.end(), &maxz[0]); + assert(maxrV.size() == TrackerTraits::nPairs); + std::copy(maxrV.begin(), maxrV.end(), &maxr[0]); } bool doClusterCut_; @@ -72,6 +82,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { float ptCut_; int phiCuts[T::nPairs]; + int minz[T::nPairs]; + int maxz[T::nPairs]; + int maxr[T::nPairs]; template ALPAKA_FN_ACC ALPAKA_FN_INLINE bool __attribute__((always_inline)) zSizeCut(const TAcc& acc, @@ -93,8 +106,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { auto dz = hh[i].zGlobal() - hh[o].zGlobal(); auto dr = hh[i].rGlobal() - hh[o].rGlobal(); - auto innerBarrel = mi < T::last_barrel_detIndex; - auto onlyBarrel = mo < T::last_barrel_detIndex; + auto innerBarrel = mi < T::last_barrel_detIndex; //|| (mi >= 1856 && mi <=3392); + auto onlyBarrel = mo < T::last_barrel_detIndex; //|| (mo >= 1856 && hh[o].detectorIndex() <=3392); if (not innerBarrel and not onlyBarrel) return false; @@ -196,8 +209,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { ALPAKA_ASSERT_ACC(j < innerLayerCumulativeSize[pairLayerId]); ALPAKA_ASSERT_ACC(0 == pairLayerId || j >= innerLayerCumulativeSize[pairLayerId - 1]); - uint8_t inner = TrackerTraits::layerPairs[2 * pairLayerId]; - uint8_t outer = TrackerTraits::layerPairs[2 * pairLayerId + 1]; + uint8_t inner = TrackerTraits::layerPairs[2 * pairLayerId]; // layer id + uint8_t outer = TrackerTraits::layerPairs[2 * pairLayerId + 1]; // layer id ALPAKA_ASSERT_ACC(outer > inner); auto hoff = PhiBinner::histOff(outer); @@ -208,9 +221,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { ALPAKA_ASSERT_ACC(i < offsets[inner + 1]); // found hit corresponding to our worker thread, now do the job - if (hh[i].detectorIndex() > pixelClustering::maxNumModules) + if ((outer < TrackerTraits::numberOfPixelLayers && hh[i].detectorIndex() > pixelClustering::maxNumModules)) continue; // invalid - /* maybe clever, not effective when zoCut is on auto bpos = (mi%8)/4; // if barrel is 1 for z>0 auto fpos = (outer>3) & (outer<7); @@ -219,10 +231,33 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { auto mez = hh[i].zGlobal(); - if (mez < TrackerTraits::minz[pairLayerId] || mez > TrackerTraits::maxz[pairLayerId]) + if (mez < cuts.minz[pairLayerId] || mez > cuts.maxz[pairLayerId]) continue; - - if (doClusterCut && outer > pixelTopology::last_barrel_layer && cuts.clusterCut(acc, hh, i)) +#ifdef GPU_DEBUG + if (doClusterCut && outer > pixelTopology::last_barrel_layer && outer > TrackerTraits::numberOfPixelLayers) { + bool innerB1 = i < T::last_bpix1_detIndex; + bool isOuterLadder = 0 == (i / 8) % 2; + bool innerB2 = (i >= T::last_bpix1_detIndex) && (i < T::last_bpix2_detIndex); + auto mes = (!innerB1) || isOuterLadder ? hh[i].clusterSizeY() : -1; + //if (hh[oi].detectorIndex() <=3392){ + if (mes > 0 && innerB1) + printf("ClusterSizeY innerB1 StripBarrel: %d", hh[i].clusterSizeY()); + if (mes > 0 && innerB2) + printf("ClusterSizeY innerB2 StripBarrel: %d", hh[i].clusterSizeY()); + } + if (doClusterCut && outer > pixelTopology::last_barrel_layer && outer < TrackerTraits::numberOfPixelLayers) { + bool innerB1 = i < T::last_bpix1_detIndex; + bool isOuterLadder = 0 == (i / 8) % 2; + bool innerB2 = (i >= T::last_bpix1_detIndex) && (i < T::last_bpix2_detIndex); + auto mes = (!innerB1) || isOuterLadder ? hh[i].clusterSizeY() : -1; + if (mes > 0 && innerB1) + printf("ClusterSizeY innerB1 FPix: %d", hh[i].clusterSizeY()); + if (mes > 0 && innerB2) + printf("ClusterSizeY innerB2 FPix: %d", hh[i].clusterSizeY()); + } +#endif + if (doClusterCut && outer > pixelTopology::last_barrel_layer && cuts.clusterCut(acc, hh, i) && + outer < TrackerTraits::numberOfPixelLayers) continue; auto mep = hh[i].iphi(); @@ -240,7 +275,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { auto zo = hh[j].zGlobal(); auto ro = hh[j].rGlobal(); auto dr = ro - mer; - return dr > TrackerTraits::maxr[pairLayerId] || dr < 0 || std::abs((mez * ro - mer * zo)) > z0cut * dr; + return dr > cuts.maxr[pairLayerId] || dr < 0 || std::abs((mez * ro - mer * zo)) > z0cut * dr; }; auto iphicut = cuts.phiCuts[pairLayerId]; @@ -273,31 +308,26 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { ALPAKA_ASSERT_ACC(oi >= offsets[outer]); ALPAKA_ASSERT_ACC(oi < offsets[outer + 1]); auto mo = hh[oi].detectorIndex(); - - // invalid - if (mo > pixelClustering::maxNumModules) + if ((outer < TrackerTraits::numberOfPixelLayers && mo > pixelClustering::maxNumModules)) continue; - if (doZ0Cut && z0cutoff(oi)) continue; auto mop = hh[oi].iphi(); uint16_t idphi = std::min(std::abs(int16_t(mop - mep)), std::abs(int16_t(mep - mop))); - if (idphi > iphicut) continue; - if (doClusterCut && cuts.zSizeCut(acc, hh, i, oi)) + if (doClusterCut && cuts.zSizeCut(acc, hh, i, oi) && oi < TrackerTraits::numberOfPixelLayers) continue; if (doPtCut && ptcut(oi, idphi)) continue; - auto ind = alpaka::atomicAdd(acc, nCells, (uint32_t)1, alpaka::hierarchy::Blocks{}); if (ind >= maxNumOfDoublets) { alpaka::atomicSub(acc, nCells, (uint32_t)1, alpaka::hierarchy::Blocks{}); break; - } // move to SimpleVector?? + } // move to SimpleVector?? */ cells[ind].init(*cellNeighbors, *cellTracks, hh, pairLayerId, i, oi); isOuterHitOfCell[oi].push_back(acc, ind); #ifdef GPU_DEBUG @@ -309,8 +339,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { } // #endif #ifdef GPU_DEBUG - if (tooMany > 0 or tot > 0) - printf("OuterHitOfCell for %d in layer %d/%d, %d,%d %d, %d %.3f %.3f %s\n", + if (tooMany > 0 or tot > 0) { + printf("i,inner,outer,nmin,tot,tooMany,iphicut,cuts.minz[pairLayerId],cuts.maxz[pairLayerId]"); + printf("OuterHitOfCell for %d in layer %d/%d, %d,%d %d, %d %.3d %.3d %s\n", i, inner, outer, @@ -318,9 +349,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets { tot, tooMany, iphicut, - TrackerTraits::minz[pairLayerId], - TrackerTraits::maxz[pairLayerId], + cuts.minz[pairLayerId], + cuts.maxz[pairLayerId], tooMany > 0 ? "FULL!!" : "not full."); + } #endif } // loop in block... } diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.cc b/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.cc index d0fe19233b225..5c0a0c8d29044 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.cc @@ -18,4 +18,5 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class HelixFit; template class HelixFit; template class HelixFit; + template class HelixFit; } // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.h b/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.h index f3e75e83106a7..a9e20f69066ba 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.h +++ b/RecoTracker/PixelSeeding/plugins/alpaka/HelixFit.h @@ -8,9 +8,11 @@ #include "DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h" #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h" #include "RecoTracker/PixelTrackFitting/interface/alpaka/FitResult.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "HeterogeneousCore/AlpakaInterface/interface/config.h" #include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforDevice.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/alpaka/FrameSoACollection.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/FrameSoALayout.h" #include "CAStructures.h" @@ -69,16 +71,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { ~HelixFit() { deallocate(); } void setBField(double bField) { bField_ = bField; } - void launchRiemannKernels(const HitConstView &hv, - ParamsOnDevice const *cpeParams, - uint32_t nhits, - uint32_t maxNumberOfTuples, - Queue &queue); - void launchBrokenLineKernels(const HitConstView &hv, - ParamsOnDevice const *cpeParams, - uint32_t nhits, - uint32_t maxNumberOfTuples, - Queue &queue); + void launchRiemannKernels( + const HitConstView &hv, const FrameSoAConstView &fr, uint32_t nhits, uint32_t maxNumberOfTuples, Queue &queue); + void launchBrokenLineKernels( + const HitConstView &hv, const FrameSoAConstView &fr, uint32_t nhits, uint32_t maxNumberOfTuples, Queue &queue); void allocate(TupleMultiplicity const *tupleMultiplicity, OutputSoAView &helix_fit_results); void deallocate(); diff --git a/RecoTracker/PixelSeeding/plugins/alpaka/RiemannFit.dev.cc b/RecoTracker/PixelSeeding/plugins/alpaka/RiemannFit.dev.cc index 3e6f6e9c8ed98..42cb96e98e4fa 100644 --- a/RecoTracker/PixelSeeding/plugins/alpaka/RiemannFit.dev.cc +++ b/RecoTracker/PixelSeeding/plugins/alpaka/RiemannFit.dev.cc @@ -27,16 +27,18 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { class Kernel_FastFit { public: template >> - ALPAKA_FN_ACC void operator()(TAcc const &acc, - Tuples const *__restrict__ foundNtuplets, - TupleMultiplicity const *__restrict__ tupleMultiplicity, - uint32_t nHits, - TrackingRecHitSoAConstView hh, - pixelCPEforDevice::ParamsOnDeviceT const *__restrict__ cpeParams, - double *__restrict__ phits, - float *__restrict__ phits_ge, - double *__restrict__ pfast_fit, - uint32_t offset) const { + ALPAKA_FN_ACC void operator()( + TAcc const &acc, + Tuples const *__restrict__ foundNtuplets, + TupleMultiplicity const *__restrict__ tupleMultiplicity, + uint32_t nHits, + TrackingRecHitSoAConstView hh, + FrameSoAConstView fr, + // pixelCPEforDevice::ParamsOnDeviceT> const *__restrict__ cpeParams, + double *__restrict__ phits, + float *__restrict__ phits_ge, + double *__restrict__ pfast_fit, + uint32_t offset) const { constexpr uint32_t hitsInFit = N; ALPAKA_ASSERT_ACC(hitsInFit <= nHits); @@ -74,7 +76,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { for (unsigned int i = 0; i < hitsInFit; ++i) { auto hit = hitId[i]; float ge[6]; - cpeParams->detParams(hh[hit].detectorIndex()).frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); + fr.detFrame(hh.detectorIndex(hit)).toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); + // cpeParams->detParams(hh[hit].detectorIndex()).frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge); hits.col(i) << hh[hit].xGlobal(), hh[hit].yGlobal(), hh[hit].zGlobal(); hits_ge.col(i) << ge[0], ge[1], ge[2], ge[3], ge[4], ge[5]; @@ -209,11 +212,15 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { }; template - void HelixFit::launchRiemannKernels(const TrackingRecHitSoAConstView &hv, - pixelCPEforDevice::ParamsOnDeviceT const *cpeParams, - uint32_t nhits, - uint32_t maxNumberOfTuples, - Queue &queue) { + // void HelixFit::launchRiemannKernels(const TrackingRecHitSoAConstView &hv, + // pixelCPEforDevice::ParamsOnDeviceT const *cpeParams, + void HelixFit::launchRiemannKernels( + const TrackingRecHitSoAConstView &hv, + const FrameSoAConstView &fr, + // pixelCPEforDevice::ParamsOnDeviceT> const *cpeParams, + uint32_t nhits, + uint32_t maxNumberOfTuples, + Queue &queue) { assert(tuples_); auto blockSize = 64; @@ -242,7 +249,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tupleMultiplicity_, 3, hv, - cpeParams, + fr, hitsDevice.data(), hits_geDevice.data(), fast_fit_resultsDevice.data(), @@ -281,7 +288,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tupleMultiplicity_, 4, hv, - cpeParams, + fr, hitsDevice.data(), hits_geDevice.data(), fast_fit_resultsDevice.data(), @@ -321,7 +328,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tupleMultiplicity_, 5, hv, - cpeParams, + fr, hitsDevice.data(), hits_geDevice.data(), fast_fit_resultsDevice.data(), @@ -360,7 +367,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { tupleMultiplicity_, 5, hv, - cpeParams, + fr, hitsDevice.data(), hits_geDevice.data(), fast_fit_resultsDevice.data(), @@ -397,5 +404,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class HelixFit; template class HelixFit; template class HelixFit; + template class HelixFit; } // namespace ALPAKA_ACCELERATOR_NAMESPACE diff --git a/RecoTracker/PixelTrackFitting/BuildFile.xml b/RecoTracker/PixelTrackFitting/BuildFile.xml index c21f4634d0308..52aac5fd20ea3 100644 --- a/RecoTracker/PixelTrackFitting/BuildFile.xml +++ b/RecoTracker/PixelTrackFitting/BuildFile.xml @@ -9,7 +9,9 @@ + + diff --git a/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc b/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc index ef65feb4f8b28..556d579cdbc42 100644 --- a/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc +++ b/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc @@ -107,7 +107,7 @@ void PixelTrackProducerFromSoAT::fillDescriptions(edm::Configurat edm::ParameterSetDescription desc; desc.add("beamSpot", edm::InputTag("offlineBeamSpot")); desc.add("trackSrc", edm::InputTag("pixelTracksSoA")); - desc.add("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplittingLegacy")); + desc.add("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplitting")); desc.add("minNumberOfHits", 0); desc.add("minQuality", "loose"); descriptions.addWithDefaultLabel(desc); diff --git a/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoAAlpaka.cc b/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoAAlpaka.cc index 5769c8c53976c..80d7b13979f34 100644 --- a/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoAAlpaka.cc +++ b/RecoTracker/PixelTrackFitting/plugins/PixelTrackProducerFromSoAAlpaka.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include "DataFormats/BeamSpot/interface/BeamSpot.h" @@ -27,7 +28,7 @@ #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" #include "FWCore/Utilities/interface/EDGetToken.h" #include "FWCore/Utilities/interface/InputTag.h" -#include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" #include "RecoTracker/PixelTrackFitting/interface/alpaka/FitUtils.h" @@ -35,6 +36,12 @@ #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h" #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h" +#include "DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h" +#include "Geometry/CommonTopologies/interface/GluedGeomDet.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoTracker/PixelTrackFitting/interface/alpaka/FitUtils.h" #include "storeTracks.h" /** @@ -62,27 +69,34 @@ class PixelTrackProducerFromSoAAlpaka : public edm::global::EDProducer<> { // Event Data tokens const edm::EDGetTokenT tBeamSpot_; + const edm::EDGetTokenT tokenTrack_; - const edm::EDGetTokenT cpuHits_; - const edm::EDGetTokenT hmsToken_; + const edm::EDGetTokenT cpuPixelHits_; + edm::EDGetTokenT cpuStripHits_; + edm::EDGetTokenT hmsToken_; + // Event Setup tokens + const edm::ESGetToken geomToken_; const edm::ESGetToken idealMagneticFieldToken_; const edm::ESGetToken ttTopoToken_; int32_t const minNumberOfHits_; pixelTrack::Quality const minQuality_; + const bool useStripHits_; }; template PixelTrackProducerFromSoAAlpaka::PixelTrackProducerFromSoAAlpaka(const edm::ParameterSet &iConfig) : tBeamSpot_(consumes(iConfig.getParameter("beamSpot"))), tokenTrack_(consumes(iConfig.getParameter("trackSrc"))), - cpuHits_(consumes(iConfig.getParameter("pixelRecHitLegacySrc"))), - hmsToken_(consumes(iConfig.getParameter("pixelRecHitLegacySrc"))), + cpuPixelHits_(consumes(iConfig.getParameter("pixelRecHitLegacySrc"))), + // hmsToken_(consumes(iConfig.getParameter("pixelRecHitLegacySrc"))), + geomToken_(esConsumes()), idealMagneticFieldToken_(esConsumes()), ttTopoToken_(esConsumes()), minNumberOfHits_(iConfig.getParameter("minNumberOfHits")), - minQuality_(pixelTrack::qualityByName(iConfig.getParameter("minQuality"))) { + minQuality_(pixelTrack::qualityByName(iConfig.getParameter("minQuality"))), + useStripHits_(iConfig.getParameter("useStripHits")) { if (minQuality_ == pixelTrack::Quality::notQuality) { throw cms::Exception("PixelTrackConfiguration") << iConfig.getParameter("minQuality") + " is not a pixelTrack::Quality"; @@ -98,6 +112,13 @@ PixelTrackProducerFromSoAAlpaka::PixelTrackProducerFromSoAAlpaka( // around a rare race condition in framework scheduling produces(); produces(); + if (useStripHits_) { + cpuStripHits_ = + consumes(iConfig.getParameter("stripRecHitLegacySrc")); + hmsToken_ = consumes(iConfig.getParameter("hitModuleStartSrc")); + } else { + hmsToken_ = consumes(iConfig.getParameter("pixelRecHitLegacySrc")); + } } template @@ -105,9 +126,14 @@ void PixelTrackProducerFromSoAAlpaka::fillDescriptions(edm::Confi edm::ParameterSetDescription desc; desc.add("beamSpot", edm::InputTag("offlineBeamSpot")); desc.add("trackSrc", edm::InputTag("pixelTracksAlpaka")); - desc.add("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplittingLegacy")); + desc.add("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplitting")); desc.add("minNumberOfHits", 0); desc.add("minQuality", "loose"); + + desc.add("hitModuleStartSrc", edm::InputTag("siPixelRecHitsPreSplittingAlpaka")); + desc.add("useStripHits", false); + desc.add("stripRecHitLegacySrc", edm::InputTag("siStripMatchedRecHits", "matchedRecHit")); + descriptions.addWithDefaultLabel(desc); } @@ -141,36 +167,94 @@ void PixelTrackProducerFromSoAAlpaka::produce(edm::StreamID strea const auto &bsh = iEvent.get(tBeamSpot_); GlobalPoint bs(bsh.x0(), bsh.y0(), bsh.z0()); - auto const &rechits = iEvent.get(cpuHits_); - std::vector hitmap; - auto const &rcs = rechits.data(); - auto const nhits = rcs.size(); + edm::ESHandle theTrackerGeometry = iSetup.getHandle(geomToken_); - hitmap.resize(nhits, nullptr); + auto const &pixelRecHitsDSV = iEvent.get(cpuPixelHits_); + auto const &pixelRechits = pixelRecHitsDSV.data(); + auto nPixelHits = pixelRechits.size(); auto const &hitsModuleStart = iEvent.get(hmsToken_); - for (auto const &hit : rcs) { - auto const &thit = static_cast(hit); - auto const detI = thit.det()->index(); - auto const &clus = thit.firstClusterRef(); - assert(clus.isPixel()); - auto const idx = hitsModuleStart[detI] + clus.pixelCluster().originalId(); - if (idx >= hitmap.size()) - hitmap.resize(idx + 256, nullptr); // only in case of hit overflow in one module - - assert(nullptr == hitmap[idx]); - hitmap[idx] = &hit; + size_t nStripHits = 0; + const edmNew::DetSetVector *stripRechitsDSV = nullptr; + + if (useStripHits_) { + stripRechitsDSV = &iEvent.get(cpuStripHits_); + nStripHits = hitsModuleStart[TrackerTraits::numberOfModules] - hitsModuleStart[TrackerTraits::numberOfPixelModules]; + } + + size_t nhits = nPixelHits + nStripHits; + + std::vector hitmap(nhits, nullptr); + // std::vector counter(nhits, 0); + + for (const auto &moduleHits : pixelRecHitsDSV) { + auto *det = theTrackerGeometry->idToDet(moduleHits.detId()); + const auto detI = det->index(); + for (const auto &hit : moduleHits) { + auto const &clus = hit.firstClusterRef(); + auto const idx = hitsModuleStart[detI] + clus.pixelCluster().originalId(); + + if (idx >= hitsModuleStart[detI + 1]) { + std::cout << "excess pixel hit" << std::endl; + continue; + } + if (idx >= hitmap.size()) + hitmap.resize(idx + 256, nullptr); /// only in case of hit overflow in one module + + assert(nullptr == hitmap[idx]); + + hitmap[idx] = &hit; + // ++counter[idx]; + } } + // for (auto const &hit : pixelRechits) { + // auto const &thit = static_cast(hit); + // auto const detI = thit.det()->index(); + // auto const &clus = thit.firstClusterRef(); + // assert(clus.isPixel()); + // auto const idx = hitsModuleStart[detI] + clus.pixelCluster().originalId(); + // if (idx >= hitmap.size()) + // hitmap.resize(idx + 256, nullptr); // only in case of hit overflow in one module + + // if (nullptr != hitmap[idx]) + // throw std::runtime_error("duplicate hit id: " + std::to_string(idx)); + // hitmap[idx] = &hit; + // ++counter[idx]; + // } + + if (useStripHits_) { + for (const auto &moduleHits : *stripRechitsDSV) { + const GluedGeomDet *theStripDet = + dynamic_cast(theTrackerGeometry->idToDet(moduleHits[0].geographicalId())); + int moduleIdx = TrackerTraits::mapIndex(theStripDet->stereoDet()->index()); + if (moduleIdx >= TrackerTraits::numberOfModules) + continue; + for (auto i = 0u; i < moduleHits.size(); ++i) { + auto j = hitsModuleStart[moduleIdx] + i; + hitmap[j] = &*(moduleHits.begin() + i); + // ++counter[j]; + } + } + } + +#ifdef GPU_DEBUG + std::cout << "hitmap nulls:" << std::count(hitmap.begin(), hitmap.end(), nullptr) << std::endl; +#endif std::vector hits; - hits.reserve(5); + hits.reserve(5); //TODO: Move to the average depending on tracker traits auto const &tsoa = iEvent.get(tokenTrack_); auto const *quality = tsoa.view().quality(); auto const &hitIndices = tsoa.view().hitIndices(); auto nTracks = tsoa.view().nTracks(); +#ifdef GPU_DEBUG + std::cout << "max hit index = " << *std::max_element(hitIndices.begin(), hitIndices.end()) << std::endl; + std::cout << "hitmap.size() = " << hitmap.size() << std::endl; +#endif + tracks.reserve(nTracks); int32_t nt = 0; @@ -190,9 +274,12 @@ void PixelTrackProducerFromSoAAlpaka::produce(edm::StreamID strea indToEdm.resize(sortIdxs.size(), -1); for (const auto &it : sortIdxs) { auto nHits = TracksHelpers::nHits(tsoa.view(), it); + assert(nHits >= 3); auto q = quality[it]; - +#ifdef GPU_DEBUG + std::cout << " nHits " << nHits << " quality: " << int(q) << std::endl; +#endif if (q < minQuality_) continue; if (nHits < minNumberOfHits_) //move to nLayers? @@ -239,7 +326,10 @@ void PixelTrackProducerFromSoAAlpaka::produce(edm::StreamID strea math::XYZVector mom(pp.x(), pp.y(), pp.z()); auto track = std::make_unique(chi2, ndof, pos, mom, gp.charge(), CurvilinearTrajectoryError(mo)); - +#ifdef GPU_DEBUG + std::cout << "chi2 " << chi2 << " ndof: " << ndof << "pos " << pos << " mom " << mom << " gp.charge() " + << gp.charge() << std::endl; +#endif // bad and edup not supported as fit not present or not reliable auto tkq = recoQuality[int(q)]; track->setQuality(tkq); @@ -266,8 +356,10 @@ void PixelTrackProducerFromSoAAlpaka::produce(edm::StreamID strea using PixelTrackProducerFromSoAAlpakaPhase1 = PixelTrackProducerFromSoAAlpaka; using PixelTrackProducerFromSoAAlpakaPhase2 = PixelTrackProducerFromSoAAlpaka; using PixelTrackProducerFromSoAAlpakaHIonPhase1 = PixelTrackProducerFromSoAAlpaka; +using PixelTrackProducerFromSoAAlpakaPhase1Strip = PixelTrackProducerFromSoAAlpaka; #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(PixelTrackProducerFromSoAAlpakaPhase1); DEFINE_FWK_MODULE(PixelTrackProducerFromSoAAlpakaPhase2); DEFINE_FWK_MODULE(PixelTrackProducerFromSoAAlpakaHIonPhase1); +DEFINE_FWK_MODULE(PixelTrackProducerFromSoAAlpakaPhase1Strip); diff --git a/RecoTracker/PixelTrackFitting/python/PixelTracks_cff.py b/RecoTracker/PixelTrackFitting/python/PixelTracks_cff.py index 79b257ab64fa3..efcb3fcb47eed 100644 --- a/RecoTracker/PixelTrackFitting/python/PixelTracks_cff.py +++ b/RecoTracker/PixelTrackFitting/python/PixelTracks_cff.py @@ -137,3 +137,40 @@ # Convert the pixel tracks from SoA to legacy format pixelTracks) ) + +# Patatrack with strip hits (alpaka-only) + +from Configuration.ProcessModifiers.stripNtupletFit_cff import stripNtupletFit +from RecoLocalTracker.Configuration.RecoLocalTracker_cff import striptrackerlocalrecoTask +from RecoLocalTracker.SiStripRecHitConverter.siStripRecHitSoAPhase1_cfi import siStripRecHitSoAPhase1 +from SimTracker.TrackAssociatorProducers.quickTrackAssociatorByHits_cfi import quickTrackAssociatorByHits, quickTrackAssociatorByHitsTrackerHitAssociator +from Validation.RecoTrack.TrackValidation_cff import quickTrackAssociatorByHitsPreSplitting +from RecoTracker.PixelTrackFitting.pixelTrackProducerFromSoAAlpakaPhase1Strip_cfi import pixelTrackProducerFromSoAAlpakaPhase1Strip as _pixelTrackProducerFromSoAAlpakaPhase1Strip + +from RecoTracker.PixelSeeding.caHitNtupletAlpakaPhase1Strip_cfi import caHitNtupletAlpakaPhase1Strip as _pixelTracksAlpakaPhase1Strip + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(pixelTracks, _pixelTrackProducerFromSoAAlpakaPhase1Strip.clone(useStripHits = cms.bool(True), minQuality = cms.string('loose'),hitModuleStartSrc = cms.InputTag("siStripRecHitSoAPhase1"))) + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(pixelTracksAlpaka, _pixelTracksAlpakaPhase1Strip.clone(pixelRecHitSrc = cms.InputTag("siStripRecHitSoAPhase1"))) + +siStripRecHitSoAPhase1Serial = makeSerialClone(siStripRecHitSoAPhase1) + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(pixelTracksAlpakaSerial,makeSerialClone(_pixelTracksAlpakaPhase1Strip, + pixelRecHitSrc = 'siStripRecHitSoAPhase1Serial' +)) + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(quickTrackAssociatorByHits, quickTrackAssociatorByHitsTrackerHitAssociator.clone(cluster2TPSrc = "tpClusterProducerPreSplitting")) + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(pixelTracksTask, cms.Task( + # build legacy strip hits + striptrackerlocalrecoTask, + # mix pixel and strip hits in a SoA + siStripRecHitSoAPhase1, + siStripRecHitSoAPhase1Serial, + # Build the pixel ntuplets and the pixel tracks in SoA format with alpaka on the device + pixelTracksAlpaka, + # # Build the pixel ntuplets and the pixel tracks in SoA format with alpaka on the cpu (if requested by the validation) + pixelTracksAlpakaSerial, + # Convert the pixel tracks from SoA to legacy format + pixelTracks +)) diff --git a/RecoTracker/TkSeedGenerator/plugins/SeedFromConsecutiveHitsCreator.cc b/RecoTracker/TkSeedGenerator/plugins/SeedFromConsecutiveHitsCreator.cc index d7bb7c4f7ae3e..cb99112d392ff 100644 --- a/RecoTracker/TkSeedGenerator/plugins/SeedFromConsecutiveHitsCreator.cc +++ b/RecoTracker/TkSeedGenerator/plugins/SeedFromConsecutiveHitsCreator.cc @@ -13,6 +13,8 @@ #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h" #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h" #include "FWCore/Utilities/interface/Likely.h" +#include "Geometry/CommonTopologies/interface/GluedGeomDet.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" namespace { @@ -170,12 +172,41 @@ void SeedFromConsecutiveHitsCreator::buildSeed(TrajectorySeedCollection& seedCol const TrackingRecHit* hit = nullptr; for (unsigned int iHit = 0; iHit < hits.size(); iHit++) { hit = hits[iHit]->hit(); + // added by Brunella + const GeomDet* baseDet = trackerGeometry_->idToDet(hit->geographicalId()); + const GluedGeomDet* gluedDet = dynamic_cast(baseDet); + //if (iHit == 0) std::cout << "Start analyzing!" << std::endl; + //std::cout << "Everything fine until here!" << std::endl; + //DetId detId = hit->geographicalId(); + //if (detId.subdetId() == PixelSubdetector::PixelBarrel) { + //std::cout << "Hit is in the Pixel Barrel." << std::endl; + //} else if (detId.subdetId() == PixelSubdetector::PixelEndcap) { + // std::cout << "Hit is in the Pixel Endcap." << std::endl; + //} else if (detId.subdetId() == StripSubdetector::TIB) { + // std::cout << "Hit is in the Tracker Inner Barrel." << std::endl; + //} else if (detId.subdetId() == StripSubdetector::TOB) { + // std::cout << "Hit is in the Tracker Outer Barrel." << std::endl; + //} else if (detId.subdetId() == StripSubdetector::TEC) { + // std::cout << "Hit is in the Tracker Endcap." << std::endl; + //} else if (detId.subdetId() == StripSubdetector::TID) { + // std::cout << "Hit is in the Tracker Inner Disc." << std::endl; + //} else { + // std::cout << "Hit is in an unknown detector." << std::endl; + //} + + //pixelTopology::Phase1Strip::mapIndex(trackerGeometry_->idToDet(hit->stereoId()->index()))) TrajectoryStateOnSurface state = - (iHit == 0) ? propagator_->propagate(fts, trackerGeometry_->idToDet(hit->geographicalId())->surface()) - : propagator_->propagate(updatedState, trackerGeometry_->idToDet(hit->geographicalId())->surface()); - if (!state.isValid()) + (gluedDet) + ? ((iHit == 0) + ? propagator_->propagate(fts, trackerGeometry_->idToDet(hit->geographicalId())->surface()) + : propagator_->propagate(updatedState, trackerGeometry_->idToDet(hit->geographicalId())->surface())) + : ((iHit == 0) + ? propagator_->propagate(fts, trackerGeometry_->idToDet(hit->geographicalId())->surface()) + : propagator_->propagate(updatedState, trackerGeometry_->idToDet(hit->geographicalId())->surface())); + if (!state.isValid()) { + //std::cerr << "Error: TrajectoryStateOnSurface is not valid!" << std::endl; return; - + } SeedingHitSet::ConstRecHitPointer tth = hits[iHit]; std::unique_ptr newtth(refitHit(tth, state)); diff --git a/RecoVertex/Configuration/python/RecoPixelVertexing_cff.py b/RecoVertex/Configuration/python/RecoPixelVertexing_cff.py index 96dec2b62dd1b..aea59c3abd64d 100644 --- a/RecoVertex/Configuration/python/RecoPixelVertexing_cff.py +++ b/RecoVertex/Configuration/python/RecoPixelVertexing_cff.py @@ -26,6 +26,12 @@ (pp_on_AA & ~phase2_tracker).toReplaceWith(pixelVerticesAlpaka,_pixelVerticesAlpakaHIonPhase1.clone(doSplitting = False, maxVertices = 32)) from RecoVertex.PixelVertexFinding.pixelVertexFromSoAAlpaka_cfi import pixelVertexFromSoAAlpaka as _pixelVertexFromSoAAlpaka +# strip tracks +from RecoVertex.PixelVertexFinding.pixelVertexProducerAlpakaPhase1Strip_cfi import pixelVertexProducerAlpakaPhase1Strip as _pixelVertexProducerAlpakaPhase1Strip +from Configuration.ProcessModifiers.stripNtupletFit_cff import stripNtupletFit + +(alpaka & stripNtupletFit & ~phase2_tracker).toReplaceWith(pixelVerticesAlpaka, _pixelVertexProducerAlpakaPhase1Strip.clone()) + alpaka.toReplaceWith(pixelVertices, _pixelVertexFromSoAAlpaka.clone()) # pixel vertex SoA producer with alpaka on the cpu, for validation diff --git a/RecoVertex/PixelVertexFinding/plugins/alpaka/PixelVertexProducerAlpaka.cc b/RecoVertex/PixelVertexFinding/plugins/alpaka/PixelVertexProducerAlpaka.cc index 8770239f132ab..4cb392ff2a8eb 100644 --- a/RecoVertex/PixelVertexFinding/plugins/alpaka/PixelVertexProducerAlpaka.cc +++ b/RecoVertex/PixelVertexFinding/plugins/alpaka/PixelVertexProducerAlpaka.cc @@ -1,6 +1,7 @@ #include #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h" +#include "Geometry/CommonTopologies/interface/SimpleSeedingLayersTopology.h" #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Utilities/interface/StreamID.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -107,9 +108,11 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { using PixelVertexProducerAlpakaPhase1 = PixelVertexProducerAlpaka; using PixelVertexProducerAlpakaPhase2 = PixelVertexProducerAlpaka; using PixelVertexProducerAlpakaHIonPhase1 = PixelVertexProducerAlpaka; + using PixelVertexProducerAlpakaPhase1Strip = PixelVertexProducerAlpaka; } // namespace ALPAKA_ACCELERATOR_NAMESPACE DEFINE_FWK_ALPAKA_MODULE(PixelVertexProducerAlpakaPhase1); DEFINE_FWK_ALPAKA_MODULE(PixelVertexProducerAlpakaPhase2); DEFINE_FWK_ALPAKA_MODULE(PixelVertexProducerAlpakaHIonPhase1); +DEFINE_FWK_ALPAKA_MODULE(PixelVertexProducerAlpakaPhase1Strip); diff --git a/RecoVertex/PixelVertexFinding/plugins/alpaka/vertexFinder.dev.cc b/RecoVertex/PixelVertexFinding/plugins/alpaka/vertexFinder.dev.cc index 2350f73a76cbb..9766976c46ddd 100644 --- a/RecoVertex/PixelVertexFinding/plugins/alpaka/vertexFinder.dev.cc +++ b/RecoVertex/PixelVertexFinding/plugins/alpaka/vertexFinder.dev.cc @@ -210,5 +210,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE { template class Producer; template class Producer; template class Producer; + template class Producer; } // namespace vertexFinder } // namespace ALPAKA_ACCELERATOR_NAMESPACE